cctbx-cvs Mailing List for Computational Crystallography Toolbox (Page 974)

Update of /cvsroot/cctbx/cctbx/cctbx/xray
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv24598/cctbx/xray

Modified Files:
	minimization.py 
Log Message:
cctbx::xray::minimization::apply_shifts(): all if (value < 0) value = 0 style constructs removed; two new cctbx.xray.minimization types: u_penalty_exp, u_penalty_singular_at_zero

Index: minimization.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/xray/minimization.py,v
retrieving revision 1.23
retrieving revision 1.24
diff -C2 -d -r1.23 -r1.24
*** minimization.py	6 Sep 2004 07:40:50 -0000	1.23
--- minimization.py	7 Feb 2005 23:31:05 -0000	1.24
***************
*** 6,14 ****
--- 6,51 ----
  import scitbx.lbfgs
  from scitbx.python_utils.misc import adopt_init_args
+ from libtbx import introspection
  from stdlib import math
  
+ class u_penalty_singular_at_zero:
+ 
+   def __init__(self,
+         penalty_factor=1,
+         penalty_scale=8*math.pi**2,
+         u_min=1.e-6,
+         u_max=1.0):
+     introspection.adopt_init_args()
+ 
+   def functional(self, u):
+     if (u > self.u_max): return 0
+     u = max(u, self.u_min)
+     s = self.penalty_scale
+     return 1/(s*u)*math.exp(-(s*u)*self.penalty_factor)
+ 
+   def gradient(self, u):
+     if (u > self.u_max): return 0
+     u = max(u, self.u_min)
+     s = self.penalty_scale
+     return -1/(math.exp(s*self.penalty_factor*u)*s*u**2) \
+            - self.penalty_factor/(math.exp(s*self.penalty_factor*u)*u)
+ 
+ class u_penalty_exp:
+ 
+   def __init__(self, penalty_factor=1, penalty_scale=10*8*math.pi**2):
+     introspection.adopt_init_args()
+ 
+   def functional(self, u):
+     s = self.penalty_scale
+     return math.exp(-s*u*self.penalty_factor)
+ 
+   def gradient(self, u):
+     s = self.penalty_scale
+     return -s*self.penalty_factor*math.exp(-s*u*self.penalty_factor)
+ 
  class lbfgs:
  
    def __init__(self, target_functor, gradient_flags, xray_structure,
+                      u_penalty=None,
                       lbfgs_termination_params=None,
                       lbfgs_core_params=None,
***************
*** 27,30 ****
--- 64,68 ----
      self.n = xray_structure.n_parameters(gradient_flags)
      self.x = flex.double(self.n, 0)
+     xray_structure.tidy_us(u_min=1.e-6)
      self._scatterers_start = xray_structure.scatterers()
      self._scattering_dict = xray_structure.scattering_dict()
***************
*** 46,55 ****
      scatterers_shifted = ext.minimization_apply_shifts(
        unit_cell=unit_cell,
-       space_group_type=self.xray_structure.space_group_info().type(),
        scatterers=self._scatterers_start,
-       scattering_dict=self._scattering_dict,
        gradient_flags=self.gradient_flags,
!       shifts=self.x,
!       d_min=self._d_min)
      site_symmetry_table = self.xray_structure.site_symmetry_table()
      for i_seq in site_symmetry_table.special_position_indices():
--- 84,90 ----
      scatterers_shifted = ext.minimization_apply_shifts(
        unit_cell=unit_cell,
        scatterers=self._scatterers_start,
        gradient_flags=self.gradient_flags,
!       shifts=self.x)
      site_symmetry_table = self.xray_structure.site_symmetry_table()
      for i_seq in site_symmetry_table.special_position_indices():
***************
*** 78,81 ****
--- 113,121 ----
      if (self.first_target_value is None):
        self.first_target_value = self.f
+     if (self.u_penalty is not None and self.gradient_flags.u_iso):
+       u_isos = self.xray_structure.scatterers().extract_u_iso(
+         unit_cell=self.xray_structure.unit_cell())
+       for u_iso in u_isos:
+         self.f += self.u_penalty.functional(u=u_iso)
      sf = self.structure_factor_gradients(
        xray_structure=self.xray_structure,
***************
*** 86,89 ****
--- 126,139 ----
        algorithm=self.structure_factor_algorithm)
      self.g = sf.packed()
+     if (self.u_penalty is not None and self.gradient_flags.u_iso):
+       g = flex.double()
+       for u_iso in u_isos:
+         g.append(self.u_penalty.gradient(u=u_iso))
+       del u_isos
+       ext.minimization_add_u_iso_gradients(
+         scatterers=self.xray_structure.scatterers(),
+         gradient_flags=self.gradient_flags,
+         xray_gradients=self.g,
+         u_iso_gradients=g)
      if (self.verbose > 1):
        print "xray.minimization line search: f,rms(g):",

Update of /cvsroot/cctbx/cctbx/xray/boost_python
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv24598/xray/boost_python

Modified Files:
	minimization.cpp tst_xray.py 
Log Message:
cctbx::xray::minimization::apply_shifts(): all if (value < 0) value = 0 style constructs removed; two new cctbx.xray.minimization types: u_penalty_exp, u_penalty_singular_at_zero

Index: minimization.cpp
===================================================================
RCS file: /cvsroot/cctbx/cctbx/xray/boost_python/minimization.cpp,v
retrieving revision 1.7
retrieving revision 1.8
diff -C2 -d -r1.7 -r1.8
*** minimization.cpp	25 Jan 2005 23:26:29 -0000	1.7
--- minimization.cpp	7 Feb 2005 23:31:07 -0000	1.8
***************
*** 14,31 ****
        (af::shared<scatterer<> >(*)(
          uctbx::unit_cell const&,
-         sgtbx::space_group_type const&,
          af::const_ref<scatterer<> > const&,
-         scattering_dictionary const&,
          gradient_flags const&,
!         af::const_ref<double> const&,
!         double const&))
            minimization::apply_shifts, (
        arg_("unit_cell"),
-       arg_("space_group_type"),
        arg_("scatterers"),
-       arg_("scattering_dict"),
        arg_("gradient_flags"),
!       arg_("shifts"),
!       arg_("d_min")));
  
      def("minimization_add_site_gradients",
--- 14,25 ----
        (af::shared<scatterer<> >(*)(
          uctbx::unit_cell const&,
          af::const_ref<scatterer<> > const&,
          gradient_flags const&,
!         af::const_ref<double> const&))
            minimization::apply_shifts, (
        arg_("unit_cell"),
        arg_("scatterers"),
        arg_("gradient_flags"),
!       arg_("shifts")));
  
      def("minimization_add_site_gradients",

Index: tst_xray.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/xray/boost_python/tst_xray.py,v
retrieving revision 1.50
retrieving revision 1.51
diff -C2 -d -r1.50 -r1.51
*** tst_xray.py	2 Feb 2005 10:39:14 -0000	1.50
--- tst_xray.py	7 Feb 2005 23:31:08 -0000	1.51
***************
*** 479,483 ****
  def exercise_minimization_apply_shifts():
    uc = uctbx.unit_cell((20, 20, 23))
-   sg = sgtbx.space_group_info("P 4")
    scatterers = flex.xray_scatterer((
      xray.scatterer("Si1", site=(0.01,0.02,0.3), fp=-1, fdp=2),
--- 479,482 ----
***************
*** 485,494 ****
                     u=adptbx.u_cart_as_u_star(uc,
                       (0.04,0.05,0.06,-.005,0.02,-0.002)))))
-   scattering_dict = xray.ext.scattering_dictionary(scatterers)
-   scattering_dict.assign_from_table("WK1995")
    f = xray.ext.gradient_flags(True, True, True, True, True, True)
    shifts = flex.double(19, 0.001)
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, sg.type(), scatterers, scattering_dict, f, shifts, 0)
    for a,b in zip(scatterers, shifted_scatterers):
      assert a.scattering_type == b.scattering_type
--- 484,491 ----
                     u=adptbx.u_cart_as_u_star(uc,
                       (0.04,0.05,0.06,-.005,0.02,-0.002)))))
    f = xray.ext.gradient_flags(True, True, True, True, True, True)
    shifts = flex.double(19, 0.001)
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, scatterers, f, shifts)
    for a,b in zip(scatterers, shifted_scatterers):
      assert a.scattering_type == b.scattering_type
***************
*** 506,515 ****
    shifted_scatterers = xray.ext.minimization_apply_shifts(
      unit_cell=uc,
-     space_group_type=sg.type(),
      scatterers=shifted_scatterers,
-     scattering_dict=scattering_dict,
      gradient_flags=f,
!     shifts=shifts,
!     d_min=0)
    for a,b in zip(scatterers, shifted_scatterers):
      assert a.scattering_type == b.scattering_type
--- 503,509 ----
    shifted_scatterers = xray.ext.minimization_apply_shifts(
      unit_cell=uc,
      scatterers=shifted_scatterers,
      gradient_flags=f,
!     shifts=shifts)
    for a,b in zip(scatterers, shifted_scatterers):
      assert a.scattering_type == b.scattering_type
***************
*** 523,527 ****
    shifts = flex.double((-1,2,-3,4,-5,-6))
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, sg.type(), scatterers, scattering_dict, f, shifts, 0)
    assert approx_equal(
      shifted_scatterers[0].site,
--- 517,521 ----
    shifts = flex.double((-1,2,-3,4,-5,-6))
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, scatterers, f, shifts)
    assert approx_equal(
      shifted_scatterers[0].site,
***************
*** 533,575 ****
    shifts = flex.double(1, -10)
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, sg.type(), scatterers, scattering_dict, f, shifts, 0)
!   assert shifted_scatterers[0].u_iso == 0
    f = xray.ext.gradient_flags(False, False, True, False, False, False)
    shifts = flex.double(6, -100)
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, sg.type(), scatterers, scattering_dict, f, shifts, 0)
!   assert not approx_equal(scatterers[1].u_star, shifted_scatterers[1].u_star)
!   shifts = flex.double(6, -1000)
!   shifted_scatterers_2 = xray.ext.minimization_apply_shifts(
!     uc, sg.type(), scatterers, scattering_dict, f, shifts, 0)
!   # due to eigenvalue filtering two extreme shifts
!   # should produce the same result
!   assert approx_equal(shifted_scatterers[1].u_star,
!                       shifted_scatterers_2[1].u_star)
    f = xray.ext.gradient_flags(False, False, False, True, False, False)
    shifts = flex.double(2, -10)
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, sg.type(), scatterers, scattering_dict, f, shifts, 0)
    for i in xrange(2):
!     assert shifted_scatterers[i].occupancy == 0
    f = xray.ext.gradient_flags(False, False, False, False, True, False)
    shifts = flex.double(2, -10)
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, sg.type(), scatterers, scattering_dict, f, shifts, 0)
    assert approx_equal(shifted_scatterers[0].fp, -11)
    assert shifted_scatterers[1].fp == -10
-   shifted_scatterers = xray.ext.minimization_apply_shifts(
-     uc, sg.type(), scatterers, scattering_dict, f, shifts, 3)
    for i in xrange(2):
-     assert approx_equal(
-       shifted_scatterers[i].fp,
-       -(scattering_dict
-           .lookup(shifted_scatterers[i].scattering_type)
-           .gaussian.at_d_star_sq(1/9.)))
      assert shifted_scatterers[i].fdp == scatterers[i].fdp
    f = xray.ext.gradient_flags(False, False, False, False, False, True)
    shifts = flex.double((2,3))
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, sg.type(), scatterers, scattering_dict, f, shifts, 3)
    assert shifted_scatterers[0].fp == -1
    assert approx_equal(shifted_scatterers[0].fdp, 4)
--- 527,556 ----
    shifts = flex.double(1, -10)
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, scatterers, f, shifts)
!   assert approx_equal(shifted_scatterers[0].u_iso, -10)
    f = xray.ext.gradient_flags(False, False, True, False, False, False)
    shifts = flex.double(6, -100)
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, scatterers, f, shifts)
!   assert not approx_equal(shifted_scatterers[1].u_star,
!     [u_ij-100 for u_ij in scatterers[1].u_star])
    f = xray.ext.gradient_flags(False, False, False, True, False, False)
    shifts = flex.double(2, -10)
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, scatterers, f, shifts)
    for i in xrange(2):
!     assert approx_equal(shifted_scatterers[i].occupancy, -9)
    f = xray.ext.gradient_flags(False, False, False, False, True, False)
    shifts = flex.double(2, -10)
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, scatterers, f, shifts)
    assert approx_equal(shifted_scatterers[0].fp, -11)
    assert shifted_scatterers[1].fp == -10
    for i in xrange(2):
      assert shifted_scatterers[i].fdp == scatterers[i].fdp
    f = xray.ext.gradient_flags(False, False, False, False, False, True)
    shifts = flex.double((2,3))
    shifted_scatterers = xray.ext.minimization_apply_shifts(
!     uc, scatterers, f, shifts)
    assert shifted_scatterers[0].fp == -1
    assert approx_equal(shifted_scatterers[0].fdp, 4)
***************
*** 578,583 ****
    shifts = flex.double(1)
    try:
!     xray.ext.minimization_apply_shifts(
!       uc, sg.type(), scatterers, scattering_dict, f, shifts, 0)
    except Exception, e:
      assert str(e) == "scitbx Error: Array of shifts is too small."
--- 559,563 ----
    shifts = flex.double(1)
    try:
!     xray.ext.minimization_apply_shifts(uc, scatterers, f, shifts)
    except Exception, e:
      assert str(e) == "scitbx Error: Array of shifts is too small."
***************
*** 586,591 ****
    shifts = flex.double(3)
    try:
!     xray.ext.minimization_apply_shifts(
!       uc, sg.type(), scatterers, scattering_dict, f, shifts, 0)
    except Exception, e:
      assert str(e) == "cctbx Error: Array of shifts is too large."
--- 566,570 ----
    shifts = flex.double(3)
    try:
!     xray.ext.minimization_apply_shifts(uc, scatterers, f, shifts)
    except Exception, e:
      assert str(e) == "cctbx Error: Array of shifts is too large."

Update of /cvsroot/cctbx/cctbx/include/cctbx/xray
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv24598/include/cctbx/xray

Modified Files:
	minimization.h 
Log Message:
cctbx::xray::minimization::apply_shifts(): all if (value < 0) value = 0 style constructs removed; two new cctbx.xray.minimization types: u_penalty_exp, u_penalty_singular_at_zero

Index: minimization.h
===================================================================
RCS file: /cvsroot/cctbx/cctbx/include/cctbx/xray/minimization.h,v
retrieving revision 1.10
retrieving revision 1.11
diff -C2 -d -r1.10 -r1.11
*** minimization.h	25 Jan 2005 23:25:27 -0000	1.10
--- minimization.h	7 Feb 2005 23:31:07 -0000	1.11
***************
*** 5,9 ****
  #include <cctbx/xray/gradient_flags.h>
  #include <cctbx/xray/packing_order.h>
- #include <cctbx/sgtbx/space_group_type.h>
  #include <scitbx/array_family/block_iterator.h>
  
