math/fast_cov and math/fast_corr: made recommended changes - remove axis variable (only works on columns of provided matrices), remove print_trace, added spearman correlation

Author: Dave Lahr

.gitignore

@@ -1,2 +1,5 @@
 __pycache__/
 *.pyc
+cmapPy.egg-info/
+.vscode
+.gitignore

cmapPy/math/fast_corr.py

@@ -2,42 +2,56 @@
 import cmapPy.pandasGEXpress.setup_GCToo_logger as setup_logger
 import numpy
 import cmapPy.math.fast_cov as fast_cov
+import pandas
 
 
 logger = logging.getLogger(setup_logger.LOGGER_NAME)
 
 
-def fast_corr(x, y=None, axis=0, do_print_trace=False):
-    """calculate the pearson correlation matrix for x, or optionally, the correlaton matrix between x and y.
-    axis parameter defines whether to treat the rows or columns as vectors that the covariance is calculated
-    between
+def fast_corr(x, y=None):
+    """calculate the pearson correlation matrix for the columns of x (MxN), or optionally, the correlaton matrix between x and y (OxP).
+    In the language of statistics the columns are the variables and the rows are the observations.
 
     Args:
-        x (numpy array-like) NxM in shape
-        y (numpy array-like) OxP in shape
-        axis (integer) indicates which axis of the arrays is the vectors.  axis=0 indicates calculate
-            the covariance between the rows, axis=1 indicates calculate the covariance between the columns
-        do_print_trace (bool) whether or not to use logger.debug to print out intermediate mathematical steps
-            (WARNING: will print out all intermediate matrices, potentially lots of output!)
-
-        returns (numpy array-like) QxR array of the covariance values
-            for defaults (y=None, axis=0), shape is NxN
-            if y=None and axis=1, shape is MxM
-            if y is provied, axis=0, shape is NxO (and P==M in other words number of columns in x & y are the same)
-            if y is provided, axis=1, shape is MxP (and N==O in other words number of rows in x & y are the same)
+        x (numpy array-like) MxN in shape
+        y (optional, numpy array-like) OxP in shape
+
+        returns (numpy array-like) array of the covariance values
+            for defaults (y=None), shape is NxN
+            if y is provied, shape is NxP
     """
     if y is None:
         y = x
 
-    cov_mat = fast_cov.fast_cov(x, y, axis, do_print_trace)
-    if do_print_trace:  logger.debug("cov_mat:\n{}".format(cov_mat))
+    cov_mat = fast_cov.fast_cov(x, y)
 
-    std_x = numpy.std(x.T, axis=axis, ddof=1)
-    if do_print_trace:  logger.debug("std_x:\n{}".format(std_x))
-    std_y = numpy.std(y.T, axis=axis, ddof=1)
-    if do_print_trace:  logger.debug("std_y:\n{}".format(std_y))
+    std_x = numpy.std(x, axis=0, ddof=1)
+    std_y = numpy.std(y, axis=0, ddof=1)
 
     std_outer = numpy.outer(std_x, std_y)
-    if do_print_trace:  logger.debug("std_outer:\n{}".format(std_outer))
 
     return cov_mat / std_outer
+
+
+def fast_spearman(x, y=None):
+    """calculate the spearnab correlation matrix for the columns of x (MxN), or optionally, the spearmancorrelaton matrix between x and y (OxP).
+    In the language of statistics the columns are the variables and the rows are the observations.
+
+    Args:
+        x (numpy array-like) MxN in shape
+        y (optional, numpy array-like) OxP in shape
+
+        returns:
+            (numpy array-like) array of the covariance values
+                for defaults (y=None), shape is NxN
+                if y is provied, shape is NxP
+    """
+    logger.debug("x.shape:  {}".format(x.shape))
+    if hasattr(y, "shape"):  
+        logger.debug("y.shape:  {}".format(y.shape))
+
+    x_ranks = pandas.DataFrame(x).rank(method="average").values
+    logger.debug("some min and max ranks of x_ranks:\n{}\n{}".format(numpy.min(x_ranks[:10], axis=0), numpy.max(x_ranks[:10], axis=0)))
+    y_ranks = pandas.DataFrame(y).rank(method="average").values if y is not None else None
+
+    return fast_corr(x_ranks, y_ranks)

cmapPy/math/fast_cov.py

@@ -6,59 +6,38 @@
 logger = logging.getLogger(setup_logger.LOGGER_NAME)
 
 
-def fast_cov(x, y=None, axis=0, do_print_trace=False):
-    """calculate the covariance matrix for x, or optionally, the covariance matrix between x and y.
-    axis parameter defines whether to treat the rows or columns as vectors that the correlation is calculated
-    between
+def fast_cov(x, y=None):
+    """calculate the covariance matrix for the columns of x (MxN), or optionally, the covariance matrix between the
+    columns of x and and the columns of y (MxP).  (In the language of statistics, the columns are variables, the rows
+    are observations).
 
     Args:
-        x (numpy array-like) NxM in shape
-        y (numpy array-like) OxP in shape
-        axis (integer) indicates which axis of the arrays is the vectors.  axis=0 indicates calculate
-            the covariance between the rows, axis=1 indicates calculate the covariance between the columns
-        do_print_trace (bool) whether or not to use logger.debug to print out intermediate mathematical steps
-            (WARNING: will print out all intermediate matrices, potentially lots of output!)
-
-        returns (numpy array-like) QxR array of the covariance values
-            for defaults (y=None, axis=0), shape is NxN
-            if y=None and axis=1, shape is MxM
-            if y is provied, axis=0, shape is NxO (and P==M in other words number of columns in x & y are the same)
-            if y is provided, axis=1, shape is MxP (and N==O in other words number of rows in x & y are the same)
+        x (numpy array-like) MxN in shape
+        y (numpy array-like) MxP in shape
+
+        returns (numpy array-like) array of the covariance values
+            for defaults (y=None), shape is NxN
+            if y is provided, shape is NxP
     """
-    validate_axis(axis)
-    validate_x_y(x, y, axis)
+    validate_x_y(x, y)
 
     if y is None:
         y = x
 
-    mean_x = numpy.mean(x.T, axis=axis)
-    if do_print_trace:  logger.debug("mean_x:\n{}".format(mean_x))
-    mean_y = numpy.mean(y.T, axis=axis)
-    if do_print_trace:  logger.debug("mean_y:\n{}".format(mean_y))
-    
-    mean_centered_x = (x.T - mean_x).T if axis == 0 else (x - mean_x).T
-    if do_print_trace:  logger.debug("mean_centered_x:\n{}".format(mean_centered_x))
-    mean_centered_y = y.T - mean_y if axis == 0 else y - mean_y
-    if do_print_trace:  logger.debug("mean_centered_y:\n{}".format(mean_centered_y))
-    
-    dotprod = numpy.dot(mean_centered_x, mean_centered_y)
-    if do_print_trace:  logger.debug("dotprod:\n{}".format(dotprod))
+    mean_x = numpy.mean(x, axis=0)    
+    mean_y = numpy.mean(y, axis=0)
+
+    mean_centered_x = x - mean_x
+    mean_centered_y = y - mean_y
 
-    other_axis = 1 if axis == 0 else 0
-    if do_print_trace:  logger.debug("other_axis:\n{}".format(other_axis))
+    dotprod = numpy.dot(mean_centered_x.T, mean_centered_y)
 
-    denom = float(x.shape[other_axis] - 1)
-    if do_print_trace:  logger.debug("denom:\n{}".format(denom))
+    denom = x.shape[0] - 1
 
     return dotprod / denom
 
 
-def validate_axis(axis):
-    if not (0 == axis or 1 == axis):
-        raise CmapPyMathFastCovInvalidInputAxis("invalid value for axis provided to fast_cov - axis:  {}".format(axis))
-
-
-def validate_x_y(x, y, axis):
+def validate_x_y(x, y):
     error_msg = ""
 
     if not hasattr(x, "shape"):
@@ -67,14 +46,8 @@ def validate_x_y(x, y, axis):
     if y is not None:
         if not hasattr(y, "shape"):
             error_msg += "y needs to be numpy array-like but it does not have \"shape\" attribute - type(y):  {}\n".format(type(y))
-        else:
-            base_msg = "with axis={} (finding the covariance between {}), the number of {} in x and y must be the same but they are not - x.shape:  {}  y.shape:  {}\n"
-            if 0 == axis:
-                if x.shape[1] != y.shape[1]:
-                    error_msg += base_msg.format(axis, "rows", "columns", x.shape, y.shape)
-            else:
-                if x.shape[0] != y.shape[0]:
-                    error_msg += base_msg.format(axis, "columns", "rows", x.shape, y.shape)
+        elif x.shape[0] != y.shape[0]:
+            error_msg += "the number of rows in the x and y matrices must be the same".format(x.shape, y.shape)
 
     if error_msg != "":
         raise CmapPyMathFastCovInvalidInputXY(error_msg)

cmapPy/math/tests/test_fast_corr.py

@@ -3,6 +3,7 @@
 import cmapPy.pandasGEXpress.setup_GCToo_logger as setup_logger
 import cmapPy.math.fast_corr as fast_corr
 import numpy
+import pandas
 
 
 logger = logging.getLogger(setup_logger.LOGGER_NAME)
