IMP logo

domino examples

custom filter

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

#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
    pc= IMP.container.ListPairContainer([(ps[p[0]], ps[p[1]]) for p in pairs],
                                         "Restrained pairs")
    pr= IMP.container.PairsRestraint(score, pc)
    m.set_maximum_score(pr, .01)
    m.add_restraint(pr)
    d= IMP.core.DistanceToSingletonScore(IMP.core.HarmonicUpperBound(2,1),
                                         IMP.algebra.Vector3D(2,0,0))
    return m.get_restraints()

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):
    pst= IMP.domino.ParticleStatesTable()
    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_is_ok(self, state):
            ret= state[self.pos]==self.value
            return ret
    def get_strength(self, s, excluded):
        # return the maximum value since it dictates the position
        return 1
    def __init__(self, p, s):
        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, 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))
    # filter states that score worse than the cutoffs in the Model
    filters.append(IMP.domino.RestraintScoreSubsetFilterTable(m, pst))
    filters[-1].set_log_level(IMP.SILENT)
    mf=MyFilterTable(ps[1], 0)
    # try with and without this line
    filters.append(mf)
    states= IMP.domino.BranchAndBoundSubsetStatesTable(pst, filters)
    #states.set_log_level(IMP.SILENT);
    s.set_subset_states_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, pst)
print "creating score function"
rs=create_scoring(m, ps)

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:
        print IMP.core.XYZ(p).get_x(), IMP.core.XYZ(p).get_y()

domino approach

Optimize six particles

#### NOT FULLY IMPLEMENTED YET!
import IMP
import IMP.domino
import IMP.core
import sys

def optimize_subsets(subsets):
    pst= IMP.domino.ParticleStatesTable()
    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)

rigid body excluded volume

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()
    h0= IMP.atom.read_pdb(IMP.domino.get_example_path('helix_0.pdb'), m, IMP.atom.CAlphaPDBSelector())
    h1= IMP.atom.read_pdb(IMP.domino.get_example_path('helix_1.pdb'), m, IMP.atom.CAlphaPDBSelector())
    for h in [h0, h1]:
        IMP.atom.create_rigid_body(h)
    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)
    cpf=IMP.core.RigidClosePairsFinder()
    cpf.set_log_level(IMP.SILENT)
    nbl= IMP.container.ClosePairContainer(lsc, 0, cpf)
    nbl.set_log_level(IMP.SILENT)
    ps= IMP.core.SoftSpherePairScore(1000)
    evr= IMP.container.PairsRestraint(ps, nbl)
    m.add_restraint(evr)
    m.set_maximum_score(evr, .01)


# creating the discrete states for domino
def  create_discrete_states(m,helices):
    pst= IMP.domino.ParticleStatesTable()
    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, 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))
    # filter states that score worse than the cutoffs in the Model
    filters.append(IMP.domino.RestraintScoreSubsetFilterTable(m, pst))
    states= IMP.domino.BranchAndBoundSubsetStatesTable(pst, filters)
    s.set_subset_states_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= IMP.display.HierarchyGeometry(h)
        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"
create_excluded_volume(m,helices)

print "creating discrete states"
pst=create_discrete_states(m,helices)

print "creating sampler"
s=create_sampler(m, 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")

six particles loopy optimization

Like the normal six particle optimization example, but use loopy optimization to compute the same answer.

import IMP
import IMP.domino
import IMP.core

#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
    pc= IMP.container.ListPairContainer([(ps[p[0]], ps[p[1]]) for p in pairs],
                                         "Restrained pairs")
    pr= IMP.container.PairsRestraint(score, pc)
    m.set_maximum_score(pr, .01)
    m.add_restraint(pr)
    d= IMP.core.DistanceToSingletonScore(IMP.core.HarmonicUpperBound(2,1),
                                         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])
    m.add_restraint(dr)
    # we are not interested in conformations which don't fit the distances
    # exactly, but using 0 is tricky
    m.set_maximum_score(dr, .01)
    return m.get_restraints()

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):
    pst= IMP.domino.ParticleStatesTable()
    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
    for p in ps[1:]:
        pst.set_particle_states(p, states)
    return pst

def create_sampler(m, pst):
    s=IMP.domino.DominoSampler(m, pst)
    s.set_log_level(IMP.VERBOSE)
    # instead of the default junction tree, use the restraint graph
    IMP.set_log_level(IMP.VERBOSE)
    sg= IMP.domino.get_restraint_graph(m.get_root_restraint_set(), pst)
    IMP.set_log_level(IMP.TERSE)
    print sg
    s.set_subset_graph(sg)
    # 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))
    # filter states that score worse than the cutoffs in the Model
    filters.append(IMP.domino.RestraintScoreSubsetFilterTable(m, pst))
    filters[-1].set_log_level(IMP.SILENT)
    states= IMP.domino.BranchAndBoundSubsetStatesTable(pst, filters);
    states.set_log_level(IMP.SILENT);
    s.set_subset_states_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, pst)

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:
        print IMP.core.XYZ(p).get_x(), IMP.core.XYZ(p).get_y()

six particles optimization

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

#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
    pc= IMP.container.ListPairContainer([(ps[p[0]], ps[p[1]]) for p in pairs],
                                         "Restrained pairs")
    pr= IMP.container.PairsRestraint(score, pc)
    m.set_maximum_score(pr, .01)
    m.add_restraint(pr)
    d= IMP.core.DistanceToSingletonScore(IMP.core.HarmonicUpperBound(2,1),
                                         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])
    m.add_restraint(dr)
    # we are not interested in conformations which don't fit the distances
    # exactly, but using 0 is tricky
    m.set_maximum_score(dr, .01)
    return m.get_restraints()

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):
    pst= IMP.domino.ParticleStatesTable()
    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
    for p in ps[1:]:
        pst.set_particle_states(p, states)
    return pst

def create_sampler(m, pst):
    s=IMP.domino.DominoSampler(m, pst)
    # 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))
    # filter states that score worse than the cutoffs in the Model
    filters.append(IMP.domino.RestraintScoreSubsetFilterTable(m, pst))
    filters[-1].set_log_level(IMP.SILENT)
    states= IMP.domino.BranchAndBoundSubsetStatesTable(pst, filters);
    states.set_log_level(IMP.SILENT);
    s.set_subset_states_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, pst)

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:
        print IMP.core.XYZ(p).get_x(), IMP.core.XYZ(p).get_y()

Generated on Thu Mar 24 2011 02:01:40 for IMP by doxygen 1.7.3