--- 5,8 ----
***************
*** 15,31 ****
    apply_shifts(
      uctbx::unit_cell const& unit_cell,
-     sgtbx::space_group_type const& space_group_type,
      af::const_ref<XrayScattererType> const& scatterers,
-     scattering_dictionary const& scattering_dict,
      xray::gradient_flags const& gradient_flags,
!     af::const_ref<FloatType> const& shifts,
!     FloatType const& d_min)
    {
-     CCTBX_ASSERT(scattering_dict.n_scatterers() == scatterers.size());
      BOOST_STATIC_ASSERT(packing_order_convention == 1);
      typedef typename XrayScattererType::float_type sc_f_t;
      af::shared<XrayScattererType> result((af::reserve(scatterers.size())));
-     sc_f_t d_star_sq_max = 0;
-     if (d_min > 0) d_star_sq_max = 1 / (d_min * d_min);
      scitbx::af::const_block_iterator<FloatType> next_shifts(
        shifts, "Array of shifts is too small.");
--- 14,24 ----
    apply_shifts(
      uctbx::unit_cell const& unit_cell,
      af::const_ref<XrayScattererType> const& scatterers,
      xray::gradient_flags const& gradient_flags,
!     af::const_ref<FloatType> const& shifts)
    {
      BOOST_STATIC_ASSERT(packing_order_convention == 1);
      typedef typename XrayScattererType::float_type sc_f_t;
      af::shared<XrayScattererType> result((af::reserve(scatterers.size())));
      scitbx::af::const_block_iterator<FloatType> next_shifts(
        shifts, "Array of shifts is too small.");
***************
*** 38,42 ****
          if (gradient_flags.u_iso) {
            sc.u_iso += next_shifts();
-           if (sc.u_iso < 0) sc.u_iso = 0;
          }
        }
--- 31,34 ----
***************
*** 46,67 ****
              unit_cell, sc.u_star);
            u_cart += scitbx::sym_mat3<sc_f_t>(next_shifts(6));
!           sc.u_star = adptbx::u_cart_as_u_star(
!             unit_cell, adptbx::eigenvalue_filtering(u_cart));
          }
        }
        if (gradient_flags.occupancy) {
          sc.occupancy += next_shifts();
-         if (sc.occupancy < 0) sc.occupancy = 0;
        }
        if (gradient_flags.fp) {
          sc.fp += next_shifts();
-         if (d_star_sq_max != 0) {
-           sc_f_t f0 = scattering_dict
-             .lookup(sc.scattering_type)
-             .gaussian.at_d_star_sq(d_star_sq_max);
-           if (f0 + sc.fp < 0) {
-             sc.fp = -f0;
-           }
-         }
        }
        if (gradient_flags.fdp) {
--- 38,49 ----
              unit_cell, sc.u_star);
            u_cart += scitbx::sym_mat3<sc_f_t>(next_shifts(6));
!           sc.u_star = adptbx::u_cart_as_u_star(unit_cell, u_cart);
          }
        }
        if (gradient_flags.occupancy) {
          sc.occupancy += next_shifts();
        }
        if (gradient_flags.fp) {
          sc.fp += next_shifts();
        }
        if (gradient_flags.fdp) {

Update of /cvsroot/cctbx/cctbx/cctbx/development
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv23468/cctbx/development

Modified Files:
	random_structure.py 
Log Message:
new random_u_iso_min parameter

Index: random_structure.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/development/random_structure.py,v
retrieving revision 1.20
retrieving revision 1.21
diff -C2 -d -r1.20 -r1.21
*** random_structure.py	2 Feb 2005 09:48:46 -0000	1.20
--- random_structure.py	7 Feb 2005 23:27:58 -0000	1.21
***************
*** 129,132 ****
--- 129,133 ----
                 random_f_double_prime_scale=0.6,
                 random_u_iso=False,
+                random_u_iso_min=0,
                 random_u_iso_scale=0.3,
                 u_iso=0,
***************
*** 200,204 ****
          u_iso = self.u_iso
          if (not u_iso and self.random_u_iso):
!           u_iso = random.random() * self.random_u_iso_scale
          scatterer.u_iso = u_iso
        else:
--- 201,206 ----
          u_iso = self.u_iso
          if (not u_iso and self.random_u_iso):
!           u_iso = random.random() * self.random_u_iso_scale \
!                 + self.random_u_iso_min
          scatterer.u_iso = u_iso
        else:

Update of /cvsroot/cctbx/cctbx/cctbx/regression
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv12619/cctbx/regression

Modified Files:
	tst_sampled_model_density.py 
Log Message:
FFT structure factors and gradients: tolerates negative occupancies, large negative fp, non-positive-definite displacement parameters (same results as direct summation)

Index: tst_sampled_model_density.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/regression/tst_sampled_model_density.py,v
retrieving revision 1.29
retrieving revision 1.30
diff -C2 -d -r1.29 -r1.30
*** tst_sampled_model_density.py	6 Sep 2004 07:40:49 -0000	1.29
--- tst_sampled_model_density.py	7 Feb 2005 18:12:13 -0000	1.30
***************
*** 1,4 ****
--- 1,6 ----
  from cctbx import xray
  from cctbx import maptbx
+ from cctbx import crystal
+ from cctbx import adptbx
  import cctbx.eltbx.xray_scattering
  from cctbx import eltbx
***************
*** 8,11 ****
--- 10,14 ----
  from cctbx.array_family import flex
  from scitbx import fftpack
+ from libtbx.test_utils import approx_equal
  import random
  import sys
***************
*** 173,176 ****
--- 176,249 ----
      verbose=verbose)
  
+ def exercise_negative_parameters(verbose=0):
+   structure_default = xray.structure(
+     crystal_symmetry = crystal.symmetry(
+       unit_cell=((10,13,17,75,80,85)),
+       space_group_symbol="P 1"),
+     scatterers=flex.xray_scatterer([
+       xray.scatterer(label="C", site=(0,0,0), u=0.25)]))
+   negative_gaussian = eltbx.xray_scattering.gaussian((1,2), (2,3), -4)
+   for i_trial in xrange(7):
+     structure = structure_default.deep_copy_scatterers()
+     scatterer = structure.scatterers()[0]
+     if (i_trial == 1):
+       scatterer.occupancy *= -1
+     elif (i_trial == 2):
+       structure.scattering_dict(custom_dict={"C": negative_gaussian})
+     elif (i_trial == 3):
+       scatterer.u_iso *= -1
+     elif (i_trial == 4):
+       u_cart = adptbx.random_u_cart(u_scale=1, u_min=-1.1)
+       assert max(adptbx.eigenvalues(u_cart)) < 0
+       u_star = adptbx.u_cart_as_u_star(structure.unit_cell(), u_cart)
+       scatterer.u_star = u_star
+       scatterer.anisotropic_flag = True
+     elif (i_trial == 5):
+       scatterer.fp = -10
+     elif (i_trial == 6):
+       scatterer.fp = -3
+     f_direct = structure.structure_factors(
+       d_min=1, algorithm="direct", cos_sin_table=False).f_calc()
+     f_fft = structure.structure_factors(
+       d_min=1, algorithm="fft",
+       quality_factor=1.e8, wing_cutoff=1.e-10).f_calc()
+     if (i_trial == 2):
+       assert negative_gaussian.at_d_star_sq(f_fft.d_star_sq().data()).all_lt(0)
+     if (i_trial in [5,6]):
+       f = structure.scattering_dict().lookup("C").gaussian.at_d_star_sq(
+         f_fft.d_star_sq().data())
+       if (i_trial == 5):
+         assert flex.max(f) + scatterer.fp < 0
+       else:
+         assert flex.max(f) + scatterer.fp > 0
+         assert flex.min(f) + scatterer.fp < 0
+     cc = flex.linear_correlation(
+       abs(f_direct).data(),
+       abs(f_fft).data()).coefficient()
+     if (cc < 0.999):
+       raise AssertionError("i_trial=%d, correlation=%.6g" % (i_trial, cc))
+     elif (0 or verbose):
+       print "correlation=%.6g" % cc
+     #
+     # very simple test of gradient calculations with negative parameters
+     structure_factor_gradients = \
+       cctbx.xray.structure_factors.gradients(
+         miller_set=f_direct,
+         cos_sin_table=False)
+     target_functor = xray.target_functors.intensity_correlation(
+       f_obs=abs(f_direct))
+     target_result = target_functor(f_fft, True)
+     gradient_flags = xray.structure_factors.gradient_flags(default=True)
+     grads = []
+     for algorithm in ["direct", "fft"]:
+       grads.append(structure_factor_gradients(
+         xray_structure=structure,
+         miller_set=f_direct,
+         d_target_d_f_calc=target_result.derivatives(),
+         gradient_flags=gradient_flags,
+         n_parameters=structure.n_parameters(gradient_flags=gradient_flags),
+         algorithm="direct").packed())
+     assert approx_equal(grads[0], grads[1]) # all very close to zero
+ 
  def run_call_back(flags, space_group_info):
    for anomalous_flag in (False, True)[:]: #SWITCH
***************
*** 187,190 ****
--- 260,264 ----
  def run():
    debug_utils.parse_options_loop_space_groups(sys.argv[1:], run_call_back)
+   exercise_negative_parameters()
  
  if (__name__ == "__main__"):

Update of /cvsroot/cctbx/cctbx/include/cctbx/xray
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv12619/include/cctbx/xray

Modified Files:
	fast_gradients.h sampled_model_density.h sampling_base.h 
Log Message:
FFT structure factors and gradients: tolerates negative occupancies, large negative fp, non-positive-definite displacement parameters (same results as direct summation)

