Merge pull request #26 from sauln/fixup-sparsedm

Fixup sparsedm

Author: ctralie

ripser/pyRips.pxd

@@ -1,10 +1,10 @@
 from libcpp.vector cimport vector
 
 cdef extern from "ripser.cpp":
-	vector[float] pythondm(float* D, int N, int modulus, 
-					int dim_max, float threshold, int do_cocycles)
+	vector[float] rips_dm(float* D, int N, int modulus, 
+					     int dim_max, float threshold, int do_cocycles)
 
 cdef extern from "ripser.cpp":
-	vector[float] pythondmsparse(int* I, int* J, float* V, int NEdges, 
-								 int N, int modulus, int dim_max, 
-								 float threshold, int do_cocycles)
+	vector[float] rips_dm_sparse(int* I, int* J, float* V, int NEdges, 
+							   int N, int modulus, int dim_max, 
+							   float threshold, int do_cocycles)

ripser/pyRipser.pyx

@@ -10,7 +10,7 @@ def doRipsFiltrationDM(np.ndarray[float,ndim=1,mode="c"] DParam not None, int ma
 
 	cdef int N = DParam.shape[0]
 
-	res = pyRips.pythondm(&DParam[0], N, coeff, maxHomDim, thresh, do_cocycles)
+	res = pyRips.rips_dm(&DParam[0], N, coeff, maxHomDim, thresh, do_cocycles)
 
 	return res
 
@@ -20,6 +20,6 @@ def doRipsFiltrationDMSparse(np.ndarray[int,ndim=1,mode="c"] I not None, np.ndar
 
 	cdef int NEdges = I.size
 
-	res = pyRips.pythondmsparse(&I[0], &J[0], &V[0], NEdges, N, coeff, maxHomDim, thresh, do_cocycles)
+	res = pyRips.rips_dm_sparse(&I[0], &J[0], &V[0], NEdges, N, coeff, maxHomDim, thresh, do_cocycles)
 
 	return res

ripser/ripser.cpp

@@ -932,7 +932,7 @@ void ripser<sparse_distance_matrix>::assemble_columns_to_reduce(
 	simplices.swap(next_simplices);
 }
 
-std::vector<value_t> pythondm(float* D, int N, int modulus, int dim_max, float threshold, int do_cocycles) {
+std::vector<value_t> rips_dm(float* D, int N, int modulus, int dim_max, float threshold, int do_cocycles) {
 	//Setup distance matrix and figure out threshold
 	std::vector<value_t> retvec;
 	std::vector<value_t> distances(D, D+N);
@@ -961,7 +961,6 @@ std::vector<value_t> pythondm(float* D, int N, int modulus, int dim_max, float t
 		if (d <= threshold) ++num_edges;
 	}
 
-	std::cout << "dim_max = " << dim_max << "\n";
 	if (threshold >= max) {
 		ripser<compressed_lower_distance_matrix> r(std::move(dist), dim_max, threshold, ratio,
 		                                         modulus, do_cocycles);
@@ -979,7 +978,7 @@ std::vector<value_t> pythondm(float* D, int N, int modulus, int dim_max, float t
 }
 
 
-std::vector<value_t> pythondmsparse(int* I, int* J, float* V, int NEdges,
+std::vector<value_t> rips_dm_sparse(int* I, int* J, float* V, int NEdges,
 								 int N, int modulus, int dim_max, float threshold, int do_cocycles) {
 	//Setup distance matrix and figure out threshold
 	std::vector<value_t> retvec;
@@ -1010,13 +1009,13 @@ int unwrapvector(std::vector<value_t> vec, float** out) {
 extern "C" {
     int cripser(float** out, float* D, int N,
                 int modulus, int dim_max, float threshold, int do_cocycles) {
-        std::vector<value_t> resvec = pythondm(D, N, modulus, dim_max, threshold, do_cocycles);
+        std::vector<value_t> resvec = rips_dm(D, N, modulus, dim_max, threshold, do_cocycles);
         return unwrapvector(resvec, out);
     }
 
     int cripsersparse(float** out, int* I, int* J, float* V, int NEdges, int N,
                       int modulus, int dim_max, float threshold, int do_cocycles) {
-        std::vector<value_t> resvec = pythondmsparse(I, J, V, NEdges, N, modulus, dim_max, threshold, do_cocycles);
+        std::vector<value_t> resvec = rips_dm_sparse(I, J, V, NEdges, N, modulus, dim_max, threshold, do_cocycles);
         return unwrapvector(resvec, out);
     }
 }

ripser/ripser.py

@@ -4,6 +4,7 @@
 """
 
 from itertools import cycle
+import warnings
 
 import matplotlib.pyplot as plt
 import matplotlib as mpl
@@ -63,18 +64,20 @@ def __init__(self, maxdim=1, thresh=np.inf, coeff=2, do_cocycles=False, verbose=
         self.do_cocycles = do_cocycles
         self.verbose = verbose
 
+        # Internal variables
         self.dgm_ = None
         self.cocycles_ = {}
         self.dm_ = None  # Distance matrix
-        self.metric_ = None
-        self.num_edges_ = None #Number of edges added
+        self.metric_ = None 
+        self.num_edges_ = None # Number of edges added
 
         if self.verbose:
-            print("Rips(maxdim={}, thres={}, coef={}, verbose={})".format(
-                maxdim, thresh, coeff, verbose))
+            print("Rips(maxdim={}, thresh={}, coeff={}, do_cocycles={}, verbose={})".format(
+                maxdim, thresh, coeff, do_cocycles, verbose))
 
     def transform(self, X, distance_matrix=False, metric='euclidean'):
-        """Compute persistence diagrams for X data array.
+        """ Compute persistence diagrams for X data array. If X is not a distance matrix,
+            it will be converted to a distance matrix using the chosen metric.
 
         Parameters
         ----------
@@ -93,17 +96,19 @@ def transform(self, X, distance_matrix=False, metric='euclidean'):
 
         if not distance_matrix:
             if X.shape[0] == X.shape[1]:
-                from warnings import warn
-                warn("The input matrix is square, but the distance_matrix flag is off.  Did you mean to indicate that this was a distance matrix?")
+                warnings.warn("The input matrix is square, but the distance_matrix flag is off.  Did you mean to indicate that this was a distance matrix?")
             elif X.shape[0] < X.shape[1]:
-                from warnings import warn
-                warn("The input point cloud has more columns than rows; did you mean to transpose?")
+                warnings.warn("The input point cloud has more columns than rows; did you mean to transpose?")
+
+            self.metric_ = metric
             X = pairwise_distances(X, metric=metric)
         elif sparse.issparse(X):
             #Sparse distance matrix
             X = sparse.csr_matrix.astype(X.tocsr(), dtype=np.float32)
+
         if not (X.shape[0] == X.shape[1]):
             raise Exception('Distance matrix is not square')
+
         self.dm_ = X
 
         dgm = self._compute_rips(X)
@@ -144,19 +149,22 @@ def _compute_rips(self, dm):
             res = DRFDMSparse(coo.row, coo.col, coo.data, n_points, \
                         self.maxdim, self.thresh, self.coeff, self.do_cocycles)
         else:
-            [I, J] = np.meshgrid(np.arange(n_points), np.arange(n_points))
+            I, J = np.meshgrid(np.arange(n_points), np.arange(n_points))
             DParam = np.array(dm[I > J], dtype=np.float32)
             res = DRFDM(DParam, self.maxdim, self.thresh,
                         self.coeff, int(self.do_cocycles))
+
         pds = []
         for dim in range(self.maxdim + 1):
             # Number of homology classes in this dimension
             n_classes = int(res[0])
+
             # First extract the persistence diagram
             res = res[1::]
             pd = np.array(res[0:n_classes*2])
             pds.append(np.reshape(pd, (n_classes, 2)))
             res = res[n_classes*2::]
+
             # Now extract the representative cocycles if they were computed
             if self.do_cocycles and dim > 0:
                 self.cocycles_[dim] = []

test/test_ripser.py

@@ -29,6 +29,27 @@ def test_instantiate(self):
         rip = Rips()
         assert rip is not None
 
+class TestTransform():
+    def test_input_warnings(self):
+        
+        rips = Rips()
+        data = np.random.random((3,10))
+
+        with pytest.warns(UserWarning, match='has more columns than rows') as w:
+            rips.transform(data)
+
+        data = np.random.random((3,3))
+        with pytest.warns(UserWarning, match='input matrix is square, but the distance_matrix') as w:
+            rips.transform(data)
+
+    def test_non_square_dist_matrix(self):
+        rips = Rips()
+        data = np.random.random((3,10))
+
+        with pytest.raises(Exception):
+            rips.transform(data, distance_matrix=True)
+
+
 class TestParams():
     def test_defaults(self):
         data = np.random.random((100,3))