 |
This example looks like the six particle optimization example except a filter is used instead of a restraint to remove the flip degree of freedom. The filter is written is python, which makes for quick prototyping, but slow run times.
import IMP
import IMP.domino
import IMP.core
import IMP.container
#set restraints
def create_scoring(m, ps):
pairs=[[0,1],[0,2],[1,3],[2,3],[3,4],[4,5],[1,5]]
score= IMP.core.HarmonicDistancePairScore(1, 1)
# the restraint will be broken apart during optimization
"Restrained pairs")
pr= IMP.container.PairsRestraint(score, pc)
pr.set_maximum_score(.01)
pr.set_model(m)
IMP.algebra.Vector3D(2,0,0))
return [pr]
def create_representation(m):
ps=[]
for i in range(0,6):
p=IMP.Particle(m)
IMP.core.XYZ.setup_particle(p,IMP.algebra.Vector3D(i,0.,0.))
ps.append(p)
return ps
def create_discrete_states(ps):
vs=[IMP.algebra.Vector3D(1,0,0),
IMP.algebra.Vector3D(0,1,0),
IMP.algebra.Vector3D(1,1,0),
IMP.algebra.Vector3D(2,1,0),
IMP.algebra.Vector3D(2,0,0)]
vs= vs+[-v for v in vs]
print len(vs), "states for each particle"
print vs[1]
states= IMP.domino.XYZStates(vs)
# special case ps[0] to remove a sliding degree of freedom
for p in ps[1:]:
pst.set_particle_states(p, states)
return pst
# force particle p to be in state s
class MyFilterTable(IMP.domino.SubsetFilterTable):
class MyFilter(IMP.domino.SubsetFilter):
def __init__(self, pos, value):
#print "Filtering with", pos, value
IMP.domino.SubsetFilter.__init__(self, "MF"+str(pos) + " " +str(value))
self.pos=pos
self.value=value
def get_next_state(self, pos, s):
# suggest that the sampler try the correct state
# this method is only called if the filter failed, so pos must be
# self.pos
if s[self.pos] > self.value:
# a very large number
return 2**29
else:
return self.value
def get_is_ok(self, state):
# it is only OK if it has the required state
ret= state[self.pos]==self.value
return ret
def get_strength(self, s, excluded):
# return the maximum value since it dictates the state for the particle
return 1
def __init__(self, p, s):
# set the name of the object to something vaguely useful
IMP.domino.SubsetFilterTable.__init__(self, "MFT"+p.get_name()+" at "+str(s))
self.p=p
self.s=s
def get_subset_filter(self, subset, excluded):
# create a filter if self.p is in subset but not in excluded
if self.p in subset and self.p not in sum([list(x) for x in excluded], []):
# pass the position of self.p and the value that it must have
return self.MyFilter(list(subset).index(self.p), self.s)
else:
return None
def create_sampler(m, ps, rs, pst):
s=IMP.domino.DominoSampler(m, pst)
#s.set_log_level(IMP.VERBOSE)
# the following lines recreate the defaults and so are optional
filters=[]
# do not allow particles with the same ParticleStates object
# to have the same state index
filters.append(IMP.domino.ExclusionSubsetFilterTable(pst))
rc= IMP.domino.RestraintCache(pst)
rc.add_restraints(rs)
# filter states that score worse than the cutoffs in the Model
filters.append(IMP.domino.RestraintScoreSubsetFilterTable(rc))
filters[-1].set_log_level(IMP.SILENT)
mf=MyFilterTable(ps[1], 0)
# try with and without this line
filters.append(mf)
states= IMP.domino.BranchAndBoundAssignmentsTable(pst, filters)
#states.set_log_level(IMP.SILENT);
s.set_assignments_table(states)
s.set_subset_filter_tables(filters)
return s
IMP.set_log_level(IMP.TERSE)
m=IMP.Model()
m.set_log_level(IMP.SILENT)
print "creating representation"
ps=create_representation(m)
print "creating discrete states"
pst=create_discrete_states(ps)
print "creating score function"
rs=create_scoring(m, ps)
print "creating sampler"
s=create_sampler(m, ps, rs, pst)
#s.set_log_level(IMP.SILENT)
print "sampling"
cs=s.get_sample()
print "found ", cs.get_number_of_configurations(), "solutions"
for i in range(cs.get_number_of_configurations()):
cs.load_configuration(i)
print "solution number:",i," is:", m.evaluate(False)
for p in ps:
|
|
Optimize six particles
#### NOT FULLY IMPLEMENTED YET!
import IMP
import IMP.domino
import IMP.core
import sys
def optimize_subsets(subsets):
for subset in subsets:
None
#TODO - do the actual sampling
#combine back to ParticleStatesTable
#there is now a write_particles_binary (multiple states into one file)
return pst
def get_subsets(ps):
mdl=ps[0].get_model()
inter_graph=IMP.domino.get_interaction_graph(mdl.get_root_restraint_set(), ps)
jt=IMP.domino.get_junction_tree(inter_graph)
return IMP.domino.get_subsets(jt)
def setup_scoring_function(ps):
m=ps[0].get_model()
pairs=[[0,1],[0,2],[1,2],[2,3],[3,4],[4,5],[3,5]]
sf = IMP.core.Harmonic(1.0, 1)
for pair in pairs:
r=IMP.core.DistanceRestraint(sf, ps[pair[0]], ps[pair[1]])
m.add_restraint(r)
#Initiate a set of states for each particle in ps
def initiate_configuration(domino_smp,ps):
#Generate a discrete set of states
vector_states=[]
for i in range(0,4):
vector_states.append(IMP.algebra.Vector3D(i,0,0))
vector_states.append(IMP.algebra.Vector3D(i,1,0))
#Generate a discrete set of states for each of the particles
states= IMP.domino.XYZStates(vector_states)
#Map states of particles
for p in ps:
domino_smp.set_particle_states(p,states)
sys.exit()
#### REPRESENTATION
#1. setting up the representation (6 particles)
mdl=IMP.Model()
mdl.set_log_level(IMP.SILENT)
ps=[]
for i in range(0,6):
p=IMP.Particle(mdl)
IMP.core.XYZ.setup_particle(p,IMP.algebra.Vector3D(0.,0.,0.))
ps.append(p)
####SCORING
#1. setting up the scoring function
setup_scoring_function(ps)
#1. get the subsets
subsets=get_subsets(ps)
#optimize each one (returning ParticleStatesTable)
pst = optimize_subsets(subsets)
sys.exit()
#subsets=[]
#jt.show()
#2. sample each subset
#3. gathering
####OPTIMIZATION
#1. load domino sampler and set required properties
domino_smp=IMP.domino.DominoSampler(m)
domino_smp.set_maximum_score(.2)
#2. initiate configuration
initiate_configuration(domino_smp,ps)
#3. construct subset of variables AND
# optimize subsets of variables subject to the discrete sampling space AND
# gather subset of variables into global solutions
cs=domino_smp.get_sample()
####ANALYSE
#4. report solutions
print "Found ", cs.get_number_of_configurations(), "solutions"
for i in range(cs.get_number_of_configurations()):
cs.load_configuration(i)
#print the configuration:
print "solution number:",i," scored:", m.evaluate(False)
|
IMP::domino::DominoSampler supports an interactive mode where by each step of the sampling process is called explicitly and all intermediate results are exposed. This usage mode can be used to help understand how domino behaves as well as to distribute a domino computation to multiple cores or multiple machines. For the latter, it is useful to use the IMP::domino::get_assignments() and IMP::domino::set_assignments() functions to save the states to a file.
import IMP.domino
import IMP.algebra
import IMP.container
import IMP
m = IMP.Model()
# create some particles
s= IMP.core.HarmonicDistancePairScore(1, 1)
print [(p[0].get_name(), p[1].get_name()) for p in lpc.get_particle_pairs()]
r= IMP.container.PairsRestraint(s, lpc)
r.set_model(m)
r.set_maximum_score(.1)
for p in ps:
pst.set_particle_states(p, space)
m.set_log_level(IMP.SILENT)
# make sure to break up the
mt= IMP.domino.get_merge_tree([r], pst)
try:
except:
print "Unable to display graph using 'dot'"
ds= IMP.domino.DominoSampler(m, pst)
# use the default setup for filters
ds.set_scoring_function([r])
ds.set_merge_tree(mt)
ds.set_log_level(IMP.SILENT)
# recurse down the tree getting the assignments and printing them
def get_assignments(vertex):
on= mt.get_out_neighbors(vertex)
if len(on)==0:
# we have a leaf
ret= ds.get_vertex_assignments(vertex)
else:
# recurse on the two children
if on[0] > on[1]:
# the get_vertex_assignment methods expects the children in sorted order
on=[on[1], on[0]]
a0= get_assignments(on[0])
a1= get_assignments(on[1])
ret= ds.get_vertex_assignments(vertex, a0, a1)
print mt.get_vertex_name(vertex), [str(r) for r in ret]
return ret
# the root is the last vetex
get_assignments(mt.get_vertices()[-1])
schedule=[]
# we could instead decompose the tree into independent sets of jobs
def schedule_job(vertex):
# first figure out the maximum phase for my children, then add me to the
# next higher phase
on= mt.get_out_neighbors(vertex)
max_child_time=-1
for n in on:
max_child_time=max(schedule_job(n), max_child_time)
my_time=max_child_time+1
while len(schedule) < my_time+1:
schedule.append([])
# add myself to the next free phase
schedule[my_time].append(vertex)
return my_time
schedule_job(mt.get_vertices()[-1])
print "The merging can be scheduled as", schedule
|
Run domino storing the intermediate and final results in an hdf5 database. This has the advantage that if you interrupt the run at any point, you have a list of everything computed so far and so can get some more idea of what is going on.
import IMP.domino
import IMP.algebra
import IMP.container
import IMP
import RMF
m = IMP.Model()
# create some particles
s= IMP.core.HarmonicDistancePairScore(1, 1)
print [(p[0].get_name(), p[1].get_name()) for p in lpc.get_particle_pairs()]
r= IMP.container.PairsRestraint(s, lpc)
r.set_model(m)
r.set_maximum_score(.1)
for p in ps:
pst.set_particle_states(p, space)
m.set_log_level(IMP.SILENT)
# make sure to break up the
mt= IMP.domino.get_merge_tree([r], pst)
try:
except:
print "Unable to display graph using 'dot'"
rc= IMP.domino.RestraintCache(pst)
rc.add_restraints([r])
filters=[IMP.domino.RestraintScoreSubsetFilterTable(rc),
leaf_table=IMP.domino.BranchAndBoundAssignmentsTable(pst, filters)
# create a database to store the results
name=IMP.create_temporary_file_name("assignments", ".hdf5")
root= RMF.create_hdf5_file(name)
# recurse down the tree getting the assignments and printing them
def get_assignments(vertex):
on= mt.get_out_neighbors(vertex)
ss= mt.get_vertex_name(vertex)
print "computing assignments for", ss
ssn= str(ss)
dataset= root.add_child_index_data_set_2d(ssn)
dataset.set_size([0, len(ss)])
mine= IMP.domino.WriteHDF5AssignmentContainer(dataset, ss, pst.get_particles(), ssn)
if len(on)==0:
# we have a leaf
IMP.domino.load_leaf_assignments(ss, leaf_table, mine)
else:
# recurse on the two children
(ss0, a0)= get_assignments(on[0])
(ss1, a1)= get_assignments(on[1])
IMP.domino.load_merged_assignments(ss0, a0, ss1, a1, filters, mine)
print ss, mine.get_number_of_assignments()
# make sure that the cache is flushed
del mine
return (ss, IMP.domino.ReadHDF5AssignmentContainer(dataset, ss, pst.get_particles(), ssn))
# the root is the last vetex
all=get_assignments(mt.get_vertices()[-1])
all[1].set_was_used(True)
print 'try: h5dump', name
|
|
The example shows how to generate and inspect a merge tree for use in Domino.
import IMP
import IMP.domino
import IMP.core
import IMP.container
import IMP.algebra
IMP.set_log_level(IMP.TERSE)
m=IMP.Model()
# don't print messages about evaluation
m.set_log_level(IMP.SILENT)
bb= IMP.algebra.BoundingBox3D((0,0,0), (10,10,10))
allc=[]
for i in range(0,7):
p= IMP.Particle(m)
d= IMP.core.XYZR.setup_particle(p, IMP.algebra.Sphere3D(IMP.algebra.get_random_vector_in(bb), 1))
allc.append(d.get_coordinates())
ss= IMP.domino.XYZStates(allc)
for p in m.get_particles():
pst.set_particle_states(p, ss)
# generate a set of restraints based on the close pairs in this randomly chosen configuration
cp.set_distance(1)
cps=cp.get_close_pairs(m.get_particles())
if len(cps)>0:
# one cannot create a container from an empty list
acp= IMP.container.ListPairContainer(cps)
else:
acp= IMP.container.ListPairContainer(m)
ps= IMP.core.SoftSpherePairScore(1)
r= IMP.container.PairsRestraint(ps, acp)
r.set_model(m)
# compute the interaction graph based on all the restraints
ig= IMP.domino.get_interaction_graph([r],
pst)
# generate a junction tree from the interaction graph
jt= IMP.domino.get_junction_tree(ig)
print dir(jt)
print type(jt)
# create a merge tree from the junction tree, this can be passed to IMP.domin.DominoSampler
mt= IMP.domino.get_merge_tree(jt)
s=pst.get_subset()
print s, type(s)
|
|
We are interested in applying domino to problems systematically in a multiscale manner. This script experiments with those approaches.
import IMP.domino
import IMP.core
# Use faster built-in 'set' type on newer Pythons; fall back to the older
# 'sets' module on older Pythons
try:
x = set
del x
except NameError:
import sets
set = sets.Set
m= IMP.Model()
m.set_log_level(IMP.SILENT)
for i,d in enumerate(ds):
d.set_radius(1)
k=1
h= IMP.core.Harmonic(0,k)
ds[0], "0 at origin")
r0.set_model(m)
IMP.core.XYZ.get_xyz_keys()[0]),
ds[1], "1 on axis")
r1.set_model(m)
rs=[r0, r1]
for pr in [(0,1), (1,2), (0,2)]:
(ds[pr[0]], ds[pr[1]]),
"R for "+str(pr))
r.set_model(m)
rs.append(r)
IMP.algebra.Vector2D(4,4))
covers=[]
for i in range(0,6):
cur=IMP.algebra.get_grid_interior_cover_by_spacing(bb, 4.0/2**i)
print cur
def setup(cover, scale):
st= IMP.domino.XYZStates(cover)
for p in ds:
pst.set_particle_states(p, st)
for r in rs:
r.set_maximum_score(.5*scale**2)
lf= IMP.domino.ListSubsetFilterTable(pst)
rc= IMP.domino.RestraintCache(pst)
rc.add_restraints(rs);
lf]
sampler= IMP.domino.DominoSampler(m, pst)
sampler.set_subset_filter_tables(fs)
sampler.set_log_level(IMP.SILENT)
return (sampler, lf, pst)
(sampler,lf, pst)= setup(covers[0], 4.0)
subset=IMP.domino.Subset(ds)
ac= sampler.get_sample_assignments(subset)
print ac
def get_mapping(cover0, cover1):
nn= IMP.algebra.NearestNeighbor3D(cover0)
ret=[[] for c in cover0]
for i,p in enumerate(cover1):
nns=nn.get_nearest_neighbor(p)
ret[nns].append(i)
return ret
mw= IMP.display.PymolWriter("mapping.pym")
def display_mapping(index, cover0, cover1, mapping):
mw.set_frame(index)
for i,c in enumerate(mapping):
for p in c:
g= IMP.display.PointGeometry(cover1[p])
g.set_color(IMP.display.get_display_color(i))
g.set_name("fine")
mw.add_geometry(g)
for i,c in enumerate(cover0):
g.set_color(IMP.display.get_display_color(i))
g.set_name("coarse")
mw.add_geometry(g)
for curi in range(1,len(covers)):
scale= 4.0/2**curi
print scale
mapping=get_mapping(covers[curi-1], covers[curi])
print mapping
display_mapping(curi-1, covers[curi-1], covers[curi], mapping)
(sampler, lf, pst)= setup(covers[curi], scale)
lac= ac
cac=[]
for a in lac:
for i,p in enumerate(subset):
s= a[i]
allowed= mapping[s]
lf.set_allowed_states(p, allowed)
ccac=sampler.get_sample_assignments(subset)
print a, ccac
cac= cac+ccac
ac= list(set(cac))
print "for scale", scale, "got", ac
for i,a in enumerate(ac):
IMP.domino.load_particle_states(subset, a, pst)
sw.set_frame(i)
for p in ds:
g= IMP.core.XYZRGeometry(p)
sw.add_geometry(g)
for c in covers[curi]:
g.set_color(IMP.display.Color(1,1,1))
g.set_name("grid")
sw.add_geometry(g)
|
|
Caching restraint scores so that restraints are not evaluated repeatedly for the same configuration is an important part of domino. Without caching, sub assignments that are shared between subsets will be rescored. The IMP::domino::RestraintCache provides a centralized place for this. To use it, one creates one and then adds the restraints you want to use for filtering and scoring to it. You can then pass the cache to the IMP::domino::RestraintScoreFilterTable and it will filter based on those restraints. You can also extract scores from the table directly, using it to manage the loading of particle states.
import IMP.domino
import IMP.algebra
import IMP.container
import IMP.atom
import IMP.rmf
import RMF
resolution=.5
m= IMP.Model()
# suppress messages about restraint evaluations
m.set_log_level(IMP.SILENT)
# create some particles and a restraint that scores based on the
# length of the chain
ps= [IMP.core.XYZR.setup_particle(IMP.Particle(m))
for i in range(0,int(2.0/resolution))]
eps[0].set_coordinates(IMP.algebra.Vector3D(-1,0,0))
eps[1].set_coordinates(IMP.algebra.Vector3D(1,0,0))
for p in ps:
p.set_radius(resolution/2.0)
for p in eps:
p.set_radius(.4*resolution)
order=[eps[0]]+ps+[eps[1]]
# create a hierarchy with the particles to aid writing them to RMF
for i,p in enumerate(order):
h.add_child(IMP.atom.Hierarchy.setup_particle(p))
color= IMP.display.get_jet_color(float(i)/(len(order)-1))
IMP.display.Colored.setup_particle(p, color)
# allow a maximum distance of 2
r.set_maximum_score(0.01)
# dumb way to generate points on a gride
g= IMP.algebra.DenseFloatGrid3D(resolution, IMP.algebra.get_unit_bounding_box_3d())
# create some random sites
vs= [g.get_center(i) for i in g.get_all_indexes()]
print "States are", vs
states=IMP.domino.XYZStates(vs)
for p in ps:
pst.set_particle_states(p, states)
# now create a restraint cache
# cache the most recently used one million scores
rc= IMP.domino.RestraintCache(pst, 1000000)
rc.add_restraints([r])
s=IMP.domino.DominoSampler(m, pst)
s.set_check_level(IMP.NONE)
# pass the cache to the restraint score based filter
# it will use all the restraints in the cache
rssft= IMP.domino.RestraintScoreSubsetFilterTable(rc)
s.set_subset_filter_tables([ef, rssft])
# create a subset with all the particles
# Subsets keep the particles in sorted order and so are different than
# a simple list of particles
alls= IMP.domino.Subset(ps)
# get all the assignments that fit
sts= s.get_sample_assignments(alls)
# create a temporary file to write things to
print "saving configurations to", rmf.get_name()
IMP.rmf.add_hierarchy(rmf, h)
# the restraint cache processes the restraints to make them more efficient
# for use with domino. So we want to get back the processed restraints before
# looking at them further
domino_restraints= rc.get_restraints()
IMP.rmf.add_restraints(rmf, domino_restraints)
# for each assignment load it, and add the configuration to an rmf file
print "found assignments", sts
# we don't care about messages during saving
for i, s in enumerate(sts):
IMP.domino.load_particle_states(alls, s, pst)
# go through the restraints and get the score
# here, we only care about its side effect of setting the last score
for r in domino_restraints:
rc.load_last_score(r, alls, s)
# save the configuration and scores to the rmf
IMP.rmf.save_frame(rmf, i)
|
|
This example shows using two rigid bodies and doing excluded volume with them.
import IMP
import IMP.core
import IMP.algebra
import IMP.domino
import IMP.atom
import math
# create a rigid body per helix
def create_representation():
m=IMP.Model()
for h in [h0, h1]:
return (m, [h0, h1])
def create_excluded_volume(m, helices):
# this is the interesting function:
# it uses a KClosePairsPair score to generate the list of close atoms on the fly
all=[]
for h in helices:
all.extend(IMP.atom.get_by_type(h, IMP.atom.ATOM_TYPE))
lsc= IMP.container.ListSingletonContainer(all)
evr= IMP.core.ExcludedVolumeRestraint(lsc, 1)
evr.set_model(m)
evr.set_maximum_score(.01)
return [evr]
# creating the discrete states for domino
def create_discrete_states(m,helices):
rot00= IMP.algebra.get_rotation_about_axis(IMP.algebra.Vector3D(0,1,0), math.pi/2.0)
trs=[]
zv=IMP.algebra.Vector3D(0,0,0)
trs.append(IMP.algebra.ReferenceFrame3D(IMP.algebra.Transformation3D(rot00,zv)))
for dx in range(0,15):
tr=IMP.algebra.Vector3D(1.0*dx,0,0)
trs.append(IMP.algebra.ReferenceFrame3D(IMP.algebra.Transformation3D(rot00,tr)))
pstate= IMP.domino.RigidBodyStates(trs)
for h in helices:
pst.set_particle_states(IMP.core.RigidMember(h).get_rigid_body(), pstate)
return pst
# setting up domino (and filters)
def create_sampler(m, rs, pst):
s=IMP.domino.DominoSampler(m, pst)
filters=[]
# do not allow particles with the same ParticleStates object
# to have the same state index
filters.append(IMP.domino.ExclusionSubsetFilterTable(pst))
rc= IMP.domino.RestraintCache(pst)
rc.add_restraints(rs)
# filter states that score worse than the cutoffs in the Model
filters.append(IMP.domino.RestraintScoreSubsetFilterTable(rc))
states= IMP.domino.BranchAndBoundAssignmentsTable(pst, filters)
s.set_assignments_table(states)
s.set_subset_filter_tables(filters)
return s
def display(m,helices,name):
m.update()
w= IMP.display.PymolWriter(name)
for i,h in enumerate(helices):
g.set_color(IMP.display.get_display_color(i))
w.add_geometry(g)
IMP.set_log_level(IMP.SILENT)
print "creating representation"
(m,helices)=create_representation()
print "creating score function"
rs=create_excluded_volume(m,helices)
print "creating discrete states"
pst=create_discrete_states(m,helices)
print "creating sampler"
s=create_sampler(m, rs, pst)
m.set_log_level(IMP.SILENT)
IMP.set_log_level(IMP.VERBOSE)
print "sampling"
cs=s.get_sample()
print "found ", cs.get_number_of_configurations(), "solutions"
score=[]
for i in range(cs.get_number_of_configurations()):
cs.load_configuration(i)
ss=m.evaluate(False)
score.append(ss)
print "** solution number:",i," is:",ss
display(m,helices,"sol_"+str(i)+".pym")
|
|
It can often be useful to save the assignments to a file. IMP provides support to do this to an a data set in an hdf5 file.
You can use h5dump to view the contents of the created file.
import IMP.domino
import RMF
# create a model and some particles, they are just used as markers here
m= IMP.Model()
# create a subset with a few of the particles
ss= IMP.domino.Subset([ps[3], ps[5], ps[7]])
file_name= IMP.create_temporary_file_name("assignments", ".hdf5")
print "File name is", file_name
# create a list of assignments
wcn=IMP.create_temporary_file_name("assignments", ".asn")
asl=IMP.domino.WriteAssignmentContainer(wcn, ss, ps, "writer")
written=[]
for i in range(0,5):
for j in range(0,5):
for k in range(0,5):
a=IMP.domino.Assignment([i,j,k])
written.append(a)
asl.add_assignment(a)
del asl
# to check, we can read it back immediately
back_asl= IMP.domino.ReadAssignmentContainer(wcn, ss, ps, "reader")
if back_asl.get_assignments()==written:
print "They match!"
else:
print "They don't match :-("
raise 1
# More interestingly, we can create a new model and read back the assignments for that
mp= IMP.Model()
# create a subset with a few of the particles
ssp= IMP.domino.Subset([psp[3], psp[5], psp[7]])
print "Note the subsets are differently ordered (most of the time): ", ss, ssp
back_aslp= IMP.domino.ReadAssignmentContainer(wcn, ssp, psp, "reader2")
# however, the states are demuxed so they match the particles
print [str(a) for a in back_aslp.get_assignments()]
|
Optimize six particles on a 2D unit grid. In order to remove translation degrees of freedom, the 0th particle is pinned at the origin by allowing it only a single conformation. To remove flips, the first particle is restrained to have a positive x coordinate.
import IMP
import IMP.domino
import IMP.core
import IMP.container
#set restraints
def create_scoring(m, ps):
pairs=[[0,1],[0,2],[1,3],[2,3],[3,4],[4,5],[1,5]]
# we will restrain various pairs to be 1 apart
score= IMP.core.HarmonicDistancePairScore(1, 1)
# the restraint will be broken apart during optimization
# map the indices above to actual particles
"Restrained pairs")
pr= IMP.container.PairsRestraint(score, pc)
pr.set_maximum_score(.01)
pr.set_model(m)
IMP.algebra.Vector3D(2,0,0))
# force ps[1] to be on the positive side to remove flip degree of freedom
dr= IMP.core.SingletonRestraint(d, ps[1])
dr.set_model(m)
# we are not interested in conformations which don't fit the distances
# exactly, but using 0 is tricky
dr.set_maximum_score(.01)
print m.get_root_restraint_set()
return [pr, dr]
def create_representation(m):
ps=[]
# create size particles, initial coordinates don't matter.
for i in range(0,6):
p=IMP.Particle(m)
IMP.core.XYZ.setup_particle(p,IMP.algebra.Vector3D(i,0.,0.))
ps.append(p)
return ps
def create_discrete_states(ps):
vs=[IMP.algebra.Vector3D(1,0,0),
IMP.algebra.Vector3D(0,1,0),
IMP.algebra.Vector3D(1,1,0),
IMP.algebra.Vector3D(2,1,0),
IMP.algebra.Vector3D(2,0,0)]
vs= vs+[-v for v in vs]
print len(vs), "states for each particle"
states= IMP.domino.XYZStates(vs)
# special case ps[0] to remove a sliding degree of freedom
# all other particles are given the same set of states
for p in ps[1:]:
pst.set_particle_states(p, states)
return pst
def create_sampler(m, r, pst):
# create the sampler and pass it the states for each patricle
s=IMP.domino.DominoSampler(m, pst)
# the following lines recreate the defaults and so are optional
filters=[]
# create a restraint cache to avoid re-evaluating restraints
rc= IMP.domino.RestraintCache(pst)
# add the list of restraints we want to use
rc.add_restraints(r)
# do not allow particles with the same ParticleStates object
# to have the same state index
filters.append(IMP.domino.ExclusionSubsetFilterTable(pst))
# filter states that score worse than the cutoffs in the Model
filters.append(IMP.domino.RestraintScoreSubsetFilterTable(rc))
filters[-1].set_log_level(IMP.SILENT)
# try to be intelligent about enumerating the states in each subset
states= IMP.domino.BranchAndBoundAssignmentsTable(pst, filters);
states.set_log_level(IMP.SILENT);
s.set_assignments_table(states)
s.set_subset_filter_tables(filters)
return s
IMP.set_log_level(IMP.TERSE)
m=IMP.Model()
# don't print information during Model.evaluate
m.set_log_level(IMP.SILENT)
print "creating representation"
ps=create_representation(m)
print "creating discrete states"
pst=create_discrete_states(ps)
print "creating score function"
rs=create_scoring(m, ps)
print "creating sampler"
s=create_sampler(m, rs, pst)
print "sampling"
# get an IMP.ConfigurationSet with the sampled states. If there are very
# many, it might be better to use s.get_sample_states() and then
# IMP.domino.load_particle_states() to handle the states as that takes
# much less memory, and time.
cs=s.get_sample()
print "found ", cs.get_number_of_configurations(), "solutions"
for i in range(cs.get_number_of_configurations()):
cs.load_configuration(i)
print "solution number:",i," is:", m.evaluate(False)
for p in ps:
|
|