Index: sampling_base.h
===================================================================
RCS file: /cvsroot/cctbx/cctbx/include/cctbx/xray/sampling_base.h,v
retrieving revision 1.14
retrieving revision 1.15
diff -C2 -d -r1.14 -r1.15
*** sampling_base.h	21 Oct 2004 17:37:24 -0000	1.14
--- sampling_base.h	7 Feb 2005 18:12:18 -0000	1.15
***************
*** 160,164 ****
      {
        FloatType d = b_all.determinant();
!       CCTBX_ASSERT(d != 0);
        scitbx::sym_mat3<FloatType> cfmt = b_all.co_factor_matrix_transposed();
        as = eight_pi_pow_3_2 * a / std::sqrt(d);
--- 160,164 ----
      {
        FloatType d = b_all.determinant();
!       CCTBX_ASSERT(d > 0);
        scitbx::sym_mat3<FloatType> cfmt = b_all.co_factor_matrix_transposed();
        as = eight_pi_pow_3_2 * a / std::sqrt(d);
***************
*** 174,178 ****
      {
        FloatType d = b_all.determinant();
!       CCTBX_ASSERT(d != 0);
        return eight_pi_pow_3_2 * a / std::sqrt(d);
      }
--- 174,178 ----
      {
        FloatType d = b_all.determinant();
!       CCTBX_ASSERT(d > 0);
        return eight_pi_pow_3_2 * a / std::sqrt(d);
      }
***************
*** 679,682 ****
--- 679,683 ----
    {
      CCTBX_ASSERT(scattering_dict.n_scatterers() == scatterers.size());
+     CCTBX_ASSERT(u_base > 0);
      scitbx::mat3<FloatType> orth_mx = unit_cell_.orthogonalization_matrix();
      if (orth_mx[3] != 0 || orth_mx[6] != 0 || orth_mx[7] != 0) {
***************
*** 693,696 ****
--- 694,698 ----
      FloatType sum_w_gaussian_a = 0;
      FloatType sum_w_u_equiv = 0;
+     bool have_u_min = false;
      for(std::size_t i_seq=0;i_seq<scatterers.size();i_seq++) {
        XrayScattererType const& scatterer = scatterers[i_seq];
***************
*** 708,715 ****
        FloatType u_min;
        if (!scatterer.anisotropic_flag) {
-         CCTBX_ASSERT(scatterer.u_iso >= 0);
          u_radius_cache_.push_back(scatterer.u_iso);
          u_min = scatterer.u_iso;
!         sum_w_u_equiv += w * scatterer.u_iso;
        }
        else {
--- 710,717 ----
        FloatType u_min;
        if (!scatterer.anisotropic_flag) {
          u_radius_cache_.push_back(scatterer.u_iso);
          u_min = scatterer.u_iso;
!         sum_w_u_equiv += w * std::max(static_cast<FloatType>(0),
!           scatterer.u_iso);
        }
        else {
***************
*** 718,730 ****
          scitbx::vec3<FloatType>
            u_cart_eigenvalues = adptbx::eigenvalues(u_cart_cache_.back());
-         CCTBX_ASSERT(adptbx::is_positive_definite(
-           u_cart_eigenvalues, tolerance_positive_definite));
          u_radius_cache_.push_back(af::max(u_cart_eigenvalues));
          u_min = af::min(u_cart_eigenvalues);
!         sum_w_u_equiv += w * adptbx::u_cart_as_u_iso(u_cart_cache_.back());
        }
!       CCTBX_ASSERT(u_radius_cache_.back() >= 0);
!       if (u_min_ < 0) {
          u_min_ = u_min;
        }
        else {
--- 720,731 ----
          scitbx::vec3<FloatType>
            u_cart_eigenvalues = adptbx::eigenvalues(u_cart_cache_.back());
          u_radius_cache_.push_back(af::max(u_cart_eigenvalues));
          u_min = af::min(u_cart_eigenvalues);
!         sum_w_u_equiv += w * std::max(static_cast<FloatType>(0),
!           adptbx::u_cart_as_u_iso(u_cart_cache_.back()));
        }
!       if (not have_u_min) {
          u_min_ = u_min;
+         have_u_min = true;
        }
        else {
***************
*** 732,736 ****
        }
      }
!     u_extra_ = u_base_ - u_min_;
      if (sum_w != 0) {
        ave_u_iso_or_equiv_ = sum_w_u_equiv / sum_w;
--- 733,739 ----
        }
      }
!     if (have_u_min) {
!       u_extra_ = u_base_ - u_min_;
!     }
      if (sum_w != 0) {
        ave_u_iso_or_equiv_ = sum_w_u_equiv / sum_w;

Index: sampled_model_density.h
===================================================================
RCS file: /cvsroot/cctbx/cctbx/include/cctbx/xray/sampled_model_density.h,v
retrieving revision 1.24
retrieving revision 1.25
diff -C2 -d -r1.24 -r1.25
*** sampled_model_density.h	16 Apr 2004 14:06:12 -0000	1.24
--- sampled_model_density.h	7 Feb 2005 18:12:18 -0000	1.25
***************
*** 115,119 ****
        eltbx::xray_scattering::gaussian const& gaussian=scattering_dict.lookup(
          scatterer.scattering_type).gaussian;
-       CCTBX_ASSERT(scatterer.weight() >= 0);
        if (scatterer.weight() == 0) continue;
        FloatType fdp = scatterer.fdp;
--- 115,118 ----

Index: fast_gradients.h
===================================================================
RCS file: /cvsroot/cctbx/cctbx/include/cctbx/xray/fast_gradients.h,v
retrieving revision 1.34
retrieving revision 1.35
diff -C2 -d -r1.34 -r1.35
*** fast_gradients.h	29 Nov 2004 18:50:38 -0000	1.34
--- fast_gradients.h	7 Feb 2005 18:12:17 -0000	1.35
***************
*** 482,486 ****
        FloatType gr_fp(0);
        FloatType gr_fdp(0);
-       CCTBX_ASSERT(scatterer.weight() >= 0);
        if (scatterer.weight() != 0) {
          FloatType fdp = scatterer.fdp;
--- 482,485 ----

Update of /cvsroot/cctbx/cctbx/cctbx/xray
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv12619/cctbx/xray

Modified Files:
	structure.py 
Log Message:
FFT structure factors and gradients: tolerates negative occupancies, large negative fp, non-positive-definite displacement parameters (same results as direct summation)

Index: structure.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/xray/structure.py,v
retrieving revision 1.49
retrieving revision 1.50
diff -C2 -d -r1.49 -r1.50
*** structure.py	4 Feb 2005 02:57:56 -0000	1.49
--- structure.py	7 Feb 2005 18:12:14 -0000	1.50
***************
*** 232,236 ****
                                quality_factor=None,
                                u_base=None,
!                               b_base=None):
      if (anomalous_flag is None):
        if (self.scatterers().count_anomalous() != 0):
--- 232,237 ----
                                quality_factor=None,
                                u_base=None,
!                               b_base=None,
!                               wing_cutoff=None):
      if (anomalous_flag is None):
        if (self.scatterers().count_anomalous() != 0):
***************
*** 250,254 ****
        quality_factor=quality_factor,
        u_base=u_base,
!       b_base=b_base)(
          xray_structure=self,
          miller_set=miller_set,
--- 251,256 ----
        quality_factor=quality_factor,
        u_base=u_base,
!       b_base=b_base,
!       wing_cutoff=wing_cutoff)(
          xray_structure=self,
          miller_set=miller_set,

Update of /cvsroot/cctbx/scitbx/include/scitbx/math
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv6841

Modified Files:
	principal_axes_of_inertia.h 
Log Message:
missing include added

Index: principal_axes_of_inertia.h
===================================================================
RCS file: /cvsroot/cctbx/scitbx/include/scitbx/math/principal_axes_of_inertia.h,v
retrieving revision 1.4
retrieving revision 1.5
diff -C2 -d -r1.4 -r1.5
*** principal_axes_of_inertia.h	4 Feb 2005 06:19:25 -0000	1.4
--- principal_axes_of_inertia.h	4 Feb 2005 17:21:44 -0000	1.5
***************
*** 3,6 ****
--- 3,7 ----
  
  #include <scitbx/math/eigensystem.h>
+ #include <cstdio>
  
  namespace scitbx { namespace math {

Update of /cvsroot/cctbx/iotbx/iotbx
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv12449

Modified Files:
	tst_reflection_file_utils.py 
Log Message:
work around strange Python 2.2.1 C assert failure

Index: tst_reflection_file_utils.py
===================================================================
RCS file: /cvsroot/cctbx/iotbx/iotbx/tst_reflection_file_utils.py,v
retrieving revision 1.5
retrieving revision 1.6
diff -C2 -d -r1.5 -r1.6
*** tst_reflection_file_utils.py	26 Jan 2005 07:47:00 -0000	1.5
--- tst_reflection_file_utils.py	4 Feb 2005 06:34:48 -0000	1.6
***************
*** 1,4 ****
- from iotbx.reflection_file_utils import reflection_file_server
  from iotbx import reflection_file_reader
  from iotbx import mtz
  from cctbx import miller
--- 1,4 ----
  from iotbx import reflection_file_reader
+ from iotbx.reflection_file_utils import reflection_file_server
  from iotbx import mtz
  from cctbx import miller

Update of /cvsroot/cctbx/cctbx/cctbx/regression
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv12413

Modified Files:
	tst_k_b_scaling.py tst_reflection_statistics.py 
Log Message:
work around strange Python 2.2.1 C assert failure

Index: tst_k_b_scaling.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/regression/tst_k_b_scaling.py,v
retrieving revision 1.14
retrieving revision 1.15
diff -C2 -d -r1.14 -r1.15
*** tst_k_b_scaling.py	1 Feb 2005 07:22:20 -0000	1.14
--- tst_k_b_scaling.py	4 Feb 2005 06:34:25 -0000	1.15
***************
*** 1,9 ****
- from cctbx.array_family import flex
- from cctbx import mintbx
  from cctbx import miller
  from cctbx import adptbx
! from scitbx import lbfgs
  from cctbx.development import debug_utils
  from cctbx.development import random_structure
  from scitbx.python_utils.misc import adopt_init_args
  import sys
--- 1,9 ----
  from cctbx import miller
  from cctbx import adptbx
! from cctbx import mintbx
! from cctbx.array_family import flex
  from cctbx.development import debug_utils
  from cctbx.development import random_structure
+ from scitbx import lbfgs
  from scitbx.python_utils.misc import adopt_init_args
  import sys

Index: tst_reflection_statistics.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/regression/tst_reflection_statistics.py,v
retrieving revision 1.4
retrieving revision 1.5
diff -C2 -d -r1.4 -r1.5
*** tst_reflection_statistics.py	28 Jan 2005 17:06:12 -0000	1.4
--- tst_reflection_statistics.py	4 Feb 2005 06:34:25 -0000	1.5
***************
*** 1,2 ****
--- 1,3 ----
+ from cctbx.development import debug_utils
  from cctbx import miller
  from cctbx import crystal
***************
*** 6,10 ****
  import cctbx.sgtbx.subgroups
  import cctbx.sgtbx.cosets
- from cctbx.development import debug_utils
  from cctbx.array_family import flex
  from libtbx.test_utils import approx_equal
--- 7,10 ----

Update of /cvsroot/cctbx/scitbx/include/scitbx/math
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv10604

Modified Files:
	principal_axes_of_inertia.h 
Log Message:
RuntimeError: weight=-0.0603412 is negative (must be >=0)

Index: principal_axes_of_inertia.h
===================================================================
RCS file: /cvsroot/cctbx/scitbx/include/scitbx/math/principal_axes_of_inertia.h,v
retrieving revision 1.3
retrieving revision 1.4
diff -C2 -d -r1.3 -r1.4
*** principal_axes_of_inertia.h	3 Feb 2005 12:51:56 -0000	1.3
--- principal_axes_of_inertia.h	4 Feb 2005 06:19:25 -0000	1.4
***************
*** 6,9 ****
--- 6,26 ----
  namespace scitbx { namespace math {
  
+   namespace detail {
+ 
+     std::string
+     report_negative_weight(
+       double weight,
+       const char* where_file_name,
+       long where_line_number)
+     {
+       char buf[256];
+       std::sprintf(buf,
+         "weight=%.6g is negative (must be >=0) (%s, line %ld)",
+         weight, where_file_name, where_line_number);
+       return std::string(buf);
+     }
+ 
+   } // namespace detail
+ 
    /*! Principal axes of inertia given discrete points in 3-dimensional
        space and optionally weights.
***************
*** 60,63 ****
--- 77,84 ----
          for(std::size_t i_p=0;i_p<points.size();i_p++) {
            FloatType w = weights[i_p];
+           if (w < 0) {
+             throw std::runtime_error(detail::report_negative_weight(
+               static_cast<double>(w), __FILE__, __LINE__));
+           }
            center_of_mass_ += w * points[i_p];
            sum_weights += w;

Update of /cvsroot/cctbx/cctbx/cctbx/development
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv5055

Modified Files:
	electron_density_sampling.py 
Log Message:
avoid division by zero time

Index: electron_density_sampling.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/development/electron_density_sampling.py,v
retrieving revision 1.3
retrieving revision 1.4
diff -C2 -d -r1.3 -r1.4
*** electron_density_sampling.py	6 Sep 2004 07:40:48 -0000	1.3
--- electron_density_sampling.py	4 Feb 2005 05:39:46 -0000	1.4
***************
*** 31,36 ****
          times.append(f_calc_object.manager().estimate_time_fft.time_sampling)
          print "  %.2f seconds," % times[-1]
!       print "d_min=%d: %.2f s / %.2f s = %.2f" % (
!         d_min, times[0], times[1], times[0]/times[1])
        sys.stdout.flush()
    print
--- 31,38 ----
          times.append(f_calc_object.manager().estimate_time_fft.time_sampling)
          print "  %.2f seconds," % times[-1]
!       print "d_min=%d: %.2f s / %.2f s" % (d_min, times[0], times[1]),
!       if (times[1] != 0):
!         print "= %.2f" % times[1],
!       print
        sys.stdout.flush()
    print

Update of /cvsroot/cctbx/iotbx/iotbx/command_line
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7326/iotbx/command_line

Modified Files:
	show_distances.py 
Log Message:
cctbx.xray.structure:
  - show_pairs -> cctbx.crystal.show_distances
  - new show_distances() method

cctbx.crystal.pair_asu_table:
  - new show_distances() method
  - new pair_counts() method

cctbx.crystal.neighbors_simple_pair_generator:
cctbx.crystal.neighbors_fast_pair_generator:
  - new max_distance_sq() method

scitbx.math.icosahedron:
  - new next_neighbors_distance() method
  - slight modification of make_icosahedron() algorithm:
      mid points are now immediately normalized; this leads to a
      slightly more uniform distribution of next-neighbor distances


Index: show_distances.py
===================================================================
RCS file: /cvsroot/cctbx/iotbx/iotbx/command_line/show_distances.py,v
retrieving revision 1.9
retrieving revision 1.10
diff -C2 -d -r1.9 -r1.10
*** show_distances.py	1 Oct 2004 13:45:16 -0000	1.9
--- show_distances.py	4 Feb 2005 02:58:51 -0000	1.10
***************
*** 16,20 ****
    xray_structure.show_summary().show_scatterers()
    print
!   pairs = xray_structure.show_pairs(
      distance_cutoff=distance_cutoff,
      show_cartesian=show_cartesian,
--- 16,20 ----
    xray_structure.show_summary().show_scatterers()
    print
!   pairs = xray_structure.show_distances(
      distance_cutoff=distance_cutoff,
      show_cartesian=show_cartesian,

Update of /cvsroot/cctbx/iotbx/iotbx/pdb
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7326/iotbx/pdb

Modified Files:
	tst_pdb.py 
Log Message:
cctbx.xray.structure:
  - show_pairs -> cctbx.crystal.show_distances
  - new show_distances() method

cctbx.crystal.pair_asu_table:
  - new show_distances() method
  - new pair_counts() method

cctbx.crystal.neighbors_simple_pair_generator:
cctbx.crystal.neighbors_fast_pair_generator:
  - new max_distance_sq() method

scitbx.math.icosahedron:
  - new next_neighbors_distance() method
  - slight modification of make_icosahedron() algorithm:
      mid points are now immediately normalized; this leads to a
      slightly more uniform distribution of next-neighbor distances


Index: tst_pdb.py
===================================================================
RCS file: /cvsroot/cctbx/iotbx/iotbx/pdb/tst_pdb.py,v
retrieving revision 1.27
retrieving revision 1.28
diff -C2 -d -r1.27 -r1.28
*** tst_pdb.py	22 Jan 2005 14:18:26 -0000	1.27
--- tst_pdb.py	4 Feb 2005 02:58:51 -0000	1.28
***************
*** 9,12 ****
--- 9,13 ----
  from cctbx.development import random_structure
  from cctbx.array_family import flex
+ import scitbx.math
  from cStringIO import StringIO
  import sys, os
***************
*** 455,458 ****
--- 456,469 ----
        assert abs(regression.y_intercept()) < 0.1
  
+ def write_icosahedron():
+   for level in xrange(3):
+     icosahedron = scitbx.math.icosahedron(level=level)
+     scale = 1.5/icosahedron.next_neighbors_distance()
+     f = open("icosahedron_%d.pdb"%level, "w")
+     for i,site in enumerate(icosahedron.sites*scale):
+       print >> f, iotbx.pdb.format_atom_record(serial=i+1, site=site)
+     print >> f, "END"
+     f.close()
+ 
  def run():
    verbose = "--verbose" in sys.argv[1:]
***************
*** 467,470 ****
--- 478,482 ----
    for anisotropic_flag in (False, True):
      exercise_xray_structure(anisotropic_flag, verbose=verbose)
+   write_icosahedron()
    print "OK"
  

Update of /cvsroot/cctbx/cctbx/cctbx/geometry_restraints
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7086/cctbx/geometry_restraints

Modified Files:
	distance_least_squares.py 
Log Message:
cctbx.xray.structure:
  - show_pairs -> cctbx.crystal.show_distances
  - new show_distances() method

cctbx.crystal.pair_asu_table:
  - new show_distances() method
  - new pair_counts() method

cctbx.crystal.neighbors_simple_pair_generator:
cctbx.crystal.neighbors_fast_pair_generator:
  - new max_distance_sq() method

scitbx.math.icosahedron:
  - new next_neighbors_distance() method
  - slight modification of make_icosahedron() algorithm:
      mid points are now immediately normalized; this leads to a
      slightly more uniform distribution of next-neighbor distances


Index: distance_least_squares.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/geometry_restraints/distance_least_squares.py,v
retrieving revision 1.7
retrieving revision 1.8
diff -C2 -d -r1.7 -r1.8
*** distance_least_squares.py	20 Sep 2004 19:39:05 -0000	1.7
--- distance_least_squares.py	4 Feb 2005 02:57:55 -0000	1.8
***************
*** 129,135 ****
      asu_mappings=si_asu_mappings)
    si_pair_asu_table.add_all_pairs(distance_cutoff=distance_cutoff)
!   si_pairs = xray.show_pairs(
!     xray_structure=si_structure,
!     pair_asu_table=si_pair_asu_table)
    if (connectivities is not None):
      assert list(si_pairs.pair_counts) == connectivities
--- 129,133 ----
      asu_mappings=si_asu_mappings)
    si_pair_asu_table.add_all_pairs(distance_cutoff=distance_cutoff)
!   si_pairs = si_structure.show_distances(pair_asu_table=si_pair_asu_table)
    if (connectivities is not None):
      assert list(si_pairs.pair_counts) == connectivities
***************
*** 145,150 ****
      asu_mappings=si_o_asu_mappings)
    si_o_bond_asu_table.add_pair_sym_table(sym_table=si_o.bond_sym_table)
!   si_o_bonds = xray.show_pairs(
!     xray_structure=si_o.structure,
      pair_asu_table=si_o_bond_asu_table)
    n_si = si_pairs.pair_counts.size()
--- 143,147 ----
      asu_mappings=si_o_asu_mappings)
    si_o_bond_asu_table.add_pair_sym_table(sym_table=si_o.bond_sym_table)
!   si_o_bonds = si_o.structure.show_distances(
      pair_asu_table=si_o_bond_asu_table)
    n_si = si_pairs.pair_counts.size()
***************
*** 156,162 ****
      si_o_structure=si_o.structure,
      si_o_bond_asu_table=si_o_bond_asu_table)
!   o_si_o_pairs = xray.show_pairs(
!     xray_structure=si_o.structure,
!     pair_asu_table=o_si_o_asu_table)
    assert o_si_o_pairs.pair_counts[:n_si].all_eq(0)
    if (si_pairs.pair_counts.count(4) == n_si):
--- 153,157 ----
      si_o_structure=si_o.structure,
      si_o_bond_asu_table=si_o_bond_asu_table)
!   o_si_o_pairs = si_o.structure.show_distances(pair_asu_table=o_si_o_asu_table)
    assert o_si_o_pairs.pair_counts[:n_si].all_eq(0)
    if (si_pairs.pair_counts.count(4) == n_si):
***************
*** 260,266 ****
        minimized[1].final_target_result.target())
    print
!   xray.show_pairs(
!     xray_structure=minimized_structure,
!     pair_asu_table=si_o_bond_asu_table)
    print
    sites_cart = minimized_structure.sites_cart()
--- 255,259 ----
        minimized[1].final_target_result.target())
    print
!   minimized_structure.show_distances(pair_asu_table=si_o_bond_asu_table)
    print
    sites_cart = minimized_structure.sites_cart()

Update of /cvsroot/cctbx/cctbx/cctbx/crystal
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7086/cctbx/crystal

Modified Files:
	__init__.py 
Log Message:
cctbx.xray.structure:
  - show_pairs -> cctbx.crystal.show_distances
  - new show_distances() method

cctbx.crystal.pair_asu_table:
  - new show_distances() method
  - new pair_counts() method

cctbx.crystal.neighbors_simple_pair_generator:
cctbx.crystal.neighbors_fast_pair_generator:
  - new max_distance_sq() method

scitbx.math.icosahedron:
  - new next_neighbors_distance() method
  - slight modification of make_icosahedron() algorithm:
      mid points are now immediately normalized; this leads to a
      slightly more uniform distribution of next-neighbor distances


Index: __init__.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/crystal/__init__.py,v
retrieving revision 1.34
retrieving revision 1.35
diff -C2 -d -r1.34 -r1.35
*** __init__.py	2 Feb 2005 09:48:46 -0000	1.34
--- __init__.py	4 Feb 2005 02:57:55 -0000	1.35
***************
*** 335,338 ****
--- 335,443 ----
              print >> f, "    j_syms:", list(j_syms)
  
+   def show_distances(self,
+         site_labels=None,
+         sites_frac=None,
+         sites_cart=None,
+         show_cartesian=False,
+         keep_pair_asu_table=False,
+         out=None):
+     return show_distances(
+       pair_asu_table=self,
+       site_labels=site_labels,
+       sites_frac=sites_frac,
+       sites_cart=sites_cart,
+       show_cartesian=show_cartesian,
+       keep_pair_asu_table=keep_pair_asu_table,
+       out=out)
+ 
+ class show_distances:
+ 
+   def __init__(self,
+         pair_asu_table,
+         site_labels=None,
+         sites_frac=None,
+         sites_cart=None,
+         show_cartesian=False,
+         keep_pair_asu_table=False,
+         out=None):
+     assert [sites_frac, sites_cart].count(None) == 1
+     if (out is None): out = sys.stdout
+     if (keep_pair_asu_table):
+       self.pair_asu_table = pair_asu_table
+     else:
+       self.pair_asu_table = None
+     self.distances = flex.double()
+     self.pair_counts = flex.size_t()
+     asu_mappings = pair_asu_table.asu_mappings()
+     unit_cell = asu_mappings.unit_cell()
+     if (sites_frac is None):
+       sites_frac = unit_cell.fractionalization_matrix() * sites_cart
+     if (site_labels is None):
+       label_len = len("%d" % (sites_frac.size()+1))
+       label_fmt = "site_%%0%dd" % label_len
+       label_len += 5
+     else:
+       label_len = 1
+       for label in site_labels:
+         label_len = max(label_len, len(label))
+       label_fmt = "%%-%ds" % (label_len+1)
+     for i_seq,asu_dict in enumerate(pair_asu_table.table()):
+       rt_mx_i_inv = asu_mappings.get_rt_mx(i_seq, 0).inverse()
+       site_frac_i = sites_frac[i_seq]
+       pair_count = 0
+       dists = flex.double()
+       j_seq_i_group = []
+       for j_seq,j_sym_groups in asu_dict.items():
+         site_frac_j = sites_frac[j_seq]
+         for i_group,j_sym_group in enumerate(j_sym_groups):
+           pair_count += j_sym_group.size()
+           j_sym = j_sym_group[0]
+           rt_mx_ji = rt_mx_i_inv.multiply(asu_mappings.get_rt_mx(j_seq, j_sym))
+           distance = unit_cell.distance(site_frac_i, rt_mx_ji * site_frac_j)
+           dists.append(distance)
+           j_seq_i_group.append((j_seq,i_group))
+       if (site_labels is None):
+         s = label_fmt % (i_seq+1)
+       else:
+         s = label_fmt % site_labels[i_seq]
+       s += " pair count: %3d" % pair_count
+       if (show_cartesian):
+         formatted_site = [" %7.2f" % x
+           for x in unit_cell.orthogonalize(site_frac_i)]
+       else:
+         formatted_site = [" %7.4f" % x for x in site_frac_i]
+       print >> out, ("%%-%ds" % (label_len+23)) % s, \
+         "<<"+",".join(formatted_site)+">>"
+       permutation = flex.sort_permutation(data=dists)
+       for j_seq,i_group in flex.select(j_seq_i_group, permutation):
+         site_frac_j = sites_frac[j_seq]
+         j_sym_groups = asu_dict[j_seq]
+         j_sym_group = j_sym_groups[i_group]
+         for i_j_sym,j_sym in enumerate(j_sym_group):
+           rt_mx_ji = rt_mx_i_inv.multiply(
+             asu_mappings.get_rt_mx(j_seq, j_sym))
+           site_frac_ji = rt_mx_ji * site_frac_j
+           distance = unit_cell.distance(site_frac_i, site_frac_ji)
+           self.distances.append(distance)
+           if (site_labels is None):
+             print >> out, " ", label_fmt % (j_seq+1) + ":",
+           else:
+             print >> out, " ", label_fmt % (site_labels[j_seq] + ":"),
+           print >> out, "%8.4f" % distance,
+           if (i_j_sym != 0):
+             s = "sym. equiv."
+           else:
+             s = "           "
+           if (show_cartesian):
+             formatted_site = [" %7.2f" % x
+               for x in unit_cell.orthogonalize(site_frac_ji)]
+           else:
+             formatted_site = [" %7.4f" % x for x in site_frac_ji]
+           s += " (" + ",".join(formatted_site) +")"
+           print >> out, s
+       if (pair_count == 0):
+         print >> out, "  no neighbors"
+       self.pair_counts.append(pair_count)
+ 
  class sym_pair:
  

Update of /cvsroot/cctbx/cctbx/include/cctbx/crystal
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7086/include/cctbx/crystal

Modified Files:
	neighbors_bpl.cpp neighbors_fast.h neighbors_simple.h 
	pair_tables.h pair_tables_bpl.cpp tst_ext.py 
Log Message:
cctbx.xray.structure:
  - show_pairs -> cctbx.crystal.show_distances
  - new show_distances() method

cctbx.crystal.pair_asu_table:
  - new show_distances() method
  - new pair_counts() method

cctbx.crystal.neighbors_simple_pair_generator:
cctbx.crystal.neighbors_fast_pair_generator:
  - new max_distance_sq() method

scitbx.math.icosahedron:
  - new next_neighbors_distance() method
  - slight modification of make_icosahedron() algorithm:
      mid points are now immediately normalized; this leads to a
      slightly more uniform distribution of next-neighbor distances


Index: neighbors_bpl.cpp
===================================================================
RCS file: /cvsroot/cctbx/cctbx/include/cctbx/crystal/neighbors_bpl.cpp,v
retrieving revision 1.4
retrieving revision 1.5
diff -C2 -d -r1.4 -r1.5
*** neighbors_bpl.cpp	25 Sep 2004 08:12:29 -0000	1.4
--- neighbors_bpl.cpp	4 Feb 2005 02:57:57 -0000	1.5
***************
*** 47,50 ****
--- 47,51 ----
          .def("restart", &w_t::restart)
          .def("count_pairs", &w_t::count_pairs)
+         .def("max_distance_sq", &w_t::max_distance_sq)
          .def("neighbors_of", &w_t::neighbors_of, (arg_("primary_selection")))
          .def("is_simple_interaction", &w_t::is_simple_interaction)
***************
*** 76,79 ****
--- 77,81 ----
          .def("restart", &w_t::restart)
          .def("count_pairs", &w_t::count_pairs)
+         .def("max_distance_sq", &w_t::max_distance_sq)
          .def("neighbors_of", &w_t::neighbors_of, (arg_("primary_selection")))
        ;

Index: neighbors_fast.h
===================================================================
RCS file: /cvsroot/cctbx/cctbx/include/cctbx/crystal/neighbors_fast.h,v
retrieving revision 1.12
retrieving revision 1.13
diff -C2 -d -r1.12 -r1.13
*** neighbors_fast.h	14 Oct 2004 01:36:11 -0000	1.12
--- neighbors_fast.h	4 Feb 2005 02:57:57 -0000	1.13
***************
*** 104,108 ****
        }
  
!       //! Count the number of pairs.
        std::size_t
        count_pairs()
--- 104,108 ----
        }
  
!       //! Counts the number of pairs.
        std::size_t
        count_pairs()
***************
*** 116,119 ****
--- 116,130 ----
        }
  
+       //! Maximum distance squared of all remaining pairs.
+       FloatType
+       max_distance_sq()
+       {
+         FloatType result = -1;
+         while (!this->at_end_) {
+           result = std::max(result, next().dist_sq);
+         }
+         return result;
+       }
+ 
        /*! \brief Selection of all neighbors within distance_cutoff
            of primary_selection.

Index: pair_tables_bpl.cpp
===================================================================
RCS file: /cvsroot/cctbx/cctbx/include/cctbx/crystal/pair_tables_bpl.cpp,v
retrieving revision 1.10
retrieving revision 1.11
diff -C2 -d -r1.10 -r1.11
*** pair_tables_bpl.cpp	10 Dec 2004 01:35:15 -0000	1.10
--- pair_tables_bpl.cpp	4 Feb 2005 02:57:57 -0000	1.11
***************
*** 82,85 ****
--- 82,86 ----
          .def("__eq__", &w_t::operator==)
          .def("__ne__", &w_t::operator!=)
+         .def("pair_counts", &w_t::pair_counts)
          .def("add_all_pairs", &w_t::add_all_pairs,
            add_all_pairs_overloads((

Index: pair_tables.h
===================================================================
RCS file: /cvsroot/cctbx/cctbx/include/cctbx/crystal/pair_tables.h,v
retrieving revision 1.14
retrieving revision 1.15
diff -C2 -d -r1.14 -r1.15
*** pair_tables.h	10 Dec 2004 01:35:14 -0000	1.14
--- pair_tables.h	4 Feb 2005 02:57:57 -0000	1.15
***************
*** 194,197 ****
--- 194,220 ----
        }
  
+       //! Pair count for each site.
+       af::shared<std::size_t>
+       pair_counts() const
+       {
+         af::const_ref<pair_asu_dict> tab = table_.const_ref();
+         af::shared<std::size_t> result((af::reserve(tab.size())));
+         for(unsigned i_seq=0;i_seq<tab.size();i_seq++) {
+           std::size_t count = 0;
+           pair_asu_dict const& dict = tab[i_seq];
+           for(pair_asu_dict::const_iterator
+                 dict_i=dict.begin();
+                 dict_i!=dict.end();
+                 dict_i++) {
+             pair_asu_j_sym_groups const& jgs = dict_i->second;
+             for(unsigned ig=0;ig<jgs.size();ig++) {
+               count += jgs[ig].size();
+             }
+           }
+           result.push_back(count);
+         }
+         return result;
+       }
+ 
        /*! \brief Uses neighbors::fast_pair_generator to add all pairs with
            distances <= distance_cutoff*(1+epsilon).

Index: neighbors_simple.h
===================================================================
RCS file: /cvsroot/cctbx/cctbx/include/cctbx/crystal/neighbors_simple.h,v
retrieving revision 1.18
retrieving revision 1.19
diff -C2 -d -r1.18 -r1.19
*** neighbors_simple.h	25 Sep 2004 08:12:29 -0000	1.18
--- neighbors_simple.h	4 Feb 2005 02:57:57 -0000	1.19
***************
*** 113,116 ****
--- 113,127 ----
        }
  
+       //! Maximum distance squared of all remaining pairs.
+       FloatType
+       max_distance_sq()
+       {
+         FloatType result = -1;
+         while (!this->at_end_) {
+           result = std::max(result, next().dist_sq);
+         }
+         return result;
+       }
+ 
        /*! \brief Selection of all neighbors within distance_cutoff
            of primary_selection.

Index: tst_ext.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/include/cctbx/crystal/tst_ext.py,v
retrieving revision 1.18
retrieving revision 1.19
diff -C2 -d -r1.18 -r1.19
*** tst_ext.py	9 Dec 2004 02:38:23 -0000	1.18
--- tst_ext.py	4 Feb 2005 02:57:57 -0000	1.19
***************
*** 9,12 ****
--- 9,13 ----
  from libtbx.itertbx import count
  from cStringIO import StringIO
+ import md5
  import sys
  
***************
*** 285,291 ****
--- 286,298 ----
          assert approx_equal(dist_sq_sorted, short_dist_sq_sorted)
        assert pair_generator.count_pairs() == 0
+       assert pair_generator.max_distance_sq() == -1
        pair_generator.restart()
        assert pair_generator.count_pairs() == len(index_pairs)
        pair_generator.restart()
+       if (len(index_pairs) == 0):
+         assert pair_generator.max_distance_sq() == -1
+       else:
+         assert approx_equal(pair_generator.max_distance_sq(), 0.1)
+       pair_generator.restart()
        primary_selection = flex.bool(
          pair_generator.asu_mappings().mappings().size(), False)
***************
*** 417,420 ****
--- 424,428 ----
    assert asu_table.table().size() == 3
    assert not asu_table.contains(i_seq=0, j_seq=1, j_sym=2)
+   assert asu_table.pair_counts().all_eq(0)
    assert asu_table.add_all_pairs(distance_cutoff=3.5) is asu_table
    assert [d.size() for d in asu_table.table()] == [2,3,2]
***************
*** 434,437 ****
--- 442,446 ----
      for pair in pair_generator:
        asu_table.add_pair(pair=pair)
+     asu_table.pair_counts().all_eq(4)
      check_pair_asu_table(asu_table, expected_asu_pairs)
    for p in expected_asu_pairs:
***************
*** 783,800 ****
      max_shell=3)
    check_pair_asu_table(shell_asu_tables[0], expected_asu_pairs)
!   if (0):
!     from cctbx.crystal import distance_ls
!     distance_ls.show_pairs(
!       structure=structure,
!       pair_asu_table=shell_asu_tables[1])
!     print
!     s1_sym_table = shell_asu_tables[1].extract_pair_sym_table(
!       skip_j_seq_less_than_i_seq=False)
!     s1_asu_table = crystal.pair_asu_table(asu_mappings=asu_mappings)
!     s1_asu_table.add_pair_sym_table(sym_table=s1_sym_table)
!     distance_ls.show_pairs(
!       structure=structure,
!       pair_asu_table=s1_asu_table)
!     print
  
  def exercise_symmetry():
--- 792,808 ----
      max_shell=3)
    check_pair_asu_table(shell_asu_tables[0], expected_asu_pairs)
!   s = StringIO()
!   structure.show_distances(pair_asu_table=shell_asu_tables[1], out=s)
!   print >> s
!   s1_sym_table = shell_asu_tables[1].extract_pair_sym_table(
!     skip_j_seq_less_than_i_seq=False)
!   s1_asu_table = crystal.pair_asu_table(asu_mappings=asu_mappings)
!   s1_asu_table.add_pair_sym_table(sym_table=s1_sym_table)
!   structure.show_distances(pair_asu_table=s1_asu_table, out=s)
!   print >> s
!   s = s.getvalue().replace("-0.0000", " 0.0000")
!   if (md5.md5(s).hexdigest() != "2644729423835c824fafc53919d2ccf4"):
!     sys.stderr.write(s)
!     raise AssertionError("Unexpected show_distances() output.")
  
  def exercise_symmetry():

Update of /cvsroot/cctbx/cctbx/cctbx/regression
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7086/cctbx/regression

Modified Files:
	tst_coordination_sequences.py tst_pair_asu_table.py 
	tst_xray.py 
Log Message:
cctbx.xray.structure:
  - show_pairs -> cctbx.crystal.show_distances
  - new show_distances() method

cctbx.crystal.pair_asu_table:
  - new show_distances() method
  - new pair_counts() method

cctbx.crystal.neighbors_simple_pair_generator:
cctbx.crystal.neighbors_fast_pair_generator:
  - new max_distance_sq() method

scitbx.math.icosahedron:
  - new next_neighbors_distance() method
  - slight modification of make_icosahedron() algorithm:
      mid points are now immediately normalized; this leads to a
      slightly more uniform distribution of next-neighbor distances


Index: tst_coordination_sequences.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/regression/tst_coordination_sequences.py,v
retrieving revision 1.9
retrieving revision 1.10
diff -C2 -d -r1.9 -r1.10
*** tst_coordination_sequences.py	22 Jan 2005 14:25:01 -0000	1.9
--- tst_coordination_sequences.py	4 Feb 2005 02:57:56 -0000	1.10
***************
*** 42,49 ****
    for shell_asu_table in shell_asu_tables:
      if (0 or verbose):
!       pairs_1 = xray.show_pairs(
!         xray_structure=structure,
!         pair_asu_table=shell_asu_table)
        print list(pairs_1.pair_counts)
        print
      sym_table = shell_asu_table.extract_pair_sym_table(
--- 42,48 ----
    for shell_asu_table in shell_asu_tables:
      if (0 or verbose):
!       pairs_1 = structure.show_distances(pair_asu_table=shell_asu_table)
        print list(pairs_1.pair_counts)
+       assert pairs_1.pair_counts == shell_asu_table.pair_counts()
        print
      sym_table = shell_asu_table.extract_pair_sym_table(
***************
*** 52,59 ****
      asu_table.add_pair_sym_table(sym_table=sym_table)
      if (0 or verbose):
!       pairs_2 = xray.show_pairs(
!         xray_structure=structure,
!         pair_asu_table=asu_table)
        print list(pairs_2.pair_counts)
        print
        assert pairs_1.pair_counts.all_eq(pairs_2.pair_counts)
--- 51,57 ----
      asu_table.add_pair_sym_table(sym_table=sym_table)
      if (0 or verbose):
!       pairs_2 = structure.show_distances(pair_asu_table=asu_table)
        print list(pairs_2.pair_counts)
+       assert pairs_2.pair_counts == asu_table.pair_counts()
        print
        assert pairs_1.pair_counts.all_eq(pairs_2.pair_counts)

Index: tst_xray.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/regression/tst_xray.py,v
retrieving revision 1.54
retrieving revision 1.55
diff -C2 -d -r1.54 -r1.55
*** tst_xray.py	2 Feb 2005 10:39:13 -0000	1.54
--- tst_xray.py	4 Feb 2005 02:57:56 -0000	1.55
***************
*** 247,251 ****
    xs2.apply_symmetry_u_stars()
    s = StringIO()
!   xs1.show_pairs(distance_cutoff=0.1, out=s)
    assert s.getvalue() == """\
  C  pair count:   2       <<  0.0200,  0.0000,  0.1000>>
--- 247,251 ----
    xs2.apply_symmetry_u_stars()
    s = StringIO()
!   xs1.show_distances(distance_cutoff=0.1, out=s)
    assert s.getvalue() == """\
  C  pair count:   2       <<  0.0200,  0.0000,  0.1000>>

Index: tst_pair_asu_table.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/regression/tst_pair_asu_table.py,v
retrieving revision 1.5
retrieving revision 1.6
diff -C2 -d -r1.5 -r1.6
*** tst_pair_asu_table.py	22 Jan 2005 14:25:01 -0000	1.5
--- tst_pair_asu_table.py	4 Feb 2005 02:57:56 -0000	1.6
***************
*** 2,9 ****
--- 2,57 ----
  from cctbx import geometry_restraints
  from cctbx import crystal
+ from cctbx.array_family import flex
  from scitbx.python_utils.misc import adopt_init_args
+ import scitbx.math
+ from libtbx.test_utils import approx_equal
+ from cStringIO import StringIO
  import math
  import sys, os
  
+ def exercise_icosahedron(max_level=2, verbose=0):
+   for level in xrange(0,max_level+1):
+     if (0 or verbose):
+       print "level:", level
+     icosahedron = scitbx.math.icosahedron(level=level)
+     try:
+       distance_cutoff = icosahedron.next_neighbors_distance()*(1+1.e-3)
+       estimated_distance_cutoff = False
+     except RuntimeError, e:
+       assert str(e) == "next_neighbors_distance not known."
+       distance_cutoff = 0.4/(2**(level-1))
+       estimated_distance_cutoff = True
+     asu_mappings = crystal.direct_space_asu.non_crystallographic_asu_mappings(
+       sites_cart=icosahedron.sites)
+     pair_asu_table = crystal.pair_asu_table(asu_mappings=asu_mappings)
+     pair_asu_table.add_all_pairs(distance_cutoff=distance_cutoff)
+     if (0 or verbose):
+       ps = pair_asu_table.show_distances(sites_cart=icosahedron.sites)
+       print "level", level, "min", flex.min(ps.distances)
+       print "     ", " ",   "max", flex.max(ps.distances)
+       assert ps.pair_counts.all_eq(pair_asu_table.pair_counts())
+       if (level == 0):
+         for d in ps.distances:
+           assert approx_equal(d, 1.0514622242382672)
+     elif (level < 2):
+       s = StringIO()
+       ps = pair_asu_table.show_distances(sites_cart=icosahedron.sites, out=s)
+       assert ps.pair_counts.all_eq(pair_asu_table.pair_counts())
+       assert len(s.getvalue().splitlines()) == [72,320][level]
+       del s
+     if (level == 0):
+       assert pair_asu_table.pair_counts().all_eq(5)
+     else:
+       assert pair_asu_table.pair_counts().all_eq(3)
+     del pair_asu_table
+     max_distance = crystal.neighbors_fast_pair_generator(
+       asu_mappings=asu_mappings,
+       distance_cutoff=distance_cutoff).max_distance_sq()**.5
+     if (0 or verbose):
+       print "max_distance:", max_distance
+     if (not estimated_distance_cutoff):
+       assert approx_equal(max_distance, icosahedron.next_neighbors_distance())
+       assert approx_equal(max_distance/icosahedron.next_neighbors_distance(),1)
+ 
  def is_sym_equiv_interaction_simple(unit_cell,
                                      i_seq,
***************
*** 150,153 ****
--- 198,202 ----
  def run():
    verbose = "--verbose" in sys.argv[1:]
+   exercise_icosahedron(verbose=verbose)
    default_distance_cutoff=3.5
    regression_misc = os.path.join(

Update of /cvsroot/cctbx/cctbx/cctbx/xray
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7086/cctbx/xray

Modified Files:
	structure.py 
Log Message:
cctbx.xray.structure:
  - show_pairs -> cctbx.crystal.show_distances
  - new show_distances() method

cctbx.crystal.pair_asu_table:
  - new show_distances() method
  - new pair_counts() method

cctbx.crystal.neighbors_simple_pair_generator:
cctbx.crystal.neighbors_fast_pair_generator:
  - new max_distance_sq() method

scitbx.math.icosahedron:
  - new next_neighbors_distance() method
  - slight modification of make_icosahedron() algorithm:
      mid points are now immediately normalized; this leads to a
      slightly more uniform distribution of next-neighbor distances


Index: structure.py
===================================================================
RCS file: /cvsroot/cctbx/cctbx/cctbx/xray/structure.py,v
retrieving revision 1.48
retrieving revision 1.49
diff -C2 -d -r1.48 -r1.49
*** structure.py	2 Feb 2005 10:39:14 -0000	1.48
--- structure.py	4 Feb 2005 02:57:56 -0000	1.49
***************
*** 455,472 ****
      return pair_asu_table
  
!   def show_pairs(self,
!         distance_cutoff,
          asu_mappings_buffer_thickness=None,
          asu_mappings_is_inside_epsilon=None,
          show_cartesian=False,
          keep_pair_asu_table=False,
          out=None):
!     assert distance_cutoff is not None
!     return show_pairs(
!       xray_structure=self,
!       pair_asu_table=self.pair_asu_table(
          distance_cutoff=distance_cutoff,
          asu_mappings_buffer_thickness=asu_mappings_buffer_thickness,
!         asu_mappings_is_inside_epsilon=asu_mappings_is_inside_epsilon),
        show_cartesian=show_cartesian,
        keep_pair_asu_table=keep_pair_asu_table,
--- 455,475 ----
      return pair_asu_table
  
!   def show_distances(self,
!         distance_cutoff=None,
          asu_mappings_buffer_thickness=None,
          asu_mappings_is_inside_epsilon=None,
+         pair_asu_table=None,
          show_cartesian=False,
          keep_pair_asu_table=False,
          out=None):
!     assert [distance_cutoff, pair_asu_table].count(None) == 1
!     if (pair_asu_table is None):
!       pair_asu_table = self.pair_asu_table(
          distance_cutoff=distance_cutoff,
          asu_mappings_buffer_thickness=asu_mappings_buffer_thickness,
!         asu_mappings_is_inside_epsilon=asu_mappings_is_inside_epsilon)
!     return pair_asu_table.show_distances(
!       site_labels=self.scatterers().extract_labels(),
!       sites_frac=self.sites_frac(),
        show_cartesian=show_cartesian,
        keep_pair_asu_table=keep_pair_asu_table,
***************
*** 478,554 ****
    def rms_difference(self, other):
      return self.sites_cart().rms_difference(other.sites_cart())
- 
- class show_pairs:
- 
-   def __init__(self,
-         xray_structure,
-         pair_asu_table,
-         show_cartesian=False,
-         keep_pair_asu_table=False,
-         out=None):
-     if (out is None): out = sys.stdout
-     if (keep_pair_asu_table):
-       self.pair_asu_table = pair_asu_table
-     else:
-       self.pair_asu_table = None
-     self.distances = flex.double()
-     self.pair_counts = flex.size_t()
-     label_len = 1
-     for scatterer in xray_structure.scatterers():
-       label_len = max(label_len, len(scatterer.label))
-     label_fmt = "%%-%ds" % (label_len+1)
-     unit_cell = xray_structure.unit_cell()
-     scatterers = xray_structure.scatterers()
-     sites_frac = xray_structure.sites_frac()
-     asu_mappings = pair_asu_table.asu_mappings()
-     for i_seq,asu_dict in enumerate(pair_asu_table.table()):
-       rt_mx_i_inv = asu_mappings.get_rt_mx(i_seq, 0).inverse()
-       site_frac_i = sites_frac[i_seq]
-       pair_count = 0
-       dists = flex.double()
-       j_seq_i_group = []
-       for j_seq,j_sym_groups in asu_dict.items():
-         site_frac_j = sites_frac[j_seq]
-         for i_group,j_sym_group in enumerate(j_sym_groups):
-           pair_count += j_sym_group.size()
-           j_sym = j_sym_group[0]
-           rt_mx_ji = rt_mx_i_inv.multiply(asu_mappings.get_rt_mx(j_seq, j_sym))
-           distance = unit_cell.distance(site_frac_i, rt_mx_ji * site_frac_j)
-           dists.append(distance)
-           j_seq_i_group.append((j_seq,i_group))
-       s = label_fmt % scatterers[i_seq].label \
-         + " pair count: %3d" % pair_count
-       if (show_cartesian):
-         formatted_site = [" %7.2f" % x
-           for x in unit_cell.orthogonalize(site_frac_i)]
-       else:
-         formatted_site = [" %7.4f" % x for x in site_frac_i]
-       print >> out, ("%%-%ds" % (label_len+23)) % s, \
-         "<<"+",".join(formatted_site)+">>"
-       permutation = flex.sort_permutation(data=dists)
-       for j_seq,i_group in flex.select(j_seq_i_group, permutation):
-         site_frac_j = sites_frac[j_seq]
-         j_sym_groups = asu_dict[j_seq]
-         j_sym_group = j_sym_groups[i_group]
-         for i_j_sym,j_sym in enumerate(j_sym_group):
-           rt_mx_ji = rt_mx_i_inv.multiply(
-             asu_mappings.get_rt_mx(j_seq, j_sym))
-           site_frac_ji = rt_mx_ji * site_frac_j
-           distance = unit_cell.distance(site_frac_i, site_frac_ji)
-           self.distances.append(distance)
-           print >> out, " ", label_fmt % (scatterers[j_seq].label + ":"),
-           print >> out, "%8.4f" % distance,
-           if (i_j_sym != 0):
-             s = "sym. equiv."
-           else:
-             s = "           "
-           if (show_cartesian):
-             formatted_site = [" %7.2f" % x
-               for x in unit_cell.orthogonalize(site_frac_ji)]
-           else:
-             formatted_site = [" %7.4f" % x for x in site_frac_ji]
-           s += " (" + ",".join(formatted_site) +")"
-           print >> out, s
-       if (pair_count == 0):
-         print >> out, "  no neighbors"
-       self.pair_counts.append(pair_count)
--- 481,482 ----

Update of /cvsroot/cctbx/scitbx/math/boost_python
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7015/math/boost_python

Modified Files:
	icosahedron.cpp tst_math.py 
Log Message:
cctbx.xray.structure:
  - show_pairs -> cctbx.crystal.show_distances
  - new show_distances() method

cctbx.crystal.pair_asu_table:
  - new show_distances() method
  - new pair_counts() method

cctbx.crystal.neighbors_simple_pair_generator:
cctbx.crystal.neighbors_fast_pair_generator:
  - new max_distance_sq() method

scitbx.math.icosahedron:
  - new next_neighbors_distance() method
  - slight modification of make_icosahedron() algorithm:
      mid points are now immediately normalized; this leads to a
      slightly more uniform distribution of next-neighbor distances


Index: icosahedron.cpp
===================================================================
RCS file: /cvsroot/cctbx/scitbx/math/boost_python/icosahedron.cpp,v
retrieving revision 1.2
retrieving revision 1.3
diff -C2 -d -r1.2 -r1.3
*** icosahedron.cpp	3 Feb 2005 14:00:05 -0000	1.2
--- icosahedron.cpp	4 Feb 2005 02:57:41 -0000	1.3
***************
*** 24,27 ****
--- 24,28 ----
          .def_readonly("level", &w_t::level)
          .add_property("sites", make_getter(&w_t::sites, rbv()))
+         .def("next_neighbors_distance", &w_t::next_neighbors_distance)
        ;
      }

Index: tst_math.py
===================================================================
RCS file: /cvsroot/cctbx/scitbx/math/boost_python/tst_math.py,v
retrieving revision 1.38
retrieving revision 1.39
diff -C2 -d -r1.38 -r1.39
*** tst_math.py	3 Feb 2005 14:00:05 -0000	1.38
--- tst_math.py	4 Feb 2005 02:57:41 -0000	1.39
***************
*** 941,1274 ****
  
  def exercise_icosahedron():
!   expected_sites_level_2=[
!         (-0.13039500000000001, -0.82034300000000004, -0.55680700000000005),
!         (-0.298564, -0.64480300000000002, -0.70362499999999994),
!         (-0.27460699999999999, -0.76751199999999997, -0.57923899999999995),
!         (-0.39410099999999998, -0.79938500000000001, -0.45350600000000002),
!         (-0.44216899999999998, -0.40295599999999998, -0.80132000000000003),
!         (-0.55441700000000005, -0.137651, -0.82077599999999995),
!         (-0.56907799999999997, -0.28447099999999997, -0.771509),
!         (-0.69738800000000001, -0.30169400000000002, -0.65010100000000004),
!         (-0.44827800000000001, -0.54707600000000001, -0.70693399999999995),
!         (-0.70693399999999995, -0.44827800000000001, -0.54707600000000001),
!         (-0.57735000000000003, -0.57735000000000003, -0.57735000000000003),
!         (-0.54707600000000001, -0.70693399999999995, -0.44827800000000001),
!         (-0.63921799999999995, -0.70696800000000004, -0.30264799999999997),
!         (-0.79803500000000005, -0.44999699999999998, -0.40080300000000002),
!         (-0.76232699999999998, -0.58046699999999996, -0.28621099999999999),
!         (-0.816797, -0.55912499999999998, -0.142204),
!         (0.0, -0.77053700000000003, -0.63739500000000004),
!         (-0.154582, -0.58980699999999997, -0.79261099999999995),
!         (0.0, -0.66262100000000002, -0.74895500000000004),
!         (0.154582, -0.58980699999999997, -0.79261099999999995),
!         (-0.310166, -0.35253499999999999, -0.88290199999999996),
!         (-0.42387200000000003, -0.087787000000000004, -0.90145799999999998),
!         (-0.30607400000000001, -0.18401300000000001, -0.93405499999999997),
!         (-0.15959100000000001, -0.096272999999999997, -0.98247799999999996),
!         (-0.159858, -0.43690899999999999, -0.88518699999999995),
!         (0.0, -0.17825299999999999, -0.983985),
!         (0.0, -0.35682199999999997, -0.934172),
!         (0.159858, -0.43690899999999999, -0.88518699999999995),
!         (0.30616599999999999, -0.345864, -0.88692800000000005),
!         (0.14734900000000001, -0.088893, -0.98508200000000001),
!         (0.29425600000000002, -0.17688799999999999, -0.93921500000000002),
!         (0.41692099999999999, -0.087887000000000007, -0.90468400000000004),
!         (-0.080588999999999994, -0.90093199999999996, -0.42641200000000001),
!         (-0.33910499999999999, -0.88837100000000002, -0.30952400000000002),
!         (-0.16971600000000001, -0.93722799999999995, -0.30463200000000001),
!         (-0.088985999999999996, -0.983908, -0.154942),
!         (-0.58344099999999999, -0.79470399999999997, -0.167458),
!         (-0.76652100000000001, -0.64220500000000003, -0.0043920000000000001),
!         (-0.65778300000000001, -0.75309099999999995, -0.013266999999999999),
!         (-0.59060100000000004, -0.79366099999999995, 0.14591999999999999),
!         (-0.43690899999999999, -0.88518699999999995, -0.159858),
!         (-0.43690899999999999, -0.88518699999999995, 0.159858),
!         (-0.35682199999999997, -0.934172, 0.0),
!         (-0.17825299999999999, -0.983985, 0.0),
!         (-0.088893, -0.98508200000000001, 0.14734900000000001),
!         (-0.345864, -0.88692800000000005, 0.30616599999999999),
!         (-0.17688799999999999, -0.93921500000000002, 0.29425600000000002),
!         (-0.087887000000000007, -0.90468400000000004, 0.41692099999999999),
!         (0.080588999999999994, -0.90093199999999996, -0.42641200000000001),
!         (0.088985999999999996, -0.983908, -0.154942),
!         (0.16971600000000001, -0.93722799999999995, -0.30463200000000001),
!         (0.33910499999999999, -0.88837100000000002, -0.30952400000000002),
!         (0.081582000000000002, -0.98639699999999997, 0.142708),
!         (0.080682000000000004, -0.90417199999999998, 0.41948000000000002),
!         (0.162546, -0.94225400000000004, 0.29280600000000001),
!         (0.33237699999999998, -0.89229400000000003, 0.30551099999999998),
!         (0.17825299999999999, -0.983985, 0.0),
!         (0.43690899999999999, -0.88518699999999995, 0.159858),
!         (0.35682199999999997, -0.934172, 0.0),
!         (0.43690899999999999, -0.88518699999999995, -0.159858),
!         (0.58427899999999999, -0.79586100000000004, -0.15881700000000001),
!         (0.58427899999999999, -0.79586100000000004, 0.15881700000000001),
!         (0.65300400000000003, -0.757355, 0.0),
!         (0.76248000000000005, -0.64701200000000003, 0.0),
!         (0.13039500000000001, -0.82034300000000004, -0.55680700000000005),
!         (0.298564, -0.64480300000000002, -0.70362499999999994),
!         (0.27460699999999999, -0.76751199999999997, -0.57923899999999995),
!         (0.39410099999999998, -0.79938500000000001, -0.45350600000000002),
!         (0.44216899999999998, -0.40295599999999998, -0.80132000000000003),
!         (0.55441700000000005, -0.137651, -0.82077599999999995),
!         (0.56907799999999997, -0.28447099999999997, -0.771509),
!         (0.69738800000000001, -0.30169400000000002, -0.65010100000000004),
!         (0.44827800000000001, -0.54707600000000001, -0.70693399999999995),
!         (0.70693399999999995, -0.44827800000000001, -0.54707600000000001),
!         (0.57735000000000003, -0.57735000000000003, -0.57735000000000003),
!         (0.54707600000000001, -0.70693399999999995, -0.44827800000000001),
!         (0.63921799999999995, -0.70696800000000004, -0.30264799999999997),
!         (0.79803500000000005, -0.44999699999999998, -0.40080300000000002),
!         (0.76232699999999998, -0.58046699999999996, -0.28621099999999999),
!         (0.816797, -0.55912499999999998, -0.142204),
!         (-0.63739500000000004, 0.0, -0.77053700000000003),
!         (-0.79261099999999995, -0.154582, -0.58980699999999997),
!         (-0.74895500000000004, 0.0, -0.66262100000000002),
!         (-0.79261099999999995, 0.154582, -0.58980699999999997),
!         (-0.88290199999999996, -0.310166, -0.35253499999999999),
!         (-0.90145799999999998, -0.42387200000000003, -0.087787000000000004),
!         (-0.93405499999999997, -0.30607400000000001, -0.18401300000000001),
!         (-0.98247799999999996, -0.15959100000000001, -0.096272999999999997),
!         (-0.88518699999999995, -0.159858, -0.43690899999999999),
!         (-0.983985, 0.0, -0.17825299999999999),
!         (-0.934172, 0.0, -0.35682199999999997),
!         (-0.88518699999999995, 0.159858, -0.43690899999999999),
!         (-0.88692800000000005, 0.30616599999999999, -0.345864),
!         (-0.98508200000000001, 0.14734900000000001, -0.088893),
!         (-0.93921500000000002, 0.29425600000000002, -0.17688799999999999),
!         (-0.90468400000000004, 0.41692099999999999, -0.087887000000000007),
!         (-0.42641200000000001, 0.080588999999999994, -0.90093199999999996),
!         (-0.154942, 0.088985999999999996, -0.983908),
!         (-0.30463200000000001, 0.16971600000000001, -0.93722799999999995),
!         (-0.30952400000000002, 0.33910499999999999, -0.88837100000000002),
!         (0.142708, 0.081582000000000002, -0.98639699999999997),
!         (0.41948000000000002, 0.080682000000000004, -0.90417199999999998),
!         (0.29280600000000001, 0.162546, -0.94225400000000004),
!         (0.30551099999999998, 0.33237699999999998, -0.89229400000000003),
!         (0.0, 0.17825299999999999, -0.983985),
!         (0.159858, 0.43690899999999999, -0.88518699999999995),
!         (0.0, 0.35682199999999997, -0.934172),
!         (-0.159858, 0.43690899999999999, -0.88518699999999995),
!         (-0.15881700000000001, 0.58427899999999999, -0.79586100000000004),
!         (0.15881700000000001, 0.58427899999999999, -0.79586100000000004),
!         (0.0, 0.65300400000000003, -0.757355),
!         (0.0, 0.76248000000000005, -0.64701200000000003),
!         (-0.55680700000000005, 0.13039500000000001, -0.82034300000000004),
!         (-0.70362499999999994, 0.298564, -0.64480300000000002),
!         (-0.57923899999999995, 0.27460699999999999, -0.76751199999999997),
!         (-0.45350600000000002, 0.39410099999999998, -0.79938500000000001),
!         (-0.80132000000000003, 0.44216899999999998, -0.40295599999999998),
!         (-0.82077599999999995, 0.55441700000000005, -0.137651),
!         (-0.771509, 0.56907799999999997, -0.28447099999999997),
!         (-0.65010100000000004, 0.69738800000000001, -0.30169400000000002),
!         (-0.70693399999999995, 0.44827800000000001, -0.54707600000000001),
!         (-0.54707600000000001, 0.70693399999999995, -0.44827800000000001),
!         (-0.57735000000000003, 0.57735000000000003, -0.57735000000000003),
!         (-0.44827800000000001, 0.54707600000000001, -0.70693399999999995),
!         (-0.30264799999999997, 0.63921799999999995, -0.70696800000000004),
!         (-0.40080300000000002, 0.79803500000000005, -0.44999699999999998),
!         (-0.28621099999999999, 0.76232699999999998, -0.58046699999999996),
!         (-0.142204, 0.816797, -0.55912499999999998),
!         (-0.82034300000000004, -0.55680700000000005, 0.13039500000000001),
!         (-0.64480300000000002, -0.70362499999999994, 0.298564),
!         (-0.76751199999999997, -0.57923899999999995, 0.27460699999999999),
!         (-0.79938500000000001, -0.45350600000000002, 0.39410099999999998),
!         (-0.40295599999999998, -0.80132000000000003, 0.44216899999999998),
!         (-0.137651, -0.82077599999999995, 0.55441700000000005),
!         (-0.28447099999999997, -0.771509, 0.56907799999999997),
!         (-0.30169400000000002, -0.65010100000000004, 0.69738800000000001),
!         (-0.54707600000000001, -0.70693399999999995, 0.44827800000000001),
!         (-0.44827800000000001, -0.54707600000000001, 0.70693399999999995),
!         (-0.57735000000000003, -0.57735000000000003, 0.57735000000000003),
!         (-0.70693399999999995, -0.44827800000000001, 0.54707600000000001),
!         (-0.70696800000000004, -0.30264799999999997, 0.63921799999999995),
!         (-0.44999699999999998, -0.40080300000000002, 0.79803500000000005),
!         (-0.58046699999999996, -0.28621099999999999, 0.76232699999999998),
!         (-0.55912499999999998, -0.142204, 0.816797),
!         (-0.90093199999999996, -0.42641200000000001, 0.080588999999999994),
!         (-0.983908, -0.154942, 0.088985999999999996),
!         (-0.93722799999999995, -0.30463200000000001, 0.16971600000000001),
!         (-0.88837100000000002, -0.30952400000000002, 0.33910499999999999),
!         (-0.98639699999999997, 0.142708, 0.081582000000000002),
!         (-0.90417199999999998, 0.41948000000000002, 0.080682000000000004),
!         (-0.94225400000000004, 0.29280600000000001, 0.162546),
!         (-0.89229400000000003, 0.30551099999999998, 0.33237699999999998),
!         (-0.983985, 0.0, 0.17825299999999999),
!         (-0.88518699999999995, 0.159858, 0.43690899999999999),
!         (-0.934172, 0.0, 0.35682199999999997),
!         (-0.88518699999999995, -0.159858, 0.43690899999999999),
!         (-0.79586100000000004, -0.15881700000000001, 0.58427899999999999),
!         (-0.79586100000000004, 0.15881700000000001, 0.58427899999999999),
!         (-0.757355, 0.0, 0.65300400000000003),
!         (-0.64701200000000003, 0.0, 0.76248000000000005),
!         (0.0, -0.77053700000000003, 0.63739500000000004),
!         (-0.154582, -0.58980699999999997, 0.79261099999999995),
!         (0.0, -0.66262100000000002, 0.74895500000000004),
!         (0.154582, -0.58980699999999997, 0.79261099999999995),
!         (-0.310166, -0.35253499999999999, 0.88290199999999996),
!         (-0.42387200000000003, -0.087787000000000004, 0.90145799999999998),
!         (-0.30607400000000001, -0.18401300000000001, 0.93405499999999997),
!         (-0.15959100000000001, -0.096272999999999997, 0.98247799999999996),
!         (-0.159858, -0.43690899999999999, 0.88518699999999995),
!         (0.0, -0.17825299999999999, 0.983985),
!         (0.0, -0.35682199999999997, 0.934172),
!         (0.159858, -0.43690899999999999, 0.88518699999999995),
!         (0.30616599999999999, -0.345864, 0.88692800000000005),
!         (0.14734900000000001, -0.088893, 0.98508200000000001),
!         (0.29425600000000002, -0.17688799999999999, 0.93921500000000002),
!         (0.41692099999999999, -0.087887000000000007, 0.90468400000000004),
!         (0.13039500000000001, -0.82034300000000004, 0.55680700000000005),
!         (0.39410099999999998, -0.79938500000000001, 0.45350600000000002),
!         (0.27460699999999999, -0.76751199999999997, 0.57923899999999995),
!         (0.298564, -0.64480300000000002, 0.70362499999999994),
!         (0.63386200000000004, -0.71312200000000003, 0.29946099999999998),
!         (0.81638500000000003, -0.56152299999999999, 0.134937),
!         (0.75824199999999997, -0.59054499999999999, 0.27627200000000002),
!         (0.79602099999999998, -0.461285, 0.39187699999999998),
!         (0.54707600000000001, -0.70693399999999995, 0.44827800000000001),
!         (0.70693399999999995, -0.44827800000000001, 0.54707600000000001),
!         (0.57735000000000003, -0.57735000000000003, 0.57735000000000003),
!         (0.44827800000000001, -0.54707600000000001, 0.70693399999999995),
!         (0.44999699999999998, -0.40080300000000002, 0.79803500000000005),
!         (0.70696800000000004, -0.30264799999999997, 0.63921799999999995),
!         (0.58046699999999996, -0.28621099999999999, 0.76232699999999998),
!         (0.55912499999999998, -0.142204, 0.816797),
!         (0.55680700000000005, 0.13039500000000001, -0.82034300000000004),
!         (0.45350600000000002, 0.39410099999999998, -0.79938500000000001),
!         (0.57923899999999995, 0.27460699999999999, -0.76751199999999997),
!         (0.70362499999999994, 0.298564, -0.64480300000000002),
!         (0.29946099999999998, 0.63386200000000004, -0.71312200000000003),
!         (0.134937, 0.81638500000000003, -0.56152299999999999),
!         (0.27627200000000002, 0.75824199999999997, -0.59054499999999999),
!         (0.39187699999999998, 0.79602099999999998, -0.461285),
!         (0.44827800000000001, 0.54707600000000001, -0.70693399999999995),
!         (0.54707600000000001, 0.70693399999999995, -0.44827800000000001),
!         (0.57735000000000003, 0.57735000000000003, -0.57735000000000003),
!         (0.70693399999999995, 0.44827800000000001, -0.54707600000000001),
!         (0.79803500000000005, 0.44999699999999998, -0.40080300000000002),
!         (0.63921799999999995, 0.70696800000000004, -0.30264799999999997),
!         (0.76232699999999998, 0.58046699999999996, -0.28621099999999999),
!         (0.816797, 0.55912499999999998, -0.142204),
!         (0.63739500000000004, 0.0, -0.77053700000000003),
!         (0.79261099999999995, -0.154582, -0.58980699999999997),
!         (0.74895500000000004, 0.0, -0.66262100000000002),
!         (0.79261099999999995, 0.154582, -0.58980699999999997),
!         (0.88290199999999996, -0.310166, -0.35253499999999999),
!         (0.90145799999999998, -0.42387200000000003, -0.087787000000000004),
!         (0.93405499999999997, -0.30607400000000001, -0.18401300000000001),
!         (0.98247799999999996, -0.15959100000000001, -0.096272999999999997),
!         (0.88518699999999995, -0.159858, -0.43690899999999999),
!         (0.983985, 0.0, -0.17825299999999999),
!         (0.934172, 0.0, -0.35682199999999997),
!         (0.88518699999999995, 0.159858, -0.43690899999999999),
!         (0.88692800000000005, 0.30616599999999999, -0.345864),
!         (0.98508200000000001, 0.14734900000000001, -0.088893),
!         (0.93921500000000002, 0.29425600000000002, -0.17688799999999999),
!         (0.90468400000000004, 0.41692099999999999, -0.087887000000000007),
!         (-0.77053700000000003, 0.63739500000000004, 0.0),
!         (-0.58980699999999997, 0.79261099999999995, -0.154582),
!         (-0.66262100000000002, 0.74895500000000004, 0.0),
!         (-0.58980699999999997, 0.79261099999999995, 0.154582),
!         (-0.35253499999999999, 0.88290199999999996, -0.310166),
!         (-0.087787000000000004, 0.90145799999999998, -0.42387200000000003),
!         (-0.18401300000000001, 0.93405499999999997, -0.30607400000000001),
!         (-0.096272999999999997, 0.98247799999999996, -0.15959100000000001),
!         (-0.43690899999999999, 0.88518699999999995, -0.159858),
!         (-0.17825299999999999, 0.983985, 0.0),
!         (-0.35682199999999997, 0.934172, 0.0),
!         (-0.43690899999999999, 0.88518699999999995, 0.159858),
!         (-0.345864, 0.88692800000000005, 0.30616599999999999),
!         (-0.088893, 0.98508200000000001, 0.14734900000000001),
!         (-0.17688799999999999, 0.93921500000000002, 0.29425600000000002),
!         (-0.087887000000000007, 0.90468400000000004, 0.41692099999999999),
!         (-0.82034300000000004, 0.55680700000000005, 0.13039500000000001),
!         (-0.79938500000000001, 0.45350600000000002, 0.39410099999999998),
!         (-0.76751199999999997, 0.57923899999999995, 0.27460699999999999),
!         (-0.64480300000000002, 0.70362499999999994, 0.298564),
!         (-0.71312200000000003, 0.29946099999999998, 0.63386200000000004),
!         (-0.56152299999999999, 0.134937, 0.81638500000000003),
!         (-0.59054499999999999, 0.27627200000000002, 0.75824199999999997),
!         (-0.461285, 0.39187699999999998, 0.79602099999999998),
!         (-0.70693399999999995, 0.44827800000000001, 0.54707600000000001),
!         (-0.44827800000000001, 0.54707600000000001, 0.70693399999999995),
!         (-0.57735000000000003, 0.57735000000000003, 0.57735000000000003),
!         (-0.54707600000000001, 0.70693399999999995, 0.44827800000000001),
!         (-0.40080300000000002, 0.79803500000000005, 0.44999699999999998),
!         (-0.30264799999999997, 0.63921799999999995, 0.70696800000000004),
!         (-0.28621099999999999, 0.76232699999999998, 0.58046699999999996),
!         (-0.142204, 0.816797, 0.55912499999999998),
!         (0.080588999999999994, 0.90093199999999996, -0.42641200000000001),
!         (0.088985999999999996, 0.983908, -0.154942),
!         (0.16971600000000001, 0.93722799999999995, -0.30463200000000001),
!         (0.33910499999999999, 0.88837100000000002, -0.30952400000000002),
!         (0.081582000000000002, 0.98639699999999997, 0.142708),
!         (0.080682000000000004, 0.90417199999999998, 0.41948000000000002),
!         (0.162546, 0.94225400000000004, 0.29280600000000001),
!         (0.33237699999999998, 0.89229400000000003, 0.30551099999999998),
!         (0.17825299999999999, 0.983985, 0.0),
!         (0.43690899999999999, 0.88518699999999995, 0.159858),
!         (0.35682199999999997, 0.934172, 0.0),
!         (0.43690899999999999, 0.88518699999999995, -0.159858),
!         (0.58427899999999999, 0.79586100000000004, -0.15881700000000001),
!         (0.58427899999999999, 0.79586100000000004, 0.15881700000000001),
!         (0.65300400000000003, 0.757355, 0.0),
!         (0.76248000000000005, 0.64701200000000003, 0.0),
!         (-0.42641200000000001, 0.080588999999999994, 0.90093199999999996),
!         (-0.30952400000000002, 0.33910499999999999, 0.88837100000000002),
!         (-0.30463200000000001, 0.16971600000000001, 0.93722799999999995),
!         (-0.154942, 0.088985999999999996, 0.983908),
!         (-0.167458, 0.58344099999999999, 0.79470399999999997),
!         (-0.0043920000000000001, 0.76652100000000001, 0.64220500000000003),
!         (-0.013266999999999999, 0.65778300000000001, 0.75309099999999995),
!         (0.14591999999999999, 0.59060100000000004, 0.79366099999999995),
!         (-0.159858, 0.43690899999999999, 0.88518699999999995),
!         (0.159858, 0.43690899999999999, 0.88518699999999995),
!         (0.0, 0.35682199999999997, 0.934172),
!         (0.0, 0.17825299999999999, 0.983985),
!         (0.14734900000000001, 0.088893, 0.98508200000000001),
!         (0.30616599999999999, 0.345864, 0.88692800000000005),
!         (0.29425600000000002, 0.17688799999999999, 0.93921500000000002),
!         (0.41692099999999999, 0.087887000000000007, 0.90468400000000004),
!         (0.90093199999999996, -0.42641200000000001, 0.080588999999999994),
!         (0.88837100000000002, -0.30952400000000002, 0.33910499999999999),
!         (0.93722799999999995, -0.30463200000000001, 0.16971600000000001),
!         (0.983908, -0.154942, 0.088985999999999996),
!         (0.79470399999999997, -0.167458, 0.58344099999999999),
!         (0.64220500000000003, -0.0043920000000000001, 0.76652100000000001),
!         (0.75309099999999995, -0.013266999999999999, 0.65778300000000001),
!         (0.79366099999999995, 0.14591999999999999, 0.59060100000000004),
!         (0.88518699999999995, -0.159858, 0.43690899999999999),
!         (0.88518699999999995, 0.159858, 0.43690899999999999),
!         (0.934172, 0.0, 0.35682199999999997),
!         (0.983985, 0.0, 0.17825299999999999),
!         (0.98508200000000001, 0.14734900000000001, 0.088893),
!         (0.88692800000000005, 0.30616599999999999, 0.345864),
!         (0.93921500000000002, 0.29425600000000002, 0.17688799999999999),
!         (0.90468400000000004, 0.41692099999999999, 0.087887000000000007),
!         (0.13039500000000001, 0.82034300000000004, 0.55680700000000005),
!         (0.298564, 0.64480300000000002, 0.70362499999999994),
!         (0.27460699999999999, 0.76751199999999997, 0.57923899999999995),
!         (0.39410099999999998, 0.79938500000000001, 0.45350600000000002),
!         (0.44216899999999998, 0.40295599999999998, 0.80132000000000003),
!         (0.55441700000000005, 0.137651, 0.82077599999999995),
!         (0.56907799999999997, 0.28447099999999997, 0.771509),
!         (0.69738800000000001, 0.30169400000000002, 0.65010100000000004),
!         (0.44827800000000001, 0.54707600000000001, 0.70693399999999995),
!         (0.70693399999999995, 0.44827800000000001, 0.54707600000000001),
!         (0.57735000000000003, 0.57735000000000003, 0.57735000000000003),
!         (0.54707600000000001, 0.70693399999999995, 0.44827800000000001),
!         (0.63921799999999995, 0.70696800000000004, 0.30264799999999997),
!         (0.79803500000000005, 0.44999699999999998, 0.40080300000000002),
!         (0.76232699999999998, 0.58046699999999996, 0.28621099999999999),
!         (0.816797, 0.55912499999999998, 0.142204)]
!   ico = icosahedron(level=2)
!   assert ico.level == 2
!   diff = ico.sites - flex.vec3_double(expected_sites_level_2)
!   for d in diff.min(): assert abs(d) < 0.02
!   for d in diff.max(): assert abs(d) < 0.02
!   for level in xrange(1,6):
      ico = icosahedron(level=level)
      assert ico.level == level
!     assert ico.sites.size() == 80 * 4**(level-1)
  
  def run():
--- 941,961 ----
  
  def exercise_icosahedron():
!   ico = icosahedron(level=0)
!   for level in xrange(6):
      ico = icosahedron(level=level)
      assert ico.level == level
!     if (level == 0):
!       assert ico.sites.size() == 12
!     else:
!       assert ico.sites.size() == 80 * 4**(level-1)
!     assert approx_equal(ico.sites.mean(), [0,0,0])
!     assert approx_equal(ico.sites.dot(), [1]*ico.sites.size())
!     d = ico.next_neighbors_distance()
!     m = flex.min((ico.sites[1:] - ico.sites[0]).dot())**0.5
!     if (level == 0):
!       assert approx_equal(d, m)
!     else:
!       assert d > m
!       assert d/2 < m
  
  def run():

Update of /cvsroot/cctbx/scitbx/include/scitbx/math
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7015/include/scitbx/math

Modified Files:
	icosahedron.h 
Log Message:
cctbx.xray.structure:
  - show_pairs -> cctbx.crystal.show_distances
  - new show_distances() method

cctbx.crystal.pair_asu_table:
  - new show_distances() method
  - new pair_counts() method

cctbx.crystal.neighbors_simple_pair_generator:
cctbx.crystal.neighbors_fast_pair_generator:
  - new max_distance_sq() method

scitbx.math.icosahedron:
  - new next_neighbors_distance() method
  - slight modification of make_icosahedron() algorithm:
      mid points are now immediately normalized; this leads to a
      slightly more uniform distribution of next-neighbor distances


Index: icosahedron.h
===================================================================
RCS file: /cvsroot/cctbx/scitbx/include/scitbx/math/icosahedron.h,v
retrieving revision 1.2
retrieving revision 1.3
diff -C2 -d -r1.2 -r1.3
*** icosahedron.h	3 Feb 2005 12:51:56 -0000	1.2
--- icosahedron.h	4 Feb 2005 02:57:40 -0000	1.3
***************
*** 3,125 ****
  
  #include <scitbx/array_family/shared.h>
! #include <math.h>
  
  namespace scitbx { namespace math {
-   template <typename FloatType = double>
-     class icosahedron
-     {
-       public:
-       int level;
-       af::shared<scitbx::vec3<FloatType> > sites;
-       private:
-       af::shared<scitbx::vec3<FloatType> > v;
-       FloatType distsq(scitbx::vec3<FloatType> s1, scitbx::vec3<FloatType> s2) {
-         return (s1-s2).length();
-       }
- 
-       scitbx::vec3<FloatType> midpt(scitbx::vec3<FloatType> s1, scitbx::vec3<FloatType> s2) {
-         return (s1+s2)/2.0;
-       }
  
!       void make_icosahedron() {
!         scitbx::vec3<FloatType> temp;
!         static const FloatType mag_val = atan((1.0+sqrt(5.0))/2.0);
!         static const FloatType s = sin(mag_val);
!         static const FloatType c = cos(mag_val);
!         FloatType a = s;
!         FloatType b = c;
  
!         for (int i = 0; i < 2; i++) {
            a = -a;
!           for (int j = 0; j < 2; j++) {
              b = -b;
! 
!             temp[0] = 0.0;
!             temp[1] = a;
!             temp[2] = b;
!             v.push_back(temp);
! 
!             temp[0] = b;
!             temp[1] = 0.0;
!             temp[2] = a;
!             v.push_back(temp);
! 
!             temp[0] = a;
!             temp[1] = b;
!             temp[2] = 0.0;
!             v.push_back(temp);
            }
          }
        }
  
!       /*-
!        * Recursively subdivides triangle ABC levl times, storing the result in
!        * the probe array, PROBE->uco.
!        *
!        * STRATEGY: if not the bottom level, find the midpoints of each of the
!        * sides of this sub_triangle and recurse, calling each of the four sub
!        * -sub_triangles in turn.  If it's the bottom level, find the center, and
!        * add that point to the output probe object.
!        -*/
!       void sub_triangle(scitbx::vec3<FloatType> a, scitbx::vec3<FloatType> b, scitbx::vec3<FloatType> c, int lvl) {
!         scitbx::vec3<FloatType> d, e, f;
  
!         if (lvl > 0) {
!           /* find new midpoints */
!           d = midpt(a, b);
!           e = midpt(b, c);
!           f = midpt(a, c);
  
!           sub_triangle(a, d, f, lvl - 1);
!           sub_triangle(d, b, e, lvl - 1);
!           sub_triangle(d, e, f, lvl - 1);
!           sub_triangle(f, e, c, lvl - 1);
!         } else {
!           /* normalize the vectors */
!           d = a.normalize();
!           e = b.normalize();
!           f = c.normalize();
!           /* find center */
!           af::shared<scitbx::vec3<FloatType> > tmp;
!           tmp.push_back(d);
!           tmp.push_back(e);
!           tmp.push_back(f);
!           scitbx::vec3<FloatType> g;
!           for(int i = 0; i < 3; i++) {
!             g[i] = (d[i]+e[i]+f[i])/3.0;
            }
!           sites.push_back(g.normalize());
          }
        }
  
!       public:
!       icosahedron() { }
! 
!       /*-
!         builds a probe sphere with uniformly distributed points, by recursively
!         sub-dividing an icosahedrom level times.
!         -*/
!       icosahedron(int level_) {
!         int i, j, k;
  
!         SCITBX_ASSERT(level_>=1);
!         level=level_;
!         sites.clear();
  
!         make_icosahedron();
!         for (i = 0; i < 10; i++) {
!           for (j = i + 1; j < 11; j++) {
!             if (distsq(v[i], v[j]) < 1.21) {
!               for (k = j + 1; k < 12; k++) {
!                 if ((distsq(v[i], v[k]) < 1.21) &&
!                     (distsq(v[j], v[k]) < 1.21))
!                   sub_triangle(v[i], v[j], v[k],level);
!               } /* k */
!             } /* endif */
!           } /* j */
!         } /* i */
        }
!     };
! }}
  
  #endif // SCITBX_MATH_ICOSAHEDRON_H
--- 3,143 ----
  
  #include <scitbx/array_family/shared.h>
! #include <cmath>
  
  namespace scitbx { namespace math {
  
!   //! Cartesian coordinates of points on a sphere.
!   /*! This code was contributed by Erik McKee, Texas A&M University.
!    */
!   template <typename FloatType=double>
!   class icosahedron
!   {
!     public:
!       unsigned level;
!       af::shared<scitbx::vec3<FloatType> > sites;
  
!     private:
!       static
!       void
!       make_icosahedron(scitbx::vec3<FloatType>* v)
!       {
!         static const FloatType phi_ratio = (1.0+std::sqrt(5.0))/2.0;
!         static const FloatType b_ = 1.0 / (std::sqrt(1.0+phi_ratio*phi_ratio));
!         static const FloatType a_ = phi_ratio * b_;
!         FloatType a = a_;
!         FloatType b = b_;
!         for (unsigned i = 0; i < 2; i++) {
            a = -a;
!           for (unsigned j = 0; j < 2; j++) {
              b = -b;
!             (*v)[0] = 0.0;
!             (*v)[1] = a;
!             (*v)[2] = b;
!             v++;
!             (*v)[0] = b;
!             (*v)[1] = 0.0;
!             (*v)[2] = a;
!             v++;
!             (*v)[0] = a;
!             (*v)[1] = b;
!             (*v)[2] = 0.0;
!             v++;
            }
          }
        }
  
!       void
!       sub_triangle(
!         scitbx::vec3<FloatType> const& a,
!         scitbx::vec3<FloatType> const& b,
!         scitbx::vec3<FloatType> const& c,
!         unsigned lvl)
!       {
!         if (lvl == 0) {
!           sites.push_back((a+b+c).normalize());
!         }
!         else {
!           /* compute midpoints */
!           scitbx::vec3<FloatType> d = (a+b).normalize();
!           scitbx::vec3<FloatType> e = (b+c).normalize();
!           scitbx::vec3<FloatType> f = (a+c).normalize();
!           sub_triangle(a, d, f, lvl-1);
!           sub_triangle(d, b, e, lvl-1);
!           sub_triangle(d, e, f, lvl-1);
!           sub_triangle(f, e, c, lvl-1);
!         }
!       }
  
!     public:
!       //! Default constructor. Some data members are not initialized!
!       icosahedron() {}
  
!       /*! Builds points on a sphere.
!        */
!       /*! level=0 builds the vertices of the base icosahedron. The
!           sites are uniformly distributed.
! 
!           level>0 recursively sub-divides the 20 triangles of the base
!           icosahedron level times. The number of sites is therefore
!           80 * 4**(level-1). The sites are approximately uniformly
!           distributed.
!        */
!       icosahedron(unsigned level_)
!       :
!         level(level_)
!       {
!         if (level == 0) {
!           sites.resize(12);
!           make_icosahedron(sites.begin());
!         }
!         else {
!           af::tiny<scitbx::vec3<FloatType>, 12> v;
!           make_icosahedron(v.begin());
!           std::size_t four_pow_level_minus_one = 1;
!           for (unsigned i=0;i<level-1;i++) four_pow_level_minus_one *= 4;
!           sites.reserve(80 * four_pow_level_minus_one);
!           for (unsigned i = 0; i < 10; i++) {
!             for (unsigned j = i + 1; j < 11; j++) {
!               if ((v[i]-v[j]).length_sq() < 1.2) {
!                 for (unsigned k = j + 1; k < 12; k++) {
!                   if (   (v[i]-v[k]).length_sq() < 1.2
!                       && (v[j]-v[k]).length_sq() < 1.2)
!                     sub_triangle(v[i], v[j], v[k], level);
!                 }
!               }
!             }
            }
!           SCITBX_ASSERT(sites.size() == 80 * four_pow_level_minus_one);
          }
        }
  
!       //! Maximum distance to the next neighbors.
!       /*! For level=0 this is the distance to the five next neighbors
!           of each icosahedron vertex.
  
!           For level>0 this is the maximum distance of any site
!           to the three next neighbors.
  
!           Distances are tabulated for levels 0 through 7.
!           An exception is thrown if level>7.
!        */
!       double
!       next_neighbors_distance() const
!       {
!         static const af::tiny<double, 8> known_distances(
!           1.05146222424,
!           0.353098248494,
!           0.185386249656,
!           0.0947464326266,
!           0.0476510500603,
!           0.0238609877705,
!           0.011934950279,
!           0.00596803292972);
!         if (level < known_distances.size()) return known_distances[level];
!         throw std::runtime_error("next_neighbors_distance not known.");
        }
!   };
! 
! }} // namespace scitbx::math
  
  #endif // SCITBX_MATH_ICOSAHEDRON_H

Update of /cvsroot/cctbx/scitbx/math/boost_python
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv1320

Modified Files:
	icosahedron.cpp tst_math.py 
Log Message:
keyword argument support; unneeded includes removed

Index: icosahedron.cpp
===================================================================
RCS file: /cvsroot/cctbx/scitbx/math/boost_python/icosahedron.cpp,v
retrieving revision 1.1
retrieving revision 1.2
diff -C2 -d -r1.1 -r1.2
*** icosahedron.cpp	3 Feb 2005 09:12:35 -0000	1.1
--- icosahedron.cpp	3 Feb 2005 14:00:05 -0000	1.2
***************
*** 1,10 ****
  #include <scitbx/array_family/boost_python/flex_fwd.h>
  #include <boost/python/class.hpp>
  #include <boost/python/return_value_policy.hpp>
  #include <boost/python/return_by_value.hpp>
- #include <boost/python/copy_const_reference.hpp>
- #include <boost/python/return_internal_reference.hpp>
- #include <boost/python/def.hpp>
- 
  
  #include <scitbx/math/icosahedron.h>
--- 1,7 ----
  #include <scitbx/array_family/boost_python/flex_fwd.h>
  #include <boost/python/class.hpp>
+ #include <boost/python/args.hpp>
  #include <boost/python/return_value_policy.hpp>
  #include <boost/python/return_by_value.hpp>
  
  #include <scitbx/math/icosahedron.h>
***************
*** 16,20 ****
    struct icosahedron_wrappers
    {
!     typedef icosahedron<double> w_t;
  
      static void
--- 13,17 ----
    struct icosahedron_wrappers
    {
!     typedef icosahedron<> w_t;
  
      static void
***************
*** 23,32 ****
        using namespace boost::python;
        typedef return_value_policy<return_by_value> rbv;
-       typedef return_value_policy<copy_const_reference> ccr;
-       typedef return_internal_reference<> rir;
        class_<w_t>("icosahedron", no_init)
!         .def(init<int>())
!         .add_property("level",make_getter(&w_t::level))
!         .add_property("sites",make_getter(&w_t::sites,rbv()))
        ;
      }
--- 20,27 ----
        using namespace boost::python;
        typedef return_value_policy<return_by_value> rbv;
        class_<w_t>("icosahedron", no_init)
!         .def(init<int>((arg_("level"))))
!         .def_readonly("level", &w_t::level)
!         .add_property("sites", make_getter(&w_t::sites, rbv()))
        ;
      }

Index: tst_math.py
===================================================================
RCS file: /cvsroot/cctbx/scitbx/math/boost_python/tst_math.py,v
retrieving revision 1.37
retrieving revision 1.38
diff -C2 -d -r1.37 -r1.38
*** tst_math.py	3 Feb 2005 12:54:15 -0000	1.37
--- tst_math.py	3 Feb 2005 14:00:05 -0000	1.38
***************
*** 1262,1266 ****
          (0.76232699999999998, 0.58046699999999996, 0.28621099999999999),
          (0.816797, 0.55912499999999998, 0.142204)]
!   ico = icosahedron(2)
    assert ico.level == 2
    diff = ico.sites - flex.vec3_double(expected_sites_level_2)
--- 1262,1266 ----
          (0.76232699999999998, 0.58046699999999996, 0.28621099999999999),
          (0.816797, 0.55912499999999998, 0.142204)]
!   ico = icosahedron(level=2)
    assert ico.level == 2
    diff = ico.sites - flex.vec3_double(expected_sites_level_2)
***************
*** 1268,1272 ****
    for d in diff.max(): assert abs(d) < 0.02
    for level in xrange(1,6):
!     ico = icosahedron(level)
      assert ico.level == level
      assert ico.sites.size() == 80 * 4**(level-1)
--- 1268,1272 ----
    for d in diff.max(): assert abs(d) < 0.02
    for level in xrange(1,6):
!     ico = icosahedron(level=level)
      assert ico.level == level
      assert ico.sites.size() == 80 * 4**(level-1)

Update of /cvsroot/cctbx/scitbx/math/boost_python
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv21605

Modified Files:
	tst_math.py 
Log Message:
test icosahedron with levels 1 trhough 5

Index: tst_math.py
===================================================================
RCS file: /cvsroot/cctbx/scitbx/math/boost_python/tst_math.py,v
retrieving revision 1.36
retrieving revision 1.37
diff -C2 -d -r1.36 -r1.37
*** tst_math.py	3 Feb 2005 09:12:35 -0000	1.36
--- tst_math.py	3 Feb 2005 12:54:15 -0000	1.37
***************
*** 941,945 ****
  
  def exercise_icosahedron():
!   tmp=[
          (-0.13039500000000001, -0.82034300000000004, -0.55680700000000005),
          (-0.298564, -0.64480300000000002, -0.70362499999999994),
--- 941,945 ----
  
  def exercise_icosahedron():
!   expected_sites_level_2=[
          (-0.13039500000000001, -0.82034300000000004, -0.55680700000000005),
          (-0.298564, -0.64480300000000002, -0.70362499999999994),
***************
*** 1262,1278 ****
          (0.76232699999999998, 0.58046699999999996, 0.28621099999999999),
          (0.816797, 0.55912499999999998, 0.142204)]
- 
    ico = icosahedron(2)
    assert ico.level == 2
!   nvals = ico.sites
!   ovals = flex.vec3_double()
!   for site in tmp:
!     ovals.append(site)
!   diff = nvals-ovals
!   dmin = diff.min()
!   dmax = diff.max()
!   for i in xrange(3):
!     assert 0.0<=math.fabs(dmin[i])<=0.02
!     assert 0.0<=math.fabs(dmax[i])<=0.02
  
  def run():
--- 1262,1274 ----
          (0.76232699999999998, 0.58046699999999996, 0.28621099999999999),
          (0.816797, 0.55912499999999998, 0.142204)]
    ico = icosahedron(2)
    assert ico.level == 2
!   diff = ico.sites - flex.vec3_double(expected_sites_level_2)
!   for d in diff.min(): assert abs(d) < 0.02
!   for d in diff.max(): assert abs(d) < 0.02
!   for level in xrange(1,6):
!     ico = icosahedron(level)
!     assert ico.level == level
!     assert ico.sites.size() == 80 * 4**(level-1)
  
  def run():

Update of /cvsroot/cctbx/scitbx/include/scitbx/math
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv21197/math

Modified Files:
	floating_point_epsilon.h golay.h icosahedron.h 
	principal_axes_of_inertia.h 
Log Message:
missing include guards added

Index: icosahedron.h
===================================================================
RCS file: /cvsroot/cctbx/scitbx/include/scitbx/math/icosahedron.h,v
retrieving revision 1.1
retrieving revision 1.2
diff -C2 -d -r1.1 -r1.2
*** icosahedron.h	3 Feb 2005 09:14:00 -0000	1.1
--- icosahedron.h	3 Feb 2005 12:51:56 -0000	1.2
***************
*** 1,2 ****
--- 1,5 ----
+ #ifndef SCITBX_MATH_ICOSAHEDRON_H
+ #define SCITBX_MATH_ICOSAHEDRON_H
+ 
  #include <scitbx/array_family/shared.h>
  #include <math.h>
***************
*** 119,120 ****
--- 122,125 ----
      };
  }}
