fixed iter mult bug

Author: alexrd

csnanalysis/csn.py

@@ -329,3 +329,13 @@ def calc_committors(self,basins,labels=None,basin_labels=None,add_basins=False,t
 
         return full_comm
 
+    def idxs_to_trim(self,idxs):
+        """
+        Converts a list of idxs to trim_idxs.
+
+        idxs -- List of states in the transition matrix. Elements should be
+                integers from 0 to nstates.
+        """
+
+        return [self.trim_indices.index(i) for i in idxs if i in self.trim_indices]
+        

csnanalysis/matrix.py

@@ -1,6 +1,7 @@
 
 import scipy
 import numpy as np
+from itertools import compress
 
 def count_to_trans(countmat):
     """
@@ -46,17 +47,20 @@ def _make_sink(transmat,sink_states):
                 sink_mat.data[i] = 1.
                 set_to_one[sink_states.index(sink_mat.col[i])] = True
 
-    # set diagonal elements to 1
-    for i in range(len(sink_states)):
-        if not set_to_one[i]:
-            # add element sink_mat[sink_states[i]][sink_states[i]] = 1
-            np.append(sink_mat.row,sink_states[i])
-            np.append(sink_mat.col,sink_states[i])
-            np.append(sink_mat.data,1.)
-
+    # set diagonal elements to 1 that haven't been set to one already
+    statelist = list(compress(sink_states, np.logical_not(set_to_one)))
+    sink_mat.row = np.append(sink_mat.row,statelist)
+    sink_mat.col = np.append(sink_mat.col,statelist)
+    sink_mat.data = np.append(sink_mat.data,[1 for i in statelist])
+    
     # remove zeros
     sink_mat.eliminate_zeros()
 
+    # check if sink_mat columns still sum to 1
+    minval = sink_mat.toarray().sum(axis=0).min()
+    if minval < 0.99999:
+        raise ValueError("Error! Columns no longer sum to one in _make_sink!")
+
     return sink_mat
 
 def eig_weights(transmat):
@@ -107,14 +111,14 @@ def _trans_mult_iter(transmat,tol,maxstep=20):
 
     var = 1
     step = 0
-    while var > tol or step > maxstep:
+    while (var > tol) and (step < maxstep):
         newmat = np.matmul(t,t)
         var = np.abs(newmat-t).max()
         t = newmat.copy()
         step += 1
 
-    if step > maxstep:
-        print("Warning: iterative multiplication not converged after",maxstep,"steps: (var = ",var)
+    if step == maxstep and var > tol:
+        print("Warning: iterative multiplication not converged after",step,"steps: (var = ",var,"), (tol = ",tol,")")
 
     return t
 
@@ -198,7 +202,7 @@ def _extend(transmat,hubstates):
     ext_mat = scipy.sparse.coo_matrix((data, (rows, cols)), shape=(2*n,2*n))
     return ext_mat
 
-def _getring(transmat,basin,eig_weights,tol,maxstep):
+def _getring(transmat,basin,wts,tol,maxstep):
     """
     Given a transition matrix, and a set of states that form a basin,
     this returns a vector describing how probability exits that basin.
@@ -215,11 +219,11 @@ def _getring(transmat,basin,eig_weights,tol,maxstep):
     for b in basin:
         for i in range(n):
             if i not in basin:
-                ringprob[i] += sink_results[i][b]
-    
-    return ringprob
+                ringprob[i] += wts[b]*sink_results[i][b]
+
+    return ringprob/wts[basin].sum()
 
-def hubscores(transmat,hubstates,basins,tol=1e-6,maxstep=20,wts=None):
+def hubscores(transmat,hubstates,basins,tol=1e-6,maxstep=30,wts=None):
     """
     This function computes hub scores, which are the probabilities that
     transitions between a set of communities will use a given community as
@@ -261,8 +265,9 @@ def hubscores(transmat,hubstates,basins,tol=1e-6,maxstep=20,wts=None):
     if wts is None:
         wts = eig_weights(transmat)
 
+
     h = np.zeros((2,2),dtype=float)
-    ring = [_getring(transmat,b,eig_weights,tol,maxstep) for b in basins]
+    ring = [_getring(transmat,b,wts,tol,maxstep) for b in basins]
 
     for source,sink in [[0,1],[1,0]]:
         for i,p in enumerate(ring[source]):