@@ -29,77 +30,25 @@ def test_fast_corr_just_x(self):
         logger.debug("*************happy path just x")
         x, _ = TestFastCorr.build_standard_x_y()
 
-        ex = numpy.corrcoef(x)
+        ex = numpy.corrcoef(x, rowvar=False)
         logger.debug("expected ex:  {}".format(ex))
 
-        r = fast_corr.fast_corr(x, do_print_trace=True)
+        r = fast_corr.fast_corr(x)
         logger.debug("r:  {}".format(r))
 
         self.assertTrue(numpy.allclose(ex, r))
 
-        #happy path just x, transposed
-        ex = numpy.corrcoef(x.T)
-        logger.debug("happy path just x, transposed, expected ex:  {}".format(ex))
-        r = fast_corr.fast_corr(x.T, do_print_trace=True)
-        logger.debug("r:  {}".format(r))
-        self.assertTrue(numpy.allclose(ex, r))
-
-    def test_fast_corr_just_x_different_axis(self):
-        logger.debug("*************happy path just x, use different axis")
-        x, _ = TestFastCorr.build_standard_x_y()
-
-        ex = numpy.corrcoef(x, rowvar=False)
-        logger.debug("happy path just x, use different axis, expected ex:  {}".format(ex))
-        r = fast_corr.fast_corr(x, axis=1, do_print_trace=True)
-        logger.debug("r:  {}".format(r))
-        self.assertTrue(numpy.allclose(ex, r))
-
-        #happy path just x, use different axis, transpose
-        ex = numpy.corrcoef(x)
-        logger.debug("happy path just x, use different axis, transpose, expected ex:  {}".format(ex))
-        r = fast_corr.fast_corr(x.T, axis=1, do_print_trace=True)
+        #happy path just x, other direction
+        ex = numpy.corrcoef(x, rowvar=True)
+        logger.debug("happy path just x, other direction, expected ex:  {}".format(ex))
+        r = fast_corr.fast_corr(x.T)
         logger.debug("r:  {}".format(r))
         self.assertTrue(numpy.allclose(ex, r))
 
     def test_fast_corr_x_and_y(self):
         logger.debug("*************happy path x and y")
         x, y = TestFastCorr.build_standard_x_y()
 
-        combined = numpy.vstack([x, y])
-        logger.debug("combined:  {}".format(combined))
-        logger.debug("combined.shape:  {}".format(combined.shape))
-
-        off_diag_ind = combined.shape[0] / 2
-
-        raw_ex = numpy.corrcoef(combined)
-        logger.debug("raw expected produced from numpy.cov on full combined - raw_ex:  {}".format(raw_ex))
-        ex = raw_ex[:off_diag_ind, off_diag_ind:]
-        logger.debug("expected ex:  {}".format(ex))
-
-        r = fast_corr.fast_corr(x, y, do_print_trace=True)
-        logger.debug("r:  {}".format(r))
-        self.assertTrue(numpy.allclose(ex, r))
-
-        #happy path x and y, transposed
-        combined = numpy.vstack([x.T, y.T])
-        logger.debug("*************happy path x and y, transposed - combined:  {}".format(combined))
-        logger.debug("combined.shape:  {}".format(combined.shape))
-
-        off_diag_ind = combined.shape[0] / 2
-
-        raw_ex = numpy.corrcoef(combined)
-        logger.debug("raw expected produced from numpy.cov on full combined - raw_ex:  {}".format(raw_ex))
-        ex = raw_ex[:off_diag_ind, off_diag_ind:]
-        logger.debug("expected ex:  {}".format(ex))
-
-        r = fast_corr.fast_corr(x.T, y.T, do_print_trace=True)
-        logger.debug("r:  {}".format(r))
-        self.assertTrue(numpy.allclose(ex, r))
-
-    def test_fast_corr_x_and_y_different_axis(self):
-        logger.debug("*************happy path, x and y using different axis")
-        x, y = TestFastCorr.build_standard_x_y()
-
         combined = numpy.hstack([x, y])
         logger.debug("combined:  {}".format(combined))
         logger.debug("combined.shape:  {}".format(combined.shape))
@@ -111,16 +60,13 @@ def test_fast_corr_x_and_y_different_axis(self):
         ex = raw_ex[:off_diag_ind, off_diag_ind:]
         logger.debug("expected ex:  {}".format(ex))
 
-        r = fast_corr.fast_corr(x, y, axis=1, do_print_trace=True)
+        r = fast_corr.fast_corr(x, y)
         logger.debug("r:  {}".format(r))
         self.assertTrue(numpy.allclose(ex, r))
 