+ 
+ #endif // SCITBX_MATH_ICOSAHEDRON_H

Index: golay.h
===================================================================
RCS file: /cvsroot/cctbx/scitbx/include/scitbx/math/golay.h,v
retrieving revision 1.1
retrieving revision 1.2
diff -C2 -d -r1.1 -r1.2
*** golay.h	23 Dec 2003 21:19:40 -0000	1.1
--- golay.h	3 Feb 2005 12:51:56 -0000	1.2
***************
*** 1,2 ****
--- 1,5 ----
+ #ifndef SCITBX_MATH_GOLAY_H
+ #define SCITBX_MATH_GOLAY_H
+ 
  #include <scitbx/array_family/loops.h>
  #include <scitbx/array_family/tiny.h>
***************
*** 78,79 ****
--- 81,84 ----
  
  }} // namespace scitbx::math
+ 
+ #endif // SCITBX_MATH_GOLAY_H

Index: floating_point_epsilon.h
===================================================================
RCS file: /cvsroot/cctbx/scitbx/include/scitbx/math/floating_point_epsilon.h,v
retrieving revision 1.1
retrieving revision 1.2
diff -C2 -d -r1.1 -r1.2
*** floating_point_epsilon.h	2 Jan 2004 21:19:11 -0000	1.1
--- floating_point_epsilon.h	3 Feb 2005 12:51:56 -0000	1.2
***************
*** 1,2 ****
--- 1,5 ----
+ #ifndef SCITBX_MATH_FLOATING_POINT_EPSILON_H
+ #define SCITBX_MATH_FLOATING_POINT_EPSILON_H
+ 
  #include <scitbx/serialization/base_256.h>
  
