hdiff output

r32034/_compute_rates_utils.py 2017-03-07 17:30:18.909423707 +0000 r32033/_compute_rates_utils.py 2017-03-07 17:30:20.169440381 +0000
109:         self.committor_dict.update(dict(((b, 1) for b in self.B)))109:         self.committor_dict.update(dict(((b, 1) for b in self.B)))
110:         self.make_matrix()110:         self.make_matrix()
111:         if self.matrix is None:111:         if self.matrix is None:
112:             return self.committor_dict112:             return self.committor_dict
113:             113:             
114:         if self.right_side.size == 1:114:         if self.right_side.size == 1:
115:             # some versions of scipy can't handle matrices of size 1115:             # some versions of scipy can't handle matrices of size 1
116:             committors = np.array([self.right_side[0] / self.matrix[0,0]])116:             committors = np.array([self.right_side[0] / self.matrix[0,0]])
117:         else:117:         else:
118:             t0 = time.clock()118:             t0 = time.clock()
119:             #committors = scipy.sparse.linalg.spsolve(self.matrix, self.right_side)119:             committors = scipy.sparse.linalg.spsolve(self.matrix, self.right_side)
120:             #committors, res = scipy.sparse.linalg.minres(self.matrix, self.right_side) 
121:             committors, res = scipy.sparse.linalg.lgmres(self.matrix, self.right_side) 
122:             self.time_solve += time.clock() - t0120:             self.time_solve += time.clock() - t0
123:         eps = 1e-5121:         eps = 1e-10
124:         if np.any(committors < -eps) or np.any(committors > 1+eps):122:         if np.any(committors < -eps) or np.any(committors > 1+eps):
125:             qmax = committors.max()123:             qmax = committors.max()
126:             qmin = committors.min()124:             qmin = committors.min()
127:             raise LinalgError("The committors are not all between 0 and 1.  max=%.18g, min=%.18g" % (qmax, qmin))125:             raise LinalgError("The committors are not all between 0 and 1.  max=%.18g, min=%.18g" % (qmax, qmin))
128:         self.committor_dict.update(dict(((node, c) for node, c in izip(self.node_list, committors))))126:         self.committor_dict.update(dict(((node, c) for node, c in izip(self.node_list, committors))))
129: #        self.committors = committors127: #        self.committors = committors
130: #        print "committors", committors128: #        print "committors", committors
131:         return self.committor_dict129:         return self.committor_dict
132:     130:     
133: 131: 
173:         node_list = list(intermediates)171:         node_list = list(intermediates)
174:         n = len(node_list)172:         n = len(node_list)
175:         matrix = scipy.sparse.dok_matrix((n,n))173:         matrix = scipy.sparse.dok_matrix((n,n))
176:         node2i = dict([(node,i) for i, node in enumerate(node_list)])174:         node2i = dict([(node,i) for i, node in enumerate(node_list)])
177:         175:         
178:         for iu, u in enumerate(node_list):176:         for iu, u in enumerate(node_list):
179:             matrix[iu,iu] = -self.sum_out_rates[u]177:             matrix[iu,iu] = -self.sum_out_rates[u]
180:         178:         
181:         for uv, rate in self.rates.iteritems():179:         for uv, rate in self.rates.iteritems():
182:             u, v = uv180:             u, v = uv
183:             if u in intermediates and v in intermediates:181:             if u in intermediates and v in intermediates: 
184:                 ui = node2i[u]182:                 ui = node2i[u]
185:                 vi = node2i[v]183:                 vi = node2i[v]
186:                 assert ui != vi184:                 assert ui != vi
187:                 matrix[ui,vi] = rate185:                 matrix[ui,vi] = rate
188:         186:         
189:         self.node_list = node_list187:         self.node_list = node_list
190:         self.node2i = node2i188:         self.node2i = node2i
191:         self.matrix =  matrix.tocsr()189:         self.matrix =  matrix.tocsr()
192:     190:     
193:     def compute_mfpt(self, use_umfpack=True):191:     def compute_mfpt(self, use_umfpack=True):
194:         if not hasattr(self, "matrix"):192:         if not hasattr(self, "matrix"):
195:             self.make_matrix(self.nodes - self.B)193:             self.make_matrix(self.nodes - self.B)
196:         t0 = time.clock()194:         t0 = time.clock()
197:         #times = scipy.sparse.linalg.spsolve(self.matrix, -np.ones(self.matrix.shape[0]),195:         times = scipy.sparse.linalg.spsolve(self.matrix, -np.ones(self.matrix.shape[0]),
198:         #                                    use_umfpack=use_umfpack)196:                                             use_umfpack=use_umfpack)
199:         #times = scipy.sparse.linalg.factorized(self.matrix)(-np.ones(self.matrix.shape[0])) 
200:         #times = scipy.sparse.linalg.bicgstab(self.matrix, -np.ones(self.matrix.shape[0])) 
201:         #times, res = scipy.sparse.linalg.minres(self.matrix, -np.ones(self.matrix.shape[0])) 
202:         times, res = scipy.sparse.linalg.lgmres(self.matrix, -np.ones(self.matrix.shape[0])) 
203:         #import pdb; pdb.set_trace() 
204:         self.time_solve += time.clock() - t0197:         self.time_solve += time.clock() - t0
205:         self.mfpt_dict = dict(((node, time) for node, time in izip(self.node_list, times)))198:         self.mfpt_dict = dict(((node, time) for node, time in izip(self.node_list, times)))
206:         if np.any(times < 0):199:         if np.any(times < 0):
207:             raise LinalgError("error the mean first passage times are not all greater than zero")200:             raise LinalgError("error the mean first passage times are not all greater than zero")
208:         return self.mfpt_dict201:         return self.mfpt_dict
209: 202: 
210:     def compute_mfpt_subgroups(self, use_umfpack=True):203:     def compute_mfpt_subgroups(self, use_umfpack=True):
211:         for group in self.subgroups:204:         for group in self.subgroups:
212:             self.make_matrix(group)205:             self.make_matrix(group)
213:             t0 = time.clock()206:             t0 = time.clock()
214:             #times = scipy.sparse.linalg.spsolve(self.matrix, -np.ones(self.matrix.shape[0]),207:             times = scipy.sparse.linalg.spsolve(self.matrix, -np.ones(self.matrix.shape[0]),
215:             #                                    use_umfpack=use_umfpack)208:                                                 use_umfpack=use_umfpack)
216:             #times, res = scipy.sparse.linalg.minres(self.matrix, -np.ones(self.matrix.shape[0])) 
217:             times, res = scipy.sparse.linalg.lgmres(self.matrix, -np.ones(self.matrix.shape[0])) 
218:             self.time_solve += time.clock() - t0209:             self.time_solve += time.clock() - t0
219:             for node, time in izip(self.node_list, times):210:             for node, time in izip(self.node_list, times):
220:                 self.mfpt_dict[node] = time211:                 self.mfpt_dict[node] = time
221: 212: 
222: 213: 
223: class TwoStateRates(object):214: class TwoStateRates(object):
224:     """compute committors and several different rates between two groups"""215:     """compute committors and several different rates between two groups"""
225:     def __init__(self, rate_constants, A, B, weights=None, check_rates=True):216:     def __init__(self, rate_constants, A, B, weights=None, check_rates=True):
226:         if check_rates:217:         if check_rates:
227:             self.rate_constants = reduce_rates(rate_constants, B, A=A)218:             self.rate_constants = reduce_rates(rate_constants, B, A=A)


