from __future__ import division from phaser.pipeline import stage from phaser.pipeline import calculation from phaser.pipeline import mr_object from phaser import output from libtbx.utils import Sorry import operator import math class SelectBest(object): """ Select best of possible candidates """ def __init__(self, expert, type): self._expert = expert self._type = type def __call__(self, board, case, context): candidates = [ c for c in board.read( type = self._type ) if self._expert not in c.processed_by ] if not candidates: context.output.debug( msg = "No candidates found" ) selected = None else: context.output.debug( msg = "%s candidates found" % len( candidates ) ) selected = max( candidates, key = lambda stage: stage.score() ) context.output.debug( msg = "Selected:" ) context.output.indent() context.output.debug( msg = str( selected ) ) context.output.dedent() return selected class SelectRange(object): """ Select best N possible candidates """ def __init__(self, expert, type, count): self._expert = expert self._type = type self._count = count def __call__(self, board, case, context): candidates = [ c for c in board.read( type = self._type ) if self._expert not in c.processed_by ] if not candidates: context.output.debug( msg = "No candidates found" ) selected = None else: context.output.debug( msg = "%s candidates found" % len( candidates ) ) candidates.sort( key = lambda s: s.score(), reverse = True ) selected = candidates[ : self._count ] context.output.debug( msg = "Selected:" ) for ( i, c ) in enumerate( selected, start = 1 ): context.output.debug( "%d: %s" % ( i, c.short() ) ) return selected class Suspendable(object): """ Suspend if certain objects are present """ def __init__(self, marker, selection): self._marker = marker self._selection = selection def __call__(self, board, case, context): candidates = [ c for c in board.read( type = self._marker ) ] if candidates: context.output.debug( msg = "Knowledge source suspended" ) return None else: return self._selection( board = board, case = case, context = context ) class FindEnsemble(object): """ Find ensembles for missing part """ SCANS = [ ( "Models:", mr_object.Ensemble, stage.Ensemble ), ( "Collections:", mr_object.ModelCollection, stage.ModelCollection ), ( "Templates:", mr_object.Template, stage.ModelTemplate ), ( "Homology search hits:", mr_object.HomologySearchHit, stage.TemplateHit ), ] def __init__(self): self.inspect = SelectBest( expert = self, type = stage.Initial ) def manipulate(self, candidate, board, case, context): candidate.processed_by.add( self ) context.output.verbose( msg = output.subtitle( text = "Composition analysis" ) ) context.output.notify( message = "status", data = "Composition analysis" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = str( candidate ) ) for ( caption, intype, outtype ) in self.SCANS: context.output.verbose( msg = output.heading( text = caption ) ) context.output.indent() self.fetch_and_identify( candidate = candidate, board = board, case = case, context = context, intype = intype, outtype = outtype ) context.output.dedent() context.output.verbose( msg = output.section_close( text = "Composition analysis finished" ) ) def fetch_and_identify(self, candidate, board, case, context, intype, outtype): available = case.problem.data.read( type = intype ) + case.data.read( type = intype ) available.sort( key = lambda e: e.score(), reverse = True ) context.output.verbose( msg = "Available: %d" % len( available ) ) context.output.verbose( msg = "" ) accepted = [] for ( index, obj ) in enumerate( available, start = 1 ): if not candidate.missing.contains( other = obj.composition ): context.output.verbose( msg = "%d: %s - not in missing composition" % ( index, obj ) ) continue context.output.verbose( msg = "%d: %s - applicable for extension" % ( index, obj ) ) accepted.append( obj ) for ( index, obj ) in enumerate( accepted, start = 1 ): st = outtype( initial = candidate, intype = obj, index = index, ) board.write( obj = st ) def __str__(self): return "EnsembleFind" class RotationSearch(object): """ Performs rotation search for ensemble stages """ def __init__(self, cutoff): self.inspect = Suspendable( marker = stage.Suspension, selection = self.select ) self.cutoff = cutoff def select(self, board, case, context): submitteds = set( [ ( prog.calculation.partial, prog.calculation.ensemble ) for prog in case.problem.data.read( type = calculation.Progress ) if isinstance( prog.calculation, calculation.RotationSearch ) ] ) candidates = [ c for c in board.read( type = stage.Ensemble ) if ( self not in c.processed_by and ( c.partial(), c.ensemble ) not in submitteds ) ] if not candidates: context.output.debug( msg = "No candidates found" ) selected = None else: context.output.debug( msg = "%s candidates found" % len( candidates ) ) selected = max( candidates, key = lambda stage: stage.score() ) context.output.debug( msg = "Selected:" ) context.output.indent() context.output.debug( msg = selected ) context.output.dedent() return selected def manipulate(self, candidate, board, case, context): candidate.processed_by.add( self ) context.output.verbose( msg = output.subtitle( text = "Rotation search" ) ) context.output.notify( message = "status", data = "Rotation search" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = candidate ) # There is no need to check for equivalent searches - this is done by # the following object calc = calculation.RotationSearch( data = case.calculation_data(), partial = candidate.partial(), ensemble = candidate.ensemble, cutoff = self.cutoff ) context.output.verbose( msg = output.section_close( text = "Rotation job submitted" ) ) case.problem.data.write( obj = calculation.Progress( calculation = calc ) ) context.queue.submit( owner = self, calculation = calc, extra = candidate ) def process(self, calculation, result, extra, board, case, context): from phaser.pipeline import calculation as calc context.output.verbose( msg = output.subtitle( text = "Rotation function processing" ) ) context.output.notify( message = "status", data = "Rotation function evaluation" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = str( extra ) ) context.output.verbose( msg = output.blank() ) progress = [ p for p in case.problem.data.read( type = calc.Progress ) if p.calculation == calculation ] assert len( progress ) == 1 case.problem.data.delete( progress[0] ) res = calculation.process( results = result ) context.output.debug( msg = "Logfile:" ) context.output.debug( msg = res.logfile ) if calculation.partial == case.problem.root: context.output.verbose( msg = "Rotation function valid for all space groups: stored" ) # Remove data calculation.cleanup() case.problem.data.write( obj = res ) self.place_peaks( result = res, ensemble = extra, board = board, case = case, context = context ) def place_peaks(self, result, ensemble, board, case, context): context.output.verbose( msg = "Peaks above threshold: %d" % len( result.data ) ) for ( i, p ) in enumerate( result.data, start = 1 ): context.output.verbose( msg = "" ) context.output.verbose( msg = "%d: %s" % ( i, p ) ) st = stage.NormalRotation( ensemble = ensemble, peak = p, index = i ) board.write( obj = st ) context.output.verbose( msg = output.section_close( text = "Rotation function processing finished" ) ) def __str__(self): return "RotationSearch" class ApplyRotationSearchResults(object): """ Reinstates existing rotation searches (only important for first molecule) """ def __init__(self, rotation_search_expert): self.rotation_search_expert = rotation_search_expert def inspect(self, board, case, context): selected = self.rotation_search_expert.inspect( board = board, case = case, context = context ) if not selected: return None calcs = ( res for res in case.problem.data.read( type = calculation.Result ) if ( isinstance( res.calculation, calculation.RotationSearch ) and res.calculation.partial == selected.partial() and res.calculation.ensemble == selected.ensemble ) ) try: search = calcs.next() except StopIteration: context.output.debug( msg = "No search has been performed for this peak" ) return None context.output.debug( msg = "Selected rotation search has been performed" ) return ( selected, search ) def manipulate(self, candidate, board, case, context): ( selected, search ) = candidate selected.processed_by.add( self.rotation_search_expert ) context.output.verbose( msg = output.subtitle( text = "Apply rotation search results" ) ) context.output.notify( message = "status", data = "Apply rotation search results for other space group" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = selected ) self.rotation_search_expert.place_peaks( result = search, ensemble = selected, board = board, case = case, context = context ) def __str__(self): return "ApplyRotationSearchResults" class TranslationSearch(object): """ Performs translation search for rotation stages """ def __init__(self): self.inspect = SelectBest( expert = self, type = stage.Rotation ) def manipulate(self, candidate, board, case, context): candidate.processed_by.add( self ) context.output.verbose( msg = output.subtitle( text = "Translation search" ) ) context.output.notify( message = "status", data = "Translation search" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = str( candidate ) ) if len( candidate.partial().peaks ) == 0 and case.space_group_number() == 1: context.output.verbose( msg = "P1 space group with empty located part: skip translation search" ) self.place_peaks( peaks = [ mr_object.TranslationPeak( ensemble = candidate.peak.ensemble, rotation = candidate.peak.rotation, translation = ( 0.0, 0.0, 0.0 ), tf = candidate.peak.rf, tfz = candidate.peak.rfz, bfactor = 0.0 ), ], rotation = candidate, board = board, case = case, context = context ) candidate.translation_done = True context.output.verbose( msg = output.section_close( text = "Translation function finished" ) ) else: calc = calculation.TranslationSearch( data = case.calculation_data(), partial = candidate.partial(), rpeak = candidate.peak ) context.output.verbose( msg = output.section_close( text = "Translation job submitted" ) ) context.queue.submit( owner = self, calculation = calc, extra = candidate ) def process(self, calculation, result, extra, board, case, context): context.output.verbose( msg = output.subtitle( text = "Translation function processing" ) ) context.output.notify( message = "status", data = "Translation function evaluation" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = str( extra ) ) context.output.verbose( msg = output.blank() ) res = calculation.process( results = result ) context.output.debug( msg = "Logfile:" ) context.output.debug( msg = res.logfile ) context.output.verbose( msg = "Peaks above threshold: %d" % len( res.data ) ) self.place_peaks( peaks = res.data, rotation = extra, board = board, case = case, context = context ) extra.translation_done = True context.output.verbose( msg = output.section_close( text = "Translation function processing finished" ) ) def place_peaks(self, peaks, rotation, board, case, context): for ( i, p ) in enumerate( peaks, start = 1 ): context.output.verbose( msg = "" ) context.output.verbose( msg = "%d: %s" % ( i, p ) ) st = stage.Translation( rotation = rotation, peak = p, index = i ) board.write( obj = st ) def __str__(self): return "TranslationSearch" class PooledPackingCheck(object): """ Performs packing check for translation stages pooling all peaks """ def __init__(self, pool): self.inspect = SelectRange( type = stage.Translation, expert = self, count = pool ) def manipulate(self, candidate, board, case, context): context.output.verbose( msg = output.subtitle( text = "Packing check" ) ) context.output.notify( message = "status", data = "Packing check" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = "Pooled: %d" % len( candidate ) ) context.output.verbose( msg = "" ) context.output.indent() for c in candidate: context.output.verbose( msg = str( c ) ) context.output.verbose( msg = "" ) c.processed_by.add( self ) calc = calculation.SinglePoolPackingCheck( data = case.calculation_data(), partials = [ c.partial() for c in candidate ], tpeaks = [ c.peak for c in candidate ] ) context.output.dedent() context.output.verbose( msg = output.section_close( text = "Packing job submitted" ) ) context.queue.submit( owner = self, calculation = calc, extra = candidate ) def process(self, calculation, result, extra, board, case, context): context.output.verbose( msg = output.subtitle( text = "Packing check processing" ) ) context.output.notify( message = "status", data = "Packing check evaluation" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) res = calculation.process( results = result ) context.output.debug( msg = "Logfile:" ) context.output.debug( msg = res.logfile ) assert len( res.data ) == len( extra ) for ( d, c ) in zip( res.data, extra ): if d: verdict = "PASS" st = stage.Packed( translation = c, peak = d ) else: verdict = "FAILS" st = stage.Clash( translation = c ) context.output.verbose( msg = "Peak: %s" % verdict ) context.output.indent() context.output.verbose( msg = str( c ) ) context.output.dedent() context.output.verbose( msg = "" ) board.write( obj = st ) c.packing_done = True context.output.verbose( msg = output.section_close( text = "Packing check processing finished" ) ) def __str__(self): return "PooledPackingCheck" class SuperposeSolve(object): """ Finds ensembles that could be superposed on a solutions to get new solutions """ def inspect(self, board, case, context): candidates = ( c for c in board.readiter( type = stage.Scored ) if c.significant() and self not in c.processed_by ) try: selected = candidates.next() except StopIteration: context.output.debug( msg = "No candidates found" ) selected = None else: context.output.debug( msg = "Selected:" ) context.output.indent() context.output.debug( msg = selected ) context.output.dedent() return selected def manipulate(self, candidate, board, case, context): context.output.verbose( msg = output.subtitle( text = "Template structure solution" ) ) context.output.notify( message = "status", data = "Superpose solution" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = str( candidate ) ) context.output.verbose( msg = "" ) candidate.processed_by.add( self ) assert candidate.significant() context.output.verbose( msg = "This peak is significant and could be used" ) context.output.verbose( msg = "" ) ensemble = candidate.peak().ensemble initial = candidate.initial() estages = [ e for e in board.read( type = stage.Ensemble ) if ( e.ensemble != ensemble and e.initial() == initial ) ] applicable = [] for e in estages: if not ensemble.composition.contains( other = e.ensemble.composition ): continue overlap = ensemble.composition.overlap( other = e.ensemble.composition ) if overlap.superposable(): applicable.append( e ) context.output.verbose( msg = "Applicable: %s cases" % len( applicable ) ) context.output.verbose( msg = "" ) for es in applicable: context.output.verbose( msg = "Replacement ensemble: %s" % es.ensemble ) assert es.ensemble != ensemble context.output.indent() try: peak = context.superposition.substitute( peak = candidate.peak(), ensemble = es.ensemble, cell = case.uctbx_unit_cell(), stream = context.output.verbose ) except RuntimeError, e: context.output.verbose( msg = "Superposition error: %s" % e ) context.output.verbose( msg = "Alternative rejected" ) continue finally: context.output.dedent() if candidate.peak().ensemble.superpose_template() != es.ensemble.superpose_template(): context.output.verbose( msg = "Different ensembles: do refinement before scoring" ) calc = calculation.RefinedEvaluation( data = case.calculation_data(), partial = initial.partial(), peak = peak, protocol = context.b_factor_refinement, ) method = self.process_refined else: context.output.verbose( msg = "Ensembles have the same template: no refinement before scoring" ) calc = calculation.Evaluation( data = case.calculation_data(), partial = initial.partial(), peak = peak, ) method = self.process_unrefined context.output.verbose( msg = "Evaluation job submitted" ) context.output.verbose( msg = "" ) context.queue.submit( owner = self, calculation = calc, extra = ( candidate, es, method ) ) context.output.verbose( msg = output.section_close( text = "Template structure solution finished" ) ) def process(self, calculation, result, extra, board, case, context): context.output.verbose( msg = output.subtitle( text = "Template structure evaluation processing" ) ) context.output.notify( message = "status", data = "Superpose solution evaluation" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) ( template, es, method ) = extra context.output.verbose( msg = str( template ) ) context.output.verbose( msg = output.heading( text = "Replacement ensemble: %s" % es.ensemble ) ) res = calculation.process( results = result ) ( is_pak, peak ) = ( res.data[0], res.data[1] ) if not is_pak: context.output.verbose( msg = "Generated peak fails packing test" ) context.output.verbose( msg = "Discarding this solution" ) return context.output.verbose( msg = "Accepting this solution" ) st = method( peak = peak, data = res.data, ensemble = es, template = template, context = context, ) board.write( obj = st ) context.output.verbose( msg = output.section_close( text = "Template structure evaluation processing finished" ) ) def process_unrefined(self, peak, data, ensemble, template, context): context.output.verbose( msg = "Peak statistics: tf = %.3f, tfz = %.2f" % ( peak.tf, peak.tfz ) ) return stage.Templated( ensemble = ensemble, template = template, peak = peak, refined = False, ) def process_refined(self, peak, data, ensemble, template, context): llg = data[2] context.output.verbose( msg = "Peak statistics: tf = %.3f, tfz = %.2f, llg = %.2f" % ( peak.tf, peak.tfz, llg ) ) return stage.Templated( ensemble = ensemble, template = template, peak = peak, refined = True, llg = llg ) def __str__(self): return "SuperposeSolve" class RotationPeakSalvage(object): """ Transfer good rotation peaks to the next round (if there is one) in order to save rotation function calculation """ def __init__(self): self.inspect = SelectBest( expert = self, type = stage.Ensemble ) def manipulate(self, candidate, board, case, context): context.output.verbose( msg = output.subtitle( text = "Previously located significant rotation peaks" ) ) context.output.notify( message = "status", data = "Reinstate significant rotation peaks found in last cycle" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = str( candidate ) ) context.output.verbose( msg = "" ) candidate.processed_by.add( self ) related = [ hr for hr in case.data.read( type = mr_object.HistoricRotation ) if ( hr.structure == candidate.partial() and hr.peak.ensemble == candidate.ensemble ) ] if not related: context.output.verbose( msg = output.section_close( text = "No previously located peaks found" ) ) return related.sort( key = lambda r: r.peak.rf ) rotations = [] for ( index, hr ) in enumerate( related, start = 1 ): rp = stage.PreviouslyEstablishedRotation( ensemble = candidate, peak = hr.peak, index = index ) board.write( obj = rp ) rotations.append( rp ) case.data.delete( obj = hr ) context.output.verbose( msg = "Reinstating: %d peaks" % len( related ) ) suspension = stage.Suspension( condition = RotationsExplored( board = board, rotations = rotations ) ) board.write( obj = suspension ) context.output.verbose( msg = output.section_close( text = "Rotation searches are suspended until these peaks are processed" ) ) def __call__(self, board, case, context, solutions, postprocessings): context.output.verbose( msg = output.subtitle( text = "Rotation peak salvage - identification" ) ) context.output.notify( message = "status", data = "Search for untested significant rotation peaks" ) useds = set( [ s.rotation for s in board.read( type = stage.Translation ) ] ) unuseds = [ c for c in board.read( type = stage.Rotation ) if 7.0 <= c.peak.rfz and c not in useds ] context.output.verbose( msg = "Unexplored significant rotation peaks: %d" % len( unuseds ) ) peaks = [] for s in solutions: initial = s.initial() found = [ r for r in unuseds if initial == r.initial() ] if not found: continue context.output.verbose( msg = output.heading( text = "Associations:" ) ) context.output.indent() context.output.verbose( msg =output.heading( text = "Structure:" ) ) context.output.indent() context.output.verbose( msg = str( s ) ) context.output.dedent() context.output.verbose( msg = output.heading( text = "Peaks:" ) ) context.output.indent() found.sort( key = lambda s: s.score(), reverse = True ) # Add new peaks for ( i, r ) in enumerate( found, start = 1 ): context.output.verbose( msg = "%d: %s" % ( i, r.peak ) ) context.output.verbose( msg = "" ) hrot = mr_object.GenuineRotation( structure = s.structure, peak = r.peak ) peaks.append( hrot ) context.output.dedent() context.output.dedent() context.output.verbose( msg = output.section_close( text = "Rotation peak salvage finished" ) ) return peaks def __str__(self): return "RotationPeakSalvage" class RotationsExplored(object): """ Suspension until packing jobs starting from rotation peak finish """ def __init__(self, board, rotations): self.board = board self.rotations = set( rotations ) def __call__(self): if not all( r.translation_done for r in self.rotations ): return False return all( s.packing_done for s in self.board.read( type = stage.Translation ) if s.rotation in self.rotations ) class PeakEvaluate(object): """ Evaluate potential peaks that have been associated with a structure """ def __init__(self): self.inspect = SelectBest( expert = self, type = stage.Ensemble ) def manipulate(self, candidate, board, case, context): context.output.verbose( msg = output.subtitle( text = "Peak evaluation" ) ) context.output.notify( message = "status", data = "Evaluate potential peaks" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = candidate ) context.output.verbose( msg = output.blank() ) candidate.processed_by.add( self ) partial = candidate.partial() data = case.calculation_data() related = [ pp for pp in case.data.read( type = mr_object.PotentialPeak ) if pp.is_related( partial = partial, ensemble = candidate.ensemble ) ] if not related: context.output.verbose( msg = output.section_close( text = "No previously located peaks found" ) ) return for ( index, peak ) in enumerate( related, start = 1 ): handler = peak.get_handler_for( partial = partial, index = index ) context.output.verbose( msg = handler ) context.output.indent() context.output.verbose( msg = handler.peak ) context.output.dedent() calc = handler.calculation( data = data, context = context ) context.output.verbose( msg = "Evaluation job (%s) submitted" % calc.__class__.__name__ ) suspension = stage.Suspension( condition = lambda: False ) board.write( obj = suspension ) context.queue.submit( owner = self, calculation = calc, extra = ( candidate, handler, suspension ) ) # Prevent symmetric jobs from being submitted case.data.delete( obj = peak ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = "Rotation searches suspended until these peaks are processed" ) context.output.verbose( msg = output.section_close( text = "Peak evaluation finished" ) ) def process(self, calculation, result, extra, board, case, context): context.output.verbose( msg = output.subtitle( text = "Peak evaluation processing" ) ) context.output.notify( message = "status", data = "Evaluate peak rescoring results" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) ( ensemble, handler, suspension ) = extra board.delete( obj = suspension ) context.output.verbose( msg = ensemble ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = handler ) context.output.verbose( msg = "%s (trial)" % handler.peak ) context.output.verbose( msg = output.blank() ) res = handler.process( result = calculation.process( results = result ) ) if not res.is_pak: context.output.verbose( msg = "Generated peak fails packing test" ) context.output.verbose( msg = "Discarding this solution" ) failure = mr_object.FillFailureRecord( structure = handler.structure, peak = handler.peak ) case.data.write( obj = failure ) else: context.output.verbose( msg = "%s (final)" % res.peak ) context.output.verbose( msg = "Peak statistics: %s" % res.statistics ) context.output.verbose( msg = "Accepting this solution" ) board.write( obj = handler.stage( parent = ensemble, result = res ) ) context.output.verbose( msg = output.section_close( text = "Peak evaluation processing finished" ) ) def __str__(self): return "PeakEvaluate" class Amalgamation(object): """ Mix solutions if coming from the same starting point """ def __call__(self, board, case, context, solutions, postprocessings): context.output.verbose( msg = output.subtitle( text = "Amalgamation - identification" ) ) context.output.notify( message = "status", data = "Find solutions that can be amalgamated" ) # Setup handling based on origin definition # Assume that we either have an undefined origin or a defined one if self.origin_undefined( solutions = solutions ): context.output.verbose( msg = "Initial structure empty - undefined origin" ) from phaser import AlterOriginSymmate origins = AlterOriginSymmate.GetAllOrigins( case.space_group_number() ) if self.is_floating_origin( origins = origins ): context.output.verbose( msg = "Floating origin - extracting rotation component only" ) collect = self.rotation_peaks_collect collate = self.rotation_peaks_collate else: context.output.verbose( msg = "Discrete origin - adding possible origin shifts" ) collect = self.associated_peaks_collect collate = self.associated_peaks_collate else: origins = [ ( 0, 0, 0 ) ] collect = self.associated_peaks_collect collate = self.associated_peaks_collate # Find new associations significants = [ s for s in solutions if s.structure.significant ] extra = [ s for s in postprocessings if s.structure.significant ] context.output.verbose( msg = "Significant peaks: %d" % len( significants ) ) context.output.verbose( msg = output.heading( text = "Processing:" ) ) results = [] from phaser import matching for st in significants: context.output.verbose( msg = "Structure:" ) context.output.indent() context.output.verbose( msg = str( st ) ) context.output.dedent() initial = st.initial() usables = [ s for s in ( significants + extra ) if st.evaluated.mixable_for_amalgamation( other = s.evaluated ) ] context.output.verbose( msg = "" ) context.output.verbose( "Potential peaks: %d" % len( usables ) ) # Filter peaks that are the same but with alternative models accepteds = [] peak = st.structure.peaks[-1] tmf = case.create_matching_mr_structure( structure = mr_object.Structure( peaks = [ peak ] ), cached = True ) context.output.indent() for u in usables: context.output.verbose( msg = "" ) context.output.verbose( msg = str( u.structure ) ) latest = u.structure.peaks[-1] if peak.ensemble == latest.ensemble: context.output.verbose( msg = "ACCEPTED: unique peak with same model" ) accepteds.append( u ) continue if peak.ensemble.composition != latest.ensemble.composition: context.output.verbose( msg = "ACCEPTED: peak with non-alternative model" ) accepteds.append( u ) continue context.output.verbose( msg = "Relating %s to %s:" % ( peak.ensemble, latest.ensemble ) ) context.output.indent() try: subst_peak = context.superposition.substitute( peak = latest, ensemble = peak.ensemble, cell = case.uctbx_unit_cell(), stream = context.output.debug, ) except RuntimeError, e: context.output.verbose( msg = "REJECTED: no superposition" ) continue finally: context.output.dedent() equiv = case.create_matching_mr_structure( structure = mr_object.Structure( peaks = [ subst_peak ] ) ) if not matching.is_equivalent( left = tmf, right = equiv ): context.output.verbose( msg = "ACCEPTED: unique peak with alternative model" ) accepteds.append( u ) continue context.output.verbose( msg = "REJECTED: identical peak with alternative model" ) context.output.dedent() context.output.verbose( msg = "" ) context.output.verbose( msg = "After filtering identical solutions with alternative models: %d" % len( accepteds ) ) if not accepteds: context.output.verbose( msg = "No mixable peaks found" ) context.output.verbose( msg = "\n" ) continue context.output.verbose( msg = "Accepted peaks: %d" % len( accepteds ) ) results.extend( collect( structure = st.structure, relateds = accepteds, context = context ) ) context.output.verbose( msg = "" ) possibilities = collate( collection = results, origins = origins, context = context ) context.output.verbose( msg = output.section_close( text = "Amalgamation finished" ) ) return possibilities @staticmethod def origin_undefined(solutions): ods = [] for s in solutions: ods.append( len( s.initial().partial().peaks ) == 0 ) assert ods == [ True ] * len( solutions ) or ods == [ False ] * len( solutions ), ods return any( ods ) @staticmethod def rotation_peaks_collect(structure, relateds, context): rotations = [] for ( i, r ) in enumerate( relateds, start = 1 ): peak = r.structure.peaks[-1] rp = mr_object.RotationPeak( ensemble = peak.ensemble, rotation = peak.rotation, rf = 0.0, rfz = 0.0, ) context.output.verbose( msg = "%d: %s" % ( i, rp ) ) context.output.verbose( msg = "" ) exrp = mr_object.ExtractedRotation( structure = structure, peak = rp ) rotations.append( exrp ) return rotations @staticmethod def rotation_peaks_collate(collection, origins, context): context.output.verbose( msg = "Number of extracted rotations: %d" % len( collection ) ) return collection @staticmethod def associated_peaks_collect(structure, relateds, context): associations = [] for ( i, r ) in enumerate( relateds, start = 1 ): current = r.structure context.output.verbose( msg = "%d: %s" % ( i, current.peaks[-1] ) ) context.output.verbose( msg = "" ) associations.append( frozenset( [ structure, current ] ) ) return associations @staticmethod def associated_peaks_collate(collection, origins, context): uniques = set( collection ) context.output.verbose( msg = "Number of possibilities found (taking symmetry into account): %d" % len( uniques ) ) return [ mr_object.Amalgamation( structures = fs, origins = origins ) for fs in uniques ] @staticmethod def is_floating_origin(origins): tokens = set( reduce( operator.add, origins, () ) ) return bool( tokens.intersection( [ "x", "y", "z" ] ) ) def __str__(self): return "Amalgamation" def substitute_peak(peak, ensemble, cell, superposition, stream): if not peak.ensemble.substituteable_with( other = ensemble ): raise RuntimeError, "different composition" try: speak = superposition.substitute( peak = peak, ensemble = ensemble, cell = cell, stream = stream ) except RuntimeError, e: raise RuntimeError, "no superposition" return speak def get_superposition(reference, moving, superposition, stream): if not reference.substituteable_with( other = moving ): raise RuntimeError, "different composition" try: sup = superposition.superpose( reference = reference, moving = moving, stream = stream, ) except RuntimeError, e: raise RuntimeError, "no superposition" return sup class GeneratorSearch(object): """ Recognizes generators (pure rotations around an axis) Only active in post-processing """ def __init__(self, tolerances): self.tolerances = tolerances def __call__(self, board, case, context): context.output.verbose( msg = output.underlined( text = "Generator search" ) ) context.output.notify( message = "status", data = "Find generators" ) significants = [ l for l in board.read( type = mr_object.ForSymmetryAnalysis ) ] significants.sort( key = lambda s: s.index ) context.output.verbose( msg = "Analysing %s structures" % len( significants ) ) for label in significants: context.output.verbose( msg = "%d. structure" % label.index ) relations = self.find_assemblies_in( structure = label.structure, case = case, context = context ) context.output.verbose( msg = "Total %s geometric relations found" % len( relations ) ) if relations: board.write( obj = mr_object.StructureRelations( structure = label.structure, relations = relations, ) ) context.output.verbose( msg = "Saved relations for further analysis" ) context.output.verbose( msg = "" ) context.output.verbose( msg = output.section_close( text = "Generator search finished" ) ) def find_assemblies_in(self, structure, case, context): cell = case.uctbx_unit_cell() relations = [] for ( i_ref, reference ) in enumerate( structure.peaks[:-1] ): context.output.verbose( msg = "Reference peak:" ) context.output.verbose( msg = "%s: %s" % ( i_ref + 1, reference ) ) context.output.verbose( msg = output.blank() ) context.output.debug( msg = "Comparing with other peaks:" ) context.output.indent() inverse = reference.as_scitbx_rt_operator( cell = cell ).inverse() einfo = reference.ensemble.matching_ensemble_info() for ( i_other, other ) in enumerate( structure.peaks[ i_ref + 1 :], start = i_ref + 1 ): context.output.verbose( msg = "%s: %s" % ( i_other + 1, other ) ) try: speak = substitute_peak( peak = other, ensemble = reference.ensemble, cell = cell, superposition = context.superposition, stream = context.output.debug ) except RuntimeError, e: context.output.debug( msg = " REJECTED: %s" % e ) context.output.debug( msg = output.blank() ) continue context.output.debug( msg = " ACCEPTED" ) members = self.potential_assembly_members_for( reference = reference, other = speak, case = case, ) context.output.debug( msg = " Testing %s symmetry equivalents for assembly membership" % len( members ) ) for ( i, m ) in enumerate( members, start = 1 ): context.output.verbose( msg = " %s. potential mate:" % i ) context.output.verbose( msg = " %s" % m ) transf = inverse * m.as_scitbx_rt_operator( cell = cell ) try: gen = self.to_generator_form( op = transf, ensemble = einfo ) except RuntimeError, e: context.output.verbose( msg = " REJECTED: %s" % e ) continue context.output.verbose( msg = " ACCEPTED" ) context.output.verbose( msg = " Generator:" ) context.output.verbose( msg = " %s" % gen ) relations.append( mr_object.ModelRelation( setting = mr_object.GeneratorSetting( ensemble = reference.ensemble, generator = gen, ), relating = ( reference, other ), ) ) context.output.dedent() context.output.verbose( msg = output.blank() ) return relations def potential_assembly_members_for(self, reference, other, case): assert reference.ensemble == other.ensemble tmf_ref = case.create_matching_mr_peak( peak = reference ).positional() tmf_other = case.create_matching_mr_peak( peak = other ).positional() potentials = [] import scitbx.matrix import itertools extent = reference.ensemble.matching_ensemble_info().extent cell = case.uctbx_unit_cell() for ( symop, delta ) in tmf_ref.orthogonal_differences_between( other = tmf_other ): transformed = scitbx.matrix.col( symop.sgop * tmf_other.centre ) shift = scitbx.matrix.col( [ round( n ) for n in ( tmf_ref.centre - transformed ) ] ) rotation = mr_object.Rotation.Matrix( elements = ( symop.orth_rotation * scitbx.matrix.sqr( other.rotation.matrix() ) ).elems ) translation = scitbx.matrix.col( symop.sgop * other.translation ) for neighbour in itertools.product( [ -1, 0, 1 ], repeat = 3 ): ishift = shift + scitbx.matrix.col( neighbour ) npos = scitbx.matrix.col( cell.orthogonalize( transformed + ishift - tmf_ref.centre ) ) if self.sufficiently_close( delta = npos, extent = extent ): potentials.append( mr_object.Peak( ensemble = reference.ensemble, rotation = rotation, translation = ( translation + ishift ).elems, bfactor = 0, ) ) return potentials def to_generator_form(self, op, ensemble): try: generator = mr_object.Generator.from_scitbx_rt_operator( op = op, max_fold = self.tolerances.max_fold, angular = self.tolerances.max_angular_error, spatial = self.tolerances.max_spatial_error ) except ValueError, e: raise RuntimeError, e for o in generator.nontrivial_series(): current = o * ensemble.centre delta = ensemble.centre - current if not self.sufficiently_close( delta = delta, extent = ensemble.extent ): raise RuntimeError, "Generated assembly larger than allowed" return generator def sufficiently_close(self, delta, extent): return delta.length() / extent <= self.tolerances.max_relative_assembly_size def __str__(self): return "Generator search" class GeneratorClassifyAverage(object): """ Averages near-identical generators Only active in post-processing """ def __init__(self, tolerances): self.tolerances = tolerances def __call__(self, board, case, context): context.output.verbose( msg = output.underlined( text = "Generator classification and averaging" ) ) context.output.notify( message = "status", data = "Classify and average generators" ) analysis_for = dict( [ ( r.structure, r.relations ) for r in board.read( type = mr_object.StructureRelations ) ] ) context.output.verbose( msg = "Analysing geometric relations data for %s structures" % len( analysis_for ) ) origin_of = dict( reduce( operator.add, [ [ ( r, s ) for r in rels ] for ( s, rels ) in analysis_for.items() ], [], ) ) context.output.verbose( msg = "%s relations found" % len( origin_of ) ) components_of = self.merge_equivalent_assemblies( relations = origin_of.keys(), context = context ) context.output.verbose( msg = "%s clusters formed" % len( components_of ) ) for ( average, members ) in components_of.items(): gclust = mr_object.GeneratorCluster( average = average ) for rel in members: gclust.add( structure = origin_of[ rel ], relation = rel ) board.write( obj = gclust ) context.output.verbose( msg = output.section_close( text = "Generator classification finished" ) ) def merge_equivalent_assemblies(self, relations, context): remaining = set( relations ) members_of = {} while remaining: reference = remaining.pop() context.output.verbose( msg = "Reference:" ) context.output.verbose( msg = reference.setting ) context.output.verbose( msg = "Transform other assemblies to same ensemble" ) context.output.indent() equiv = mr_object.EquivalenceRelation( relation = lambda left, right: ( left[0].approx_equal_spatially( other = right[0], max_misalignment = self.tolerances.min_axis_alignment, max_spatial_error = self.tolerances.max_spatial_error ) ) ) equiv.add( value = ( reference.setting.generator, reference ) ) for trial in list( remaining ): context.output.verbose( msg = "Trial:" ) context.output.verbose( msg = trial.setting ) try: sup = get_superposition( reference = reference.setting.ensemble, moving = trial.setting.ensemble, superposition = context.superposition, stream = context.output.debug ) except RuntimeError, e: context.output.debug( msg = " REJECTED: %s" % e ) context.output.debug( msg = output.blank() ) continue tr = trial.setting.generator.transformed( operation = sup.transformation ) context.output.verbose( msg = "Transformed generator: %s" % tr ) equiv.add( value = ( tr, trial ) ) remaining.remove( trial ) context.output.dedent() for ( index, cluster ) in enumerate( equiv.classes(), start = 1 ): context.output.verbose( msg = "Cluster %s: %s elements" % ( index, len( cluster ) ) ) assert 0 < len( cluster ) axis = reduce( operator.add, [ op.axis for ( op, rel ) in cluster ] ).normalize() centre = 1.0 / len( cluster ) * reduce( operator.add, [ op.centre for ( op, rel ) in cluster ] ) unmerged_folds = set( [ op.fold for ( op, rel ) in cluster ] ) while unmerged_folds: highest = max( unmerged_folds ) unmerged_folds.remove( highest ) generators_with = {} generators_with[ highest ] = [ rel for ( op, rel ) in cluster if op.fold == highest ] for of in list( unmerged_folds ): assert of < highest if highest % of == 0: unmerged_folds.remove( of ) generators_with[ of ] = [ rel for ( op, rel ) in cluster if op.fold == of ] aver = mr_object.GeneratorSetting( ensemble = reference.setting.ensemble, generator = mr_object.Generator.new( fold = self.determine_generator_fold( generators_with = generators_with ), axis = axis, centre = centre, ), ) members = reduce( operator.add, generators_with.values() ) context.output.verbose( msg = "Averaged operator: %s (from %s operators)" % ( aver.generator, len( members ) ) ) members_of[ aver ] = members return members_of def determine_generator_fold(self, generators_with): assert generators_with highest_count = max( [ len( v ) for v in generators_with.values() ] ) for key in sorted( generators_with, reverse = True ): if self.tolerances.min_fold_fraction <= len( generators_with[ key ] ) / highest_count: return key raise RuntimeError, "This statement should not be reached" def __str__(self): return "Generator classification and averaging" class PointGroupDetermine(object): """ Analyzes averaged generators and turn them into point groups Only active in post-processing """ def __init__(self, tolerances): self.tolerances = tolerances def __call__(self, board, case, context): context.output.verbose( msg = output.underlined( text = "Point group determination" ) ) context.output.notify( message = "status", data = "Determine point group symmetry" ) averageds = board.read( type = mr_object.GeneratorCluster ) context.output.verbose( msg = "%s generators found" % len( averageds ) ) analyses = [] for ( name, method ) in self.METHOD_FOR.items(): symmetries = method( self, clusters = averageds ) context.output.verbose( msg = "%s %s found" % ( len( symmetries ), name ) ) analyses.extend( symmetries ) context.output.verbose( msg = "%s point groups found" % len( analyses ) ) context.output.verbose( msg = "Merging subgroups into supergroups" ) highests = self.find_supergroups( analyses = analyses, context = context ) for obj in highests: context.output.verbose( msg = "Symmetry id: %s" % id( obj.setting.point_group ) ) board.write( obj = obj ) context.output.verbose( msg = output.section_close( text = "Point group determination finished" ) ) def find_rotation_axes(self, clusters): symmetries = [] for clust in clusters: symmetries.append( mr_object.SymmetryAnalysis( setting = mr_object.SymmetrySetting( ensemble = clust.average.ensemble, point_group = mr_object.RotationAxis( axis = clust.average.generator ), ), relations_from = clust.relations_from, formed_by = [ clust ], ) ) return symmetries def find_supergroups(self, analyses, context): remaining = set( analyses ) supergroups = [] while remaining: highest = max( remaining, key = lambda s: len( s.setting.point_group.series() ) ) context.output.verbose( msg = "Highest symmetry: %s" % highest.setting ) supergroups.append( highest ) remaining.remove( highest ) for trial in list( remaining ): context.output.verbose( msg = "Testing: %s" % trial.setting ) if highest.contains( other = trial ): context.output.verbose( msg = " CONTAINED" ) remaining.remove( trial ) else: context.output.verbose( msg = " NOT contained" ) return supergroups def __str__(self): return "Point group determination" METHOD_FOR = { "rotation axes": find_rotation_axes, } def find_symmetry_related_assemblies(relations): """ Determines which molecules are related by a specific symnetry """ ec = mr_object.EquivalenceRelationChain() for r in relations: ec.add( relation = r.relating ) return ec.classes class SymmetryCompare(object): """ Compares new symmetries with that has been located before """ def __init__(self, tolerances): self.tolerances = tolerances def __call__(self, board, case, context): context.output.verbose( msg = output.underlined( text = "Symmetry comparison" ) ) context.output.notify( message = "status", data = "Compare found and existing point group symmetries" ) known = case.data.read( type = mr_object.SymmetrySetting ) context.output.verbose( msg = "%s symmetries known" % len( known ) ) analyses = board.read( type = mr_object.SymmetryAnalysis ) for ( index, analysis ) in enumerate( analyses, start = 1 ): context.output.verbose( msg = "%s. symmetry" % index ) context.output.verbose( msg = "Compare with existing point groups:" ) context.output.indent() trial_contains = [] trial_contained_by = [] for existing in known: context.output.verbose( msg = "Compare: %s" % existing ) try: ( contained_by, contains ) = self.determine_containment_relationship( existing = existing, trial = analysis.setting, context = context, ) except RuntimeError, e: context.output.verbose( msg = " REJECTED: %s" % e ) if contained_by: context.output.verbose( msg = " KNOWN CONTAINS NOVEL" ) trial_contained_by.append( existing ) if contains: context.output.verbose( msg = " NOVEL CONTAINS KNOWN" ) trial_contains.append( existing ) context.output.dedent() relations = mr_object.SymmetrySetRelationAnalysis( analysis = analysis, contained_by = trial_contained_by, contains = trial_contains, ) board.write( obj = relations ) context.output.verbose( msg = output.section_close( text = "Symmetry comparison finished" ) ) def determine_containment_relationship(self, existing, trial, context): sup = get_superposition( reference = existing.ensemble, moving = trial.ensemble, superposition = context.superposition, stream = context.output.debug ) tr = trial.point_group.transformed( operation = sup.transformation ) contained_by = existing.point_group.approx_contains( other = tr, min_cos_angle = self.tolerances.min_axis_alignment, max_spatial_disp = self.tolerances.max_spatial_error, ) contains = tr.approx_contains( other = existing.point_group, min_cos_angle = self.tolerances.min_axis_alignment, max_spatial_disp = self.tolerances.max_spatial_error, ) return ( contained_by, contains ) class SymmetryCreationPolicy(object): """ Creates novel point groups """ def __call__(self, board, case, context): context.output.verbose( msg = output.underlined( text = "Symmetry creation" ) ) context.output.notify( message = "status", data = "Rationalise point groups" ) context.output.verbose( msg = "Create novel point groups" ) supergroups = [ sr for sr in board.read( type = mr_object.SymmetrySetRelationAnalysis ) if not sr.contained_by ] context.output.verbose( msg = "%s novel point groups identified" % len( supergroups ) ) for sra in supergroups: case.data.write( obj = sra.analysis.setting ) assemblies_in = dict( [ ( struct, find_symmetry_related_assemblies( relations = rels ) ) for ( struct, rels ) in sra.analysis.relations_from.items() ] ) board.write( obj = mr_object.NewSymmetry( analysis = sra, assemblies_in = assemblies_in, ) ) context.output.verbose( msg = "%s: SAVED" % sra ) context.output.verbose( msg = output.section_close( text = "Symmetry creation finished" ) ) class SymmetrySubgroupPolicy(object): """ Controls what to do with point groups that are not the highest available """ def __call__(self, board, case, context): context.output.verbose( msg = output.underlined( text = "Symmetry update" ) ) context.output.notify( message = "status", data = "Rationalise point groups" ) context.output.verbose( msg = "Remove subgroups of new symmetry" ) supergroups = board.read( type = mr_object.NewSymmetry ) context.output.verbose( msg = "%s new symmetries created" % len( supergroups ) ) for ( index, ns ) in enumerate( supergroups, start = 1 ): context.output.verbose( msg = "%s. %s" % ( index, ns ) ) context.output.indent() self.remove( symmetry = ns, board = board, case = case, context = context ) context.output.dedent() context.output.verbose( msg = output.section_close( text = "Symmetry update finished" ) ) @staticmethod def remove(symmetry, board, case, context): moribunds = set( symmetry.analysis.contains ) presents = set( [ ss for ss in case.data.read( type = mr_object.SymmetrySetting ) if ss in moribunds ] ) for ss in symmetry.analysis.contains: if ss in presents: case.data.delete( obj = ss ) board.write( obj = mr_object.ObsoletedBySuperSymmetry( setting = ss, supersymmetry = symmetry.analysis.analysis.setting, ) ) context.output.verbose( msg = "DELETED: %s" % ss ) class SymmetryUnsupportedPolicy(object): """ Controls what to do with point groups that are too high, i.e. not found in current best models """ def __call__(self, board, case, context): context.output.verbose( msg = output.underlined( text = "Symmetry removal - unsupported" ) ) context.output.notify( message = "status", data = "Remove point groups with symmetry not found in current solutions" ) context.output.verbose( msg = "Remove supergroups not supported by data" ) ssras = board.read( type = mr_object.SymmetrySetRelationAnalysis ) subgroups = set( reduce( operator.add, [ s.contains for s in ssras ], [] ) ) supergroup_of = {} for rel in ssras: for symm in rel.contained_by: supergroup_of.setdefault( symm, set() ).add( rel ) supergroups = [ ss for ss in case.data.read( type = mr_object.SymmetrySetting ) if ss not in subgroups ] context.output.verbose( msg = "%s symmetries not detected" % len( supergroups ) ) for ( index, ss ) in enumerate( supergroups, start = 1 ): case.data.delete( obj = ss ) context.output.verbose( msg = "%s: REMOVED" % ss ) if ss in supergroup_of: highest = max( [ ra.analysis.setting for ra in supergroup_of[ ss ] ], key = lambda s: len( s.point_group.series() ) ) context.output.verbose( msg = "Highest subgroup symmetry: %s" % highest ) case.data.write( obj = highest ) else: highest = None context.output.verbose( msg = "No subgroup symmetry detected" ) board.write( obj = mr_object.SymmetryDowngrade( setting = ss, highest = highest ) ) context.output.verbose( msg = output.section_close( text = "Symmetry removal finished" ) ) class SymmetryUpdate(object): """ Handles updates to existing symmetries """ def __init__(self, tolerances): self.tolerances = tolerances self.max_spatial_error = 1.0 self.max_misalignment = 0.99 def __call__(self, board, case, context): context.output.verbose( msg = output.underlined( text = "Symmetry update" ) ) context.output.notify( message = "status", data = "Update existing symmetries" ) analyses_for = {} for ssra in board.read( type = mr_object.SymmetrySetRelationAnalysis ): for symm in ssra.contained_by + ssra.contains: analyses_for.setdefault( symm, set() ).add( ssra.analysis ) context.output.verbose( msg = "Up to %s symmetries could be affected" % len( analyses_for ) ) for ( index, ss ) in enumerate( case.data.read( type = mr_object.SymmetrySetting ), start = 1 ): context.output.verbose( "%s. %s" % ( index, ss ) ) # This skips update for novel symmetries identified in this round as well if ss not in analyses_for: context.output.verbose( msg = "No updateable elements" ) continue context.output.indent() relateds = reduce( operator.or_, [ a.formed_by for a in analyses_for[ ss ] ], set(), ) updated = self.updated_point_group( symmetry = ss, relateds = relateds, context = context, ) context.output.dedent() if updated == ss.point_group: context.output.verbose( msg = "Skipping this symmetry" ) continue context.output.verbose( msg = "Replacing original with updated symmetry" ) relations_from = {} for analysis in analyses_for[ ss ]: for ( structure, relations ) in analysis.relations_from.items(): relations_from.setdefault( structure, [] ).extend( relations ) assemblies_in = dict( [ ( struct, find_symmetry_related_assemblies( relations = rels ) ) for ( struct, rels ) in relations_from.items() ] ) previous = ss.point_group board.write( obj = mr_object.SymmetryUpdateRecord( setting = ss, previous = ss.point_group, assemblies_in = assemblies_in ) ) ss.point_group = updated context.output.verbose( msg = output.section_close( text = "Symmetry update finished" ) ) def updated_point_group(self, symmetry, relateds, context): context.output.verbose( msg = "%s related generators found" % len( relateds ) ) context.output.verbose( msg = "Transforming to common reference" ) transformeds = [] for r in relateds: # Not sure whether this can go wrong try: sup = get_superposition( reference = symmetry.ensemble, moving = r.average.ensemble, superposition = context.superposition, stream = context.output.debug ) except RuntimeError, e: continue transformeds.append( r.average.generator.transformed( operation = sup.transformation ) ) context.output.verbose( msg = "%s generators transformed" % len( transformeds ) ) context.output.verbose( msg = "Matching with point group generator classes" ) bases = [] for ( index, gc ) in enumerate( symmetry.point_group.generator_classes(), start = 1 ): context.output.verbose( msg = "%s. class" % index ) current = [] for gen in transformeds: searchgen = ( op for ( g, op ) in zip( gc.generators, gc.operations ) if gen.approx_equal_spatially( other = g, max_misalignment = self.max_misalignment, max_spatial_error = self.max_spatial_error, ) ) try: stdop = searchgen.next() except StopIteration: continue current.append( gen.transformed( operation = stdop ) ) context.output.verbose( msg = "%s generators match" % len( current ) ) if current: context.output.verbose( msg = "Averaging spatial components" ) average = mr_object.Generator( fold = gc.standard.fold, axis = reduce( operator.add, [ op.axis for op in current ] ).normalize(), centre = 1.0 / len( current ) * reduce( operator.add, [ op.centre for op in current ] ), ) else: context.output.verbose( msg = "Using original operator" ) average = gc.standard bases.append( average ) if bases == [ gc.standard for gc in symmetry.point_group.generator_classes() ]: context.output.verbose( msg = "No generator updated, keeping original point group" ) return symmetry.point_group context.output.verbose( msg = "Forming point group from bases" ) return symmetry.point_group.from_generators( generators = bases ) class SpatiallyEquivalentMatchingPeak(object): """ A matching-compatible class that takes care of superpose transformations """ def __init__( self, peak, superposition, stream, crystal, dmin, multiplier, cached = False, ): from phaser import matching self.peak = peak self.crystal = crystal self.dmin = dmin self.multiplier = multiplier if cached: detail_orientation = matching.cached_rotation_differences detail_position = matching.cached_position_differences else: detail_orientation = matching.simple_rotation_differences detail_position = matching.simple_position_differences self.spatial = matching.MRPeakPO.phaser_style( rotation = self.peak.rotation.matrix(), translation = self.peak.translation, crystal = self.crystal, ensemble = self.peak.ensemble.matching_ensemble_info(), dmin = dmin, multiplier = multiplier, detail_position = detail_position, detail_orientation = detail_orientation, ) self.superposition = superposition self.stream = stream def operations_between(self, other): try: sup = get_superposition( reference = self.peak.ensemble, moving = other.peak.ensemble, superposition = self.superposition, stream = self.stream ) except RuntimeError, e: raise StopIteration from phaser import matching superposed = ( other.peak.as_scitbx_rt_operator( cell = self.crystal.cell ) * sup.transformation ) transformed = matching.MRPeakPO.phaser_style( rotation = superposed.r.elems, translation = self.crystal.cell.fractionalize( superposed.t.elems ), crystal = self.crystal, ensemble = self.peak.ensemble.matching_ensemble_info(), dmin = self.dmin, multiplier = self.multiplier ) for op in self.spatial.operations_between( other = transformed ): yield op class SymmetryTolerances(object): """ Small container class to hold all tolerances """ def __init__( self, max_fold = 12, max_angular_error = 0.15, max_spatial_error = 1.0, min_axis_alignment = 0.99, max_relative_assembly_size = 3.0, min_fold_fraction = 0.3, ): self.max_fold = max_fold self.max_angular_error = max_angular_error self.max_spatial_error = max_spatial_error self.min_axis_alignment = min_axis_alignment self.max_relative_assembly_size = max_relative_assembly_size self.min_fold_fraction = min_fold_fraction class LocalSymmetryExpert(object): """ A wrapper class to contain all symmetry-related experts for scheduling """ def __init__( self, add_new_symmetry = True, remove_subgroup_symmetry = True, remove_unsupported_symmetry = True, update_known_symmetry = True, ): self.tolerances = SymmetryTolerances() self.actions = [ GeneratorSearch( tolerances = self.tolerances ), GeneratorClassifyAverage( tolerances = self.tolerances ), PointGroupDetermine( tolerances = self.tolerances ), SymmetryCompare( tolerances = self.tolerances ), ] if remove_unsupported_symmetry: self.actions.append( SymmetryUnsupportedPolicy() ) if update_known_symmetry: self.actions.append( SymmetryUpdate( tolerances = self.tolerances ) ) if add_new_symmetry: self.actions.append( SymmetryCreationPolicy() ) if remove_subgroup_symmetry: self.actions.append( SymmetrySubgroupPolicy() ) def __call__(self, board, case, context, solutions, postprocessings): context.output.verbose( msg = output.subtitle( text = "Local symmetry analysis" ) ) context.output.notify( message = "status", data = "Analyse solutions to identify local symmetries" ) context.output.verbose( msg = output.underlined( text = "Structure selection" ) ) selected = [ s for s in ( solutions + postprocessings ) if s.structure.significant ] for ( index, sol ) in enumerate( selected, start = 1 ): context.output.verbose( msg = "%d. Structure:" % index ) context.output.indent() context.output.verbose( msg = str( sol ) ) context.output.dedent() board.write( obj = mr_object.ForSymmetryAnalysis( structure = sol.structure, index = index, ) ) for action in self.actions: action( board = board, case = case, context = context ) context.output.verbose( msg = output.section_close( text = "Local symmetry analysis finished" ) ) return [] def partition_structure_according_to_membership_in( structure, assembly, case, context, multiplier, ): from phaser import matching cell = case.uctbx_unit_cell() crystal = case.matching_crystal_info() tmfs = [ SpatiallyEquivalentMatchingPeak( peak = peak, superposition = context.superposition, stream = context.output.debug, crystal = crystal, dmin = case.problem.xray_data.resolution, multiplier = multiplier, cached = True, ) for peak in structure.peaks ] remaining = set( structure.peaks ) peak_for = dict( zip( tmfs, structure.peaks ) ) tmf_for = dict( zip( structure.peaks, tmfs ) ) associations = [] while remaining: head = remaining.pop() alternatives = [ [] ] for ( index, anchor ) in enumerate( assembly.anchors, start = 1 ): context.output.verbose( msg = "Anchor %s" % index ) try: sup = get_superposition( reference = head.ensemble, moving = anchor.ensemble, superposition = context.superposition, stream = context.output.debug, ) except RuntimeError, e: context.output.verbose( msg = " REJECTED: %s" % e ) continue transformation = ( head.as_scitbx_rt_operator( cell = cell ) * sup.transformation ) trials = [] for ( ense, op ) in anchor.peaks: final = transformation * op peak = mr_object.Peak( ensemble = ense, rotation = mr_object.Rotation.Matrix( elements = final.r.elems ), translation = cell.fractionalize( final.t ), bfactor = peak.bfactor ) ptmf = SpatiallyEquivalentMatchingPeak( peak = peak, superposition = context.superposition, stream = context.output.debug, crystal = crystal, dmin = case.problem.xray_data.resolution, multiplier = multiplier, ) trials.append( ptmf ) overlaps = matching.overlap( left = [ tmf_for[ p ] for p in remaining ], right = trials, ) assoc = [ l for ( l, r ) in matching.pairs( overlaps = overlaps ) ] context.output.verbose( msg = " members: %s" % ( len( assoc ) + 1 ) ) alternatives.append( assoc ) largest = [ peak_for[ t ] for t in max( alternatives, key = len ) ] context.output.verbose( msg = "Largest association has %s members" % ( len( largest ) + 1 ) ) for peak in largest: remaining.remove( peak ) associations.append( [ head ] + largest ) context.output.verbose( msg = "Number of associations found: %s" % len( associations ) ) return associations def AcceptAssemblyStrategy(assembly, structures, case, context): """ Accepts any assembly without any checks """ context.output.verbose( msg = "Accepting assembly (no condition specified)" ) return True def RejectAssemblyStrategy(assembly, structures, case, context): """ Rejects any assembly without any checks """ context.output.verbose( msg = "Rejecting assembly (no condition specified)" ) return False class AcceptObservedAssemblyStrategy(object): """ Accepts if the assembly has been observed in any structures """ def __init__(self, multiplier = 3.0): self.multiplier = multiplier def __call__(self, assembly, structures, case, context): context.output.verbose( msg = "Accepting symmetry if full assembly present" ) member_count = len( assembly.description ) context.output.verbose( msg = "Full assembly has %s members" % member_count ) context.output.verbose( msg = "Testing membership in %s structures" % len( structures ) ) for ( index, structure ) in enumerate( structures, start = 1 ): associations = partition_structure_according_to_membership_in( structure = structure, assembly = assembly, case = case, context = context, multiplier = self.multiplier, ) largest = max( [ len( a ) for a in associations ] ) context.output.verbose( msg = "Structure %s: largest assembly has %s members" % ( index, largest, ) ) if member_count <= largest: context.output.verbose( msg = "Full assembly found" ) return True context.output.verbose( msg = "Full assembly NOT found" ) return False class AssemblyCreate(object): """ Creates models with symmetries """ def __init__(self, is_valid): self.is_valid = is_valid def __call__(self, board, case, context, solutions, postprocessings): context.output.verbose( msg = output.subtitle( text = "Assembly creation" ) ) context.output.notify( message = "status", data = "Create assemblies from identified symmetries" ) useds = set( [ a.symmetry for a in case.data.read( type = mr_object.Assembly ) ] ) context.output.verbose( msg = "Current assemblies use %s symmetries" % len( useds ) ) unuseds = [ s for s in case.data.read( type = mr_object.SymmetrySetting ) if s.point_group not in useds ] context.output.verbose( msg = "%s symmetries not used in assemblies are present" % len( unuseds ) ) structures = [ s.structure for s in ( solutions + postprocessings ) if s.structure.significant ] for ( index, setting ) in enumerate( unuseds, start = 1 ): context.output.verbose( msg = "%s. %s" % ( index, setting ) ) assembly = mr_object.Assembly.Simple( ensemble = setting.ensemble, symmetry = setting.point_group, strategy = mr_object.RefineStrategy(), ) if self.is_valid( assembly = assembly, structures = structures, case = case, context = context, ): context.output.verbose( msg = " ACCEPTED" ) case.data.write( obj = assembly ) board.write( obj = mr_object.NewAssemblyRecord( assembly = assembly ) ) else: context.output.verbose( msg = " REJECTED" ) context.output.verbose( msg = output.section_close( text = "Assembly creation finished" ) ) return [] @classmethod def Always(cls): return cls( is_valid = AcceptAssemblyStrategy ) @classmethod def Observed(cls): return cls( is_valid = AcceptObservedAssemblyStrategy() ) @classmethod def Never(cls): return cls( is_valid = RejectAssemblyStrategy ) class AssemblyRemove(object): """ Removes models with obsoleted symmetries """ def __call__(self, board, case, context, solutions, postprocessings): context.output.verbose( msg = output.subtitle( text = "Assembly removal" ) ) context.output.notify( message = "status", data = "Remove assemblies with obsoleted symmetries" ) context.output.verbose( msg = output.underlined( text = "Subgroup symmetries" ) ) self.remove_subgroup_symmetries( board = board, case = case, context = context, ) context.output.verbose( msg = output.underlined( text = "Supergroup symmetries" ) ) self.remove_supergroup_symmetries( board = board, case = case, context = context, ) context.output.verbose( msg = output.section_close( text = "Assembly removal finished" ) ) return [] def remove_subgroup_symmetries(self, board, case, context): obsoletes = set( [ n.setting.point_group for n in board.read( type = mr_object.ObsoletedBySuperSymmetry ) ] ) context.output.verbose( msg = "%s symmetries obsoleted" % len( obsoletes ) ) for assem in case.data.read( type = mr_object.Assembly ): if assem.symmetry in obsoletes: context.output.verbose( msg = "%s: symmetry obsoleted, remove" % assem ) case.data.delete( obj = assem ) board.write( obj = mr_object.AssemblyRemovalRecord( assembly = assem ) ) def remove_supergroup_symmetries(self, board, case, context): downgraded_to = dict( [ ( dg.setting.point_group, dg.highest ) for dg in board.read( type = mr_object.SymmetryDowngrade ) ] ) context.output.verbose( msg = "%s symmetries downgraded" % len( downgraded_to ) ) for assem in case.data.read( type = mr_object.Assembly ): if assem.symmetry in downgraded_to: context.output.verbose( msg = "%s: symmetry obsoleted, remove" % assem ) case.data.delete( obj = assem ) board.write( obj = mr_object.AssemblyRemovalRecord( assembly = assem ) ) replacement = downgraded_to[ assem.symmetry ] if replacement is not None: context.output.verbose( msg = "Reinstating subgroup: %s" % replacement ) r_assem = mr_object.Assembly( definitions = assem.definitions, symmetry = replacement.point_group, strategy = assem.strategy, ) case.data.write( obj = r_assem ) board.write( obj = mr_object.NewAssemblyRecord( assembly = r_assem ) ) class AssemblyUpdate(object): """ Updates models with symmetries """ def __call__(self, board, case, context, solutions, postprocessings): context.output.verbose( msg = output.subtitle( text = "Assembly update" ) ) context.output.notify( message = "status", data = "Update symmetry for assemblies" ) context.output.verbose( msg = "Update assemblies according to changes made to symmetries" ) update_for = dict( [ ( u.previous, u.setting.point_group ) for u in board.read( type = mr_object.SymmetryUpdateRecord ) ] ) context.output.verbose( msg = "%s updates available" % len( update_for ) ) for assem in case.data.read( type = mr_object.Assembly ): if assem.symmetry in update_for: context.output.verbose( msg = "%s: UPDATED" % assem ) assem.symmetry = update_for[ assem.symmetry ] context.output.verbose( msg = output.section_close( text = "Assembly update finished" ) ) return [] class SearchModelMaintenance(object): """ Creates/removes AssemblyEnsembles as dictated by Assemblies """ def __init__(self, is_valid): self.is_valid = is_valid def __call__(self, board, case, context, solutions, postprocessings): context.output.verbose( msg = output.subtitle( text = "Search model assembly and maintenance" ) ) context.output.notify( message = "status", data = "Creates and search models from assemblies" ) context.output.verbose( msg = output.underlined( text = "New assemblies" ) ) structures = [ s.structure for s in ( solutions + postprocessings ) if s.structure.significant ] self.add_new_assemblies( structures = structures, board = board, case = case, context = context, ) context.output.verbose( msg = output.underlined( text = "Obsoleted assemblies" ) ) self.remove_deleted_assemblies( board = board, case = case, context = context, ) context.output.verbose( msg = output.section_close( text = "Search model maintenance finished" ) ) return [] def add_new_assemblies(self, structures, board, case, context): useds = set( [ e.assembly for e in case.data.read( type = mr_object.AssemblyEnsemble ) ] ) context.output.verbose( msg = "Current models use %s assemblies" % len( useds ) ) unuseds = [ a for a in case.data.read( type = mr_object.Assembly ) if a not in useds ] context.output.verbose( msg = "%s assemblies not used as models are present" % len( unuseds ) ) for ( index, assembly ) in enumerate( unuseds, start = 1 ): context.output.verbose( msg = "%s. %s" % ( index, assembly ) ) if self.is_valid( assembly = assembly, structures = structures, case = case, context = context, ): context.output.verbose( msg = " ACCEPTED" ) case.data.write( obj = mr_object.AssemblyEnsemble( assembly = assembly ) ) else: context.output.verbose( msg = " REJECTED" ) def remove_deleted_assemblies(self, board, case, context): deleteds = set( r.assembly for r in board.read( type = mr_object.AssemblyRemovalRecord ) ) context.output.verbose( msg = "%s assemblies deleted" % len( deleteds ) ) for ense in case.data.read( type = mr_object.AssemblyEnsemble ): if ense.assembly in deleteds: context.output.verbose( msg = "%s: removed search model" % ense ) case.data.delete( obj = ense ) @classmethod def Always(cls): return cls( is_valid = AcceptAssemblyStrategy ) @classmethod def Observed(cls): return cls( is_valid = AcceptObservedAssemblyStrategy() ) @classmethod def Never(cls): return cls( is_valid = RejectAssemblyStrategy ) class AssemblyModelCreate(object): """ Creates aggregate models from assemblies Only active in post-processing """ def __init__(self, observe_first): if observe_first: self.accept = self.accept_if_observed else: self.accept = self.accept_anyway self.min_cos_angle = 0.99 self.max_spatial_disp = 1.0 def __call__(self, board, case, context, solutions, postprocessings): context.output.verbose( msg = output.subtitle( text = "Model assembly" ) ) context.output.notify( message = "status", data = "Create models according to assemblies" ) known = case.data.read( type = mr_object.RegularEnsemble ) assemblies = case.data.read( type = mr_object.Assembly ) context.output.verbose( msg = "Assemblies: %s" % len( assemblies ) ) for assem in assemblies: context.output.verbose( msg = "Processing: %s" % assem ) context.output.verbose( msg = output.underlined( text = "Comparing with existing models" ) ) for reg in known: context.output.verbose( msg = "Trial:" ) context.output.verbose( msg = reg ) if not self.approx_contains( left = assem, right = reg.assembly, context = context ): continue context.output.verbose( msg = "There is a model for this assembly: UPDATE" ) reg.assembly = assem break else: if self.accept( assembly = assem, context = context ): reg_ense = mr_object.RegularEnsemble( assembly = assem ) case.data.write( obj = reg_ense ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = output.section_close( text = "Model assembly finished" ) ) return [] def accept_anyway(self, assembly, context): context.output.verbose( msg = "Novel assembly: CREATE model" ) return True def accept_if_observed(self, assembly, context): if assembly.present: context.output.verbose( msg = "Novel assembly (observed): CREATE model" ) return True else: context.output.verbose( msg = "Novel assembly (unobserved): DISCARD" ) return False def approx_contains(self, left, right, context): left_desc = left.description right_desc = right.description if len( left_desc ) < len( right_desc ): context.output.verbose( msg = "Containing assembly smaller" ) return False context.output.verbose( msg = "Finding seed:" ) ( l_ense, l_op ) = left_desc[0] for ( r_index, ( r_ense, r_op ) ) in enumerate( right_desc, start = 1 ): try: sup = self.ensemble_relation( left = l_ense, right = r_ense, superposition = context.superposition, stream = context.output.debug, ) except RuntimeError, e: continue context.output.verbose( msg = "Seed: 1 vs %s" % r_index ) match = ( ( ( l_ense, l_op ), ( r_ense, r_op ) ), left_desc[1:], [ ri for ri in right_desc if ri != ( r_ense, r_op ) ], sup.transformation, ) for ( l, r , s ) in self.assign( match = match, context = context ): context.output.verbose( msg = "Match: %s vs %s" % ( left_desc.index( l ) + 1 if l else "n/a", right_desc.index( r ) + 1 if r else "n/a", ) ) if l is None: context.output.verbose( msg = "Unmatched component in larger assembly, dropping seed" ) break else: context.output.verbose( msg = "All components match" ) return True return False def approx_identity_operation(self, op): return mr_object.PointGroup.approx_identity_operation( operation = op, min_cos_angle = self.min_cos_angle, max_spatial_disp = self.max_spatial_disp, ) def ensemble_relation(self, left, right, superposition, stream): if not left.substituteable_with( other = right ): raise RuntimeError, "different composition" return superposition.superpose( reference = left, moving = right, stream = stream, ) def assign(self, match, context): ( ( ( l_ense, l_op ), ( r_ense, r_op ) ), left, right, sup ) = match yield ( ( l_ense, l_op ), ( r_ense, r_op ), sup ) transf = l_op * ( r_op * sup ).inverse() r_transf = [ ( e, transf * op ) for ( e, op ) in right ] for ( l, r, s_op ) in self.overlap( left = left, right = r_transf, context = context ): if r is not None: assert r in r_transf r = right[ r_transf.index( r ) ] yield ( l, r, s_op ) def overlap(self, left, right, context): copy = set( right ) for ( ref_ense, ref_op ) in left: inverse = ref_op.inverse() for ( trial_ense, trial_op ) in copy: try: sup = self.ensemble_relation( left = ref_ense, right = trial_ense, superposition = context.superposition, stream = context.output.debug, ) except RuntimeError, e: continue if self.approx_identity_operation( op = inverse * trial_op * sup.transformation ): copy.remove( ( trial_ense, trial_op ) ) yield ( ( ref_ense, ref_op ), ( trial_ense, trial_op ), sup.transformation ) break else: yield ( ( ref_ense, ref_op ), None, None ) for r in copy: yield ( None, r, None ) def __str__(self): return "Model assembly" class ObjectCompletion(object): """ Adds molecules to incomplete assemblies """ def __init__(self): self.multiplier = 3.0 def __call__(self, board, case, context, solutions, postprocessings): context.output.verbose( msg = output.subtitle( text = "Object completion - identification" ) ) context.output.notify( message = "status", data = "Find assemblies that can be completed" ) assemblies = ( case.data.read( type = mr_object.Assembly ) + case.problem.data.read( type = mr_object.Assembly ) ) significants = [ s for s in solutions if s.structure.significant ] context.output.verbose( msg = "Significant peaks: %d" % len( significants ) ) context.output.verbose( msg = output.heading( text = "Processing:" ) ) failures = [ ( rec.peak.ensemble, case.create_matching_mr_peak( peak = rec.peak, multiplier = self.multiplier, ), case.create_matching_mr_structure( structure = rec.structure ), ) for rec in case.data.read( type = mr_object.FillFailureRecord ) ] results = [] for sol in significants: generateds = self.complete_assemblies_in( structure = sol.structure, assemblies = assemblies, case = case, context = context ) # Check whether these peaks have not been tried before trial_structure_tmf = case.create_matching_mr_structure( structure = sol.structure, cached = True ) accepteds = [] for ( peak, strat ) in generateds: if not self.is_known_failure( structure = trial_structure_tmf, peak = peak, failures = failures, case = case, context = context, ): accepteds.append( ( peak, strat ) ) else: context.output.verbose( msg = peak ) context.output.verbose( msg = " REJECT: know failure" ) context.output.verbose( msg = "%s missing molecules filled in" % len( accepteds ) ) for ( peak, strategy ) in accepteds: roc = mr_object.ObjectCompletion( structure = sol.structure, peak = peak, strategy = strategy, ) context.output.verbose( msg = peak ) results.append( roc ) context.output.verbose( msg = "" ) context.output.verbose( msg = output.section_close( text = "Object completion finished" ) ) return results def process_assembly(self, peak, assembly, cell, context): generateds = [] for anchor in assembly.anchors: context.output.verbose( msg = anchor ) if not peak.ensemble.substituteable_with( other = anchor.ensemble ): context.output.verbose( msg = " REJECTED: different composition" ) context.output.verbose( msg = output.blank() ) continue try: sup = context.superposition.superpose( reference = peak.ensemble, moving = anchor.ensemble, stream = context.output.debug ) except RuntimeError, e: context.output.verbose( msg = " REJECTED: no superposition" ) context.output.verbose( msg = output.blank() ) continue context.output.verbose( msg = "Generating possible molecules using current anchor" ) transformation = ( peak.as_scitbx_rt_operator( cell = cell ) * sup.transformation ) currents = [] for ( ense, op ) in anchor.peaks: final = transformation * op currents.append( ( mr_object.Peak( ensemble = ense, rotation = mr_object.Rotation.Matrix( elements = final.r.elems ), translation = cell.fractionalize( final.t ), bfactor = peak.bfactor ), assembly.strategy, ) ) context.output.verbose( msg = " %s peak generated" % len( currents ) ) context.output.verbose( msg = output.blank() ) generateds.extend( currents ) context.output.verbose( msg = "Total %s peaks generated from this assembly" % len( generateds ) ) return generateds def complete_assemblies_in(self, structure, assemblies, case, context): context.output.verbose( msg = "Structure:" ) context.output.indent() context.output.verbose( msg = structure ) context.output.dedent() context.output.verbose( msg = output.blank() ) equiv = mr_object.EquivalenceRelation( relation = lambda left, right: any( left[1].operations_between( other = right[1] ) ) ) for peak in structure.peaks: context.output.verbose( msg = "Trial peak:" ) context.output.indent() context.output.verbose( msg = peak ) context.output.dedent() context.output.verbose( msg = output.blank() ) context.output.verbose( msg = "Processing assemblies:" ) context.output.verbose( msg = output.blank() ) currents = [] for ( index, assem ) in enumerate( assemblies, start = 1 ): context.output.verbose( msg = "%d: %s" % ( index, assem ) ) context.output.indent() currents.extend( self.process_assembly( peak = peak, assembly = assem, cell = case.uctbx_unit_cell(), context = context ) ) context.output.verbose( msg = output.blank() ) context.output.dedent() context.output.verbose( msg = "Total %s peaks generated for trial peak" % len( currents ) ) context.output.verbose( msg = output.blank() ) for ( p, s ) in currents: tmf = case.create_matching_mr_peak( peak = p, multiplier = self.multiplier ) equiv.add( value = ( ( p, s ), tmf ) ) context.output.verbose( msg = "Total %s peaks generated for this structure" % len( equiv.values ) ) context.output.verbose( msg = output.blank() ) averageds = [] presents = case.create_matching_mr_structure( structure = structure, cached = True, multiplier = self.multiplier ) structure_ensembles = set( [ p.ensemble for p in structure.peaks ] ) import scitbx.math import scitbx.matrix from phaser import matching context.output.verbose( msg = "Merging quasi-equivalent peaks" ) for ( index, cluster ) in enumerate( equiv.classes(), start = 1 ): assert cluster context.output.verbose( msg = "Cluster %d: %d elements" % ( index, len( cluster ) ) ) # Reduce equivalents to the same asymmetric unit ( ( ref, strat ), rtmf ) = cluster.pop() ref_posi = rtmf.positional() ref_ori = rtmf.spatial.orientation rots = [ ref.rotation.matrix() ] tras = [ scitbx.matrix.col( ref.translation ) ] bfacs = [ ref.bfactor ] strats = [ strat ] for ( ( p, s ), tmf ) in cluster: # Equivalence does not mean these are equal! This does not work: # ( symop, pgop ) = rtmf.operations_between( other = tmf ).next() # Has to be separated: ( symop, delta ) = min( ref_posi.orthogonal_differences_between( other = tmf.positional() ), key = lambda p: scitbx.matrix.col( p[1] ).length_sq() ) ( pgop, rotmat ) = max( ref_ori.rotation_differences_between( other = tmf.spatial.orientation ), key = lambda p: scitbx.math.r3_rotation_cos_rotation_angle_from_matrix( p[1] ) ) rots.append( symop.orth_rotation * scitbx.matrix.sqr( p.rotation.matrix() ) * pgop.r ) t = scitbx.matrix.col( symop.sgop * p.translation ) tras.append( t + scitbx.matrix.col( [ round( n ) for n in ( tras[0] - t ) ] ) ) bfacs.append( p.bfactor ) strats.append( s ) rot = scitbx.math.r3_rotation_average_rotation_matrix_from_matrices( *rots ) tra = 1.0 / len( tras ) * reduce( operator.add, tras ) bfac = 1.0 / len( bfacs ) * sum( bfacs ) strategy = max( strats, key = lambda s: s.coverage() ) trial = mr_object.Peak( ensemble = ref.ensemble, rotation = mr_object.Rotation.Matrix( elements = rot ), translation = tra.elems, bfactor = bfac ) context.output.verbose( msg = "Averaged peak:" ) context.output.verbose( msg = trial ) for ense in structure_ensembles: if not trial.ensemble.substituteable_with( other = ense ): continue try: sub_peak = context.superposition.substitute( peak = trial, ensemble = ense, cell = case.uctbx_unit_cell(), stream = context.output.debug ) except RuntimeError: continue overlaps = matching.overlap( left = presents, right = [ case.create_matching_mr_peak( peak = sub_peak, multiplier = self.multiplier ) ] ) if any( matching.pairs( overlaps = overlaps ) ): context.output.verbose( msg = "REJECT: present in known structure" ) context.output.verbose( msg = output.blank() ) break else: context.output.verbose( msg = "ACCEPT" ) context.output.verbose( msg = output.blank() ) averageds.append( ( trial, strategy ) ) context.output.verbose( msg = "Total %s peaks generated for this structure" % len( averageds ) ) return averageds def is_known_failure(self, structure, peak, failures, case, context): from phaser import matching for ( ense, p_tmf, s_tmf ) in failures: if not peak.ensemble.substituteable_with( other = ense ): continue try: sub_peak = context.superposition.substitute( peak = peak, ensemble = ense, cell = case.uctbx_unit_cell(), stream = context.output.debug ) except RuntimeError: continue tmf = case.create_matching_mr_peak( peak = sub_peak, multiplier = self.multiplier, ) if not any( tmf.operations_between( other = p_tmf ) ): continue overlap = matching.pairs( overlaps = matching.overlap( left = structure, right = s_tmf ) ) if len( list( overlap ) ) == len( s_tmf ): return True return False def __str__(self): return "Object completion" class Coordinator(object): """ Manages suspensions """ def inspect(self, board, case, context): return [ s for s in board.read( type = stage.Suspension ) if s.condition() ] def manipulate(self, candidate, board, case, context): context.output.verbose( msg = output.subtitle( text = "Process coordination" ) ) context.output.notify( message = "status", data = "Process coordination" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = output.heading( text = "Expired suspensions:" ) ) for ( index, suspension ) in enumerate( candidate, start = 1 ): context.output.verbose( "%s: %s" % ( index, suspension ) ) board.delete( obj = suspension ) context.output.verbose( msg = output.section_close( text = "Process coordination finished" ) ) def __str__(self): return "Coordinator" class Ensembler(object): """ Runs ensembler """ def __init__(self): self.inspect = SelectBest( expert = self, type = stage.ModelCollection ) from phaser import ensembler self.params = ensembler.PHIL_MASTER.extract() def manipulate(self, candidate, board, case, context): candidate.processed_by.add( self ) context.output.verbose( msg = output.subtitle( text = "Model generation - ensembler" ) ) context.output.notify( message = "status", data = "Model generation - ensembler" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = str( candidate ) ) context.output.verbose( msg = "" ) ensemble = self.process( collection = candidate.collection, case = case, context = context ) candidate.collection.ensemble = ensemble context.output.verbose( msg = "Ensemble #%d generated" % candidate.collection.index ) initial = candidate.initial() existing = [ es for es in board.read( type = stage.Ensemble ) if es.initial() == initial ] es = stage.Ensemble( initial = initial, intype = ensemble, index = len( existing ) + 1 ) board.write( obj = es ) context.output.verbose( msg = output.section_close( text = "Ensembler finished" ) ) def process(self, collection, case, context): if collection.ensemble is not None: context.output.verbose( msg = "This entry has been processed" ) return collection.ensemble case.problem.data.delete( obj = collection ) logger = output.SingleStream( stream = context.output.debug ) pdbs = collection.pdb_objects() self.params.configuration.trim = collection.trim from phaser import ensembler try: infos = ensembler.process( pdbs = pdbs, alignments = [], params = self.params, logger = logger ) except RuntimeError, e: context.output.verbose( msg = "Ensembler error: %s" % e ) raise Sorry, "Ensembler: %s" % e # Copy chains chains = [ info.chain.detached_copy() for info in infos ] # Transform chains for ( info, chain ) in zip( infos, chains ): sites = ( info.ensembler_rotation.elems * chain.atoms().extract_xyz() + info.ensembler_translation ) chain.atoms().set_xyz( sites ) if self.params.configuration.trim: assert all( hasattr( info, "ensembler_trimming_mask" ) for info in infos ) logger.info( msg = "Trimming chains" ) logger.info( msg = output.blank() ) ensembler.trim_chains( infos = infos, chains = chains, logger = logger ) # Write out files try: file_names = ensembler.write_chain_and_model_sorted_output_files( root = "ensembler_%d" % collection.index, infos = infos, chains = chains, logger = logger ) except IOError, e: raise Sorry, e assert len( file_names ) == len( infos ) rmsd_for = dict( zip( collection.pdbs, collection.rmsds ) ) assert all( info.pdb.name in rmsd_for for info in infos ) chain_rmsds = [ rmsd_for[ info.pdb.name ] for info in infos ] ensemble = mr_object.MultifileEnsemble( pdbs = file_names, rmsds = chain_rmsds, composition = collection.composition ) case.problem.data.write( obj = ensemble ) return ensemble def __str__(self): return "ModelGeneration (ensembler)" class Sculptor(object): """ Runs sculptor """ def __init__(self, protocols = [ "all" ]): self.inspect = SelectBest( expert = self, type = stage.ModelTemplate ) from phaser import sculptor self.params = sculptor.PHIL_MULTIPROTOCOL_MASTER.extract() self.params.macromolecule.protocols = protocols self.sculptors = sculptor.process_multiprotocol_phil( params = self.params, logger = output.SingleStream() ) self.discarder = sculptor.get_discarder_object( params = None ) def manipulate(self, candidate, board, case, context): candidate.processed_by.add( self ) context.output.verbose( msg = output.subtitle( text = "Model generation - sculptor" ) ) context.output.notify( message = "status", data = "Model generation - sculptor" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = str( candidate ) ) context.output.verbose( msg = "" ) context.output.verbose( msg = "Number of active protocols: %s" % len( self.sculptors ) ) ensembles = self.process_and_uniqueify( template = candidate.template, case = case, context = context ) candidate.template.ensembles = ensembles context.output.verbose( msg = "Unique models generated:" ) context.output.indent() initial = candidate.initial() existing = [ es for es in board.read( type = stage.Ensemble ) if es.initial() == initial ] for ( i, ensemble ) in enumerate( ensembles, start = len( existing ) + 1 ): context.output.verbose( msg = ensemble ) es = stage.Ensemble( initial = initial, intype = ensemble, index = i ) board.write( obj = es ) context.output.dedent() context.output.verbose( msg = output.section_close( text = "Sculptor finished" ) ) def process_and_uniqueify(self, template, case, context): if template.ensembles is not None: context.output.verbose( msg = "This entry has been processed" ) return template.ensembles case.problem.data.delete( obj = template ) logger = output.SingleStream( stream = context.output.debug ) from phaser import sculptor from phaser import chisel pdb = template.pdb_object() root = pdb.object sculptor.sanitize_occupancies( root = root ) ( known, unknown ) = sculptor.chain_samples( factory = sculptor.aligned_chain_sample, chain_infos = pdb.get_chain_infos(), alignments = [ template.alignment ], sequence_for = {}, min_length = self.params.min_hssp_length, min_fraction = self.params.min_matching_fraction, logger = logger ) indexer = chisel.ForwardIndexer( root = root ) identities = sculptor.identity_process( indexer = indexer, known = known, logger = logger ) bp_unknown = sculptor.sample_process( indexer = indexer, samples = unknown, processor = self.discarder, logger = logger ) unique_mis = self.unique_modification_instructions( indexer = indexer, known = known, logger = logger ) outensembles = [] import iotbx.pdb for ( tag, bp ) in unique_mis: for d in bp_unknown.data: bp.add( datum = d ) outfile = sculptor.write_pdb_output( root = root, blueprint = bp, identities = identities, file_root = "sculpt", suffix = "%s%s" % ( template.suffix(), tag ), logger = logger ) # Read back this file and check whether valid ### Would be much better checking before writing it out ### outfile_root = iotbx.pdb.input( outfile ).construct_hierarchy() if not outfile_root.models() or not outfile_root.atoms(): context.output.verbose( msg = "%s: model empty" % outfile ) context.output.verbose( msg = "%s: discarded" % outfile ) continue context.output.verbose( msg = "%s: model generated" % outfile ) length = len( list( outfile_root.models()[0].residue_groups() ) ) rmsd = context.identity_to_rms( identity = template.identity(), num_residues = length, ) context.output.verbose( msg = "Identity=%.2f%%, length=%d => rmsd=%.3f" % ( template.identity() * 100, length, rmsd, ) ) ensemble = mr_object.SculptorEnsemble( pdbs = [ outfile ], rmsds = [ rmsd ], composition = template.composition, template = template ) case.problem.data.write( obj = ensemble ) outensembles.append( ensemble ) return outensembles def unique_modification_instructions(self, indexer, known, logger): from phaser import sculptor blueprints = [] uniques = set() for ( name, tag, s ) in self.sculptors: bp = sculptor.sample_process( indexer = indexer, samples = known, processor = s, logger = logger ) mods = frozenset( [ ( i, d ) for ( i, d ) in bp.data if i.is_structural_change() ] ) if mods in uniques: continue uniques.add( mods ) blueprints.append( ( tag, bp ) ) return blueprints def __str__(self): return "ModelGeneration (sculptor)" class DataFetch(object): """ Fetches models """ def __init__(self, site): self.inspect = SelectBest( expert = self, type = stage.TemplateHit ) self.site = site def manipulate(self, candidate, board, case, context): context.output.verbose( msg = output.subtitle( text = "Homology search - fetch" ) ) context.output.notify( message = "status", data = "Data fetch" ) context.output.verbose( msg = "Space group %s" % case.space_group_symbol() ) context.output.verbose( msg = output.blank() ) context.output.verbose( msg = str( candidate ) ) context.output.verbose( msg = "" ) candidate.processed_by.add( self ) template = self.fetch( hss = candidate.hss, case = case, context = context ) candidate.hss.fetched = template if isinstance( template, mr_object.BadHomologySearchEntry ): context.output.verbose( msg = "This entry cannot be used: %s" % template.reason ) context.output.verbose( msg = output.section_close( text = "Fetch finished" ) ) return initial = candidate.initial() existing = [ ts for ts in board.read( type = stage.ModelTemplate ) if ts.initial() == initial ] ts = stage.ModelTemplate( initial = initial, intype = template, index = len( existing ) + 1 ) board.write( obj = ts ) context.output.verbose( msg = output.section_close( text = "Fetch finished" ) ) def fetch(self, hss, case, context): if hss.fetched: context.output.verbose( msg = "This entry has been fetched" ) return hss.fetched case.problem.data.delete( obj = hss ) hit = hss.hit import iotbx.pdb.fetch import iotbx.pdb try: stream = iotbx.pdb.fetch.fetch( id = hit.identifier, mirror = self.site ) root = iotbx.pdb.input( lines = stream.read(), source_info = self.site ) \ .construct_hierarchy() except RuntimeError, e: context.output.verbose( msg = "Error: %s" % e ) return mr_object.BadHomologySearchEntry( hss = hss, reason = e ) if not root.atoms(): context.output.verbose( msg = "Error: PDB file empty" ) return mr_object.BadHomologySearchEntry( hss = hss, reason = "Empty file" ) context.output.verbose( msg = "Selecting reference chain %s" % hit.chain ) asc = root.atom_selection_cache().selection( "chain '%s'" % hit.chain ) croot = root.select( asc, True ) if not croot.atoms(): context.output.verbose( msg = "Error: No atoms after chain selection" ) return mr_object.BadHomologySearchEntry( hss = hss, reason = "No chain '%s'" % hit.chain ) file_name = "%s_%s.pdb" % ( hit.identifier, hit.chain ) context.output.verbose( msg = "Writing out file: %s" % file_name ) croot.write_pdb_file( file_name ) template = mr_object.Template( pdb_file = file_name, alignment = hss.alignment_object(), composition = hss.composition ) context.output.verbose( msg = "Saving template: %s" % template ) case.problem.data.write( obj = template ) return template def __str__(self): return "DataFetch (%s)" % self.site def download_homology_hit_model(hit, pdb, redirections, stream): if redirections.retracted( identifier = hit.identifier ): stream.write( "Entry %s has been retracted\n" % hit.identifier ) raise RuntimeError, "retracted" if redirections.obsoleted( identifier = hit.identifier ): stream.write( "This entry has been obsoleted\n" ) repl = redirections.replacement_for( identifier = hit.identifier ) assert repl stream.write( "Replacement entry: %s\n" % repl ) hit.identifier = repl stream.write( "Retrieving file...\n" ) root = pdb.fetch_entry( identifier = hit.identifier ) if not root.atoms(): stream.write( "Error: PDB file empty\n" ) raise RuntimeError, "empty" stream.write( "Selecting reference chain %s\n" % hit.chain ) asc = root.atom_selection_cache().selection( "chain '%s'" % hit.chain ) croot = root.select( asc, True ) if not croot.atoms(): stream.write( "Error: No atoms after chain selection\n" ) raise RuntimeError, "no chain '%s'" % hit.chain return croot