***************
*** 117,118 ****
--- 120,123 ----
  
  }} // namespace scitbx::math
+ 
+ #endif // SCITBX_MATH_FLOATING_POINT_EPSILON_H

Index: principal_axes_of_inertia.h
===================================================================
RCS file: /cvsroot/cctbx/scitbx/include/scitbx/math/principal_axes_of_inertia.h,v
retrieving revision 1.2
retrieving revision 1.3
diff -C2 -d -r1.2 -r1.3
*** principal_axes_of_inertia.h	11 Apr 2004 06:49:55 -0000	1.2
--- principal_axes_of_inertia.h	3 Feb 2005 12:51:56 -0000	1.3
***************
*** 1,2 ****
--- 1,5 ----
+ #ifndef SCITBX_MATH_PRINCIPAL_AXES_OF_INERTIA_H
+ #define SCITBX_MATH_PRINCIPAL_AXES_OF_INERTIA_H
+ 
  #include <scitbx/math/eigensystem.h>
  
***************
*** 99,100 ****
--- 102,105 ----
  
  }} // namespace scitbx::math
+ 
+ #endif // SCITBX_MATH_PRINCIPAL_AXES_OF_INERTIA_H

Update of /cvsroot/cctbx/scitbx/include/scitbx/math
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv15846

