Merge pull request #27 from meyer-lab/missing

Missing matrix imputation with TruncatedSVD; add and edit test cases

Author: cyrillustan

tensorpack/decomposition.py

@@ -1,10 +1,35 @@
 import pickle
 import numpy as np
 import pandas as pd
+from numpy.linalg import norm
 from sklearn.decomposition import TruncatedSVD
 from .cmtf import perform_CP, calcR2X
 
 
+def impute_missing_mat(dat):
+    miss_idx = np.where(~np.isfinite(dat))
+    if len(miss_idx[0]) <= 0:
+        return dat
+    assert np.all(np.any(np.isfinite(dat), axis=0)), "Cannot impute if an entire column is empty"
+    assert np.all(np.any(np.isfinite(dat), axis=1)), "Cannot impute if an entire row is empty"
+
+    imp = np.copy(dat)
+    col_mean = np.nanmean(dat, axis=0, keepdims=True)
+    imp[miss_idx] = np.take(col_mean, miss_idx[1])
+
+    diff = 1.0
+    while diff > 1e-3:
+        tsvd = TruncatedSVD(n_components=min(dat.shape)-1)
+        scores = tsvd.fit_transform(imp)
+        loadings = tsvd.components_
+        recon = scores @ loadings
+        new_diff = norm(imp[miss_idx] - recon[miss_idx]) / norm(recon[miss_idx])
+        assert new_diff < diff, "Matrix imputation difference is not decreasing"
+        diff = new_diff
+        imp[miss_idx] = recon[miss_idx]
+    return imp
+
+
 class Decomposition():
     def __init__(self, data, max_rr=6):
         self.data = data
@@ -20,6 +45,8 @@ def perform_tfac(self):
     def perform_PCA(self, flattenon=0):
         dataShape = self.data.shape
         flatData = np.reshape(np.moveaxis(self.data, flattenon, 0), (dataShape[flattenon], -1))
+        if not np.all(np.isfinite(flatData)):
+            flatData = impute_missing_mat(flatData)
 
         tsvd = TruncatedSVD(n_components=max(self.rrs))
         scores = tsvd.fit_transform(flatData)

tensorpack/test/test_decomposition.py

@@ -2,9 +2,13 @@
 Testing Decomposition
 """
 
+import numpy as np
+import os
+import tensorly as tl
+from tensorly.random import random_cp
 from ..decomposition import Decomposition
 from .atyeo import createCube
-import os
+from ..cmtf import perform_CP, calcR2X
 
 def test_decomp_obj():
     a = Decomposition(createCube())
@@ -21,3 +25,32 @@ def test_decomp_obj():
     assert len(b.PCAR2X) == len(b.sizePCA)
     assert len(b.TR2X) == len(b.sizeT)
     os.remove(fname)
+
+
+def test_missing_obj():
+    for miss_rate in [0.1, 0.2]:
+        dat = createCube()
+        filter = np.random.rand(*dat.shape) > 1-miss_rate
+        dat[filter] = np.nan
+        a = Decomposition(dat)
+        a.perform_tfac()
+        a.perform_PCA()
+        assert len(a.PCAR2X) == len(a.sizePCA)
+        assert len(a.TR2X) == len(a.sizeT)
+
+
+def test_known_rank():
+    shape = (100, 80, 60)
+    tFacOrig = random_cp(shape, 10, full=False)
+    tOrig = tl.cp_to_tensor(tFacOrig)
+    assert calcR2X(tFacOrig, tOrig) >= 1.0
+
+    newtFac = [calcR2X(perform_CP(tOrig, r=rr), tOrig) for rr in [1,3,5,7,9]]
+    assert np.all([newtFac[ii+1] > newtFac[ii]  for ii in range(len(newtFac)-1)])
+    assert newtFac[0] > 0.0
+    assert newtFac[-1] < 1.0
+
+    filter = np.random.rand(*shape) > 0.8
+    missT = np.copy(tOrig)
+    missT[filter] = np.nan
+    pass