hdiff output

r27480/_compute_rates_utils.py 2016-07-06 15:37:25.376709463 +0100 r27479/_compute_rates_utils.py 2016-07-06 15:37:25.744714386 +0100
  1: import time  1: import time
  2: from itertools import izip  2: from itertools import izip
  3: from collections import defaultdict  3: from collections import defaultdict
  4:   4: 
  5: import numpy as np  5: import numpy as np
  6: import scipy.sparse  6: import scipy.sparse
  7: import scipy.sparse.linalg  7: import scipy.sparse.linalg
  8: import networkx as nx  8: import networkx as nx
  9:   9: 
 10: class LinalgError(Exception): 
 11:     """raised when a the linear algebra solver fails""" 
 12:  10: 
 13: def reduce_rates(rates, B, A=None): 11: def reduce_rates(rates, B, A=None):
 14:     B = set(B) 12:     B = set(B)
 15:     if A is not None: 13:     if A is not None:
 16:         A = set(A) 14:         A = set(A)
 17:         if A.intersection(B): 15:         if A.intersection(B):
 18:             raise Exception("A and B share", len(A.intersection(B)), "nodes") 16:             raise Exception("A and B share", len(A.intersection(B)), "nodes")
 19:     graph = nx.Graph() 17:     graph = nx.Graph()
 20:     graph.add_edges_from(rates.iterkeys()) 18:     graph.add_edges_from(rates.iterkeys())
 21:      19:     
 67:         self.B = set(B) 65:         self.B = set(B)
 68:         self.nodes = set() 66:         self.nodes = set()
 69:         for u, v in self.rates.iterkeys(): 67:         for u, v in self.rates.iterkeys():
 70:             self.nodes.add(u) 68:             self.nodes.add(u)
 71:             self.nodes.add(v) 69:             self.nodes.add(v)
 72:          70:         
 73:         self.time_solve = 0. 71:         self.time_solve = 0.
 74:          72:         
 75:     def make_matrix(self): 73:     def make_matrix(self):
 76:         intermediates = self.nodes - self.A - self.B 74:         intermediates = self.nodes - self.A - self.B
 77:         if len(intermediates) == 0: 
 78:             self.matrix = None 
 79:             return 
 80:          75:         
 81:         node_list = list(intermediates) 76:         node_list = list(intermediates)
 82:         n = len(node_list) 77:         n = len(node_list)
 83:         matrix = scipy.sparse.dok_matrix((n,n)) 78:         matrix = scipy.sparse.dok_matrix((n,n))
 84:         right_side = np.zeros(n) 79:         right_side = np.zeros(n)
 85:         node2i = dict([(node,i) for i, node in enumerate(node_list)]) 80:         node2i = dict([(node,i) for i, node in enumerate(node_list)])
 86:  81:         
 87:          82:         
 88:         for uv, rate in self.rates.iteritems(): 83:         for uv, rate in self.rates.iteritems():
 89:             u, v = uv 84:             u, v = uv
 90: #            v, u = uv 85: #            v, u = uv
 91:  86: 
 92:             if u in intermediates: 87:             if u in intermediates:
 93:                 iu = node2i[u] 88:                 iu = node2i[u]
 94:                 matrix[iu,iu] -= rate 89:                 matrix[iu,iu] -= rate
 95:  90: 
 96:                 if v in intermediates: 91:                 if v in intermediates:
 98:          93:         
 99:                 if v in self.B: 94:                 if v in self.B:
