hdiff output

r27403/compute_rates.py 2016-03-16 18:34:01.739350585 +0000 r27402/compute_rates.py 2016-03-16 18:34:02.155354856 +0000
177:     177:     
178:     return generator.rate_constants, generator.Peq, generator.rate_norm178:     return generator.rate_constants, generator.Peq, generator.rate_norm
179: 179: 
180: 180: 
181: description="""Compute rates from a pathsample database.  The transition181: description="""Compute rates from a pathsample database.  The transition
182: states and minima data are read from min.data and ts.data.  The product182: states and minima data are read from min.data and ts.data.  The product
183: and reactant states are read from min.A and min.B.  The rates are computed183: and reactant states are read from min.A and min.B.  The rates are computed
184: by solving a set of linear equations using a sparse solver.184: by solving a set of linear equations using a sparse solver.
185: """185: """
186: 186: 
187: def run(directory, T, A, B, out_prefix, tstart):187: def main():
188:     rate_constants, Peq, knorm = make_rates(directory, T)188:     """the main loop"""
 189:     tstart =  time.clock()
 190:     parser = argparse.ArgumentParser(description=description)
 191: 
 192:     parser.add_argument("-d", type=str, default=".",
 193:                         help="directory with the min.data, ts.data files")
 194:     parser.add_argument("-T", type=float, default=1., help="Temperature")
 195:     args = parser.parse_args()
 196:     print "temperature", args.T
 197:     directory = args.d
 198:     print "reading from directory:", os.path.abspath(directory)
 199:     
 200:     A = read_minA(directory+"/min.A")
 201:     B = read_minA(directory+"/min.B")
 202:     rate_constants, Peq, knorm = make_rates(args.d, args.T)
189:             203:             
190: 204: 
191:     if True:205:     if True:
192:         fname = "{}.rate_consts".format(out_prefix)206:         fname = "out.rate_consts"
193:         print "saving rate constants to:", fname207:         print "saving rate constants to", fname
194:         with open(fname, "w") as fout:208:         with open(fname, "w") as fout:
195:             for (u,v), k in sorted(rate_constants.iteritems()):209:             for (u,v), k in sorted(rate_constants.iteritems()):
196:                 fout.write("%6d %6d %s\n" % (u,v,k/knorm))210:                 fout.write("%6d %6d %s\n" % (u,v,k/knorm))
197: 211: 
198:     print "checking and reducing the graph structure"212:     print "checking and reducing the graph structure"
199:     try:213:     try:
200:         rate_constants = reduce_rates(rate_constants, B, A=A)214:         rate_constants = reduce_rates(rate_constants, B, A=A)
201:     except Exception:215:     except Exception:
202:         analyze_graph_error(rate_constants, A, B)216:         analyze_graph_error(rate_constants, A, B)
203:         raise217:         raise
204:     218:     
205:     print "computing mean first passage times"219:     print "computing mean first passage times"
206:     calculator = TwoStateRates(rate_constants, A, B, weights=Peq, check_rates=False)220:     calculator = TwoStateRates(rate_constants, A, B, weights=Peq, check_rates=False)
207:     calculator.compute_rates(use_umfpack=True)221:     calculator.compute_rates(use_umfpack=True)
208:     print "computing rates"222:     print "computing rates"
209:     kAB = calculator.get_rate_AB()223:     kAB = calculator.get_rate_AB()
210:     print "k(B<-A)", kAB / knorm224:     print "k(B<-A)", kAB / knorm
211:     225:     
212:     if True:226:     if True:
213:         fname = "{}.mfpt".format(out_prefix)227:         fname = "out.mfpt"
214:         print "saving mean first passage times for all minima to reach B to file", fname228:         print "saving mean first passage times for all minima to reach B to file", fname
215:         mfpt = sorted([(m, t) for m, t in 229:         mfpt = sorted([(m, t) for m, t in 
216:                        calculator.mfpt_computer.mfpt_dict.iteritems()])230:                        calculator.mfpt_computer.mfpt_dict.iteritems()])
217:         with open(fname, "w") as fout:231:         with open(fname, "w") as fout:
218:             for i, t in mfpt:232:             for i, t in mfpt:
219:                 fout.write("%d %.12g\n" % (i, t * knorm))233:                 fout.write("%d %.12g\n" % (i, t * knorm))
220: 234: 
221:     print "computing committor probabilities"235:     print "computing committor probabilities"
222:     sys.stdout.flush()236:     sys.stdout.flush()
223:     calculator.compute_committors()237:     calculator.compute_committors()
224:     print "computing steady state rates"238:     print "computing steady state rates"
225:     kSS = calculator.get_rate_AB_SS() / knorm239:     kSS = calculator.get_rate_AB_SS() / knorm
226:     print "kSS(B<-A)", kSS240:     print "kSS(B<-A)", kSS
227:     241:     
228:     if True:242:     if True:
229:         # get the alternate definition of committors for the nodes in A and B243:         fname = "out.committors"
230:         Acomm = calculator.get_alternate_committors(A, B) 
231:         Bcomm = calculator.get_alternate_committors(B, A) 
232:         all_committors = calculator.committor_computer.committor_dict.copy() 
233:         all_committors.update(Acomm) 
234:         all_committors.update(Bcomm) 
235:          
236:         fname = "{}.committors".format(out_prefix) 
237:         print "saving committor probabilities for all minima to file", fname244:         print "saving committor probabilities for all minima to file", fname
238:         coms = sorted([(m, c) for m, c in all_committors.iteritems()])245:         coms = sorted([(m, c) for m, c in 
 246:                        calculator.committor_computer.committor_dict.iteritems()])
