[Cctbx-cvs] SF.net SVN: cctbx: [6299] trunk/smtbx/ab_initio

[<luc...@us...>]

Revision: 6299
          http://cctbx.svn.sourceforge.net/cctbx/?rev=6299&view=rev
Author:   luc_j_bourhis
Date:     2007-09-25 14:07:02 -0700 (Tue, 25 Sep 2007)

Log Message:
-----------
- remove a non-ascii character which caused a warning
- adjust delta as done in SUPERFLIP
- still does not seem to cut it though!

Modified Paths:
--------------
    trunk/smtbx/ab_initio/charge_flipping.py
    trunk/smtbx/ab_initio/tests/tst_charge_flipping.py

Modified: trunk/smtbx/ab_initio/charge_flipping.py
===================================================================
--- trunk/smtbx/ab_initio/charge_flipping.py	2007-09-25 14:44:51 UTC (rev 6298)
+++ trunk/smtbx/ab_initio/charge_flipping.py	2007-09-25 21:07:02 UTC (rev 6299)
@@ -12,7 +12,7 @@
 
   [1] G. Oszl{\'a}nyi and A. S{\"u}t{\H o}.
   Ab initio structure solution by charge flipping.
-  Acta Cryst. A, 60:134\xD0141, 2003.
+  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.
@@ -21,7 +21,7 @@
   [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
+  J. Appl. Cryst., 40:786--790, 2007
 
   Notes.
 
@@ -37,14 +37,20 @@
                       .merge_equivalents().array()\
                       .discard_sigmas()
 
+    self.set_rho_to_randomized_f_obs()
+
+    # 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()
+
+  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()
 
-    # delta being None means we need to find a starting value
-    if self.delta is None: self.initialize_delta()
-
   def __iter__(self):
     return self
 
@@ -56,7 +62,7 @@
     flipped *= -1
     rho_1d_view.set_selected(flipped_selection, flipped)
 
-    #FFT (it seems most natural to scale by the number of grid points here)
+    #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)
     self.g *= scale
@@ -73,39 +79,22 @@
   def transfer_phases_from_g_to_f_obs(self):
     return self.f_obs.phase_transfer(self.g)
 
-  def initialize_delta(self):
+  def set_delta_based_on_flipped_fraction(self):
     ## possible optimisation here: partial_sort
     rho = self.rho.as_1d()
-    perm = flex.sort_permutation(rho)
-    sorted_rho = rho.select(perm)
+    p = flex.sort_permutation(rho)
+    sorted_rho = rho.select(p)
     self.delta = sorted_rho[int(sorted_rho.size()*0.8)]
 
-  def adjust_delta(self):
+  def set_delta_based_on_c_tot_over_c_flip(self):
     rho = self.rho.as_1d()
     c_tot = flex.sum(rho)
-    abs_rho = flex.abs(rho)
-    perm = flex.sort_permutation(abs_rho)
-    increasing_abs_rho = abs_rho.select(perm)
-    i,s = flex.find_partial_sum_greater_than(increasing_abs_rho, c_tot)
-    j,t = flex.find_partial_sum_greater_than(increasing_abs_rho,
-                                             c_tot/0.8 - s,
-                                             first_index=i)
-    if j is not None:
-      k = (i+j)//2
-    else:
-      k = -1
-    self.delta = rho[perm[k]]
-    assert 0.8 < self.c_tot_over_c_flip() < 1
-    assert self.delta > 0
+    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
 
-  def c_tot_over_c_flip(self):
-    rho = self.rho.as_1d()
-    c_tot = flex.sum(rho)
-    flipped = rho < self.delta
-    c_flip = flex.sum(flex.abs(rho.select(flipped)))
-    return c_tot / c_flip
 
-
 class charge_flipping_iterator_tweaked_for_weak_reflections(
   charge_flipping_iterator):
 

Modified: trunk/smtbx/ab_initio/tests/tst_charge_flipping.py
===================================================================
--- trunk/smtbx/ab_initio/tests/tst_charge_flipping.py	2007-09-25 14:44:51 UTC (rev 6298)
+++ trunk/smtbx/ab_initio/tests/tst_charge_flipping.py	2007-09-25 21:07:02 UTC (rev 6299)
@@ -28,29 +28,35 @@
     xray_structure=structure,
     algorithm="direct").f_calc()
   f_obs = abs(f_correct)
-  if special_handling_of_weak_reflections:
-    flipped = cflip.charge_flipping_iterator_tweaked_for_weak_reflections(f_obs)
-  else:
-    flipped = cflip.charge_flipping_iterator(f_obs)
-  print "start:"
-  print "fraction of positive pixels=%.3f" % (
-    (flipped.rho > 0).count(True)/flipped.rho.size())
-  for i, state in enumerate(flipped):
-    if i == 10: state.adjust_delta()
-    r1 = state.f_obs.r1_factor(state.g, assume_index_matching=True)
-    pos = (state.rho > 0).count(True) / state.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
+  while True:
+    if special_handling_of_weak_reflections:
+      flipped = cflip.charge_flipping_iterator_tweaked_for_weak_reflections(f_obs)
+    else:
+      flipped = cflip.charge_flipping_iterator(f_obs)
+    print "start:"
+    print "fraction of positive pixels=%.3f" % (
+      (flipped.rho > 0).count(True)/flipped.rho.size())
+    next_delta_check = 10
+    for i,state in enumerate(flipped):
+      isinstance(state, cflip.charge_flipping_iterator)
+      if i == next_delta_check:
+        state.set_delta_based_on_c_tot_over_c_flip()
+        state.set_rho_to_randomized_f_obs()
+        next_delta_check += 10
+      r1 = state.f_obs.r1_factor(state.g, assume_index_matching=True)
+      pos = (state.rho > 0).count(True) / state.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
 
 
 def run():
   exercise_randomly(space_group_info=sgtbx.space_group_info("P31"),
                     elements=["C"]*5 + ["O"]*2 + ["N"],
                     anomalous_flag=False,
-                    special_handling_of_weak_reflections=True
+                    special_handling_of_weak_reflections=False
                   )
 
 if __name__ == '__main__':


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.