#!/usr/bin/env python # LabelitIndex.py # Copyright (C) 2006 CCLRC, Graeme Winter # # This code is distributed under the BSD license, a copy of which is # included in the root directory of this package. # # 2nd June 2006 # # A wrapper for labelit.index - this will provide functionality to: # # Decide the beam centre. # Index the lattce. # import os import sys import copy import shutil import math if not os.environ.has_key('XIA2CORE_ROOT'): raise RuntimeError, 'XIA2CORE_ROOT not defined' if not os.environ.has_key('XIA2_ROOT'): raise RuntimeError, 'XIA2_ROOT not defined' if not os.path.join(os.environ['XIA2CORE_ROOT'], 'Python') in sys.path: sys.path.append(os.path.join(os.environ['XIA2CORE_ROOT'], 'Python')) if not os.environ['XIA2_ROOT'] in sys.path: sys.path.append(os.environ['XIA2_ROOT']) from Driver.DriverFactory import DriverFactory from Handlers.Syminfo import Syminfo # interfaces that this inherits from ... from Schema.Interfaces.FrameProcessor import FrameProcessor from Schema.Interfaces.Indexer import Indexer # other labelit things that this uses from Wrappers.Labelit.LabelitMosflmScript import LabelitMosflmScript from Wrappers.Labelit.LabelitStats_distl import LabelitStats_distl from Modules.Indexer.IndexerSelectImages import index_select_images_lone, \ index_select_images_user from lib.bits import auto_logfiler from Handlers.Streams import Chatter, Debug, Journal from Handlers.Citations import Citations from Handlers.Flags import Flags from Handlers.Files import FileHandler from Modules.Indexer.MosflmCheckIndexerSolution import \ mosflm_check_indexer_solution def LabelitIndex(DriverType = None, indxr_print = True): '''Factory for LabelitIndex wrapper classes, with the specified Driver type.''' DriverInstance = DriverFactory.Driver(DriverType) # FIXME this class should not inherit directly from these interfaces - # would be better to have a higher level module which did that as a # bridge... class LabelitIndexWrapper(DriverInstance.__class__, FrameProcessor, Indexer): '''A wrapper for the program labelit.index - which will provide functionality for deciding the beam centre and indexing the diffraction pattern.''' def __init__(self): DriverInstance.__class__.__init__(self) # interface constructor calls FrameProcessor.__init__(self) Indexer.__init__(self) self.set_executable('labelit.index') # control over the behaviour self._refine_beam = True # this is linked to the above! self._beam_search_scope = 0.0 self._solutions = { } self._solution = None self._indxr_print = indxr_print return def set_refine_beam(self, refine_beam): self._refine_beam = refine_beam return def _write_dataset_preferences(self, n_images): '''Write the dataset_preferences.py file in the working directory - this will include the beam centres etc.''' out = open(os.path.join(self.get_working_directory(), 'dataset_preferences.py'), 'w') # only write things out if they have been overridden # from what is in the header... if self.get_distance_prov() == 'user': out.write('autoindex_override_distance = %f\n' % self.get_distance()) if self.get_wavelength_prov() == 'user': out.write('autoindex_override_wavelength = %f\n' % self.get_wavelength()) if self.get_beam_prov() == 'user': out.write('autoindex_override_beam = (%f, %f)\n' % \ self.get_beam_centre()) if self._refine_beam is False: out.write('beam_search_scope = 0.0\n') else: out.write('beam_search_scope = %f\n' % \ self._beam_search_scope) # check to see if this is an image plate *or* the # wavelength corresponds to Cu KA (1.54A) or Cr KA (2.29 A). # numbers from rigaku americas web page. if math.fabs(self.get_wavelength() - 1.54) < 0.01: out.write('distl_force_binning = True\n') out.write('distl_profile_bumpiness = 10\n') out.write('distl_binned_image_spot_size = 10\n') if math.fabs(self.get_wavelength() - 2.29) < 0.01: out.write('distl_force_binning = True\n') out.write('distl_profile_bumpiness = 10\n') out.write('distl_binned_image_spot_size = 10\n') out.write('wedgelimit = %d\n' % n_images) # new feature - index on the spot centre of mass, not the # highest pixel (should improve the RMS deviation reports.) out.write('distl_spotfinder_algorithm = "maximum_pixel"\n') out.close() return def check_labelit_errors(self): '''Check through the standard output for error reports.''' output = self.get_all_output() for o in output: if 'No_Indexing_Solution' in o: raise RuntimeError, 'indexing failed: %s' % \ o.split(':')[-1].strip() if 'InputFileError' in o: raise RuntimeError, 'indexing failed: %s' % \ o.split(':')[-1].strip() if 'INDEXING UNRELIABLE' in o: raise RuntimeError, 'indexing failed: %s' % \ o.split(':')[-1].strip() return def _index_prepare(self): # prepare to do some autoindexing if self._indxr_images == []: self._index_select_images() return def _index_select_images(self): '''Select correct images based on image headers.''' phi_width = self.get_phi_width() images = self.get_matching_images() if Flags.get_interactive(): selected_images = index_select_images_user(phi_width, images, Chatter) else: selected_images = index_select_images_lone(phi_width, images) for image in selected_images: Debug.write('Selected image %s' % image) self.add_indexer_image_wedge(image) return def _compare_cell(self, c_ref, c_test): '''Compare two sets of unit cell constants: if they differ by less than 5% / 5 degrees return True, else False.''' for j in range(3): if math.fabs((c_test[j] - c_ref[j]) / c_ref[j]) > 0.05: return False for j in range(3, 6): if math.fabs(c_test[j] - c_ref[j]) > 5: return False return True def _index(self): '''Actually index the diffraction pattern. Note well that this is not going to compute the matrix...''' # acknowledge this program Citations.cite('labelit') Citations.cite('distl') self.reset() _images = [] for i in self._indxr_images: for j in i: if not j in _images: _images.append(j) _images.sort() images_str = '%d' % _images[0] for i in _images[1:]: images_str += ', %d' % i cell_str = None if self._indxr_input_cell: cell_str = '%.2f %.2f %.2f %.2f %.2f %.2f' % \ self._indxr_input_cell if self._indxr_sweep_name: # then this is a proper autoindexing run - describe this # to the journal entry if len(self._fp_directory) <= 50: dirname = self._fp_directory else: dirname = '...%s' % self._fp_directory[-46:] Journal.block( 'autoindexing', self._indxr_sweep_name, 'labelit', {'images':images_str, 'target cell':cell_str, 'target lattice':self._indxr_input_lattice, 'template':self._fp_template, 'directory':dirname}) if len(_images) > 4: raise RuntimeError, 'cannot use more than 4 images' auto_logfiler(self) task = 'Autoindex from images:' for i in _images: task += ' %s' % self.get_image_name(i) self.set_task(task) self.add_command_line('--index_only') Debug.write('Indexing from images:') for i in _images: self.add_command_line(self.get_image_name(i)) Debug.write('%s' % self.get_image_name(i)) if self._indxr_input_lattice and False: lattice_to_spacegroup = {'aP':1, 'mP':3, 'mC':5, 'oP':16, 'oC':20, 'oF':22, 'oI':23, 'tP':75, 'tI':79, 'hP':143, 'hR':146, 'cP':195, 'cF':196, 'cI':197} self.add_command_line( 'known_symmetry=%d' % \ lattice_to_spacegroup[self._indxr_input_lattice]) if self._indxr_input_cell and False: self.add_command_line('known_cell=%f,%f,%f,%f,%f,%f' % \ self._indxr_input_cell) self._write_dataset_preferences(len(_images)) shutil.copyfile(os.path.join(self.get_working_directory(), 'dataset_preferences.py'), os.path.join(self.get_working_directory(), '%d_dataset_preferences.py' % \ self.get_xpid())) self.start() self.close_wait() sweep = self.get_indexer_sweep_name() FileHandler.record_log_file( '%s INDEX' % (sweep), self.get_log_file()) # check for errors self.check_for_errors() # check for labelit errors - if something went wrong, then # try to address it by e.g. extending the beam search area... try: self.check_labelit_errors() except RuntimeError, e: if self._refine_beam is False: raise e # can we improve the situation? if self._beam_search_scope < 4.0: self._beam_search_scope += 4.0 # try repeating the indexing! self.set_indexer_done(False) return 'failed' # otherwise this is beyond redemption raise e # ok now we're done, let's look through for some useful stuff output = self.get_all_output() counter = 0 # FIXME 03/NOV/06 something to do with the new centre search... # example output: # Beam center is not immediately clear; rigorously retesting \ # 2 solutions # Beam x 109.0 y 105.1, initial score 538; refined rmsd: 0.1969 # Beam x 108.8 y 106.1, initial score 354; refined rmsd: 0.1792 # in here want to parse the beam centre search if it was done, # and check that the highest scoring solution was declared # the "best" - though should also have a check on the # R.M.S. deviation of that solution... # do this first! for j in range(len(output)): o = output[j] if 'Beam centre is not immediately clear' in o: # read the solutions that it has found and parse the # information centres = [] scores = [] rmsds = [] num_solutions = int(o.split()[-2]) for n in range(num_solutions): record = output[j + n + 1].replace(',', ' ').replace( ';', ' ').split() x, y = float(record[2]), \ float(record[4]) centres.append((x, y)) scores.append(int(record[7])) rmsds.append(float(record[-1])) # next perform some analysis and perhaps assert the # correct solution - for the moment just raise a warning # if it looks like wrong solution may have been picked best_beam_score = (0.0, 0.0, 0) best_beam_rms = (0.0, 0.0, 1.0e8) for n in range(num_solutions): beam = centres[n] score = scores[n] rmsd = rmsds[n] if score > best_beam_score[2]: best_beam_score = (beam[0], beam[1], score) if rmsd < best_beam_rmsd[2]: best_beam_rmsd = (beam[0], beam[1], rmsd) # allow a difference of 0.1mm in either direction... if math.fabs( best_beam_score[0] - best_beam_rmsd[0]) > 0.1 or \ math.fabs(best_beam_score[1] - best_beam_rmsd[1]) > 0.1: Chatter.write( 'Labelit may have picked the wrong beam centre') # FIXME as soon as I get the indexing loop # structure set up, this should reset the # indexing done flag, set the search range to # 0, correct beam and then return... # should also allow for the possibility that # labelit has selected the best solution - so this # will need to remember the stats for this solution, # then compare them against the stats (one day) from # running with the other solution - eventually the # correct solution will result... for o in output: l = o.split() if l[:3] == ['Beam', 'center', 'x']: x = float(l[3].replace('mm,', '')) y = float(l[5].replace('mm,', '')) mosflm_beam_centre = (x, y) mosflm_detector_distance = float(l[7].replace('mm', '')) #self.set_indexer_beam_centre((x, y)) #self.set_indexer_distance(float(l[7].replace('mm', ''))) self._mosaic = float(l[10].replace('mosaicity=', '')) if l[:3] == ['Solution', 'Metric', 'fit']: break counter += 1 # if we've just broken out (counter < len(output)) then # we need to gather the output if counter >= len(output): raise RuntimeError, 'error in indexing' # FIXME this needs to check the smilie status e.g. # ":)" or ";(" or " ". # FIXME need to check the value of the RMSD and raise an # exception if the P1 solution has an RMSD > 1.0... # Change 27/FEB/08 to support user assigned spacegroups # (euugh!) have to "ignore" solutions with higher symmetry # otherwise the rest of xia will override us. Bummer. lattice_to_spacegroup = {'aP':1, 'mP':3, 'mC':5, 'oP':16, 'oC':20, 'oF':22, 'oI':23, 'tP':75, 'tI':79, 'hP':143, 'hR':146, 'cP':195, 'cF':196, 'cI':197} for i in range(counter + 1, len(output)): o = output[i][3:] smiley = output[i][:3] l = o.split() if l: if self._indxr_user_input_lattice: if lattice_to_spacegroup[l[6]] > \ lattice_to_spacegroup[self._indxr_input_lattice]: Debug.write('Ignoring solution: %s' % l[6]) continue self._solutions[int(l[0])] = {'number':int(l[0]), 'mosaic':self._mosaic, 'metric':float(l[1]), 'rmsd':float(l[3]), 'nspots':int(l[4]), 'lattice':l[6], 'cell':map(float, l[7:13]), 'volume':int(l[-1]), 'smiley':smiley} # check the RMSD from the triclinic unit cell if self._solutions[1]['rmsd'] > 1.0 and False: # don't know when this is useful - but I know when it is not! raise RuntimeError, 'high RMSD for triclinic solution' # configure the "right" solution self._solution = self.get_solution() # now store also all of the other solutions... keyed by the # lattice - however these should only be added if they # have a smiley in the appropriate record, perhaps? for solution in self._solutions.keys(): lattice = self._solutions[solution]['lattice'] if self._indxr_other_lattice_cell.has_key(lattice): if self._indxr_other_lattice_cell[lattice]['goodness'] < \ self._solutions[solution]['metric']: continue self._indxr_other_lattice_cell[lattice] = { 'goodness':self._solutions[solution]['metric'], 'cell':self._solutions[solution]['cell']} self._indxr_lattice = self._solution['lattice'] self._indxr_cell = tuple(self._solution['cell']) self._indxr_mosaic = self._solution['mosaic'] lms = LabelitMosflmScript() lms.set_working_directory(self.get_working_directory()) lms.set_solution(self._solution['number']) self._indxr_payload['mosflm_orientation_matrix'] = lms.calculate() # get the beam centre from the mosflm script - mosflm # may have inverted the beam centre and labelit will know # this! mosflm_beam_centre = lms.get_mosflm_beam() if mosflm_beam_centre: self._indxr_payload['mosflm_beam_centre'] = tuple(mosflm_beam_centre) import copy detector = copy.deepcopy(self.get_detector()) beam = copy.deepcopy(self.get_beam_obj()) from Wrappers.Mosflm.AutoindexHelpers import set_mosflm_beam_centre set_mosflm_beam_centre(detector, beam, mosflm_beam_centre) from Experts.SymmetryExpert import lattice_to_spacegroup_number from scitbx import matrix from cctbx import sgtbx, uctbx from dxtbx.model.crystal import crystal_model_from_mosflm_matrix mosflm_matrix = matrix.sqr( [float(i) for line in lms.calculate() for i in line.replace("-", " -").split() ][:9]) space_group = sgtbx.space_group_info(lattice_to_spacegroup_number( self._solution['lattice'])).group() crystal_model = crystal_model_from_mosflm_matrix( mosflm_matrix, unit_cell=uctbx.unit_cell( tuple(self._solution['cell'])), space_group=space_group) from dxtbx.model.experiment.experiment_list import Experiment, ExperimentList experiment = Experiment(beam=beam, detector=detector, goniometer=self.get_goniometer(), scan=self.get_scan(), crystal=crystal_model, ) experiment_list = ExperimentList([experiment]) self.set_indexer_experiment_list(experiment_list) # also get an estimate of the resolution limit from the # labelit.stats_distl output... FIXME the name is wrong! lsd = LabelitStats_distl() lsd.set_working_directory(self.get_working_directory()) lsd.stats_distl() resolution = 1.0e6 for i in _images: stats = lsd.get_statistics(self.get_image_name(i)) resol = 0.5 * (stats['resol_one'] + stats['resol_two']) if resol < resolution: resolution = resol self._indxr_resolution_estimate = resolution return 'ok' def _index_finish(self): '''Check that the autoindexing gave a convincing result, and if not (i.e. it gave a centred lattice where a primitive one would be correct) pick up the correct solution.''' # strictly speaking, given the right input there should be # no need to test... if self._indxr_input_lattice: return if self.get_indexer_sweep(): if self.get_indexer_sweep().get_user_lattice(): return try: raise RuntimeError, 'no, lets not' status, lattice, matrix, cell = mosflm_check_indexer_solution( self) except: status = None if status is None: # basis is primitive return if status is False: # basis is centred, and passes test return # ok need to update internals... self._indxr_lattice = lattice self._indxr_cell = cell Debug.write('Inserting solution: %s ' % lattice + '%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f' % cell) self._indxr_replace(lattice, cell, indxr_print = self._indxr_print) self._indxr_payload['mosflm_orientation_matrix'] = matrix return # things to get results from the indexing def get_solutions(self): return self._solutions def get_solution(self): '''Get the best solution from autoindexing.''' if self._indxr_input_lattice is None: # FIXME in here I need to check that there is a # "good" smiley return copy.deepcopy( self._solutions[max(self._solutions.keys())]) else: # look through for a solution for this lattice - # FIXME should it delete all other solutions? # c/f eliminate. # FIXME should also include a check for the indxr_input_cell if self._indxr_input_cell: for s in self._solutions.keys(): if self._solutions[s]['lattice'] == \ self._indxr_input_lattice: if self._compare_cell( self._indxr_input_cell, self._solutions[s]['cell']): return copy.deepcopy(self._solutions[s]) else: del(self._solutions[s]) else: del(self._solutions[s]) raise RuntimeError, \ 'no solution for lattice %s with given cell' % \ self._indxr_input_lattice else: for s in self._solutions.keys(): if self._solutions[s]['lattice'] == \ self._indxr_input_lattice: return copy.deepcopy(self._solutions[s]) else: del(self._solutions[s]) raise RuntimeError, 'no solution for lattice %s' % \ self._indxr_input_lattice return LabelitIndexWrapper()