[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [IMP-users] CG and MC.set_local_optimizer(CG)



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: