Subject: Re: [IMP-users] CG and MC.set_local_optimizer(CG)
From: Benjamin SCHWARZ <>
Date: Thu, 12 Aug 2010 18:21:46 +0200
Reply-to: Help and discussion for users of IMP <>
Hello Keren,
As you can see in the attached script (lines 168-210) I tested several configurations : FitRestraint alone or with an ExcludedVolumeRestraint; with or without weights (through a RestraintSet as suggested by Ben and Daniel).
#!/usr/bin/python
'''
Created on 22 juil. 2010
@author: schwarz
A first serious attempt to use IMP to model a complex
- a complex is forst loaded (1RDT : 2 domains and 2 peptides)
- a noisy version of the same complex is also loaded
- a map of the complex is sampled at 10A with apix 2.5 A
- the chains are extracted and rigidified
- two restraits are added to the model :
- fit sampled map with all atomic particles
- rigid bodies don't clash
- optimize m
'''
import IMP
import IMP.atom
import IMP.em
import IMP.core
import IMP.container
import IMP.algebra
import math,random
import datetime
#import time
import os
#pdbFilePrefix = "1E1Q"
pdbFilePrefix = "1RDT"
noisyFilePrefix = "1RDT-noisy-swaped-rotated-0"
#inputDataDirPath = "/Users/schwarz/Dev/IMPstuff/src/tests/fit_with_restraints_test/input"
#inputNoisyDataDirPath = "/Users/schwarz/Dev/IMPstuff/src/tests/fit_with_restraints_test/input/noisy"
#outputDataDirPath = "/Users/schwarz/Dev/IMPstuff/src/tests/fit_with_restraints_test/output"
inputDataDirPath = "input"
inputNoisyDataDirPath = "input/noisy"
outputDataDirPath = "output"
pdbFilePath = os.path.join(inputDataDirPath,pdbFilePrefix+".pdb")
noisyFilePath = os.path.join(inputNoisyDataDirPath,noisyFilePrefix+".pdb")
#
# Application parameters
#
params = {
# optimizer control
"optimizer" : {
"type":"MC",
# "type":"MCGC",
# "type":"CG",
"maxSteps": 20
},
"MC" : { "maxTranslation": 5, "maxRotation": 10},
# density map parameters
"densityMap" : { "resolution": 10, "apix": 2.5},
# restraint control
"restraints" : {
"EMFit":1000
,"VolExcl":1
},
"useRestraintWeight":True,
# "useRestraintWeight":False,
# output control
"dumpStructMap" : False
}
def forgeDateFilePrefix():
'''Forges a "unique" file prefix based on the actual date and time
usefull, for instance to identify every output files generated during one experiment of a single script '''
d=datetime.datetime.now()
genericFilePrefix = "FitRestraint--" + "-".join([str(d.year),str(d.month),str(d.day),str(d.hour),str(d.minute)]) + "--"
return genericFilePrefix
class linearMap():
''' an object of this class acts as a functor that linearly maps an interval to an other
'''
def __init__(self,min,max,MIN,MAX):
'''
'''
self._a = (MAX-MIN)/(max-min)
self._b = MAX - self._a * MIN
def map(self,value):
return self._a * value + self._b
class SigmoidRestraint(IMP.Restraint):
def __init__(self,restraint):
self._r=restraint
self.name = "pipo"
def evaluate(self,value):
''' '''
return ( 1 / (1 + math.exp(self._r.evaluate(value)) ) )
def get_is_part_of_model (self): return self._r.get_is_part_of_model()
def get_model (self): return self._r.get_model()
def get_output_objects (self): return self._r.get_output_objects()
def get_output_particles (self): return self._r.get_output_particles()
def set_model (self, model): self._r.set_model(model)
IMPversion = IMP.get_module_version_info().get_version()
if __name__ == "__main__":
# Generation of a unique fileName prefix
genericFilePrefix = forgeDateFilePrefix()
print "*** STARTING experiments "+ genericFilePrefix +" on structure "+ pdbFilePrefix + "and noisy structure " + noisyFilePrefix
print " parameters are :"
for k,v in params.iteritems():
print k+":",v
IMP.set_log_level(IMP.SILENT)
m = IMP.Model()
# load true RX structure and noisy structure
print "** read pdb files "+pdbFilePath+" and "+noisyFilePath
# load and store a structure to work on
sel = IMP.atom.ATOMPDBSelector()
mhref = IMP.atom.read_pdb(pdbFilePath,m,sel)
mh = IMP.atom.read_pdb(noisyFilePath,m,sel)
IMP.atom.add_radii(mhref)
IMP.atom.add_radii(mh)
print "** assessing complex map"
print " - compute map"
resolution,voxelSize = params["densityMap"]["resolution"],params["densityMap"]["apix"]
print " - map for config resolution:"+str(resolution)+" apix:"+str(voxelSize)
dmap = IMP.em.SampledDensityMap(mhref.get_leaves(),resolution,voxelSize)
if params["dumpStructMap"] :
mapName = genericFilePrefix + pdbFilePrefix+"-r"+str(resolution)+"-a"+str(voxelSize)+".mrc"
mapFilePath = os.path.join(outputDataDirPath,mapName)
print " - writing map to "+mapFilePath
IMP.em.write_map(dmap,mapFilePath,IMP.em.MRCReaderWriter())
print "** extract and rigidify protein chains from structure"
chains = {} # Will store chains on which I work
rigidChains = {} # Will store rigidBody for chains on which I work
refChains = {} # Will store reference chains
for chain in IMP.atom.Chains(IMP.atom.get_by_type(mh, IMP.atom.CHAIN_TYPE)) :
rigidChains[chain.get_id()] = IMP.atom.setup_as_rigid_body(chain)
chains[chain.get_id()] = chain
for chain in IMP.atom.Chains(IMP.atom.get_by_type(mhref, IMP.atom.CHAIN_TYPE)) :
refChains[chain.get_id()] = chain
print " (",str(len(refChains))," chains)"
print "** add restraints"
print " - building restraints"
myRestraints = {} # I will store my restraints and weights in this dict, and insert them in the model later
if "EMFit" in params["restraints"].keys():
print " - FitRestraint between map and all particles"
if IMPversion == "1.0" :
fitRestraint = IMP.em.FitRestraint(IMP.Particles(IMP.core.get_leaves(mh)),dmap)
else :
particles = IMP.Particles(rigidChains.values())
fitRestraint = IMP.em.FitRestraint( particles ,
dmap,
IMP.core.LeavesRefiner(chains.values()[0].get_traits()) )
fitRestraint.set_name("fitRestraint")
myRestraints["EMFit"] = fitRestraint
if "VolExcl" in params["restraints"].keys():
print " - Volume exclusion restraint"
if IMPversion == "1.0" :
lsc = IMP.container.ListSingletonContainer()
# lsc = IMP.container.ListSingletonContainer(rigidChains.values())
else :
lsc = IMP.container.ListSingletonContainer(m,"rigidChainsContainer")
for rb in rigidChains.values() :
lsc.add_particle(rb.get_particle())
exVolRestraint = IMP.core.ExcludedVolumeRestraint(
lsc,
IMP.core.LeavesRefiner(chains.values()[0].get_traits()))
# myRestraints[exVolRestraint] = params["restraints"]["VolExcl"]
exVolRestraint.set_name("excludedVolRestraint")
myRestraints["VolExcl"] = exVolRestraint
print " - inserting restraints in model"
if params["useRestraintWeight"] :
weightedRestraintSet = IMP.RestraintSet(1) # my global restraint
for restraintName,restraint in myRestraints.iteritems() :
r = IMP.RestraintSet( params["restraints"][restraintName])
r.add_restraint(restraint)
weightedRestraintSet.add_restraint(r)
weightedRestraintSet.set_name("weightedRestraintSet")
m.add_restraint(weightedRestraintSet)
else :
for restraint in myRestraints.values() :
m.add_restraint(restraint)
print " - initial model evaluation : "+str(m.evaluate(False))
print "** Optimization"
print " - create an optimizer"
#opti = IMP.core.SteepestDescent(m)
if params["optimizer"]["type"] == "CG":
opti = IMP.core.ConjugateGradients(m)
elif params["optimizer"]["type"] in ["MC","MCGC"]:
opti = IMP.core.MonteCarlo(m)
opti.set_return_best(True)
print " - add movers to MC optimizer for all distinct chains"
maxTranslation = params["MC"]["maxTranslation"]
maxRotation = params["MC"]["maxRotation"]
for chain in rigidChains.values():
rb_mover = IMP.core.RigidBodyMover(chain, maxTranslation, maxRotation )
# rb_mover.set_log_level(IMP.WARNING)
rb_mover.set_name("rbMover-chain-"+IMP.atom.Chain(chain.get_particle()).get_id())
opti.add_mover(rb_mover)
if params["optimizer"]["type"] == "MCGC":
cglo = IMP.core.ConjugateGradients(m)
opti.set_local_optimizer( cglo )
# opti.set_log_level(IMP.TERSE)
opti.set_name("optimizer")
opti.set_log_level(IMP.WARNING)
print " - optimize"
opti.optimize(params["optimizer"]["maxSteps"])
print "** Analysis"
print " - optimized model evaluation : "
print " - global score is "+str(m.evaluate(False))
for restraintName,restraint in myRestraints.iteritems() :
# m.add_restraint(restraint)
print " - "+restraintName+" score is "+str(m.evaluate(IMP.RestraintsTemp([restraint]),False))
print " - save optimized model"
optimizedModelFilePath = os.path.join(outputDataDirPath,genericFilePrefix + pdbFilePrefix+"-optimized.pdb")
IMP.atom.write_pdb(mh,optimizedModelFilePath)
print " - measure RMSD discrepency between native and optimized chains"
for chainID in refChains.keys() :
rmsd=IMP.atom.get_rmsd(
IMP.core.XYZs( chains[chainID].get_leaves() ) ,
IMP.core.XYZs( refChains[chainID].get_leaves() )
)
print " "+chainID+" : "+str(rmsd)
print " - save optimized complex map"
dmapOptimized = IMP.em.SampledDensityMap(mh.get_leaves(),resolution,voxelSize)
mapName = genericFilePrefix + pdbFilePrefix+"-optimized-r"+str(resolution)+"-a"+str(voxelSize)+".mrc"
mapFilePath = os.path.join(outputDataDirPath,mapName)
IMP.em.write_map(dmapOptimized,mapFilePath,IMP.em.MRCReaderWriter())
Le 12 août 2010 à 17:40, Keren Lasker a écrit :
> Ben -
> from the warning messages I see that you are using the FitRestraint.
> Are you using any other restraints as well?
> It will be useful to get your optimization script so we can take a look.
> Keren.
> On Aug 12, 2010, at 7:00 AM, Benjamin SCHWARZ wrote: