From: <luc...@us...> - 2007-10-01 23:08:12
|
Revision: 6317 http://cctbx.svn.sourceforge.net/cctbx/?rev=6317&view=rev Author: luc_j_bourhis Date: 2007-10-01 16:08:09 -0700 (Mon, 01 Oct 2007) Log Message: ----------- - move several features out of the charge_flipping_iterator into the miller.fft_map and miller.array - correct the implementation of the variant in ref. [2] - test cases for that variant too (all passed) - a command-line driven program to solve a structure by charge flipping (work in progress) Modified Paths: -------------- trunk/smtbx/ab_initio/charge_flipping.py trunk/smtbx/ab_initio/tests/tst_charge_flipping.py Added Paths: ----------- trunk/smtbx/ab_initio/command_line/ trunk/smtbx/ab_initio/command_line/__init__.py trunk/smtbx/ab_initio/command_line/charge_flipping.py Modified: trunk/smtbx/ab_initio/charge_flipping.py =================================================================== --- trunk/smtbx/ab_initio/charge_flipping.py 2007-10-01 23:04:35 UTC (rev 6316) +++ trunk/smtbx/ab_initio/charge_flipping.py 2007-10-01 23:08:09 UTC (rev 6317) @@ -1,35 +1,98 @@ -""" Charge flipping related tools """ +""" Charge flipping algorithm(s) and related data structures -from __future__ import generators, division +References. +[1] G. Oszl{\'a}nyi and A. S{\"u}t{\H o}. +Ab initio structure solution by charge flipping. +Acta Cryst. A, 60:134--141, 2003. + +[2] G. Oszl{\'a}nyi and A. S{\"u}t{\H o}. +Ab initio structure solution by charge flipping. +II. use of weak reflections. Acta Cryst. A, 61:147, 2004. + +[3] L. Palatinus and G. Chapuis +SUPERFLIP -- a computer program for the solution of crystal structures +by charge flipping in arbitry dimensions +J. Appl. Cryst., 40:786--790, 2007 + +""" + +from __future__ import division + +from libtbx import object_oriented_patterns as oop from cctbx.array_family import flex from cctbx import miller +from cctbx import maptbx import math -class charge_flipping_iterator(object): - """ - References. +class _array_extension(oop.extends, miller.array): - [1] G. Oszl{\'a}nyi and A. S{\"u}t{\H o}. - Ab initio structure solution by charge flipping. - Acta Cryst. A, 60:134--141, 2003. + def oszlanyi_suto_phase_transfer(self, + source, + delta_varphi=math.pi/2, + weak_reflection_fraction=0.2, + need_sorting=True): + """ As per ref. [2] """ + cut = int(weak_reflection_fraction * source.size()) + if need_sorting: + p = self.sort_permutation(by_value="data") + target = self.select(p) + source = source.select(p) + else: + target = self + source_phases = flex.arg(source.data()) + # weak reflections + phases = source_phases[:cut] + delta_varphi + moduli = flex.abs(source.data()[:cut]) + # strong ones + phases.extend(source_phases[cut:]) + moduli.extend(self.data()[cut:]) + return miller.array(self, moduli).phase_transfer(phases) - [2] G. Oszl{\'a}nyi and A. S{\"u}t{\H o}. - Ab initio structure solution by charge flipping. - II. use of weak reflections. Acta Cryst. A, 61:147, 2004. - [3] L. Palatinus and G. Chapuis - SUPERFLIP -- a computer program for the solution of crystal structures - by charge flipping in arbitry dimensions - J. Appl. Cryst., 40:786--790, 2007 +class _fft_extension(oop.extends, miller.fft_map): + """ We add those methods to fft_map so that they can be easily reused and + tested independently of the class charge_flipping_iterator and co. """ + def flipped_fraction_as_delta(self, fraction): + rho = self.real_map_unpadded().as_1d() + p = flex.sort_permutation(rho) + sorted_rho = rho.select(p) + return sorted_rho[int(fraction * sorted_rho.size())] + + def c_flip(self, delta): + rho = self.real_map_unpadded().as_1d() + flipped = rho < delta + return flex.sum(flex.abs(rho.select(flipped))) + + def c_tot(self): + return flex.sum(self.real_map()) + + def c_tot_over_c_flip_as_delta(self, low, mid, high): + pass + + def adjusted_delta_based_on_c_tot_over_c_flip(delta): + r = self.c_tot() / r.c_flip(delta) + if r < 0.8: return delta * 0.9 # magic numbers from SUPERFLIP + elif r > 1: return delta * 1.07 + else: return delta + + +class basic_charge_flipping_iterator(object): + """ An iterator over the sequence of electron densities and structure + factors obtained by repeateadly applying the basic charge flipping + described in ref. [1]. + Notes. self.rho is actually V*rho where V is the unit cell volume """ - def __init__(self, f_obs, delta=None): + def __init__(self, f_obs, delta=None, relative_delta=None): assert f_obs.data() is not None + assert delta is None or relative_delta is None + if relative_delta is not None: + delta = relative_delta * flex.max(f_obs.data()) self.delta = delta self.f_obs = f_obs.expand_to_p1()\ @@ -37,92 +100,72 @@ .merge_equivalents().array()\ .discard_sigmas() - self.set_rho_to_randomized_f_obs() + # Initial rho is the alter ego of f_obs with f_000 = 0 and random phases + f = self.f_obs.randomize_phases() + self.symmetry_flags = maptbx.use_space_group_symmetry + self.rho_map = f.fft_map(symmetry_flags=self.symmetry_flags) # delta being None means we need to find a starting value if self.delta is None: self.initialise_delta() def initialise_delta(self): - self.set_delta_based_on_flipped_fraction() + self.delta = self.rho_map.real_map_unpadded().flipped_fraction_as_delta() - def set_rho_to_randomized_f_obs(self): - # Initial rho is the alter ego of f_obs with f_000 = 0 and random phases - random_phases = (2*math.pi)*flex.random_double(self.f_obs.size()) - f = self.f_obs.phase_transfer(random_phases) - self.rho = f.fft_map().real_map_unpadded() - def __iter__(self): return self def next(self): - # flip - rho_1d_view = self.rho.as_1d() + rho = self.rho_map.real_map_unpadded() + + # flip + rho_1d_view = rho.as_1d() flipped_selection = rho_1d_view < self.delta flipped = rho_1d_view.select(flipped_selection) flipped *= -1 rho_1d_view.set_selected(flipped_selection, flipped) #FFT (it seems the natural place to scale by the number of grid points) - scale = 1/self.rho.size() - self.g = self.f_obs.structure_factors_from_map(self.rho) + scale = 1/rho.size() + self.g = self.f_obs.structure_factors_from_map(rho) self.g *= scale - self.g_000 = flex.sum(self.rho) * scale + self.g_000 = flex.sum(rho) * scale - f = self.transfer_phases_from_g_to_f_obs() + # Transfer the phases of g onto f_obs + f = self.transfer_phase_from_g_to_f_obs() # inverse FFT - self.rho = f.fft_map(f_000=self.g_000).real_map_unpadded() + self.rho_map = f.fft_map(f_000=self.g_000, + symmetry_flags=maptbx.use_space_group_symmetry) # iterator-is-its-own-state trick return self - def transfer_phases_from_g_to_f_obs(self): + def transfer_phase_from_g_to_f_obs(self): return self.f_obs.phase_transfer(self.g) - def set_delta_based_on_flipped_fraction(self): - ## possible optimisation here: partial_sort - rho = self.rho.as_1d() - p = flex.sort_permutation(rho) - sorted_rho = rho.select(p) - self.delta = sorted_rho[int(sorted_rho.size()*0.8)] + def r1_factor(self): + return self.f_obs.r1_factor(self.g, assume_index_matching=True) - def set_delta_based_on_c_tot_over_c_flip(self): - rho = self.rho.as_1d() - c_tot = flex.sum(rho) - c_flip = flex.sum(flex.abs(rho).select(rho < self.delta)) - r = c_tot/c_flip - if r < 0.8: self.delta *= 0.9 # magic numbers from SUPERFLIP - elif r > 1: self.delta *= 1.07 +class weak_reflection_improved_charge_flipping_iterator( + basic_charge_flipping_iterator): + """ The variation described in ref. [2] """ -class charge_flipping_iterator_tweaked_for_weak_reflections( - charge_flipping_iterator): - - def __init__(self, f_obs, - delta_varphi=math.pi/2, - weak_reflection_fraction=0.2 - ): - super(charge_flipping_iterator_tweaked_for_weak_reflections, - self).__init__(f_obs) - - assert ( None not in (delta_varphi, weak_reflection_fraction) - or (delta_varphi, weak_reflection_fraction) == (None, None) ) + def __init__(self, f_obs, delta=None, + delta_varphi=math.pi/2, + weak_reflection_fraction=0.2): + super(weak_reflection_improved_charge_flipping_iterator, + self).__init__(f_obs, delta) + self.delta_varphi = delta_varphi self.weak_reflection_fraction = weak_reflection_fraction - self.delta_varphi = delta_varphi - self.sorting_permutation = flex.sort_permutation(self.f_obs.data()) - self.weak_cut = int(self.f_obs.size()*self.weak_reflection_fraction) + # sort f_obs by increasing amplitudes once and for all + p = self.f_obs.sort_permutation(by_value="data", reverse=True) + self.f_obs = self.f_obs.select(p) - def transfer_phases_from_g_to_f_obs(self): - phases = self.g.phases().data() - p = self.sorting_permutation - i = self.weak_cut - - phi_weaks = phases[:i] + self.delta_varphi - phases.set_selected(p[:i], phi_weaks) - - moduli = flex.double(self.g.size()) - moduli.set_selected(p[:i], self.g.amplitudes().data()[:i]) - moduli.set_selected(p[i:], self.f_obs.data()[i:]) - - return miller.array(self.f_obs, moduli).phase_transfer(phases) + def transfer_phase_from_g_to_f_obs(self): + return self.f_obs.oszlanyi_suto_phase_transfer( + self.g, + self.delta_varphi, + self.weak_reflection_fraction, + need_sorting=False) Added: trunk/smtbx/ab_initio/command_line/__init__.py =================================================================== Added: trunk/smtbx/ab_initio/command_line/charge_flipping.py =================================================================== --- trunk/smtbx/ab_initio/command_line/charge_flipping.py (rev 0) +++ trunk/smtbx/ab_initio/command_line/charge_flipping.py 2007-10-01 23:08:09 UTC (rev 6317) @@ -0,0 +1,73 @@ +from iotbx.option_parser import option_parser +from iotbx import reflection_file_reader +from iotbx.shelx import crystal_symmetry_from_ins +from libtbx.utils import Sorry +from smtbx.ab_initio import charge_flipping + +def run(args, command_name='smtbx.charge_flipping'): + command_line = (option_parser( + usage="%s [options] reflection_file" % command_name, + description="Examples:\n" + " %s data.mtz\n" + " %s hklf4=data.hkl (implicit --symmetry=data.ins)\n" + " %s --symmetry=file.ins data.hkl=hklf4\n" + " %s --unit_cell='1 2 3 90. 105.1 90.' " + "--space_group=P21/n hklf3=data.hkl\n" + % ((command_name,)*4) + ) + .enable_symmetry_comprehensive() + .option(None, "--relative_delta", + action="store", + type="float", + default=None, + help="Parameter delta normalised to the highest reflection amplitude", + metavar="FLOAT") + .option(None, "--max_iterations", + action="store", + type="int", + default=100, + help="Number of iterations to perform", + metavar="FLOAT") + ).process(args) + if not command_line.args: + command_line.parser.show_help() + return 1 + for file_name in command_line.args: + reflection_file = reflection_file_reader.any_reflection_file(file_name) + miller_arrays = None + if reflection_file.file_type() is not None: + try: + miller_arrays = reflection_file.as_miller_arrays( + crystal_symmetry=command_line.symmetry, + force_symmetry=True, + merge_equivalents=False) + except Sorry, KeyboardInterrupt: raise + except: pass + if miller_arrays is None: + print >> sys.stderr, "Error: unknown file format:", file_name + print >> sys.stderr + return 1 + for mi in miller_arrays: + if mi.is_xray_intensity_array(): + print >> sys.stderr, "Warning: converting intensities to amplitudes" + mi = mi.f_sq_as_f() + elif not mi.is_xray_amplitude_array(): + continue + solve(mi, command_line) + +def solve(miller_array, command_line): + flipped = charge_flipping.basic_charge_flipping_iterator( + f_obs=miller_array, + relative_delta=command_line.options.relative_delta) + for i,state in enumerate(flipped): + if i == command_line.options.max_iterations: + print >> sys.stderr, "Maximum number of iterations reached: terminating" + break + r1 = state.r1_factor() + print "R_1 = %.3f" % r1 + + +if __name__ == '__main__': + import sys + import os.path + sys.exit(run(args=sys.argv[1:])) Modified: trunk/smtbx/ab_initio/tests/tst_charge_flipping.py =================================================================== --- trunk/smtbx/ab_initio/tests/tst_charge_flipping.py 2007-10-01 23:04:35 UTC (rev 6316) +++ trunk/smtbx/ab_initio/tests/tst_charge_flipping.py 2007-10-01 23:08:09 UTC (rev 6317) @@ -44,6 +44,7 @@ return f_target, f_000 + def exercise_fixed_delta(space_group_info, elements, anomalous_flag, d_min=0.8, grid_resolution_factor=1./2, @@ -57,7 +58,8 @@ f_target_in_p1 = f_target.expand_to_p1()\ .as_non_anomalous_array()\ .merge_equivalents().array() - rho_target = f_target_in_p1.fft_map(f_000=f_000).real_map_unpadded() + rho_target_map = f_target_in_p1.fft_map(f_000=f_000) + rho_target = rho_target_map.real_map_unpadded() max_rho_target = flex.max(rho_target) delta = 0.01 * max_rho_target # got that by running the debug branch just # below @@ -79,33 +81,41 @@ print "Enter delta:", delta = float(sys.stdin.readline().strip()) * max_rho_target + if verbose: + print "c_tot / c_flip of target = %.3f" % ( + rho_target_map.c_tot() / rho_target_map.c_flip(delta) ) + f_obs = abs(f_target) - flipped = cflip.charge_flipping_iterator(f_obs, delta=delta) - r1s = flex.double() - total_charges = flex.double() - for i,state in enumerate(flipped): - if i == 10: break - isinstance(state, cflip.charge_flipping_iterator) - r1 = f_target_in_p1.r1_factor(state.g, assume_index_matching=True) - r1s.append(r1) - total_charges.append(state.g_000) - pos = (state.rho > 0).count(True) / state.rho.size() - if debug: - print "%i: delta=%.5f" % (i,state.delta) - indent = " "*(len("%i" % i) + 1) - print indent, "R1=%.5f" % r1 - print indent, "fraction of positive pixels=%.3f" % pos - print indent, "total charge=%.3f" % state.g_000 - # r1 decreasing convex function of the iteration number - assert list(flex.sort_permutation(r1s, reverse=True))\ - == list(xrange(r1s.size())) - diffs = r1s[:-1] - r1s[1:] - assert list(flex.sort_permutation(diffs, reverse=True))\ - == list(xrange(diffs.size())) - assert r1s[-1] < 0.1 - assert diffs[-3:].all_lt(0.01) + basic_iter = cflip.basic_charge_flipping_iterator + advanced_iter = cflip.weak_reflection_improved_charge_flipping_iterator + final_r1 = {} + for flipped in (basic_iter(f_obs, delta), + advanced_iter(f_obs, delta), + ): + r1s = flex.double() + for i,state in enumerate(flipped): + if i == 10: break + r1 = f_target_in_p1.r1_factor(state.g) + r1s.append(r1) + if debug: + rho = state.rho_map.real_map() + pos = (rho > 0).count(True) / rho.size() + print "%i: delta=%.5f" % (i,state.delta) + indent = " "*(len("%i" % i) + 1) + print indent, "R1=%.5f" % r1 + print indent, "fraction of positive pixels=%.3f" % pos + print indent, "total charge=%.3f" % state.g_000 + # r1 decreasing function of the iteration number + assert list(flex.sort_permutation(r1s, reverse=True))\ + == list(xrange(r1s.size())) + diffs = r1s[:-1] - r1s[1:] + assert r1s[-1] < 0.1 + assert diffs[-3:].all_lt(0.01) + final_r1[type(flipped)] = r1s[-1] + assert approx_equal(final_r1[basic_iter], final_r1[advanced_iter], 0.01) + def exercise(flags, space_group_info): for i in xrange(10): n_C = random.randint(10,20) @@ -113,13 +123,14 @@ n_N = random.randint(1,5) if flags.Verbose: print "C%i O%i N%i" % (n_C, n_O, n_N) - exercise_fixed_delta(space_group_info=space_group_info, - elements=["C"]*n_C + ["O"]*n_O + ["N"]*n_N, - anomalous_flag=False, - debug=flags.Debug, - verbose=flags.Verbose - ) + exercise_fixed_delta( + space_group_info=space_group_info, + elements=["C"]*n_C + ["O"]*n_O + ["N"]*n_N, + anomalous_flag=False, + debug=flags.Debug, + verbose=flags.Verbose + ) def run(): import sys This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |