More cleaning (#34)

Remove pandas, switch over to SVD impute, add progress bars.

Author: aarmey

poetry.lock

@@ -10,7 +10,7 @@ python = ">=3.9,<3.11"
 numpy = "^1.21"
 tensorly = "^0.7"
 matplotlib = "^3.5.0"
-pandas = "^1.3.5"
+tqdm = "^4.62.3"
 
 [tool.poetry.dev-dependencies]
 pytest = "^6.2"

pyproject.toml

@@ -13,38 +13,70 @@
 from tensorly import partial_svd
 import numpy as np
 
-from .soft_impute import Solver
 
-F32PREC = np.finfo(np.float32).eps
-
-def masked_mae(X_true, X_pred, mask):
-    masked_diff = X_true[mask] - X_pred[mask]
-    return np.mean(np.abs(masked_diff))
-
-
-class IterativeSVD(Solver):
+class IterativeSVD(object):
     def __init__(
             self,
             rank,
-            convergence_threshold=0.00001,
-            max_iters=200,
-            init_fill_method="zero",
+            convergence_threshold=1e-7,
+            max_iters=500,
             random_state=None,
             min_value=None,
             max_value=None,
             verbose=False):
-        Solver.__init__(
-            self,
-            fill_method=init_fill_method,
-            min_value=min_value,
-            max_value=max_value)
+        self.min_value=min_value
+        self.max_value=max_value
         self.rank = rank
         self.max_iters = max_iters
         self.convergence_threshold = convergence_threshold
         self.verbose = verbose
         self.random_state = random_state
 
+    def clip(self, X):
+        """
+        Clip values to fall within any global or column-wise min/max constraints
+        """
+        X = np.asarray(X)
+        if self.min_value is not None:
+            X[X < self.min_value] = self.min_value
+        if self.max_value is not None:
+            X[X > self.max_value] = self.max_value
+        return X
+
+    def prepare_input_data(self, X):
+        """
+        Check to make sure that the input matrix and its mask of missing
+        values are valid. Returns X and missing mask.
+        """
+        if X.dtype != "f" and X.dtype != "d":
+            X = X.astype(float)
+
+        assert X.ndim == 2
+        missing_mask = np.isnan(X)
+        assert not missing_mask.all()
+        return X, missing_mask
+
+    def fit_transform(self, X, y=None):
+        """
+        Fit the imputer and then transform input `X`
+        Note: all imputations should have a `fit_transform` method,
+        but only some (like IterativeImputer in sklearn) also support inductive
+        mode using `fit` or `fit_transform` on `X_train` and then `transform`
+        on new `X_test`.
+        """
+        X_original, missing_mask = self.prepare_input_data(X)
+        observed_mask = ~missing_mask
+        X_filled = X_original.copy()
+        X_filled[missing_mask] = 0.0
+        assert isinstance(X_filled, np.ndarray)
+        X_result = self.solve(X_filled, missing_mask)
+        assert isinstance(X_result, np.ndarray)
+        X_result = self.clip(np.asarray(X_result))
+        X_result[observed_mask] = X_original[observed_mask]
+        return X_result
+
     def _converged(self, X_old, X_new, missing_mask):
+        F32PREC = np.finfo(np.float32).eps
         # check for convergence
         old_missing_values = X_old[missing_mask]
         new_missing_values = X_new[missing_mask]
@@ -59,19 +91,17 @@ def _converged(self, X_old, X_new, missing_mask):
             return (ssd / old_norm_squared) < self.convergence_threshold
 
     def solve(self, X, missing_mask):
-        # X = check_array(X, force_all_finite=False)
-
         observed_mask = ~missing_mask
         X_filled = X
         for i in range(self.max_iters):
             curr_rank = self.rank
-            U, S, V = partial_svd(X_filled, curr_rank, random_state=self.random_state)
-            X_reconstructed = U @ np.diag(S) @ V
+            self.U, S, V = partial_svd(X_filled, curr_rank, random_state=self.random_state)
+            X_reconstructed = self.U @ np.diag(S) @ V
             X_reconstructed = self.clip(X_reconstructed)
-            mae = masked_mae(
-                X_true=X,
-                X_pred=X_reconstructed,
-                mask=observed_mask)
+
+            # Masked mae
+            mae = np.mean(np.abs(X[observed_mask] - X_reconstructed[observed_mask]))
+
             if self.verbose:
                 print(
                     "[IterativeSVD] Iter %d: observed MAE=%0.6f" % (

tensorpack/SVD_impute.py

@@ -8,8 +8,8 @@
 from tensorly.tenalg import khatri_rao
 from copy import deepcopy
 from tensorly.decomposition._cp import initialize_cp, parafac
+from tqdm import tqdm
 from .SVD_impute import IterativeSVD
-from .soft_impute import SoftImpute
 
 
 tl.set_backend('numpy')
@@ -176,8 +176,9 @@ def initialize_cmtf(tensor: np.ndarray, matrix: np.ndarray, rank: int):
     if np.sum(~np.isfinite(unfold)) > 0:
         si = IterativeSVD(rank=rank, random_state=1)
         unfold = si.fit_transform(unfold)
-
-    factors[0] = np.linalg.svd(unfold)[0][:, :rank]
+        factors[0] = si.U
+    else:
+        factors[0] = np.linalg.svd(unfold)[0][:, :rank]
 
     unfold = tl.unfold(tensor, 1)
     unfold = unfold[:, np.all(np.isfinite(unfold), axis=0)]
@@ -205,7 +206,7 @@ def initialize_cp(tensor: np.ndarray, rank: int):
         if tensor.shape[mode] >= rank:
             unfold = tl.unfold(tensor, mode)
             if contain_missing:
-                si = SoftImpute(max_rank=rank)
+                si = IterativeSVD(rank)
                 unfold = si.fit_transform(unfold)
 
             factors[mode] = partial_svd(unfold, rank, flip=True)[0]
@@ -250,7 +251,7 @@ def perform_CP(tOrig, r=6, tol=1e-6):
     return tFac
 
 
-def perform_CMTF(tOrig, mOrig, r=9, tol=1e-6, maxiter=50, qr=False):
+def perform_CMTF(tOrig, mOrig, r=9, tol=1e-6, maxiter=50, progress=True):
     """ Perform CMTF decomposition. """
     assert tOrig.dtype == float
     assert mOrig.dtype == float
@@ -265,9 +266,10 @@ def perform_CMTF(tOrig, mOrig, r=9, tol=1e-6, maxiter=50, qr=False):
     # Precalculate the missingness patterns
     uniqueInfo = np.unique(np.isfinite(unfolded.T), axis=1, return_inverse=True)
 
-    for _ in range(maxiter):
+    tq = tqdm(range(maxiter), disable=(not progress))
+    for _ in tq:
         tensor = np.nan_to_num(tOrig) + tl.cp_to_tensor(tFac) * np.isnan(tOrig)
-        tFac = parafac(tensor, r, 200, init=tFac, verbose=False, fixed_modes=[0], mask=np.isfinite(tOrig), linesearch=True, tol=1e-9)
+        tFac = parafac(tensor, r, 2000, init=tFac, verbose=False, fixed_modes=[0], mask=np.isfinite(tOrig), linesearch=True, tol=1e-9)
 
         # Solve for the glycan matrix fit
         tFac.mFactor = np.linalg.lstsq(tFac.factors[0][missingM, :], mOrig[missingM, :], rcond=-1)[0].T
@@ -277,11 +279,9 @@ def perform_CMTF(tOrig, mOrig, r=9, tol=1e-6, maxiter=50, qr=False):
         kr = np.vstack((kr, tFac.mFactor))
         tFac.factors[0] = censored_lstsq(kr, unfolded.T, uniqueInfo)
 
-        if qr:
-            tFac.factors[0] = np.linalg.qr(tFac.factors[0])[0]
-
         R2X_last = R2X
         R2X = calcR2X(tFac, tOrig, mOrig)
+        tq.set_postfix(R2X=R2X, delta=R2X - R2X_last, refresh=False)
         assert R2X > 0.0
 
         if R2X - R2X_last < tol:
@@ -293,6 +293,4 @@ def perform_CMTF(tOrig, mOrig, r=9, tol=1e-6, maxiter=50, qr=False):
     tFac = sort_factors(tFac)
     tFac.R2X = R2X
 
-    print("R2X: " + str(tFac.R2X))
-
-    return tFac
+    return tFac

tensorpack/cmtf.py

@@ -1,6 +1,5 @@
 import pickle
 import numpy as np
-import pandas as pd
 from numpy.linalg import norm
 from tensorly import partial_svd
 from .cmtf import perform_CP, calcR2X
@@ -19,7 +18,7 @@ def impute_missing_mat(dat):
 
     diff = 1.0
     while diff > 1e-3:
-        U, S, V = partial_svd(imp, min(dat.shape)-1)
+        U, S, V = partial_svd(imp, min(dat.shape) - 1)
         scores = U @ np.diag(S)
         loadings = V
         recon = scores @ loadings
@@ -52,12 +51,11 @@ def perform_PCA(self, flattenon=0):
         scores = U @ np.diag(S)
         loadings = V
         recon = [scores[:, :rr] @ loadings[:rr, :] for rr in self.rrs]
-        self.PCAR2X = [calcR2X(c, mIn = flatData) for c in recon]
+        self.PCAR2X = [calcR2X(c, mIn=flatData) for c in recon]
         self.sizePCA = [sum(flatData.shape) * rr for rr in self.rrs]
 
-
     def Q2X_chord(self, drop=10, repeat=10):
-        self.chordQ2X = None # df
+        self.chordQ2X = None  # df
         pass
 
     def Q2X_entry(self, drop=10, repeat=10):
@@ -73,4 +71,4 @@ def load(self, pfile):
             tmp_dict = pickle.load(input_file)
             self.__dict__.update(tmp_dict)
 
-    pass
+    pass