100:                     right_side[iu] -= rate 95:                     right_side[iu] -= rate
101:          96:         
102:         self.node_list = node_list 97:         self.node_list = node_list
103:         self.node2i = node2i 98:         self.node2i = node2i
104:         self.matrix =  matrix.tocsr() 99:         self.matrix =  matrix.tocsr()
105:         self.right_side = right_side100:         self.right_side = right_side
106:         101:         
107:     def compute_committors(self):102:     def compute_committors(self):
108:         self.committor_dict = dict(((a, 0) for a in self.A)) 
109:         self.committor_dict.update(dict(((b, 1) for b in self.B))) 
110:         self.make_matrix()103:         self.make_matrix()
111:         if self.matrix is None: 
112:             return self.committor_dict 
113:              
114:         if self.right_side.size == 1:104:         if self.right_side.size == 1:
115:             # some versions of scipy can't handle matrices of size 1105:             # some versions of scipy can't handle matrices of size 1
116:             committors = np.array([self.right_side[0] / self.matrix[0,0]])106:             committors = np.array([self.right_side[0] / self.matrix[0,0]])
117:         else:107:         else:
118:             t0 = time.clock()108:             t0 = time.clock()
119:             committors = scipy.sparse.linalg.spsolve(self.matrix, self.right_side)109:             committors = scipy.sparse.linalg.spsolve(self.matrix, self.right_side)
120:             self.time_solve += time.clock() - t0110:             self.time_solve += time.clock() - t0
121:         eps = 1e-10111:         self.committor_dict = dict(((node, c) for node, c in izip(self.node_list, committors)))
122:         if np.any(committors < -eps) or np.any(committors > 1+eps): 
123:             qmax = committors.max() 
124:             qmin = committors.min() 
125:             raise LinalgError("The committors are not all between 0 and 1.  max=%.18g, min=%.18g" % (qmax, qmin)) 
126:         self.committor_dict.update(dict(((node, c) for node, c in izip(self.node_list, committors)))) 
127: #        self.committors = committors112: #        self.committors = committors
128: #        print "committors", committors113: #        print "committors", committors
129:         return self.committor_dict114:         return self.committor_dict
130:     115:     
131: 116: 
132: class MfptLinalgSparse(object):117: class MfptLinalgSparse(object):
133:     """compute mean first passage times using sparse linear algebra"""118:     """compute mean first passage times using sparse linear algebra"""
134:     def __init__(self, rates, B, sum_out_rates=None, check_graph=True):119:     def __init__(self, rates, B, sum_out_rates=None, check_graph=True):
135:         self.rates = rates120:         self.rates = rates
136:         self.B = set(B)121:         self.B = set(B)
190:     175:     
191:     def compute_mfpt(self, use_umfpack=True):176:     def compute_mfpt(self, use_umfpack=True):
192:         if not hasattr(self, "matrix"):177:         if not hasattr(self, "matrix"):
193:             self.make_matrix(self.nodes - self.B)178:             self.make_matrix(self.nodes - self.B)
194:         t0 = time.clock()179:         t0 = time.clock()
195:         times = scipy.sparse.linalg.spsolve(self.matrix, -np.ones(self.matrix.shape[0]),180:         times = scipy.sparse.linalg.spsolve(self.matrix, -np.ones(self.matrix.shape[0]),
196:                                             use_umfpack=use_umfpack)181:                                             use_umfpack=use_umfpack)
197:         self.time_solve += time.clock() - t0182:         self.time_solve += time.clock() - t0
198:         self.mfpt_dict = dict(((node, time) for node, time in izip(self.node_list, times)))183:         self.mfpt_dict = dict(((node, time) for node, time in izip(self.node_list, times)))
199:         if np.any(times < 0):184:         if np.any(times < 0):
200:             raise LinalgError("error the mean first passage times are not all greater than zero")185:             raise RuntimeError("error the mean first passage times are not all greater than zero")
201:         return self.mfpt_dict186:         return self.mfpt_dict
202: 187: 
203:     def compute_mfpt_subgroups(self, use_umfpack=True):188:     def compute_mfpt_subgroups(self, use_umfpack=True):
204:         for group in self.subgroups:189:         for group in self.subgroups:
205:             self.make_matrix(group)190:             self.make_matrix(group)
206:             t0 = time.clock()191:             t0 = time.clock()
207:             times = scipy.sparse.linalg.spsolve(self.matrix, -np.ones(self.matrix.shape[0]),192:             times = scipy.sparse.linalg.spsolve(self.matrix, -np.ones(self.matrix.shape[0]),
208:                                                 use_umfpack=use_umfpack)193:                                                 use_umfpack=use_umfpack)
209:             self.time_solve += time.clock() - t0194:             self.time_solve += time.clock() - t0
210:             for node, time in izip(self.node_list, times):195:             for node, time in izip(self.node_list, times):
284:         """compute the mean first passage times."""269:         """compute the mean first passage times."""
285:         if subgroups:270:         if subgroups:
286:             self.mfptimes = self.mfpt_computer.compute_mfpt_subgroups(use_umfpack=use_umfpack)271:             self.mfptimes = self.mfpt_computer.compute_mfpt_subgroups(use_umfpack=use_umfpack)
287:         else:272:         else:
288:             self.mfptimes = self.mfpt_computer.compute_mfpt(use_umfpack=use_umfpack)273:             self.mfptimes = self.mfpt_computer.compute_mfpt(use_umfpack=use_umfpack)
289:     274:     
290:     def compute_committors(self):275:     def compute_committors(self):
291:         """compute the committors""" 276:         """compute the committors""" 
292:         self.committor_computer = CommittorLinalg(self.rate_constants, self.A, self.B)277:         self.committor_computer = CommittorLinalg(self.rate_constants, self.A, self.B)
293:         self.committor_dict = self.committor_computer.compute_committors()278:         self.committor_dict = self.committor_computer.compute_committors()
294:         return self.committor_dict279:         
295:  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0