-    def test_fast_corr_x_and_y_different_axis_transpose(self):
-        logger.debug("*************happy path x and y different axis, transposed")
-        x, y = TestFastCorr.build_standard_x_y()
-
+        #happy path x and y, other direction
         combined = numpy.hstack([x.T, y.T])
-        logger.debug("combined:  {}".format(combined))
+        logger.debug("*************happy path x and y, other direction - combined:  {}".format(combined))
         logger.debug("combined.shape:  {}".format(combined.shape))
 
         off_diag_ind = combined.shape[1] / 2
@@ -130,7 +76,7 @@ def test_fast_corr_x_and_y_different_axis_transpose(self):
         ex = raw_ex[:off_diag_ind, off_diag_ind:]
         logger.debug("expected ex:  {}".format(ex))
 
-        r = fast_corr.fast_corr(x.T, y.T, axis=1, do_print_trace=True)
+        r = fast_corr.fast_corr(x.T, y.T)
         logger.debug("r:  {}".format(r))
         self.assertTrue(numpy.allclose(ex, r))
 
@@ -153,56 +99,53 @@ def test_fast_corr_x_and_y_different_shapes(self):
         logger.debug("expected ex:  {}".format(ex))
         logger.debug("ex.shape:  {}".format(ex.shape))
 
-        r = fast_corr.fast_corr(x, y, axis=1, do_print_trace=True)
+        r = fast_corr.fast_corr(x, y)
         logger.debug("r:  {}".format(r))
         self.assertTrue(numpy.allclose(ex, r))
 
     def test_fast_corr_functional(self):
         logger.debug("*************happy path functional test using randomly generated matrices")
 
-        axes = numpy.random.choice([0, 1], size=num_iterations_functional_tests)
-
-        for cur_axis in axes:
-            logger.debug("*******cur_axis:  {}".format(cur_axis))
-
+        for i in xrange(num_iterations_functional_tests):
             #the dimension containing the observations must have at least size 2
-            x_shape = [numpy.random.randint(1, max_dimension_functional_tests),
-                numpy.random.randint(2, max_dimension_functional_tests)]
-            if cur_axis == 1:
-                x_shape.reverse()
+            x_shape = [numpy.random.randint(2, max_dimension_functional_tests),
+                numpy.random.randint(1, max_dimension_functional_tests)]
             logger.debug("x_shape:  {}".format(x_shape))
 
             x = numpy.random.rand(x_shape[0], x_shape[1]) * numpy.random.randint(1, multiplier_max_functional_tests, size=1)
             logger.debug("x:\n{}".format(x))
 
             y_other_shape = numpy.random.randint(1, max_dimension_functional_tests, size=1)[0]
-            y_shape = (y_other_shape, x_shape[1]) if cur_axis == 0 else (x_shape[0], y_other_shape)
+            y_shape = (x_shape[0], y_other_shape)
             logger.debug("y_shape:  {}".format(y_shape))
             y = numpy.random.rand(y_shape[0], y_shape[1]) * numpy.random.randint(1, multiplier_max_functional_tests, size=1)
             logger.debug("y:\n{}".format(y))
 
-            if cur_axis == 0:
-                row_var = True
-                combined = numpy.vstack([x, y])
-                off_diag_ind = (x.shape[0], -y.shape[0])
-            else:
-                row_var = False
-                combined = numpy.hstack([x, y])
-                off_diag_ind = (x.shape[1], -y.shape[1])
+            combined = numpy.hstack([x, y])
 
-            raw_ex = numpy.corrcoef(combined, rowvar=row_var)
+            raw_ex = numpy.corrcoef(combined, rowvar=False)
             logger.debug("raw_ex.shape:  {}".format(raw_ex.shape))
 
-            ex = raw_ex[:off_diag_ind[0], off_diag_ind[1]:]
+            ex = raw_ex[:x.shape[1], -y.shape[1]:]
             logger.debug("ex:\n{}".format(ex))
             logger.debug("ex.shape:  {}".format(ex.shape))
 
-            r = fast_corr.fast_corr(x, y, axis=cur_axis)
+            r = fast_corr.fast_corr(x, y)
             logger.debug("r:\n{}".format(r))
             logger.debug("r.shape:  {}".format(r.shape))
 
             self.assertTrue(numpy.allclose(ex, r))
 
+    def test_fast_spearman(self):
+        x, y = TestFastCorr.build_standard_x_y()
+
+        ex = numpy.array([[1.0, 1.0, 1.0], [-1.0, -1.0, -1.0], [1.0, 1.0, 1.0]])
+
+        r = fast_corr.fast_spearman(x, y)
+        logger.debug("r:  {}".format(r))
+
+        self.assertTrue(numpy.allclose(ex, r))
+
 
 if __name__ == "__main__":
     setup_logger.setup(verbose=True)