from phaser.pipeline import bba from phaser.pipeline import mr_object_new as mr_object from phaser.pipeline import process_model_new as process_model from phaser.pipeline import mr_calculation from phaser.pipeline.proxies import OCAHTTP, DbfetchHTTP, ClustalWSynchronous, SuperposePdbs from phaser import chisel_new as chisel from phaser import sculptor from phaser import tbx_utils from phaser import rsam from phaser import tbx_utils from phaser import ncs from iotbx import bioinformatics import iotbx.pdb from libtbx.utils import Sorry import scitbx.matrix import os import sys import time import operator import math import fractions import tempfile from StringIO import StringIO class Controller(object): """ Coordinates knowledge sources """ def __init__(self, pool_class, max_jobs = 1): self.pool = pool_class( processes = max_jobs ) self.max_jobs = max_jobs self.jobs = set() def run(self, blackboard, experts): while True: action = False if len( self.jobs ) < self.max_jobs: blackboard.debug( src = self, msg = "available cpus - querying knowledge sources" ) for ks in experts: blackboard.debug( src = ks, msg = "evaluation phase" ) candidate = ks.inspect( blackboard = blackboard ) blackboard.debug( src = ks, msg = "evaluation finished" ) if candidate: blackboard.info( src = ks, msg = "action phase" ) blackboard.info( src = ks, msg = "candidate %s" % candidate ) ks.manipulate( context = self, candidate = candidate, blackboard = blackboard, ) blackboard.info( src = ks, msg = "action complete" ) action = True break else: blackboard.info( src = self, msg = "no knowledge sources can act" ) finished_jobs = self.finished_jobs( blackboard = blackboard ) for job in finished_jobs: blackboard.info( src = self, msg = "processing finished job" ) job.owner.process( calculation = job.calculation, initial = job.initial, results = job.results, blackboard = blackboard ) blackboard.info( src = self, msg = "job processing finished" ) self.jobs.remove( job ) if not action and not finished_jobs: if self.jobs: blackboard.debug( src = self, msg = "suspending to wait for jobs to finish..." ) time.sleep( 1 ) else: blackboard.debug( src = self, msg = "exiting as no further action possible" ) break def submit(self, owner, initial, calculation): job = Job( owner = owner, initial = initial, calculation = calculation, results = self.pool.apply_async( calculation.run ) ) self.jobs.add( job ) def finished_jobs(self): return [ j for j in self.jobs if j.results.ready() ] def __str__(self): return "Controller" class PhaserAI(object): """ Coordinates MR process """ def __init__(self, solvers, controller): self.initial_threshold = 0.75 self.solvers = solvers self.refinement_threshold = 0.75 self.refine = RefinementExpert() self.commit_threshold = 0.75 self.controller = controller def solve(self, blackboard): blackboard.logger.info( "AI: starting solution process..." ) while True: self.mr_round( blackboard = blackboard ) if not blackboard.read( type = stage.Refined ): break blackboard.info( src = self, msg = "solution process finished" ) blackboard.info( src = self, msg = "evaluating results" ) best = blackboard.problem.score.value blackboard.info( src = self, msg = "best score %.3f" % best ) initials = sorted( blackboard.read( type = stage.Initial ), key = lambda i: i.score(), reverse = True ) processed = set() potentially_smaller = False potentially_larger = len blackboard.info( "%s peaks above threshold" % len( initials ) ) for ( index, current ) in enumerate( initials ): blackboard.info( src = self, msg = "peak %s" % current ) score = current.score() blackboard.info( src = self, msg = "score %.3f" % score ) ancestors = current.partial.ancestors() missing_composition_for = dict( [ ( c.initial(), c.missing ) for c in blackboard.read( type = stage.Composition ) ] ) best_scorer = [ i for i in initials if best == i.score() ] blackboard.info( src = self, msg = "best scorer is %s" % best_scorer ) if missing_composition_for[ best_scorer ].is_empty(): blackboard.info( src = self, msg = "corresponds to target composition" ) else: blackboard.info( src = self, msg = "missing %s" % missing_composition_for[ best_scorer ] ) most_complete = min( initials, key = lambda i: missing_composition_for[ i ].mw() ) if most_complete == best_scorer: blackboard.info( src = self, msg = "this is also the most complete solution" ) else: blackboard.info( src = self, msg = "most complete is %s" % most_complete ) for current in initials: blackboard.info( src = self, msg = "peak %s" % current ) score = current.score() blackboard.info( src = self, msg = "score %.3f" % score ) if best == score: blackboard.info( ) compositions = [ c for c in blackboard.read( type = stage.Composition ) if c.initial() == initial ] assert len( compositions ) == 1 if compositions[0].missing.is_empty(): blackboard.info( src = self, msg = "corresponds to specified composition" ) continue self._logger.info( "%s were outperformed and rejected" % len( self._bad_cases ) ) self._logger.info( "%s has yielded complete solutions" % len( self._solved_cases ) ) threshold = self._ps_select.threshold_score( cases = self._finished_cases ) for c in self._finished_cases: self._logger.info( str( c ) ) partials = self._ps_select.promising_partial_structures( case = c, threshold = threshold ) self._logger.info( "\tThere are %s partials above the threshold" % len( partials ) ) self._logger.info( "\tBest score: %s" % c.statistics().best_score() ) def mr_round(self, blackboard): blackboard.info( src = self, msg = "starting MR round..." ) blackboard.info( src = self, msg = "clearing blackboard..." ) for obj in blackboard.read( type = stage.Stage ): blackboard.debug( src = self, msg = "deleting %s" % obj ) blackboard.data.remove( obj ) # Select significant peaks self.initial_select( blackboard = blackboard ) # Do MR self.controller.run( blackboard = blackboard, experts = self.solvers ) # Assess peaks self.refinement_select( blackboard = blackboard ) # Refine peaks self.controller.run( blackboard = blackboard, experts = [ self.refine ] ) self.record_search_coverage( blackboard = blackboard ) # Load refined stages back to case refs = [ st for st in blackboard.read( type = stage.Refined ) ] blackboard.info( src = self, msg = "%s peaks have been refined" % len( refs ) ) for r in refs: equiv = r.case.equivalents_for( structure = r.structure ) if not equiv: r.case.insert( structure = r.structure, score = r.llg ) else: blackboard.info( src = self, msg = "%s: not unique" % r ) for e in equiv: e.predecessors.add( r.initial().partial ) def initial_select(self, blackboard): cases = blackboard.read( type = mr_object.Case ) if not cases: blackboard.info( src = self, msg = "no active cases available" ) return blackboard.info( src = self, msg = "%s cases active" % len( cases ) ) best = blackboard.problem.score.value threshold = best - ( 1.0 - self.initial_threshold ) * abs( best ) blackboard.info( src = self, msg = "threshold: %.3f" % threshold ) for case in cases: blackboard.info( src = self, msg = "case %s" ) promising = [ p for p in case.partials if threshold <= case.score_for( partial = p ) ] if not promising: blackboard.info( "No above threshold peaks found" ) continue blackboard.info( src = self, msg = "%d found" % len( promising ) ) selected = [ promising[0] ] ancestors = promising[0].ancestors() for structure in promising[1:]: if structure not in ancestors: selected.append( structure ) ancestors.update( structure.ancestors ) blackboard.info( src = self, msg = "%d after removing ancestors with lower score" % len( selected ) ) for structure in selected: st = stage.Initial( case = case, partial = structure ) blackboard.data.add( st ) def refinement_select(self, blackboard): candidates = blackboard.read( type = stage.Packing ) if not stages: blackboard.info( src = self, msg = "no refinement candidates" ) return blackboard.info( src = self, msg = "%s refinement candidates" % len( candidates ) ) best = max( [ p.score() for p in stages ] ) threshold = best - ( 1.0 - self.refinement_threshold ) * abs( best ) blackboard.info( src = self, msg = "threshold: %.3f" % threshold ) for st in stages: if threshold <= st.score(): blackboard.info( src = self, msg = "Selected %s" % st ) ref = stage.Significant( parent = st ) blackboard.data.add( ref ) def record_search_coverage(self, blackboard): # Check for failed calculations if blackboard.read( type = stage.Failure ): blackboard.info( src = self, msg = "Calculations have failed" ) blackboard.info( src = self, msg = "Coverage calculations will be incorrect" ) # Assume that each translation peak has been fully processed # (i.e. packing and refinement) # Create records for new stages for t in blackboard.read( type = stage.Translation ): if t.processed: blackboard.info( src = self, msg = "%s has been processed" % t ) rec = record.TranslationPeak( peak = t.peak ) blackboard.data.add( rec ) # Check whether complete TFs have been fully processed processed_in_trafun = {} for t in blackboard.read( type = record.TranslationPeak ): processed_in.setdefault( t.peak.function, set() ).add( t.peak ) # Replace individual translation peak records with one rotation peak record # if trafun is fully processed for ( trafun, processed ) in processed_in_trafun.items(): if set( trafun.peaks ) == processed: blackboard.info( src = self, msg = "All peaks in %s has been processed" % trafun ) rec = record.RotationPeak( peak = trafun.calculation.rpeak ) blackboard.data.add( rec ) for r in blackboard.read( type = record.TranslationPeak ): if r.peak.function == trafun: blackboard.data.remove( r ) # Check whether complete RFs have been fully processed processed_in_rotfun = {} for r in blackboard.read( type = record.RotationPeak ): processed_in.setdefault( r.peak.function, set() ).add( r.peak ) # Replace individual rotation peak records with one ensemble record # if rotfun is fully processed for ( rotfun, processed ) in processed_in_rotfun.items(): if set( rotfun.peaks ) == processed: blackboard.info( src = self, msg = "All peaks in %s has been processed" % rotfun ) rec = record.Ensemble( case = rotfun.calculation.case, partial = rotfun.calculation.partial, ensemble = rotfun.calculation.ensemble ) blackboard.data.add( rec ) for r in blackboard.read( type = record.RotationPeak ): if r.peak.function == rotfun: blackboard.data.remove( r ) def __str__(self): return "PhaserAI" @classmethod def MRMode(cls, problem, composition, symmetries, outstream = None): solvers = [ ScoringExpert(), TemplateSolveSuperposeExpert(), TemplateSolveInitiateExpert(), EvaluationExpert(), PackingExpert(), TranslationSearchExpert(), RotationSearchExpert(), EnsembleApplicabilityExpert(), ] ai = cls( problem = problem, composition = composition, symmetries = symmetries, solvers = solvers, modellers = [], outstream = outstream ) return ai # Knowledge sources - inspect whether they are active class CalculateMissing(object): """ Calculate missing composition """ def inspect(self, blackboard): candidates = [ c for c in blackboard.read( type = stage.Initial ) if not c.processed ] if not candidates: blackboard.debug( src = self, msg = "no candidates found" ) selected = None else: blackboard.debug( src = self, msg = "%s candidates found" % len( candidates ) ) selected = max( candidates, key = lambda stage: stage.score() ) blackboard.debug( src = self, msg = "selected %s" % selected ) return selected def manipulate(self, context, candidate, blackboard): candidate.processed = True missing = candidate.case.composition.copy() for e in candidate.partial.ensembles(): missing.subtract( e.component ) blackboard.debug( src = self, msg = "calculated missing part" ) st = stage.Composition( parent = candidate, missing = missing ) blackboard.data.add( st ) class FindEnsemble(object): """ Find ensembles for missing part """ def inspect(self, blackboard): candidates = [ c for c in blackboard.read( type = stage.Composition ) if not c.processed ] if not candidates: blackboard.debug( src = self, msg = "no candidates found" ) selected = None else: blackboard.debug( src = self, msg = "%s candidates found" % len( candidates ) ) selected = max( candidates, key = lambda stage: stage.score() ) blackboard.debug( src = self, msg = "selected %s" % selected ) return selected def manipulate(self, context, candidate, blackboard): candidate.processed = True records = set( [ ( r.case, r.partial, r.ensemble ) for r in blackboard.read( type = record.Ensemble ) ] ) initial = candidate.initial() for ensemble in blackboard.read( type = process_model.Ensemble ): blackboard.info( src = self, msg = "Testing %s" % ensemble ) if ( initial.case, initial.partial, ensemble ) in records: blackboard.info( src = self, msg = "previously explored" ) continue if not missing.contains( other = ensemble.component ): blackboard.info( src = self, msg = "not in composition" ) continue blackboard.info( src = self, msg = "applicable for extension" ) st = stage.Ensemble( parent = candidate, ensemble = ensemble ) stage.attach_to( blackboard = blackboard ) class RotationSearch(object): """ Performs rotation search for ensemble stages """ def inspect(self, blackboard): candidates = [ c for c in blackboard.read( type = stage.Ensemble ) if not c.processed ] if not candidates: blackboard.debug( src = self, msg = "no candidates found" ) selected = None else: blackboard.debug( src = self, msg = "%s candidates found" % len( candidates ) ) selected = max( candidates, key = lambda stage: stage.score() ) blackboard.debug( src = self, msg = "selected %s" % selected ) return selected def manipulate(self, context, candidate, blackboard): candidate.processed = True blackboard.info( src = self, msg = "Starting rotation search" ) calc = calculation.RotationSearch( case = candidate.case, partial = candidate.partial, ensemble = candidate.ensemble ) blackboard.data.add( calc ) context.submit( owner = self, initial = candidate, calculation = calc ) def process(self, initial, calculation, results, blackboard): st = calculation.process( results = result ) blackboard.data.remove( calculation ) blackboard.data.add( st ) st.announce( initial = initial, blackboard = blackboard ) class TranslationSearch(object): """ Performs translation search for rotation stages """ def inspect(self, blackboard): candidates = [ c for c in blackboard.read( type = stage.Rotation ) if not c.processed ] if not candidates: blackboard.debug( src = self, msg = "no candidates found" ) selected = None else: blackboard.debug( src = self, msg = "%s candidates found" % len( candidates ) ) selected = max( candidates, key = lambda stage: stage.score() ) blackboard.debug( src = self, msg = "selected %s" % selected ) return selected def manipulate(self, context, candidate, blackboard): candidate.processed = True blackboard.info( src = self, msg = "Starting translation search" ) calc = calculation.TranslationSearch( case = candidate.case, partial = candidate.partial, rpeak = candidate.peak ) blackboard.data.add( calc ) context.submit( owner = self, initial = candidate, calculation = calc ) def process(self, initial, calculation, results, blackboard): st = calculation.process( results = result ) blackboard.data.remove( calculation ) blackboard.data.add( st ) st.announce( initial = initial, blackboard = blackboard ) class PackingCheck(object): """ Performs packing check for translation stages """ def inspect(self, blackboard): candidates = [ c for c in blackboard.read( type = stage.Translation ) if not c.processed ] if not candidates: blackboard.debug( src = self, msg = "no candidates found" ) selected = None else: blackboard.debug( src = self, msg = "%s candidates found" % len( candidates ) ) selected = max( candidates, key = lambda stage: stage.score() ) blackboard.debug( src = self, msg = "selected %s" % selected ) return selected def manipulate(self, context, candidate, blackboard): candidate.processed = True blackboard.info( src = self, msg = "Starting packing check" ) calc = stage.PackingCalculation( parent = candidate ) blackboard.data.add( calc ) context.submit( owner = self, calculation = calc ) def process(self, calculation, results, blackboard): st = calculation.process( results = result ) blackboard.data.remove( calculation ) blackboard.data.add( st ) st.announce( initial = initial, blackboard = blackboard ) class Refinement(object): """ Refines peaks """ def inspect(self, blackboard): candidates = [ c for c in blackboard.read( type = stage.Significant ) if not c.processed ] if not candidates: blackboard.debug( src = self, msg = "no candidates found" ) selected = None else: blackboard.debug( src = self, msg = "%s candidates found" % len( candidates ) ) selected = max( candidates, key = lambda stage: stage.score() ) blackboard.debug( src = self, msg = "selected %s" % selected ) return selected def manipulate(self, context, candidate, blackboard): candidate.processed = True blackboard.info( src = self, msg = "Starting refinement" ) calc = stage.RefinementCalculation( parent = candidate ) blackboard.data.add( calc ) context.submit( owner = self, calculation = calc ) def process(self, calculation, results, blackboard): st = calculation.process( results = result ) blackboard.data.remove( calculation ) blackboard.data.add( st ) ### Old class AI(object): """ MR solution controller """ def __init__(self): self.cpus = 1 self.experts =[] self._queue = [] self._running = set() self._engine = Sequential def solve(self, context): while True: suspend = True if 0 < self.cpus: for ks in self.experts: candidate = ks.inspect( context = context ) if candidate: context.logger.info( "%s - action phase" % ks ) context.logger.info( "Candidate: %s" % candidate ) ks.manipulate( context = context, candidate = candidate ) context.logger.info( "Action complete\n" ) suspend = False break else: if not self._running: break new_results = self.poll( context = context ) if suspend and not new_results: time.sleep( 1 ) def submit(self, owner, protocol, execution, processing): job = self._engine( owner = owner, protocol = protocol, execution = execution, processing = processing, ) self._queue.append( job ) self.spool() def poll(self, context): new_results = False for job in list( self._jobs ): if not job.finished(): continue logger.info( "Poll - job finished" ) new_results = True logger.debug( job.logfile ) if job.success: self.logger.info( "%s - processing successful job results" % job.owner ) job.owner.process_success( results = job.results, context = context, extra = job.processing ) else: logger.info( "%s - processing failed job results" % job.owner ) job.owner.process_failure( context = context, extra = job.processing ) self._jobs.remove( job ) self.cpus += 1 logger.info( "Processing complete\n" ) return new_results def spool(self): while 0 < self.cpus and self._queue: job = self._queue.pop( 0 ) job.run() self.cpus -= 1 self._running.add( job ) class Sequential(object): """ Simple engine that executes jobs on the main thread """ def __init__(self, owner, protocol, execution, processing): self.owner = owner self.processing = processing self.protocol = protocol self.processing = processing self._finished = False self.success = None self.logfile = "" self.results = None def run(self): ( self.success, self.logfile, self.results ) = self.protocol( data = self.execution ) self._finished = True def poll(self): pass def finished(self): return self._finished class PhaserAI(object): """ AI object for MR """ MESSAGES = { "ALL DONE": ( "All MRCases have been explored and decisively evaluated. Finishing" ), "REPORT - CAPTION": "There were altogether %s MRCases considered", "EVAL - FULL SOLUTION": ( "\t\tThis partial structure is complete wrt current composition and models" ), "EVAL - EXHAUSTED": ( "\tAll current models have been explored with this partial structure" ), "EVAL - CASE SOLVED": ( "\tAll selected partial structures are complete. This case can be considered solved" ), "EVAL - CASE EXHAUSTED": ( """\tAll possible extensions have been tried, but either not applicable or yield \tworse results than intermediate solutions. This may indicate that the solvent \tcontent is higher than expected and there are less molecules to find""" ), } def __init__( self, problem, composition, symmetries, solvers, modellers, outstream = None, ): if not os.path.exists( "pipeline" ): os.mkdir( "pipeline" ) self._bb = bba.Blackboard( directory = os.path.abspath( "pipeline" ) ) self._problem = problem for symm in symmetries: case = process_model.Case( composition = composition, problem = problem, space_group = symm ) case.attach_to( blackboard = self._bb ) self._ps_select = PSSelect( fraction = 0.75 ) self._solvers = solvers self._modellers = modellers self._tp_select = TPSelect( fraction = 0.75 ) self._refine = RefinementExpert() self._controller = Controller() self._logger = logging.getLogger( name = "PhaserAI.%s" % id( self ) ) self._logger.setLevel( logging.INFO ) if outstream: handler = logging.StreamHandler( outstream ) else: handler = tbx_utils.NullHandler() handler.setFormatter( logging.Formatter( "%(message)s" ) ) self._logger.addHandler( handler ) self._solved_cases = set() self._bad_cases = set() self._finished_cases = set() self._all_models_present = False def add(self, bbobj): bbobj.attach_to( blackboard = self._bb ) def logger(self): return self._logger def solve(self): self._logger.info( "Starting solution process..." ) while True: self.mr_round() if not self._bb.read_instances_of( type = process_model.Case ): self._logger.info( self.MESSAGES[ "ALL DONE" ] ) break self._logger.info( "Solution process finished" ) self._logger.info( self.MESSAGES[ "REPORT - CAPTION" ] % ( len( self._bad_cases ) + len( self._finished_cases ) ) ) self._logger.info( "%s were outperformed and rejected" % len( self._bad_cases ) ) self._logger.info( "%s has yielded complete solutions" % len( self._solved_cases ) ) threshold = self._ps_select.threshold_score( cases = self._finished_cases ) for c in self._finished_cases: self._logger.info( str( c ) ) partials = self._ps_select.promising_partial_structures( case = c, threshold = threshold ) self._logger.info( "\tThere are %s partials above the threshold" % len( partials ) ) self._logger.info( "\tBest score: %s" % c.statistics().best_score() ) def mr_round(self): self._logger.info( "Starting MR round..." ) self._logger.debug( "Clearing blackboard..." ) for obj in self._bb.read_instances_of( type = process_model.RoundObject ): obj.detach_from_blackboard() # Select significant peaks self._logger.info( "%s - analysis starts" % self._ps_select ) self._ps_select.analyse( blackboard = self._bb, logger = self._logger ) self._logger.info( "Analysis finished\n" ) # Do MR self._controller.run( blackboard = self._bb, experts = self._solvers, logger = self._logger ) # Assess translation peaks self._logger.info( "%s - analysis starts" % self._tp_select ) self._tp_select.analyse( blackboard = self._bb, logger = self._logger ) self._logger.info( "Analysis finished\n" ) # Evaluate self._logger.info( "Evaluating results..." ) if not self._all_models_present: self._all_models_present = not self._controller.any_active( experts = self._modellers, blackboard = self._bb ) if self._all_models_present: self._logger.info( "All basic models are now present" ) cases = self._bb.read_instances_of( type = process_model.Case ) self._logger.info( "There are still %s MRCases active" % len( cases ) ) exhaustions = set( [ rec.identifier() for rec in self._bb.read_instances_of( type = process_model.ExhaustionRecord ) ] ) for case in cases: self._logger.info( "Trial case: %s" % case ) initial_stages = [ stage for stage in self._bb.read_instances_of( type = process_model.InitialStage ) if stage.case() == case ] rejected_stages = [ stage for stage in self._bb.read_instances_of( type = process_model.RejectedStage ) if stage.case() == case ] descendant_stages = [ stage for stage in self._bb.read_instances_of( type = process_model.DescendantNode ) if stage.case() == case ] if rejected_stages: self._logger.info( "\tPartial structures have been selected that had already been tested" ) self._logger.info( "\tThese were not pursued any further" ) if not initial_stages: self._logger.info( "\tThis case cannot be advanced any further" ) if rejected_stages: self._logger.info( "\tAll selected partial structures have been tested" ) self._finished_cases.add( case ) else: self._logger.info( "\tNo partial structure was above threshold" ) self._bad_cases.add( case ) case.detach_from_blackboard() continue self._logger.info( "\t%s initial partial structure has been identified" % len( initial_stages ) ) terminal_stages = [] for ( index, ini ) in enumerate( initial_stages ): self._logger.info( "\tPartial structure %s - Statistics" % ( index + 1 ) ) self._logger.info( "\t\t%s ensembles considered" % len( ini.ensembles() ) ) self._logger.info( "\t\t%s tested for suitability" % ( len( ini.ensembles() ) - len( ini.untested() ), ) ) self._logger.info( "\t\t%s remain untested" % len( ini.untested() ) ) ensemble_stages = set( [ st for st in descendant_stages if ( st.parent() == ini and isinstance( st, process_model.EnsembleStage ) ) ] ) immediate_successor_stages = [ st for st in descendant_stages if st.parent() in ensemble_stages ] all_successor_stages = [ st for st in descendant_stages if ( st.predecessor( type = process_model.EnsembleStage ) in ensemble_stages ) ] self._logger.info( "\t\t%s calculations have been done from %s tested ensembles" % ( len( immediate_successor_stages ), len( ensemble_stages ), ) ) self._logger.info( "\t\t%s calculations altogether" % len( all_successor_stages ) ) if not ini.untested() and not ensemble_stages: idents = set( [ process_model.ExhaustionRecord.descriptor( case = ini.case(), partial = ini.partial(), ensemble = e ) for e in ini.ensembles() ] ) if self._all_models_present and exhaustions.issuperset( idents ): self._logger.info( self.MESSAGES[ "EVAL - FULL SOLUTION" ] ) self._solved_cases.add( case ) terminal_stages.append( ini ) if len( terminal_stages ) == len( initial_stages ): self._logger.info( self.MESSAGES[ "EVAL - CASE SOLVED" ] ) self._finished_cases.add( case ) case.detach_from_blackboard() continue # Refine peaks self._controller.run( blackboard = self._bb, experts = [ self._refine ], logger = self._logger ) # Load refined stages back to case refs = [ stage for stage in self._bb.read_instances_of( type = process_model.RefinedStage ) ] self._logger.info( "\t%s refined peaks have been found" % len( refs ) ) for r in refs: r.case().insert( partial = r.final().partial(), score = r.final().llg() ) def record_extensions(self, blackboard): for stage in blackboard.read_instances_of( type = process_model.EnsembleStage ): rot_search_exh = [ s for s in stage.searches if s.exhaustive ] if not rot_search_exh: continue rot_st_all = set( reduce( operator.add, [ s.rotations for s in rot_search_exh ] ) ) rot_st_act = [ rs for rs in stage.children if rs.rotation in rot_st_all ] assert rot_st_act for rs in rot_st_act: tra_search_exh = [ s for s in rs.searches if s.exhaustive ] if not tra_search_exh: break tra_st_all = set( reduce( operator.add, [ s.translations for s in tra_search_exh ] ) ) tra_st_act = [ ts for ts in rs.children if ts.translation in tra_st_all ] assert tra_st_act leaves = [ [ n for n in ts.descendants() if not n.children ] for ts in tra_st_act ] assert all( [ len( l ) == 1 for l in leaves ] ) if not all( [ isinstance( l[0], self.LEAVES ) for l in leaves ] ): break else: rec = process_model.ExtensionRecord( case = stage.case, partial = stage.partial, ensemble = stage.ensemble ) rec.attach_to( blackboard = blackboard ) @classmethod def MRMode(cls, problem, composition, symmetries, outstream = None): solvers = [ ScoringExpert(), TemplateSolveSuperposeExpert(), TemplateSolveInitiateExpert(), EvaluationExpert(), PackingExpert(), TranslationSearchExpert(), RotationSearchExpert(), EnsembleApplicabilityExpert(), ] ai = cls( problem = problem, composition = composition, symmetries = symmetries, solvers = solvers, modellers = [], outstream = outstream ) return ai class Controller(object): """ Coordinates knowledge sources """ def __init__(self, max_jobs = 1): self._jobs = set() self._max_jobs = 1 self._logger = logging.getLogger( name = "Controller.%s" % id( self ) ) self._logger.addHandler( tbx_utils.NullHandler() ) def run(self, blackboard, experts, logger): while True: suspend = True if self.available_cpus(): for ks in experts: candidate = ks.inspect( blackboard = blackboard, logger = logger ) if candidate: logger.info( "%s - action phase" % ks ) logger.info( "Candidate: %s" % candidate ) ks.manipulate( context = self, candidate = candidate, blackboard = blackboard, logger = logger, ) logger.info( "Action complete\n" ) suspend = False break else: if not self.running_calculations(): break new_results = self.poll( blackboard = blackboard, logger = logger ) if suspend and not new_results: time.sleep( 1 ) def available_cpus(self): return len( self._jobs ) < self._max_jobs def submit(self, owner, protocol, execution, processing): assert len( self._jobs ) < self._max_jobs job = Engine( owner = owner, protocol = protocol, execution = execution, processing = processing, ) assert job not in self._jobs self._jobs.add( job ) def running_calculations(self): return self._jobs def poll(self, blackboard, logger): new_results = False for job in list( self._jobs ): if not job.finished(): continue logger.info( "Poll - job finished" ) new_results = True logger.debug( job.logfile ) if job.success: logger.info( "%s - processing successful job results" % job.owner ) job.owner.process_success( results = job.results, blackboard = blackboard, logger = logger, **job.processing ) else: logger.info( "%s - processing failed job results" % job.owner ) job.owner.process_failure( blackboard = blackboard, logger = logger, **job.processing ) self._jobs.remove( job ) logger.info( "Processing complete\n" ) return new_results def any_active(self, experts, blackboard): return any( ks.inspect( blackboard = blackboard, logger = self._logger ) for ks in experts ) @classmethod def Autofetch(cls): return cls( knowledge_sources = [ SolutionIdentify(), TranslationSearchExpert(), RotationSearchExpert(), EnsembleApplicabilityExpert(), SculptorExpert(), AlignmentExpert(), StructureFetchExpert(), SequenceFetchExpert(), ] ) class PSSelect(object): """ Initiates the MR process: 1. Thresholds peaks 2. Identifies ensembles that have not been tried before """ def __init__(self, fraction): self._fraction = fraction def analyse(self, blackboard, logger): cases = blackboard.read_instances_of( type = process_model.Case ) if not cases: logger.info( "No active cases available" ) return logger.info( "There are %s cases active" % len( cases ) ) threshold = self.threshold_score( cases = cases ) logger.info( "Selection threshold for partial structures: %s" % threshold ) ensembles = blackboard.read_instances_of( type = process_model.Ensemble ) logger.info( "There are %s ensembles available" % len( ensembles ) ) past_selections = set( [ rec.identifier() for rec in blackboard.read_instances_of( type = process_model.SelectionRecord ) ] ) for case in cases: logger.info( "Processing %s" % case ) promising = self.promising_partial_structures( case = case, threshold = threshold ) logger.info( "There are %s partial structures above threshold" % len( promising ) ) for p in promising: logger.info( "Processing %s" % p ) rec = process_model.SelectionRecord.descriptor( case = case, partial = p ) if rec in past_selections: logger.info( "This partial has already been tested" ) stage = process_model.RejectedStage( case = case, partial = p ) else: logger.info( "Identified as potential starting point" ) stage = process_model.InitialStage( case = case, partial = p, ensembles = ensembles ) stage.attach_to( blackboard = blackboard ) def threshold_score(self, cases): best = max( [ c.statistics().best_score() for c in cases] ) return best - ( 1.0 - self._fraction ) * abs( best ) def promising_partial_structures(self, case, threshold): return [ p for p in case.partials() if threshold <= case.score_for( partial = p ) ] def __str__(self): return "PartialStructureSelect" class TPSelect(object): """ Selects above threshold peaks """ def __init__(self, fraction): self._fraction = fraction def analyse(self, blackboard, logger): stages = blackboard.read_instances_of( type = process_model.EvaluationStage ) if not stages: logger.info( "No EvaluationStages available" ) return logger.info( "There are %s EvaluationStages active" % len( stages ) ) best = max( [ p.peak().tf() for p in stages ] ) threshold = best - ( 1.0 - self._fraction ) * abs( best ) logger.info( "Selection threshold for refinement: %s" % threshold ) for st in stages: if st.is_solution() or threshold <= st.peak().tf(): logger.info( "Selected %s" % st ) ref_stage = process_model.RefinementStage( parent = st ) ref_stage.attach_to( blackboard = blackboard ) def __str__(self): return "MRPeakEvaluate" # Knowledge sources - inspect and manipulate class KnowledgeSource(object): """ Knowledge source: 1. functionality to disable if solutions have been found """ def __init__(self, selector, exhaustive): self._exhaustive = exhaustive self._selector = selector def inspect(self, blackboard, logger): logger.debug( "%s - evaluation phase" % self ) if self._exhaustive: candidates = self.fetch( bb = blackboard ) else: solved_cases = set( [ es.case() for es in self.solutions( bb = blackboard ) ] ) logger.debug( "Omitting stages belonging to %s solved cases" % len( solved_cases ) ) candidates = [ st for st in self.fetch( bb = blackboard ) if st.case() not in solved_cases ] if not candidates: logger.debug( "No candidates have been found" ) selected = None else: logger.debug( "%s candidates has been found" % len( candidates ) ) selected = self._selector( choices = candidates ) logger.debug( "Selected %s" % selected ) logger.debug( "Evaluation finished\n" ) return selected def fetch(self, bb): return [ obj for obj in bb.read_instances_of( type = self.BBOBJ ) if self.acceptable( bbobj = obj ) ] def solutions(self, bb): return [ st for st in bb.read_instances_of( type = process_model.EvaluationStage ) if st.is_solution() ] def __str__(self): return self.__class__.__name__ class EnsembleApplicabilityExpert(KnowledgeSource): """ Analyse available ensembles for each trial structure """ BBOBJ = process_model.InitialStage def __init__(self): super( EnsembleApplicabilityExpert, self ).__init__( selector = SelectBest( key = lambda s: s.partial_solution_score() ), exhaustive = False ) def manipulate(self, context, candidate, blackboard, logger): # Process all ensembles to enable sorting on a criterion missing = candidate.missing() for ensemble in list( candidate.untested() ): logger.info( "Testing %s" % ensemble ) if missing.contains( other = ensemble.composition() ): logger.info( "It is applicable for extension" ) stage = process_model.EnsembleStage( parent = candidate, search = ensemble ) stage.attach_to( blackboard = blackboard ) else: logger.info( "It does not fit into composition" ) record = process_model.ExhaustionRecord( case = candidate.case(), partial = candidate.partial(), ensemble = ensemble ) record.attach_to( blackboard = blackboard ) candidate.remove( ensemble = ensemble ) logger.info( "Adding processing record for the candidate" ) rec = process_model.SelectionRecord( case = candidate.case(), partial = candidate.partial() ) rec.attach_to( blackboard = blackboard ) def acceptable(self, bbobj): return bool( bbobj.untested() ) class RotationSearchExpert(KnowledgeSource): """ Performs rotation search for ensemble stages """ BBOBJ = process_model.EnsembleStage def __init__(self): super( RotationSearchExpert, self ).__init__( selector = SelectBest( key = lambda s: s.search().score() ), exhaustive = False ) self._fraction = 0.0 def manipulate(self, context, candidate, blackboard, logger): candidate.processed_by.add( self ) solved_ensemble_stages = set( [ st.predecessor( type = process_model.EnsembleStage ) for st in blackboard.read_instances_of( type = process_model.EvaluationStage ) if st.is_solution() ] ) if candidate in solved_ensemble_stages: logger.info( "Skipping search as a significant peak has already been found" ) return search_for = dict( [ ( s.identifier(), s ) for s in blackboard.read_instances_of( type = process_model.RotationSearch ) if s.expert == self ] ) search_id = process_model.RotationSearch.descriptor( bbobj = candidate ) if search_id in search_for: logger.info( "An equivalent search has been performed with a different case" ) logger.info( "Copying results" ) search = process_model.RotationSearch( expert = self, stage = candidate, peaks = search_for[ search_id ].peaks, exhaustive = True ) logger.info( "Selecting peaks from rotation function" ) self.evaluate_search_results( search = search, parent = candidate, blackboard = blackboard, logger = logger ) else: logger.info( "Starting rotation search" ) context.submit( owner = self, protocol = mr_calculation.rotation_function, execution = { "case": candidate.case(), "partial": candidate.partial(), "search": candidate.search(), }, processing = { "parent": candidate }, ) def process_success(self, parent, results, blackboard, logger): search = process_model.RotationSearch.from_definitions( expert = self, stage = parent, definitions = results, exhaustive = True ) logger.info( "There are %s peaks in the rotation function" % len( search.peaks ) ) self.evaluate_search_results( search = search, parent = parent, blackboard = blackboard, ) def process_failure(self, parent, blackboard, logger): pass def evaluate_search_results(self, search, parent, blackboard): search.attach_to( blackboard = blackboard ) stages = [ process_model.RotationStage( parent = parent, peak = p ) for p in search.peaks ] for s in stages: s.attach_to( blackboard = blackboard ) def acceptable(self, bbobj): return self not in bbobj.processed_by class TranslationSearchExpert(KnowledgeSource): """ Performs translation search for rotation stages """ BBOBJ = process_model.RotationStage def __init__(self): super( self.__class__, self ).__init__( selector = SelectBest( key = lambda s: s.peak().rf() ), exhaustive = False ) self._fraction = 0.75 def manipulate(self, context, candidate, blackboard, logger): candidate.processed_by.add( self ) solved_ensemble_stages = set( [ st.predecessor( type = process_model.RotationStage ) for st in blackboard.read_instances_of( type = process_model.EvaluationStage ) if st.is_solution() ] ) if candidate in solved_ensemble_stages: logger.info( "Skipping search as a significant peak has already been found" ) return logger.info( "Starting translation search" ) context.submit( owner = self, protocol = mr_calculation.translation_function, execution = { "case": candidate.case(), "partial": candidate.partial(), "peak": candidate.peak(), }, processing = { "parent": candidate }, ) def process_success(self, parent, results, blackboard, logger): search = process_model.TranslationSearch.from_definitions( expert = self, parent = parent, definitions = results, exhaustive = True, ) logger.info( "There are %s peaks in the translation function" % len( search.peaks ) ) self.evaluate_search_results( search = search, parent = parent, blackboard = blackboard, ) def process_failure(self, parent, blackboard, logger): pass def evaluate_search_results(self, search, parent, blackboard): search.attach_to( blackboard = blackboard ) stages = [ process_model.TranslationStage( parent = parent, peak = p ) for p in search.peaks ] for s in stages: s.attach_to( blackboard = blackboard ) def acceptable(self, bbobj): return self not in bbobj.processed_by class PackingExpert(KnowledgeSource): """ Performs packing analysis for translation peaks """ BBOBJ = process_model.TranslationStage def __init__(self): super( self.__class__, self ).__init__( selector = SelectBest( key = lambda s: s.peak().tf() ), exhaustive = True ) def manipulate(self, context, candidate, blackboard, logger): candidate.processed_by.add( self ) logger.info( "Starting packing test" ) context.submit( owner = self, protocol = mr_calculation.packing_function, execution = { "case": candidate.case(), "partial": candidate.partial(), "peak": candidate.peak(), }, processing = { "parent": candidate }, ) def process_success(self, parent, results, blackboard, logger): ( packs, data ) = results if packs: logger.info( "Packing test passed" ) tpeak = mr_object.TranslationPeak( ensemble = parent.peak().ensemble(), rotation = data[0], translation = data[1], tf = parent.peak().tf(), tfz = parent.peak().tfz(), ) pstage = process_model.PackingStage( parent = parent, peak = tpeak ) pstage.attach_to( blackboard = blackboard ) else: logger.info( "Packing test fail" ) def process_failure(self, parent, blackboard, logger): pass def acceptable(self, bbobj): return self not in bbobj.processed_by class RefinementExpert(KnowledgeSource): """ Refines structures """ BBOBJ = process_model.RefinementStage def __init__(self): super( self.__class__, self ).__init__( selector = SelectBest( key = lambda s: s.peak().tf() ), exhaustive = True ) def manipulate(self, context, candidate, blackboard, logger): candidate.processed_by.add( self ) logger.info( "Starting refinement" ) context.submit( owner = self, protocol = mr_calculation.refinement_function, execution = { "case": candidate.case(), "partial": candidate.partial(), "peak": candidate.peak(), }, processing = { "parent": candidate }, ) def process_success(self, parent, results, blackboard, logger): ( llg, defs ) = results logger.info( "Refined LLG: %.2f" % llg ) enses = parent.partial().ensembles() + [ parent.peak().ensemble() ] rot_stage = parent.predecessor( type = process_model.RotationStage ) if rot_stage: stats = mr_object.SearchStatistics( is_solution = parent.is_solution(), rf = rot_stage.peak().rf(), rfz = rot_stage.peak().rfz(), tf = parent.peak().tf(), tfz = parent.peak().tfz(), ) else: stats = mr_object.SearchStatistics( is_solution = parent.is_solution(), tf = parent.peak().tf(), tfz = parent.peak().tfz(), ) partial = mr_object.PartialStructure( definitions = [ ( e, r, t ) for ( e, ( r, t ) ) in zip( enses, defs ) ], statistics = stats ) refined = mr_object.RefinedStructure( partial = partial, llg = llg ) ref_stage = process_model.RefinedStage( parent = parent, final = refined, ) ref_stage.attach_to( blackboard = blackboard ) def process_failure(self, parent, blackboard, logger): pass def acceptable(self, bbobj): return self not in bbobj.processed_by class ScoringExpert(KnowledgeSource): """ Scores structures """ BBOBJ = process_model.ScoringStage def __init__(self): super( self.__class__, self ).__init__( selector = lambda choices: choices[0], exhaustive = True ) def manipulate(self, context, candidate, blackboard, logger): candidate.processed_by.add( self ) logger.info( "Starting LLG scoring" ) context.submit( owner = self, protocol = mr_calculation.scoring_function, execution = { "case": candidate.case(), "partial": candidate.partial(), "peak": candidate.peak(), }, processing = { "parent": candidate }, ) def process_success(self, parent, results, blackboard, logger): llg = results logger.info( "LLG score: %.2f" % llg ) if parent.accept( score = llg ): logger.info( "This peak conforms to expectations, accepted" ) old = parent.peak() peak = mr_object.TranslationPeak( ensemble = old.ensemble(), rotation = old.rotation(), translation = old.translation(), tf = llg, tfz = old.tfz(), ) stage = process_model.SpeculatedStage( parent = parent, peak = peak ) stage.attach_to( blackboard = blackboard ) else: logger.info( "Peak score is worse than expected, rejected" ) def process_failure(self, parent, blackboard, logger): pass def acceptable(self, bbobj): return self not in bbobj.processed_by class EvaluationExpert(KnowledgeSource): """ Identifies solutions based on Z-score """ BBOBJ = process_model.PackingStage def __init__(self, z_threshold = 7.5): super( self.__class__, self ).__init__( selector = lambda choices: choices[0], exhaustive = True ) self._z_threshold = z_threshold def manipulate(self, context, candidate, blackboard, logger): candidate.processed_by.add( self ) translation = candidate.peak() if self._z_threshold < translation.tfz(): is_solution = True logger.info( "This is an outstanding peak" ) else: is_solution = False logger.info( "This is an ordinary peak" ) stage = process_model.IdentifiedStage( parent = candidate, is_solution = is_solution, ) stage.attach_to( blackboard = blackboard ) def acceptable(self, bbobj): return self not in bbobj.processed_by class TemplateSolveInitiateExpert(KnowledgeSource): """ Superposes ensembles on solutions and creates new solutions """ BBOBJ = process_model.EvaluationStage def __init__(self): super( self.__class__, self ).__init__( selector = SelectBest( key = lambda s: s.peak().tf() ), exhaustive = True ) def manipulate(self, context, candidate, blackboard, logger): candidate.processed_by.add( self ) if isinstance( candidate, process_model.SpeculatedStage ): logger.info( "This candidate is not applicable (speculative) " ) return ensembles = set( blackboard.read_instances_of( type = process_model.Ensemble ) ) assert ensembles used = candidate.peak().ensemble() assert used in ensembles ensembles.remove( used ) repo = process_model.TemplateSolveRepository( template = candidate, ensembles = ensembles, ) repo.attach_to( blackboard = blackboard ) def acceptable(self, bbobj): return bbobj.is_solution() and self not in bbobj.processed_by class TemplateSolveSuperposeExpert(KnowledgeSource): """ Superposes ensembles on solutions and creates new solutions """ BBOBJ = process_model.TemplateSolveRepository def __init__(self): super( self.__class__, self ).__init__( selector = lambda choices: choices[0], exhaustive = True ) self._proxy = SuperposePdbs() def manipulate(self, context, candidate, blackboard, logger): untested = list( candidate.untested() ) assert untested untested.sort( key = lambda s: s.score(), reverse = True ) selected = untested[0] candidate.remove( ensemble = selected ) logger.info( "Trial ensemble: %s" % selected ) template = candidate.template() peak = template.peak() ensemble = peak.ensemble() # Check whether this ensemble is applicable if not ensemble.composition().contains( selected.composition() ): logger.info( "Template ensemble does not match trial ensemble" ) return # Check whether we have a transformation, if not, superpose superposition_for = dict( [ ( sp.identifier(), sp ) for sp in blackboard.read_instances_of( type = process_model.Superposition ) ] ) forw_id = process_model.Superposition.descriptor( fixed = ensemble, moving = selected, ) rev_id = process_model.Superposition.descriptor( fixed = selected, moving = ensemble, ) if forw_id in superposition_for: logger.info( "Using pre-determined transformation" ) transformation = superposition_for[ forw_id ].transformation() elif rev_id in superposition_for: logger.info( "Using inverted pre-determined transformation" ) transformation = superposition_for[ rev_id ].transformation().inverse() else: logger.info( "Superposing ensembles" ) try: ( transf, rmsd ) = self._proxy.run( fixed = ensemble.pdb_file(), moving = selected.pdb_file(), ) except Sorry, e: logger.info( "Superposition error: %s" % e ) logger.info( "Discarding this trial" ) logger.info( "Action complete" ) sup = process_model.Superposition( fixed = ensemble, moving = selected, transformation = transf.rt(), rmsd = rmsd, ) logger.info( "RMSD between ensembles: %.3f" % rmsd ) sup.attach_to( blackboard = blackboard ) transformation = sup.transformation() case = template.case() cell = case.problem().unit_cell() solution = scitbx.matrix.rt( ( peak.rotation().matrix(), cell.orthogonalize( peak.translation().fractional() ), ) ) superposed = solution * transformation transformed = mr_object.TranslationPeak( ensemble = selected, rotation = mr_object.Matrix( elements = superposed.r.elems ), translation = mr_object.Fractional( coordinates = cell.fractionalize( superposed.t.elems ) ), tf = 0.0, tfz = peak.tfz() ) ense_stages = [ st for st in blackboard.read_instances_of( type = process_model.EnsembleStage ) if st.case() == case and st.search() == selected ] assert len( ense_stages ) == 1 logger.info( "Superposed solution generated" ) logger.info( "Starting packing test" ) context.submit( owner = self, protocol = mr_calculation.packing_function, execution = { "case": case, "partial": template.partial(), "peak": transformed, }, processing = { "parent": ense_stages[0], "peak": transformed, "accept": lambda score: 0 <= score }, ) def process_success(self, parent, peak, accept, results, blackboard, logger): ( packs, data ) = results if packs: logger.info( "Packing test passed" ) tpeak = mr_object.TranslationPeak( ensemble = peak.ensemble(), rotation = data[0], translation = data[1], tf = peak.tf(), tfz = peak.tfz(), ) stage = process_model.ScoringStage( parent = parent, peak = tpeak, accept = accept, ) stage.attach_to( blackboard = blackboard ) logger.info( "Sending peak for scoring" ) else: logger.info( "Packing test fail" ) logger.info( "Not pursuing this peak any further" ) def process_failure(self, parent, peak, accept, blackboard, logger): pass def acceptable(self, bbobj): return bool( bbobj.untested() ) class RegularObjectCompletionInitiateExpert(KnowledgeSource): """ Finds incomplete regular assemblies in solutions completes them """ BBOBJ = process_model.EvaluationStage TWO_PI = 2.0 * math.pi EPSILON = 1E-5 def __init__(self): super( self.__class__, self ).__init__( selector = SelectBest( key = lambda s: s.peak().tf() ), exhaustive = True ) self._max_fold = 12 self._spatial_sq = 3.5 self._angular = 0.15 self._axial = 0.999 self._max_fold = 12 def manipulate(self, context, candidate, blackboard, logger): candidate.processed_by.add( self ) if isinstance( candidate, process_model.SpeculatedStage ): logger.info( "This candidate is not applicable (speculative) " ) return if len( candidate.partial().ensembles() ) == 0: logger.info( "Cannot identify operators from a single molecule" ) return logger.info( "Sorting candidate on ensembles" ) description_with = {} p = candidate.partial() cell = candidate.case().problem().unit_cell() for ( e, r, t ) in zip( p.ensembles(), p.rotations(), p.translations() ): description_with.setdefault( e, [] ).append( phaser.ncs.Operator( rotation = r.matrix(), translation = cell.orthogonalize( t.fractional() ) ) ) peak = candidate.peak() description_with.setdefault( peak.ensemble(), [] ).append( phaser.ncs.Operator( rotation = peak.rotation().matrix(), translation = cell.orthogonalize( p.translation().fractional() ) ) ) possibilities = [] for ( ense, copies ) in description_with.items(): if len( copies ) < 2: logger.info( "Cannot identify operators for %s" % ense ) continue identified = [ mr_object.TranslationPeak( ensemble = ense, rotation = mr_object.Matrix( elements = op.get_rotation() ), translation = mr_object.Fractional( coordinates = cell.fractionalize( op.get_translation() ) ), tf = 0.0, tfz = peak.tfz() ) for op in self.possible_further_molecules( existing = copies ) ] logger.info( "%s possible new positions found" % len( identified ) ) if not possibilities: logger.info( "No new positions found" ) else: bbobj = process_model.RegularObjectCompletionRepository( base = candidate, extra = possibilities, ) bbobj.attach_to( blackboard = blackboard ) def possible_further_molecules(self, existing, centre, extent): unused = existing centres = [ op.transform( point = centre ) for op in existing ] solutions = [] clash_distance_square = ( extent / 4.0 ) ** 2 for ( index, current ) in enumerate( existing ): inv = current.inverse() ops = [ op * inv for op in existing[ index + 1 : ] ] print "Considering %s operators" % len( ops ) for ( fold, gen ) in self.generators( ops = ops ): print "Processing generator: %s-fold\n%s" % ( fold, gen ) trials = [ op * current for op in self.power_series( op = gen, length = fold ) ] print "%s trial operators will be checked" % len( trials ) points = [ op.transform( point = centre ) for op in trials ] non_clashing_ops = [ op for ( op, p ) in zip( trials, points ) if not any( d_sq < clash_distance_square for d_sq in self.distance_squares( point = p, others = centres ) ) ] print "%s non-clashing" % len( non_clashing_ops ) solutions.extend( non_clashing_ops ) centres.extend( [ op.transform( point = centre ) for op in non_clashing_ops ] ) return solutions def find_generators(self, ops): found = [] for op in ops: an = self.operator_analysis( op = op ) if not an: continue for group in found: if any( self.operation_overlap( left = an, right = m ) for m in group ): group.append( an ) return [ self.merge_operators( group = g ) for g in found ] generators.sort( key = lambda s: s[0], reverse = True ) uniques = [] accounted = [] for ( fold, op ) in generators: inv = op.inverse() if any( self.approx_inverse( op1 = inv, op2 = a ) for a in accounted ): continue uniques.append( ( fold, op ) ) accounted.extend( self.power_series( op = op, length = fold - 1 ) ) return uniques def operator_analysis(self, op): aaf = scitbx.math.r3_rotation_axis_and_angle_from_matrix( op.get_rotation() ) ( fold, power ) = self.get_fold_and_power( angle = aaf.angle() ) if fold in [ 0, 1 ]: return None assert 0 < power multiples = sorted( [ ( ( m * power ) % fold , op ) for ( m, op ) in zip( range( 1, fold ), self.get_power_series( op = op, length = fold ), ) ], key = lambda p: p[0] ) assert multiples[0][0] == 1 translation = scitbx.matrix.rec( op.get_translation(), ( 3, 1 ) ) t_par_unit = scitbx.matrix.rec( aaf.axis, ( 3, 1 ) ) displacement = t_par_unit.dot( translation ) * t_par_unit return ( fold, multiples, t_par_unit, displacement, power ) def get_fold_and_power(self, angle): rotation = 0 for fold in range( 1, self._max_fold + 1 ): rotation += angle multiple = rotation / self.TWO_PI rounded = round( multiple ) if abs( multiple - rounded ) < self._angular: return ( fold, int( rounded ) % fold ) return ( 0, 0 ) def get_power_series(self, op, length): series = [ op ] for index in range( 1, length ): series.append( series[-1] * op ) return series def operation_overlap(self, left, right): ( fold1, mult1, axis1, disp1, p1 ) = left ( fold2, mult2, axis2, disp2, p2 ) = right gcd = fractions.gcd( fold1, fold2 ) if gcd == 1: print "Wrong fold" return False if axis1.dot( axis2 ) < self._axial: print "Wrong axis" return False f1 = float( fold1 ) f2 = float( fold2 ) diff = f1 / p1 / gcd * disp1 - f2 / p2/ gcd * disp2 print "Diff: %s" % diff if self._spatial_sq < diff.length_sq(): print "Large difference in parallel transl" return False op1 = mult1[ fold1 / gcd - 1 ][1] op2 = mult2[ fold2 / gcd - 1 ][1] trans1 = scitbx.matrix.rec( op1.get_translation(), ( 3, 1 ) ) print trans1 trans2 = scitbx.matrix.rec( op2.get_translation(), ( 3, 1 ) ) print trans2 if self._spatial_sq < ( trans1 - trans2 ).length_sq(): print "large difference in perp transl" return False return True def merge_operators(self, group): count = len( group ) assert 0 < count fold = max( [ g[0] for g in group ] ) assert all( [ fold % g[0] == 0 for g in group ] ) axis = 1.0 / count * reduce( operator.add, [ g[2] for g in group ] ) perpendicular = 1.0 / count * reduce( operator.add, [ scitbx.matrix.rec( g[1][0][1].get_translation(), ( 3, 1 ) ) for g in group ] ) pnorm = float( reduce( operator.add, [ g[4] for g in group ] ) ) parallel = reduce( operator.add, [ g[3] for g in group ] ) / pnorm return ncs.Operator( rotation = scitbx.math.r3_rotation_axis_and_angle_as_matrix( axis, self.TWO_PI / fold ), translation = perpendicular + parallel ) """ def possible_further_molecules(self, existing, centre, extent): centres = [ op.transform( point = centre ) for op in existing ] solutions = [] clash_distance_square = ( extent / 4.0 ) ** 2 for ( index, current ) in enumerate( existing ): inv = current.inverse() ops = [ op * inv for op in existing[ index + 1 : ] ] print "Considering %s operators" % len( ops ) for ( fold, gen ) in self.generators( ops = ops ): print "Processing generator: %s-fold\n%s" % ( fold, gen ) trials = [ op * current for op in self.power_series( op = gen, length = fold ) ] print "%s trial operators will be checked" % len( trials ) points = [ op.transform( point = centre ) for op in trials ] non_clashing_ops = [ op for ( op, p ) in zip( trials, points ) if not any( d_sq < clash_distance_square for d_sq in self.distance_squares( point = p, others = centres ) ) ] print "%s non-clashing" % len( non_clashing_ops ) solutions.extend( non_clashing_ops ) centres.extend( [ op.transform( point = centre ) for op in non_clashing_ops ] ) return solutions def generators(self, ops): generators = [] for op in ops: aaf = scitbx.math.r3_rotation_axis_and_angle_from_matrix( op.get_rotation() ) if self._spatial < abs( self.parallel_component( vec1 = op.get_translation(), vec2 = aaf.axis ) ): continue ( fold, reduced ) = self.reduce( op = op ) if fold in [ 0, 1 ]: continue translation = scitbx.matrix.rec( reduced.get_translation(), ( 3, 1 ) ) t_par_unit = scitbx.matrix.rec( aaf.axis, ( 3, 1 ) ) t_parallel = t_par_unit.dot( translation ) * t_par_unit t_perpendicular = translation - t_parallel t_perp_unit1 = ( reduced_form.t - t_parallel ).normalize() t_perp_unit2 = t_par_unit.cross( t_perp_unit1 ) scale = ( reduced_form.t - t_parallel ).length() angle_sin = math.sqrt( 1 - angle_cos ** 2 ) rotmat = scitbx.matrix.sqr( ( angle_cos, angle_sin, -angle_sin, angle_cos ) ) solu = scitbx.matrix.identity( 2 ) - rotmat assert 1E-5 < solu.determinant() centre_plane_coords = ( solu.inverse() * scitbx.matrix.rec( ( scale, 0.0 ), ( 2, 1 ) ) ) centre_coords = ( centre_plane_coords[0] * t_perp_unit1 + centre_plane_coords[1] * t_perp_unit2 + t_parallel ) idealized = ncs.Operator( rotation = reduced.get_rotation(), translation = translation.elems ) generators.append( ( fold, reduced ) ) generators.sort( key = lambda s: s[0], reverse = True ) uniques = [] accounted = [] for ( fold, op ) in generators: inv = op.inverse() if any( self.approx_inverse( op1 = inv, op2 = a ) for a in accounted ): continue uniques.append( ( fold, op ) ) accounted.extend( self.power_series( op = op, length = fold - 1 ) ) return uniques def reduce(self, op, max_fold = 12): n = 1 prod = op reduced = op while n <= max_fold: angle = prod.rotation_angle_radian() if angle <= self._angular: return ( n, reduced ) elif angle < reduced.rotation_angle_radian(): reduced = prod n = n + 1 prod *= op return ( 0, op ) """ def approx_inverse(self, op1, op2): return ( op1.approx_inverse( other = op2, angular = self._angular, spatial = self._spatial ) or op2.approx_inverse( other = op1, angular = self._angular, spatial = self._spatial ) ) def distance_squares(self, point, others): return ( ( point[0] - o[0] ) ** 2 + ( point[1] - o[1] ) ** 2 + ( point[1] - o[1] ) ** 2 for o in others ) def acceptable(self, bbobj): return bbobj.is_solution() and self not in bbobj.processed_by class Modeller(object): """ Knowledge sources responsible for model generation Disabled if a solution exists """ def __init__(self, selector): self._selector = selector def inspect(self, context): if context.quick and context.statistics.solved: context.logger.debug( "Disabled because solutions have been found" ) return None candidates = self.candidates( bb = context.blackboard ) if not candidates: context.logger.debug( "No candidates have been found" ) return None context.logger.debug( "%s candidates has been found" % len( candidates ) ) selected = self._selector( choices = candidates ) context.logger.debug( "Selected %s" % selected ) return selected def candidates(self, bb): return [ obj for obj in bb.read_instances_of( type = self.BBOBJ ) if self.acceptable( bbobj = obj ) ] def __str__(self): return self.__class__.__name__ class FetchExpert(Modeller): """ Fetches data """ def __init__(self): super( FetchExpert, self ).__init__( selector = SelectBest( key = lambda s: s.score() ) ) def manipulate(self, context, candidate): candidate.processed_by.add( self ) context.logger.info( "Fetching data for %s..." % candidate ) try: raw = self.fetch( context = context, candidate = candidate ) except ValueError, e: context.logger.warning( "Error: %s" % e ) return context.logger.info( "Parsing data for %s..." % candidate ) try: data = self.parse( context = context, raw = raw, candidate = candidate ) except ValueError, e: context.logger.warning( "Error: %s" % e ) return self.publish( candidate = candidate, data = data, context = context ) context.logger.info( "Data fetch completed" ) class SequenceFetchExpert(FetchExpert): """ Obtain sequence data """ BBOBJ = process_model.BlastHit def __init__(self): super( SequenceFetchExpert, self ).__init__() self._proxy = OCAHTTP() def fetch(self, context, candidate): return self._proxy.fetch( identifier = candidate.identifier ) def parse(self, context, raw, candidate): ( seqs, bad ) = bioinformatics.seq_sequence_parse( raw ) if bad: context.logger.warning( "Uninterpretable text encountered: %s" % bad ) if not seqs: raise ValueError, "No sequence found" if len( seqs ) != 1: context.logger.warning( "Multiple reference sequences, keeping top one" ) seq = seqs[0] seq.name = candidate.identifier return seq def publish(self, candidate, data, context): candidate.sequence = data def acceptable(self, bbobj): return not bbobj.sequence and self not in bbobj.processed_by class StructureFetchExpert(FetchExpert): """ Obtain structural data """ BBOBJ = process_model.BlastHit def __init__(self): super( SequenceFetchExpert, self ).__init__() self._proxy = DbfetchHTTP() def fetch(self, context, candidate): return self._proxy.fetch( identifier = candidate.identifier[:4] ) def parse(self, context, raw, candidate): identifier = candidate.identifier input = iotbx.pdb.input( lines = raw, source_info = "PDB" ) root = input.construct_hierarchy() selection = "chain '%s'" % identifier[5:] asc = root.atom_selection_cache().selection( selection ) new_root = root.select( asc, True ) if not new_root.atoms(): raise ValueError, "No atoms selected for %s, chain '%s'" % ( pdb_id, chain_id ) file_name = os.path.join( context.data_folder, "%s.pdb" % identifier ) context.logger.info( "Writing out selected atoms as %s..." % file_name ) new_root.write_pdb_file( file_name ) return mr_object.Structure( file_name = file_name ) def publish(self, candidate, data, context): candidate.structure = data def acceptable(self, bbobj): return not bbobj.structure and bbobj.sequence and self not in bbobj.processed_by class AlignmentExpert(FetchExpert): """ Calculates an alignment for a sequence """ BBOBJ = process_model.BlastHit TARGET = "target" MODEL = "model" def __init__(self): super( AlignmentExpert, self ).__init__() self._proxy = ClustalWSynchronous() self._aligner = "clustalw" def fetch(self, context, candidate): seq1 = bioinformatics.sequence( name = self.TARGET, sequence = candidate.component.sequence.sequence ) seq2 = bioinformatics.sequence( name = self.MODEL, sequence = candidate.sequence.sequence ) ( aln, log, err ) = self._proxy.run( sequences = "\n".join( [ str( seq1 ), str( seq2 ) ] ) ) if err: context.logger.warning( "Warning: %s" % err ) if not aln: raise ValueError, "No alignment returned" return aln def parse(self, context, raw, candidate): ( aln, bad ) = bioinformatics.clustal_alignment_parse( raw ) if not aln: raise ValueError, "Unable to interpret: %s" % raw assert aln.multiplicity() == 2 if self.TARGET == aln.names[1]: assert self.MODEL == aln.names[0] context.logger.info( "Rearranging alignment..." ) aln.names = [ self.TARGET, self.MODEL ] aln.alignments = [ aln.alignments[1], aln.alignments[0] ] assert seq1.name == self.TARGET aln.names = [ candidate.component.sequence.name, candidate.sequence.name, ] return aln def publish(self, candidate, data, context): # Check that there is no identical alignment alignments = [ a for a in context.blackboard.read_instances_of( type = process_model.Alignment ) if a.blast_hit == candidate ] for ali in alignments: if data.alignments == ali.alignment.alignments: context.logger.info( "Alignment identical to %s" % ali ) return bbobj = process_model.Alignment( blast_hit = candidate, aligner = self._aligner, alignment = data, ) bbobj.attach_to( blackboard = context.blackboard ) def acceptable(self, bbobj): return bbobj.structure and bbobj.sequence and self not in bbobj.processed_by class SculptorExpert(Modeller): """ Runs sculptor """ BBOBJ = process_model.Alignment def __init__(self): super( SculptorExpert, self ).__init__( selector = SelectBest( key = lambda s: s.blast_hit.score() ) ) mms = ( rsam.MacromoleculeSpecifics.PROTEIN, ) self._workflow = [ chisel.MainchainDelete( deletion = [ chisel.MDAGap() ], polishing = [], applicable = mms ), chisel.SidechainPrune( pruning = [ chisel.SPASchwarzenbacher( level = 2 ) ], applicable = mms ), chisel.ResidueRenumber.Target( start = 1, applicable = mms ), chisel.ResidueRename( completion = [ chisel.AGACBetaUniform() ], applicable = mms ), chisel.CleanUnknown( known = [], applicable = ( rsam.MacromoleculeSpecifics.UNKNOWN, ) ), ] def manipulate(self, context, candidate): candidate.processed_by.add( self ) context.logger.info( "Running sculptor with:" ) context.logger.info( "\tAlignment:\n%s" % candidate.alignment ) context.logger.info( "\tStructure: %s, %s" % ( candidate.blast_hit.identifier, candidate.blast_hit.structure.file_name() ) ) context.submit( owner = self, protocol = self.run, execution = ( candidate.blast_hit.structure.file_name(), candidate.alignment, candidate.aligner, ), processing = ( candidate.blast_hit.structure, candidate.alignment.identity_fraction(), candidate.blast_hit.component, ), ) def acceptable(self, bb): return bbobj.blast_hit.structure and self not in bbobj.processed_by def run(self, data): ( file_name, alignment, aligner ) = data root = iotbx.pdb.input( file_name ).construct_hierarchy() logger = logging.Logger( name = "", level = logging.INFO ) logstream = StringIO() handler = logging.StreamHandler( logstream ) handler.setFormatter( logging.Formatter( "%(message)s" ) ) logger.addHandler( handler ) vr = sculptor.process( root = root, alignments = [ alignment ], workflow = self._workflow, logger = logger ) logger.removeHandler( handler ) handler.flush() handler.close() logfile = logstream.getvalue() fingerprint = set() for command in vr.modifications(): if isinstance( command, ( chisel.DeleteInstruction, chisel.AttachInstruction ) ): fingerprint.add( command.blueprint() ) for command in vr.modifications(): command.execute() vr.regularize() outfile = tbx_utils.autoname( file_name = file_name, label = "%s_%s_sculptor_%s" % ( file_name, aligner, id( self ) ), extension = "pdb" ) vr.root().write_pdb_file( outfile ) return ( True, logfile, ( outfile, fingerprint ) ) def process_success(self, results, context, extra): ( outfile, modifications ) = results ( origin, identity, component ) = extra # Check that there is no such ensemble ensembles = [ se for a in context.blackboard.read_instances_of( type = process_model.SculptorEnsemble ) if se.origin == origin and se.composition.components == [ component ] ] for ense in ensembles: if ense.modifications == modifications: context.logger.info( "Ensemble identical to %s" % ense ) return model = process_model.SculptorEnsemble( pdb_file = outfile, identity = identity, component = component, origin = origin, modifications = modifications ) model.announce( blackboard = context.blackboard ) # Other objects class SelectBest(object): """ Selects best object """ def __init__(self, key): self._key = key def __call__(self, choices): return max( choices, key = self._key )