r32034/FA2CG 2017-03-07 17:30:17.949411004 +0000 r32033/FA2CG 2017-03-07 17:30:19.157426990 +0000
  1: svn: warning: W195007: URL 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/OPEP/FA2CG' refers to a directory  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/OPEP/FA2CG' in revision 32033
  2: svn: E200009: Could not cat all targets because some targets are directories 
  3: svn: E200009: Illegal target for the requested operation 


r32034/FA2CG.py 2017-03-07 17:30:18.153413703 +0000 r32033/FA2CG.py 2017-03-07 17:30:19.409430324 +0000
  1: # !/usr/bin/env python  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/OPEP/FA2CG/FA2CG.py' in revision 32033
  2: # -*- coding: UTF8 -*- 
  3:  
  4: import sys 
  5:  
  6: import pdbparser 
  7: import RNA 
  8:  
  9: def filtbase(base, type): 
 10:         B = base.filter(RNA.heavy[type]) 
 11:         # mean along first coordinate 
 12:         B.coords[0] = B.coords.mean(0) 
 13:         B.atomnames[0] = " "+type+" " 
 14:         return B.pdbline(0) 
 15:  
 16: def Base_(FA_pdb, resid): 
 17: #        import pdb; pdb.set_trace() 
 18:         base = FA_pdb.filter_res([resid]) 
 19:         lines = [] 
 20:         for b in RNA.bases: 
 21:                 if b in base.resnames[0]: 
 22:                         if b in RNA.purines: 
 23:                                 lines.append(filtbase(base, b+"2")) 
 24:                         lines.append(filtbase(base, b+"1")) 
 25:                         break 
 26:         return lines 
 27:  
 28: def Base(FA_pdb): 
 29:         for i in xrange(Natoms): 
 30:                 if "O5'" in FA_pdb.atomnames[i]: 
 31:                         Base_(FA_pdb, i) 
 32:  
 33: def Coarse(FA_pdb): 
 34:         CG_backbone = FA_pdb.filter(RNA.FA2CG_backbone) 
 35:         for i in xrange(CG_backbone.Natoms): 
 36:                 CG_backbone.atomnames[i] = CG_backbone.atomnames[i].replace("'", "*") 
 37:         CG_full = CG_backbone 
 38:         for i in xrange(CG_backbone.Natoms): 
 39:                 if "C1*" in CG_backbone.atomnames[i]: 
 40:                         for l in Base_(FA_pdb, CG_backbone.resnum[i]): 
 41:                                 CG_full.insertatom(l, i+1) 
 42:          
 43:         return CG_full 
 44:  
 45: def main(filename, out): 
 46:         #--------------------- Main Routine ------------------------# 
 47:          
 48:         famodels = pdbparser.readpdb(filename) 
 49:         CG = Coarse(famodels[0]) 
 50:         CG.writetofile(out) 
 51:  
 52:  
 53: if __name__ == "__main__": 
 54:         if len(sys.argv) < 2: 
 55:                 fname = "6TNA_fit.pdb" 
 56:         else: 
 57:                 fname = sys.argv[1] 
 58:         if len(sys.argv) < 3: 
 59:                 outname = None 
 60:         else: 
 61:                 outname = sys.argv[2] 
 62:          
 63:         main(fname, outname) 


r32034/pdbparser.py 2017-03-07 17:30:18.657420374 +0000 r32033/pdbparser.py 2017-03-07 17:30:19.913436993 +0000
  1: #! /usr/bin/env python  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/OPEP/FA2CG/pdbparser.py' in revision 32033
  2:  
  3: import numpy 
  4: import sys 
  5:  
  6: __author__ = [ 
  7:                 "Tristan Cranolini <tristan.cranolini@gmail.com>", 
  8:                 "Jonathan Barnoud <jonathan@barnoud.net>", 
  9: ] 
 10:  
 11: class PDB(): 
 12:         def __init__(self, coords_=None, resize=2000): 
 13:                 self.atomnames = [] 
 14:                 self.resnames = [] 
 15:                 self.resnum = [] 
 16:                 self.chainname = "" 
 17:                  
 18:                 self.resize = resize 
 19:                 # coords_ allows us to resize the coordinates array less often 
 20:                 if coords_ is None: 
 21:                         self.coords_ = numpy.zeros((resize, 3)) 
 22:                         self.Natoms = 0 
 23:                 else: 
 24:                         self.coords_ = numpy.array(coords_, dtype=float).reshape(-1,3) 
 25:                         self.Natoms = self.coords_.shape[0] 
 26:                 self.Bfactors = [] 
 27:          
 28:         def get_coords(self): 
 29:                 return self.coords_[0:self.Natoms,] 
 30:          
 31:         coords = property(get_coords) 
 32:          
 33:         def addatom(self, line):                 
 34:                 self.atomnames.append(line[12:16]) 
 35:                 self.resnames.append(line[17:20]) 
 36:                 self.resnum.append(int(line[22:26])) 
 37:  
 38:                 if self.Natoms >= self.coords_.shape[0]: 
 39:                         self.coords_ = numpy.row_stack([self.coords_, numpy.zeros((self.resize,3))]) 
 40:  
 41:                 self.coords_[self.Natoms,:] = [line[30:38], line[38:46], line[46:54]] 
 42:                 self.Natoms += 1 
 43:                 try: 
 44:                         self.Bfactors.append(float(line[60:66])) 
 45:                 except (ValueError): 
 46:                         self.Bfactors.append(0.0) 
 47:          
 48:         def insertatom(self, line, idx=None): 
 49:                 if idx is None: 
 50:                         idx = self.Natoms 
 51:                  
 52:                 self.atomnames.insert(idx, line[12:16]) 
 53:                 self.resnames.insert(idx, line[17:20]) 
 54:                 self.resnum.insert(idx, int(line[22:26])) 
 55:  
 56:                 if self.Natoms >= self.coords_.shape[0]: 
 57:                         self.coords_ = numpy.row_stack([self.coords_, numpy.zeros((self.resize,3))]) 
 58:  
 59:                 self.coords_[idx+1:self.Natoms+1,:] = self.coords_[idx:self.Natoms,:] 
 60:                 self.coords_[idx,:] = [float(line[30:38]), float(line[38:46]), float(line[46:54])] 
 61:                 self.Natoms += 1 
 62:                 try: 
 63:                         self.Bfactors.insert(idx, float(line[60:66])) 
 64:                 except (ValueError): 
 65:                         self.Bfactors.insert(idx, 0.0) 
 66:  
 67:  
 68:         def center(self): 
 69:                 return sum(self.coords) / self.coords.shape[0] 
 70:  
 71:         def movetocenter(self): 
 72:                 self.coords -= self.center() 
 73:          
 74:         def pdbline(self, idx): 
 75:                 return "%-6s%5d %4s%1s%3s %1s%4d%1s                %8.3f%8.3f%8.3f%6.2f%6.2f                                         %2s%2s\n" \ 
 76:                         % ("ATOM", idx+1, self.atomnames[idx], " ", self.resnames[idx], 
 77:                         self.chainname, self.resnum[idx], " ", 
 78:                         self.coords[idx][0], self.coords[idx][1], self.coords[idx][2], 
 79:                         0, self.Bfactors[idx], " ", " ") 
 80:                          
 81:          
 82:         def writetofile(self, filename=None, filemode="w"): 
 83:                 if filename is None: 
 84:                         pdbfile = sys.stdout 
 85:                 else: 
 86:                         pdbfile = open(filename, filemode) 
 87:                  
 88:                 #pdbfile.write("MODEL\n") 
 89:                 for idx in range(len(self.atomnames)): 
 90:                         pdbfile.write(self.pdbline(idx)) 
 91:                  
 92:                 #pdbfile.write("TER\nENDMDL\n") 
 93:                  
 94:                 pdbfile.close() 
 95:  
 96:         def _make_copy(self, filter) : 
 97:                 """ 
 98:                 Create a filtered copy of the structure. 
 99:  
100:                 Parameters : 
101:                 ------------ 
102:                 - filter : an 1D array of boolean which will be used as mask. The lenght of 
103:                 this array have to fit the number of atoms of the current structure. 
104:  
105:                 Returns : 
106:                 --------- 
107:                 - copy : a filtered copy as a PDB object. 
108:                 """ 
109:                 structure = PDB(numpy.array(self.coords[filter, 0:3])) 
110:                 structure.atomnames = list(numpy.array(self.atomnames)[filter]) 
111:                 structure.resnames = list(numpy.array(self.resnames)[filter]) 
112:                 structure.resnum = list(numpy.array(self.resnum)[filter]) 
113:                 structure.chainname = self.chainname 
114:                  
115:                 structure.resize = self.resize 
116:                 #structure.Natoms = len(structure.atomnames) 
117:                 structure.Bfactors = self.Bfactors 
118:                 return structure 
119:  
120:         def filter(self, values, exclude=False, filtprop=lambda x: x.atomnames, filtcrit=lambda x, y: x in y) : 
121:                 """ 
122:                 Return a copy filtered by atom names. 
123:  
124:                 Parameters : 
125:                 ------------ 
126:                 - values : values used to do the filtering 
127:                 - filtcrit : criterion for the filter 
128:                 - exclude : if False, only the atoms with the specified atom names will be 
129:                 kept. If True they will be exclude and only the other atoms will appier in 
130:                 the copy. (False by default) 
131:                 -  
132:  
133:                 Returns : 
134:                 --------- 
135:                 - copy : a filtered copy of the structure. 
136:                 """ 
137:                 filter = [filtcrit(i, values) for i in filtprop(self)] 
138:                 filter = numpy.array(filter, dtype=bool) != exclude 
139:                 return self._make_copy(filter) 
140:  
141:         def section(self, begin, end, exclude=False) : 
142:                 """ 
143:                 Return a copy filtered with by resid between two values. 
144:  
145:                 Parameters : 
146:                 ------------ 
147:                 - begin, end : the first and last index of the selection. 
148:                 - exclude : if True the selection will be exclude from the copy, else only 
149:                 the selection will be kept in the copy. (False by delfault) 
150:  
151:                 Returns : 
152:                 --------- 
153:                 - copy : a filtered copy of the structure. 
154:                 """ 
155:                 filter_min = (numpy.array(self.resnum) >= begin) 
156:                 filter_max = (numpy.array(self.resnum) <= end) 
157:                 filter = (filter_min == filter_max) 
158:                 filter = (filter != exclude) 
159:                 return self._make_copy(filter) 
160:  
161:         def filter_res(self, residues, exclude=False) : 
162:                 """ 
163:                 Return a copy filtered by a resid list. 
164:  
165:                 Parameters : 
166:                 ------------ 
167:                 - residues : a list of resid to put in the selection 
168:                 - exclude : if True, the selection will be excluded from the copy. (False 
169:                 by default) 
170:  
171:                 Returns : 
172:                 --------- 
173:                 -copy : a filtered copy of the structure 
174:                 """ 
175:                 filter = numpy.zeros(len(self.resnum)) 
176:                 for index, resid in enumerate(self.resnum) : 
177:                         if resid in residues : 
178:                                 filter[index] = 1 
179:                 filter = numpy.array(filter, dtype=bool) 
180:                 filter = (filter != exclude) 
181:                 return self._make_copy(filter) 
182:  
183:         def check_integrity(self) : 
184:                 """ 
185:                 Check if there is no residue missing and if the residues are well sorted. 
186:  
187:                 Returns : 
188:                 --------- 
189:                 - answer : False if something is wrong. 
190:                 """ 
191:                 last_id = self.resnum[0] 
192:                 for id in self.resnum : 
193:                         if not last_id <= id <= (last_id + 1) : 
194:                                 return False 
195:                         last_id = id 
196:                 return True 
197:  
198: def readpdb(filename): 
199:         pdbfile = open(filename, "r") 
200:         models = [] 
201:         models.append(readmodel(pdbfile)) 
202:         for line in pdbfile: 
203:                 if line.startswith("MODEL") or line.startswith("ATOM"): 
204:                         #~ models.append(readmodel(pdbfile)) 
205:                         readmodel2(pdbfile, models) 
206:          
207:         return models 
208:  
209: def isatom(line): 
210:                 return (line.startswith("ATOM") or line.startswith("HETATM"))# and line[17:20] in PDB.standard_aa_names: 
211:  
212: def readmodel(pdbfile): 
213:         model = PDB() 
214:         for line in pdbfile: 
215:                 if isatom(line):# and line[17:20] in PDB.standard_aa_names: 
216:                         model.addatom(line) 
217:                 elif line.startswith("ENDMDL"): 
218:                         break 
219:          
220:         model.coords = numpy.array(model.coords) 
221:         return model 
222:  
223: def readmodel2(pdbfile, models): 
224:         model = PDB(numpy.empty_like(models[0].coords)) 
225:         i = 0 
226:         for line in pdbfile: 
227:                 if isatom(line): 
228:                         model.insertatom(line, i) 
229:                         i += 1 
230:                 elif line.startswith("ENDMDL"): 
231:                         break 
232:          
233:         models.append(model) 
234:  
235:  
236: if __name__ == "__main__": 
237:         models = readpdb(sys.argv[1]) 
238:         print models[0].coords.shape 
239:         print len(models[0].coords) 