239:         with open(fname, "w") as fout:247:         with open(fname, "w") as fout:
240:             for i, c in coms:248:             for i, c in coms:
241:                 fout.write("%d %.12g\n" % (i, c))249:                 fout.write("%d %.12g\n" % (i, c))
242:     250:     
243:     time_solve = calculator.mfpt_computer.time_solve + calculator.committor_computer.time_solve251:     time_solve = calculator.mfpt_computer.time_solve + calculator.committor_computer.time_solve
244:     print "time spent solving linear equations", time_solve, "seconds"252:     print "time spent solving linear equations", time_solve, "seconds"
245:     print "total time", time.clock() - tstart253:     print "total time", time.clock() - tstart
246:  
247:  
248: def main(): 
249:     """the main loop""" 
250:     tstart =  time.clock() 
251:     parser = argparse.ArgumentParser(description=description) 
252:  
253:     parser.add_argument("-d", type=str, default=".", 
254:                         help="directory with the min.data, ts.data files") 
255:     parser.add_argument("-T", type=float, default=1., help="Temperature") 
256:     parser.add_argument("--prefix", type=str, default="rates.out", help="output prefix") 
257:     args = parser.parse_args() 
258:      
259:     print "temperature", args.T 
260:     directory = args.d 
261:     print "reading from directory:", os.path.abspath(directory) 
262:     out_prefix = args.prefix 
263:      
264:     A = read_minA(directory+"/min.A") 
265:     B = read_minA(directory+"/min.B") 
266:      
267:     print "\ncomputing rates from A to B" 
268:     run(directory, args.T, A, B, out_prefix + "AtoB", tstart) 
269:     print "\ncomputing rates from B to A" 
270:     run(directory, args.T, B, A, out_prefix + "BtoA", tstart) 
271:      
272:      
273:      254:      
274: 255: 
275: 256: 
276: if __name__ == "__main__":257: if __name__ == "__main__":
277:     main()258:     main()


r27403/_compute_rates_utils.py 2016-03-16 18:34:01.547348608 +0000 r27402/_compute_rates_utils.py 2016-03-16 18:34:01.959352849 +0000
219:     def get_rate_AB(self):219:     def get_rate_AB(self):
220:         """return the rate from A to B220:         """return the rate from A to B
221:         221:         
222:         the rate is the inverse mean first passage time averaged over the nodes in A222:         the rate is the inverse mean first passage time averaged over the nodes in A
223:         """223:         """
224:         rate = sum((self.weights[a] / self.mfptimes[a] for a in self.A))224:         rate = sum((self.weights[a] / self.mfptimes[a] for a in self.A))
225:         norm = sum((self.weights[a] for a in self.A))225:         norm = sum((self.weights[a] for a in self.A))
226:         226:         
227:         return rate / norm227:         return rate / norm
228:     228:     
229:     def get_alternate_committors(self, Agroup, Bgroup):229:     def get_rate_AB_SS(self):
230:         """for each node a in A, compute the probability that it ends up in 230:         """
231:         B before coming back to itself or reaching another node in A.231:         return the steady state rate from A to B
232:         """232:         """
233:         a_committors = dict([(a, 0.) for a in Agroup])233:         # for each node a in A, compute the probability that it ends up in
 234:         # B before coming back to itself or reaching another node in A.
 235:         a_committors = dict([(a, 0.) for a in self.A])
234:         for uv, rate in self.rate_constants.iteritems():236:         for uv, rate in self.rate_constants.iteritems():
235:             u, v = uv237:             u, v = uv
236:             if u in Agroup and v not in Agroup:238:             if u in self.A and v not in self.A:
237:                 if v in Bgroup:239:                 if v in self.B:
238:                     vcomm = 1.240:                     vcomm = 1.
239:                 else:241:                 else:
240:                     vcomm = self.committor_dict[v]242:                     vcomm = self.committor_dict[v]
241:                 a_committors[u] += rate * vcomm243:                 a_committors[u] += rate * vcomm
242:         for a in Agroup:244:         for a, c in a_committors.iteritems():
243:             a_committors[a] /= self.sum_out_rates[a]245:             a_committors[a] /= self.sum_out_rates[a]
244:         return a_committors246:         # the sum_out_rates cancels here, we can remove it
245:  
246:     def get_rate_AB_SS(self): 
247:         """ 
248:         return the steady state rate from A to B 
249:         """ 
250:         # for each node a in A, compute the probability that it ends up in 
251:         # B before coming back to itself or reaching another node in A. 
252:         a_committors = self.get_alternate_committors(self.A, self.B) 
253:         247:         
254:         rate = sum((self.weights[a] * a_committors[a] * self.sum_out_rates[a] for a in self.A))248:         rate = sum((self.weights[a] * a_committors[a] * self.sum_out_rates[a] for a in self.A))
255:         norm = sum((self.weights[a] for a in self.A))249:         norm = sum((self.weights[a] for a in self.A))
256:         250:         
257:         return rate / norm251:         return rate / norm
258: 252: 
259:     def get_committor(self, x):253:     def get_committor(self, x):
260:         """return the probability that a trajectory starting from x reaches B before A"""254:         """return the probability that a trajectory starting from x reaches B before A"""
261:         if x in self.A:255:         if x in self.A:
262:             return 0.256:             return 0.


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0