IMP logo
IMP Reference Guide  2.20.0
The Integrative Modeling Platform
movers.py
1 import IMP
2 import IMP.core
3 import IMP.algebra
4 import itertools
5 import math
6 import random
7 import numpy as np
8 import scipy.spatial
9 import networkx as nx
10 from IMP.pmi.tools import get_restraint_set
11 
12 
13 class MoverScheme(object):
14  def __init__(self, length):
15  self.length = length
16  self.mover_scheme = None
17  self.score_scheme = {}
18  self.score_initial_weights = {}
19  self.position = 0
20  self.numbermc = 0
21 
22  def add_mover_scheme(self, mover_scheme):
23  if len(mover_scheme) != self.length:
24  raise ("Error")
25  else:
26  mc_scheme = []
27  for mvs in mover_scheme:
28  smv = IMP.core.SerialMover(mvs)
29  mc = IMP.core.MonteCarlo(mvs[0].get_model())
30  mc.set_scoring_function(get_restraint_set(mvs[0].get_model()))
31  mc.set_return_best(False)
32  mc.set_kt(1.0)
33  mc.add_mover(smv)
34  mc_scheme.append(mc)
35  mc.set_name(str(self.numbermc))
36  self.numbermc += 1
37  self.mover_scheme = mc_scheme
38 
39  def add_montecarlo_temperatures(self, temps):
40  self.temps = temps
41 
42  def add_montecarlo_steps(self, steps):
43  self.steps = steps
44 
45  def add_score_scheme(self, score, weigth_scheme):
46  if len(weigth_scheme) != self.length:
47  raise ("Error")
48  else:
49  self.score_scheme[score] = weigth_scheme
50  self.score_initial_weights[score] = score.get_weight()
51 
52  def get_length(self):
53  return self.length
54 
55  def get_position(self, position):
56  movers = self.mover_scheme[position]
57  score_tmp_dict = {}
58  for score in self.score_scheme:
59  score_tmp_dict[score] = self.score_scheme[score][position]
60  return movers, score_tmp_dict
61 
62  def get_and_update(self):
63  mc, score_tmp_dict = self.get_position(self.position)
64  temp = self.temps[self.position]
65  mc.set_kt(temp)
66  steps = self.steps[self.position]
67  self.position += 1
68  if self.position > self.length - 1: self.position = 0
69  return mc, score_tmp_dict, steps
70 
71 
72 class HybridMonteCarloMover(IMP.core.MonteCarloMover):
73  def __init__(self, m, hier, MoverScheme):
74  IMP.core.MonteCarloMover.__init__(
75  self,
76  m,
77  "HybridMonteCarloMover")
78 
79  self.mols = IMP.pmi.tools.get_molecules(hier)
80  self.refframes = {}
81  self.xyzs = {}
82  self.ps = []
83  self.m = m
84  self.rejects = []
85  self.lenrejects = 20
86 
87  self.rbs, self.beads = IMP.pmi.tools.get_rbs_and_beads([hier])
88  for rb in self.rbs:
89  self.ps.append(rb.get_particle())
90  for bead in self.beads:
91  self.ps.append(bead.get_particle())
92 
93  self.indexes = [p.get_index() for p in self.ps]
94 
95  self.do_store_current_coordinates()
96 
97  self.mover_scheme = MoverScheme
98 
99  def do_store_current_coordinates(self):
100  for rb in self.rbs:
101  self.refframes[rb] = rb.get_reference_frame()
102  for bead in self.beads:
103  self.xyzs[bead] = IMP.core.XYZ(bead).get_coordinates()
104 
105  def do_propose(self):
106 
107  self.rejects.append(0)
108  if len(self.rejects) > self.lenrejects:
109  self.rejects.pop(0)
110  print
111  self.rejects, float(sum(self.rejects)) / len(self.rejects)
112 
113  self.do_store_current_coordinates()
114 
115  rs = get_restraint_set(self.m)
116  # print "Initial",rs.evaluate(False)
117 
118  for i in range(self.mover_scheme.get_length()):
119  mc, score_tmp_dict, steps = self.mover_scheme.get_and_update()
120 
121  mc.reset_statistics()
122 
123  for score in score_tmp_dict:
124  score.set_weight(score_tmp_dict[score])
125 
126  mc.optimize(steps)
127 
128  print
129  mc.get_name(), float(mc.get_number_of_accepted_steps()) / mc.get_number_of_proposed_steps()
130  # print i,steps,rs.evaluate(False)
131 
132  # print "Final",rs.evaluate(False)
133 
134  return IMP.core.MonteCarloMoverResult(self.indexes, 1.0)
135 
136  def do_reject(self):
137  print
138  "REJECT"
139  self.rejects[-1] = 1
140  for rb in self.refframes:
141  rb.set_reference_frame(self.refframes[rb])
142  for bead in self.xyzs:
143  IMP.core.XYZ(bead).set_coordinates(self.xyzs[bead])
144 
145  def do_get_inputs(self):
146  return self.ps
147 
148 
149 class PCAMover(IMP.core.MonteCarloMover):
150  def __init__(self, m, hier, resolution=10, depth=4):
151 
152  IMP.core.MonteCarloMover.__init__(
153  self,
154  m,
155  "PCAMover")
156  self.depth = depth
157  self.propose = 0
158  self.reject = 0
159  self.rejects = []
160  self.lenrejects = 20
161  self.mols = IMP.pmi.tools.get_molecules(hier)
162  self.resolution = resolution
163  self.rbsd = {}
164  self.beadsd = {}
165  self.refframes = {}
166  self.xyzs = {}
167  self.ps = []
168  self.particle_blocks = []
169  for mol in self.mols:
170  rbs, beads = IMP.pmi.tools.get_rbs_and_beads([mol])
171  self.rbsd[mol] = rbs
172  self.beadsd[mol] = beads
173  for rb in self.rbsd[mol]:
174  self.ps.append(rb.get_particle())
175  for bead in self.beadsd[mol]:
176  self.ps.append(bead.get_particle())
177  ps = IMP.atom.Selection(mol, resolution=self.resolution).get_selected_particles()
178  self.particle_blocks.append(ps)
179 
180  self.do_store_current_coordinates()
181 
182  def do_store_current_coordinates(self):
183  for mol in self.mols:
184  for rb in self.rbsd[mol]:
185  self.refframes[rb] = rb.get_reference_frame()
186  for bead in self.beadsd[mol]:
187  self.xyzs[bead] = IMP.core.XYZ(bead).get_coordinates()
188 
189  def get_full_graph(self):
190  '''
191  get the full graph of distances between every particle pair
192  '''
193 
194  pdist_array = np.array(
196  pdist_mat = scipy.spatial.distance.squareform(pdist_array)
197  pdist_mat[pdist_mat <= 10] = 1
198  pdist_mat[pdist_mat > 10] = 0
199  graph = nx.Graph(pdist_mat)
200  return graph
201 
202  def get_list_of_interacting_molecules(self, succ_dict, pointer):
203  if pointer in succ_dict:
204  paths = [path for p in succ_dict[pointer] for path in self.get_list_of_interacting_molecules(succ_dict, p)]
205  ret = [[pointer] + path for path in paths]
206  # print pointer,paths,ret
207  return ret
208  else:
209  return [[pointer]]
210 
211  def get_minimum_spanning_tree(self):
212  """
213  return the minimum spanning tree
214  """
215  graph = self.get_full_graph()
216  graph = nx.minimum_spanning_tree(graph)
217  return graph
218 
219  def get_sublists(self, lst, size):
220  result = [sublist for sublist in
221  (lst[x:x + size] for x in range(len(lst) - size + 1))]
222  return result
223 
224  def get_list_of_connected_mols(self, depth):
225  gr = self.get_minimum_spanning_tree()
226  all_groups = set()
227  if len(gr.edges()) > 0:
228 
229  for node in gr.nodes():
230  succ_dict = nx.dfs_successors(gr, node)
231  paths = [path for p in succ_dict[node] for path in self.get_list_of_interacting_molecules(succ_dict, p)]
232  ret = [[node] + path for path in paths]
233  for path in ret:
234  for s in range(2, min([depth + 1, len(path) + 1, len(self.mols)])):
235  subls = self.get_sublists(path, s)
236  for subl in subls: all_groups.add(frozenset(subl))
237  return all_groups
238  # succ=nx.bfs_successors(gr,0)
239 
240  # succmol=dict([(self.mols[i].get_name()+str(IMP.atom.Copy(self.mols[i]).get_copy_index()),[self.mols[k].get_name()+str(IMP.atom.Copy(self.mols[i]).get_copy_index()) for k in ks]) for i,ks in succ.iteritems()])
241  # print(succmol)
242 
243  def do_compute_pcas(self):
244  self.singletons = {}
245  for i, ps in enumerate(self.particle_blocks):
246  mol = self.mols[i]
247  pca = IMP.algebra.get_principal_components([IMP.core.XYZ(p).get_coordinates() for p in ps])
248  ## MB normalize
249  pcan = IMP.bayesianem.NormalizePCA(pca, ps)
250  self.singletons[(mol,)] = pcan
251 
252  if self.depth > 1:
253  all_groups = self.get_list_of_connected_mols(self.depth)
254 
255  for group in all_groups:
256  mols = [self.mols[n] for n in group]
257  ps = IMP.atom.Selection(mols, resolution=self.resolution).get_selected_particles()
258  pca = IMP.algebra.get_principal_components([IMP.core.XYZ(p).get_coordinates() for p in ps])
259  pcan = IMP.bayesianem.NormalizePCA(pca, ps)
260  self.singletons[tuple(mols)] = pcan
261 
262  def do_get_pairs(self, temp=20):
263  sk = self.singletons.keys()
264 
265  pairs = []
266  for s1, s2 in itertools.combinations(sk, 2):
267  rms = IMP.algebra.get_distance(self.singletons[s1].get_principal_values(),
268  self.singletons[s2].get_principal_values())
269  pairs.append(((s1, s2), math.exp(-rms / temp)))
270 
271  return pairs
272 
273  def do_compute_transforms(self, pair):
274  amols = pair[0]
275  bmols = pair[1]
276 
277  arbs = []
278  abeads = []
279  for mol in amols:
280  arbs += self.rbsd[mol]
281  abeads += self.beadsd[mol]
282  apca = self.singletons[amols]
283 
284  brbs = []
285  bbeads = []
286  for mol in bmols:
287  brbs += self.rbsd[mol]
288  bbeads += self.beadsd[mol]
289  bpca = self.singletons[bmols]
290 
291  # transatob=IMP.algebra.get_alignments_from_first_to_second(apca,bpca)
292  # MB get single transformation
293  transatob = IMP.bayesianem.PCAalign(apca, bpca)[0]
294  transbtoa = transatob.get_inverse()
295 
296  return (transatob, transbtoa, arbs, abeads, brbs, bbeads)
297 
298  def weighted_choice(self, choices):
299  total = sum(w for c, w in choices)
300  r = random.uniform(0, total)
301  upto = 0
302  for c, w in choices:
303  if upto + w >= r:
304  return c
305  upto += w
306  assert False, "Shouldn't get here"
307 
308  def do_propose(self):
309  self.rejects.append(0)
310  if len(self.rejects) > self.lenrejects:
311  self.rejects.pop(0)
312 
313  self.propose += 1
314  self.do_store_current_coordinates()
315  self.do_compute_pcas()
316  pairs = self.do_get_pairs()
317  pair = self.weighted_choice(pairs)
318  # print pair,self.reject,self.propose,float(self.propose-self.reject)/self.propose,float(sum(self.rejects))/len(self.rejects)
319  # print self.rejects
320  t = self.do_compute_transforms(pair)
321  # (ta2b,tb2a)=random.choice(zip(t[0],t[1]))
322  # MB no need to randomly select
323  (ta2b, tb2a) = (t[0], t[1])
324 
325  moved_indexes = []
326  for arb in t[2]:
327  IMP.core.transform(arb, ta2b)
328  moved_indexes.append(arb.get_particle_index())
329  for abead in t[3]:
330  IMP.core.transform(IMP.core.XYZ(abead), ta2b)
331  moved_indexes.append(arb.get_particle_index())
332  for brb in t[4]:
333  IMP.core.transform(brb, tb2a)
334  for bbead in t[5]:
335  IMP.core.transform(IMP.core.XYZ(bbead), tb2a)
336  return IMP.core.MonteCarloMoverResult(moved_indexes, 1.0)
337 
338  def do_reject(self):
339  self.rejects[-1] = 1
340 
341  self.reject += 1
342  for rb in self.refframes:
343  rb.set_reference_frame(self.refframes[rb])
344  for bead in self.xyzs:
345  IMP.core.XYZ(bead).set_coordinates(self.xyzs[bead])
346 
347  def do_get_inputs(self):
348  return self.ps
A Monte Carlo optimizer.
Definition: MonteCarlo.h:44
IMP::algebra::Transformation3Ds PCAalign(const IMP::algebra::PrincipalComponentAnalysis &pca1, const IMP::algebra::PrincipalComponentAnalysis &pca2)
def get_molecules
This function returns the parent molecule hierarchies of given objects.
Definition: tools.py:1175
Miscellaneous utilities.
Definition: tools.py:1
Return value of the MonteCarloMover::propose() function.
IMP::algebra::PrincipalComponentAnalysis NormalizePCA(const IMP::algebra::PrincipalComponentAnalysis &pca, const IMP::Particles &ps)
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
A base class for classes which perturb particles.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
PrincipalComponentAnalysisD< D > get_principal_components(const Vector< VectorD< D > > &ps)
Perform principal components analysis on a set of vectors.
Basic functionality that is expected to be used by a wide variety of IMP users.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
double get_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the distance between two vectors.
Definition: VectorD.h:213
Floats get_list_of_bipartite_minimum_sphere_distance(const ParticlesTemps &pss)
Definition: pmi/utilities.h:71
def get_rbs_and_beads
Returns unique objects in original order.
Definition: tools.py:1151
Select hierarchy particles identified by the biological name.
Definition: Selection.h:70
Applies a list of movers one at a time.
Definition: SerialMover.h:26