r32034/RNA.py 2017-03-07 17:30:18.405417038 +0000 r32033/RNA.py 2017-03-07 17:30:19.661433659 +0000
  1: # !/usr/bin/env python  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/OPEP/FA2CG/RNA.py' in revision 32033
  2: # -*- coding: UTF8 -*- 
  3:  
  4: purines = ["A", "G"] 
  5: pyrimidines = ["C", "U"] 
  6: bases = purines+pyrimidines 
  7: lenbases_AA = {'A': 33, 'G': 34, 'C': 31, 'U': 30} 
  8:  
  9:  
 10: phosphate = ["P", "OP1", "OP2"] 
 11: sugar = ["C5'", "C4'", "C3'", "C2'", "C1'", "O4'", "O3'", "O2'"] 
 12:  
 13: heavy = dict() 
 14: heavy["A1"] = [" N7 ", " N9 ", " C4 ", " C5 ", " C8 "] 
 15: heavy["A2"] = [" N1 ", " N3 ", " N6 ", " C2 ", " C4 ", " C5 ", " C6 "] 
 16: heavy["A"]= heavy["A1"] + heavy["A2"] 
 17: heavy["C1"] = [" N1 ", " N3 ", " N4 ", " C2 ", " C4 ", " C5 ", " C6 ", " O2 "] 
 18: heavy["G1"] = heavy["A1"] 
 19: heavy["G2"] = [" N1 ", " N2 ", " N3 ", " C2 ", " C4 ", " C5 ", " C6 ", " O6"] 
 20: heavy["G"]= heavy["G1"] + heavy["G2"] 
 21: heavy["U1"] = [" N1 ", " N3 ", " C2 ", " C4 ", " C5 ", " C6 ", " O2 ", " O4 "] 
 22:  
 23:  
 24: FA2CG_backbone = [" P  ", " O5'", " C5'", " C4'", " C1'"] 
 25: FA_backbone = FA2CG_backbone + ["O1P", "O2P"] 
 26: CG_backbone = [" P  ", " O5* ", " C5*", " CA", " CY"] 
 27:  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0