Added Files:
	icosahedron.h 
Log Message:
Adding icosahedron, which generates a unit sphere of evenly spaced points

--- NEW FILE: icosahedron.h ---
#include <scitbx/array_family/shared.h>
#include <math.h>

namespace scitbx { namespace math {
  template <typename FloatType = double>
    class icosahedron
    {
      public:
      int level;
      af::shared<scitbx::vec3<FloatType> > sites;
      private:
      af::shared<scitbx::vec3<FloatType> > v;
      FloatType distsq(scitbx::vec3<FloatType> s1, scitbx::vec3<FloatType> s2) {
        return (s1-s2).length();
      }

      scitbx::vec3<FloatType> midpt(scitbx::vec3<FloatType> s1, scitbx::vec3<FloatType> s2) {
        return (s1+s2)/2.0;
      }

      void make_icosahedron() {
        scitbx::vec3<FloatType> temp;
        static const FloatType mag_val = atan((1.0+sqrt(5.0))/2.0);
        static const FloatType s = sin(mag_val);
        static const FloatType c = cos(mag_val);
        FloatType a = s;
        FloatType b = c;

        for (int i = 0; i < 2; i++) {
          a = -a;
          for (int j = 0; j < 2; j++) {
            b = -b;

            temp[0] = 0.0;
            temp[1] = a;
            temp[2] = b;
            v.push_back(temp);

            temp[0] = b;
            temp[1] = 0.0;
            temp[2] = a;
            v.push_back(temp);

            temp[0] = a;
            temp[1] = b;
            temp[2] = 0.0;
            v.push_back(temp);
          }
        }
      }

