Merge pull request #52 from NelleV/fix_empty_profile

FIX LOIC now works fine when counts_profile contains 0

iced/normalization/__init__.py

@@ -70,6 +70,17 @@ def ICE_normalization(X, SS=None, max_iter=3000, eps=1e-4, copy=True,
     if verbose:
         print("Matrix is triangular superior")
 
+    if counts_profile is not None:
+        rows_to_remove = counts_profile == 0
+        if sparse.issparse(X):
+            rows_to_remove = np.where(rows_to_remove)[0]
+            X.data[np.isin(X.row, rows_to_remove)] = 0
+            X.data[np.isin(X.col, rows_to_remove)] = 0
+            X = X.eliminate_zeros()
+        else:
+            X[rows_to_remove] = 0
+            X[:, rows_to_remove] = 0
+
     if total_counts is None:
         total_counts = X.sum()
     for it in np.arange(max_iter):
@@ -93,6 +104,7 @@ def ICE_normalization(X, SS=None, max_iter=3000, eps=1e-4, copy=True,
         dbias = sum_ds.reshape((m, 1))
         if counts_profile is not None:
             dbias /= counts_profile[:, np.newaxis]
+            dbias[counts_profile == 0] = 0
         # To avoid numerical instabilities
         dbias /= dbias[dbias != 0].mean()
 

iced/normalization/tests/test_normalization.py

@@ -23,6 +23,22 @@ def test_ICE_normalization():
     assert_array_almost_equal(normed_X, X / (bias.T * bias), 6)
 
 
+def test_ICE_normalization_cancer():
+    n = 100
+    random_state = np.random.RandomState(seed=42)
+    X = random_state.randint(0, 100, size=(n, n))
+    X = X + X.T
+    profile = np.ones(n)
+    profile[:10] = 0
+    profile[50:] = 2
+    normed_X = ICE_normalization(X, eps=1e-10, counts_profile=profile)
+    assert not np.all(np.isnan(normed_X))
+    normed_X[np.isnan(normed_X)] = 0
+    inferred_profile = normed_X.sum(axis=0)
+    inferred_profile /= inferred_profile.max()
+    assert_array_almost_equal(inferred_profile, profile / profile.max())
+
+
 def test_sparse_ICE_normalization():
     n = 100
     random_state = np.random.RandomState(seed=42)