IMP  2.0.0
The Integrative Modeling Platform
multiscale.py
1 ## \example domino/multiscale.py
2 ## We are interested in applying domino to problems systematically in a multiscale manner. This script experiments with those approaches.
3 
4 import IMP.domino
5 import IMP.core
6 
7 
8 # Use faster built-in 'set' type on newer Pythons; fall back to the older
9 # 'sets' module on older Pythons
10 try:
11  x = set
12  del x
13 except NameError:
14  import sets
15  set = sets.Set
16 
17 m= IMP.Model()
18 m.set_log_level(IMP.base.SILENT)
19 ds= [IMP.core.XYZR.setup_particle(IMP.Particle(m)) for i in range(0,3)]
20 for i,d in enumerate(ds):
21  d.set_radius(1)
22  IMP.display.Colored.setup_particle(d, IMP.display.get_display_color(i))
23 
24 k=1
25 h= IMP.core.Harmonic(0,k)
26 r0= IMP.core.SingletonRestraint(IMP.core.DistanceToSingletonScore(h, IMP.algebra.Vector3D(0,0,0)),
27  ds[0], "0 at origin")
28 r0.set_model(m)
29 
30 r1= IMP.core.SingletonRestraint(IMP.core.AttributeSingletonScore(h,
31  IMP.core.XYZ.get_xyz_keys()[0]),
32  ds[1], "1 on axis")
33 r1.set_model(m)
34 
35 rs=[r0, r1]
36 for pr in [(0,1), (1,2), (0,2)]:
37  r= IMP.core.PairRestraint(IMP.core.HarmonicSphereDistancePairScore(0,k),
38  (ds[pr[0]], ds[pr[1]]),
39  "R for "+str(pr))
40  r.set_model(m)
41  rs.append(r)
42 
43 
44 bb= IMP.algebra.BoundingBox2D(IMP.algebra.Vector2D(0,0),
45  IMP.algebra.Vector2D(4,4))
46 
47 covers=[]
48 for i in range(0,6):
49  cur=IMP.algebra.get_grid_interior_cover_by_spacing(bb, 4.0/2**i)
50  print cur
51  covers.append([IMP.algebra.Vector3D(x[0], x[1], 0) for x in cur])
52 
53 
54 def setup(cover, scale):
55  pst= IMP.domino.ParticleStatesTable()
56  st= IMP.domino.XYZStates(cover)
57  for p in ds:
58  pst.set_particle_states(p, st)
59  for r in rs:
60  r.set_maximum_score(.5*scale**2)
61  lf= IMP.domino.ListSubsetFilterTable(pst)
62  rc= IMP.domino.RestraintCache(pst)
63  rc.add_restraints(rs);
64  fs=[IMP.domino.RestraintScoreSubsetFilterTable(rc),
65  lf]
66  sampler= IMP.domino.DominoSampler(m, pst)
67  sampler.set_subset_filter_tables(fs)
68  sampler.set_log_level(IMP.base.SILENT)
69  return (sampler, lf, pst)
70 
71 (sampler,lf, pst)= setup(covers[0], 4.0)
72 
73 subset=IMP.domino.Subset(ds)
74 ac= sampler.get_sample_assignments(subset)
75 
76 print ac
77 
78 def get_mapping(cover0, cover1):
79  nn= IMP.algebra.NearestNeighbor3D(cover0)
80  ret=[[] for c in cover0]
81  for i,p in enumerate(cover1):
82  nns=nn.get_nearest_neighbor(p)
83  ret[nns].append(i)
84  return ret
85 
86 mw= IMP.display.PymolWriter("mapping.pym")
87 def display_mapping(index, cover0, cover1, mapping):
88  mw.set_frame(index)
89  for i,c in enumerate(mapping):
90  for p in c:
91  g= IMP.display.PointGeometry(cover1[p])
92  g.set_color(IMP.display.get_display_color(i))
93  g.set_name("fine")
94  mw.add_geometry(g)
95  for i,c in enumerate(cover0):
96  g= IMP.display.PointGeometry(c)
97  g.set_color(IMP.display.get_display_color(i))
98  g.set_name("coarse")
99  mw.add_geometry(g)
100 
101 
102 
103 
104 
105 for curi in range(1,len(covers)):
106  scale= 4.0/2**curi
107  print scale
108  mapping=get_mapping(covers[curi-1], covers[curi])
109  print mapping
110  display_mapping(curi-1, covers[curi-1], covers[curi], mapping)
111  (sampler, lf, pst)= setup(covers[curi], scale)
112  lac= ac
113  cac=[]
114  for a in lac:
115  for i,p in enumerate(subset):
116  s= a[i]
117  allowed= mapping[s]
118  lf.set_allowed_states(p, allowed)
119  ccac=sampler.get_sample_assignments(subset)
120  print a, ccac
121  cac= cac+ccac
122  ac= list(set(cac))
123  print "for scale", scale, "got", ac
124  sw= IMP.display.PymolWriter("solutions."+str(curi)+".pym")
125  for i,a in enumerate(ac):
126  IMP.domino.load_particle_states(subset, a, pst)
127  sw.set_frame(i)
128  for p in ds:
129  g= IMP.core.XYZRGeometry(p)
130  sw.add_geometry(g)
131  for c in covers[curi]:
132  g= IMP.display.PointGeometry(c)
133  g.set_color(IMP.display.Color(1,1,1))
134  g.set_name("grid")
135  sw.add_geometry(g)