      /*-
       * Recursively subdivides triangle ABC levl times, storing the result in
       * the probe array, PROBE->uco.
       *
       * STRATEGY: if not the bottom level, find the midpoints of each of the
       * sides of this sub_triangle and recurse, calling each of the four sub
       * -sub_triangles in turn.  If it's the bottom level, find the center, and
       * add that point to the output probe object.
       -*/
      void sub_triangle(scitbx::vec3<FloatType> a, scitbx::vec3<FloatType> b, scitbx::vec3<FloatType> c, int lvl) {
        scitbx::vec3<FloatType> d, e, f;

        if (lvl > 0) {
          /* find new midpoints */
          d = midpt(a, b);
          e = midpt(b, c);
          f = midpt(a, c);

          sub_triangle(a, d, f, lvl - 1);
          sub_triangle(d, b, e, lvl - 1);
          sub_triangle(d, e, f, lvl - 1);
          sub_triangle(f, e, c, lvl - 1);
        } else {
          /* normalize the vectors */
          d = a.normalize();
          e = b.normalize();
          f = c.normalize();
          /* find center */
          af::shared<scitbx::vec3<FloatType> > tmp;
          tmp.push_back(d);
          tmp.push_back(e);
          tmp.push_back(f);
          scitbx::vec3<FloatType> g;
          for(int i = 0; i < 3; i++) {
            g[i] = (d[i]+e[i]+f[i])/3.0;
          }
          sites.push_back(g.normalize());
        }
      }

      public:
      icosahedron() { }

      /*-
        builds a probe sphere with uniformly distributed points, by recursively
        sub-dividing an icosahedrom level times.
        -*/
      icosahedron(int level_) {
        int i, j, k;

        SCITBX_ASSERT(level_>=1);
        level=level_;
        sites.clear();

        make_icosahedron();
        for (i = 0; i < 10; i++) {
          for (j = i + 1; j < 11; j++) {
            if (distsq(v[i], v[j]) < 1.21) {
              for (k = j + 1; k < 12; k++) {
                if ((distsq(v[i], v[k]) < 1.21) &&
                    (distsq(v[j], v[k]) < 1.21))
                  sub_triangle(v[i], v[j], v[k],level);
              } /* k */
            } /* endif */
          } /* j */
        } /* i */
      }
    };
}}