text: from __future__ import absolute_import, division, print_function
  
  import sys
  import random
  
  from cctbx.development import random_structure
  from cctbx.development import debug_utils
  from cctbx import adptbx
  from cctbx.array_family import flex
  
  from iotbx import csv_utils
  
  import libtbx
  import libtbx.utils
  from libtbx.test_utils import approx_equal
  
  import scitbx.random
  from scitbx.random import variate, normal_distribution, gamma_distribution
  from scitbx.math import distributions
  
  from smtbx import absolute_structure
  
  try:
    distributions.students_t_distribution(1)
  except RuntimeError as e:
    # XXX Student's t distribution is not supported with GCC 3.2 builds
    if str(e).startswith("Implementation not available in this build."):
      students_t_available = False
      print("Skipping exercise_hooft_analysis() with Student's t distribution.")
    else:
      raise RuntimeError(e)
  else:
    students_t_available = True
  
  class test_case(object):
  
    d_min=1
    use_students_t_errors=False
    elements = ("N", "C", "C", "S") * 5
  
    def __init__(self, space_group_info, **kwds):
      libtbx.adopt_optional_init_args(self, kwds)
      self.space_group_info = space_group_info
      self.structure = random_structure.xray_structure(
        space_group_info,
        elements=self.elements,
        volume_per_atom=20.,
        min_distance=1.5,
        general_positions_only=True,
        use_u_aniso=False,
        u_iso=adptbx.b_as_u(10),
      )
      self.structure.set_inelastic_form_factors(1.54, "sasaki")
      self.scale_factor = 0.05 + 10 * flex.random_double()
      fc = self.structure.structure_factors(
            anomalous_flag=True, d_min=self.d_min, algorithm="direct").f_calc()
      fo = fc.as_amplitude_array()
      fo.set_observation_type_xray_amplitude()
      if self.use_students_t_errors:
        nu = random.uniform(1, 10)
        normal_g = variate(normal_distribution())
        gamma_g = variate(gamma_distribution(0.5*nu, 2))
        errors = normal_g(fc.size())/flex.sqrt(2*gamma_g(fc.size()))
      else:
        # use gaussian errors
        g = variate(normal_distribution())
        errors = g(fc.size())
      fo2 = fo.as_intensity_array()
      self.fo2 = fo2.customized_copy(
        data=(fo2.data()+errors)*self.scale_factor,
        sigmas=flex.double(fc.size(), 1),
      )
      self.fc = fc
      xs_i = self.structure.inverse_hand()
      self.fc_i = xs_i.structure_factors(
        anomalous_flag=True, d_min=self.d_min, algorithm="direct").f_calc()
      fo2_twin = self.fc.customized_copy(
        data=self.fc.data()+self.fc_i.data()).as_intensity_array()
      self.fo2_twin = fo2_twin.customized_copy(
        data=(errors + fo2_twin.data()) * self.scale_factor,
        sigmas=self.fo2.sigmas())
  
  class hooft_analysis_test_case(test_case):
  
    def exercise(self, debug=False):
      if debug:
        distribution = distributions.normal_distribution()
        observed_deviations = (
          self.fo2.data() - self.scale_factor*self.fc.as_intensity_array().data())
        observed_deviations = observed_deviations.select(
          flex.sort_permutation(observed_deviations))
        expected_deviations = distribution.quantiles(observed_deviations.size())
        csv_utils.writer(
          open('delta_F_npp.csv', 'wb'), (expected_deviations, observed_deviations))
      # first with the correct absolute structure
      gaussian = absolute_structure.hooft_analysis(self.fo2, self.fc)
      analyses = [gaussian]
      NPP = absolute_structure.bijvoet_differences_probability_plot(gaussian)
      if self.use_students_t_errors:
        nu_calc = absolute_structure.maximise_students_t_correlation_coefficient(
          NPP.y, min_nu=1, max_nu=200)
        t_analysis = absolute_structure.students_t_hooft_analysis(
          self.fo2, self.fc, nu_calc, probability_plot_slope=NPP.fit.slope())
        analyses.append(gaussian)
        tPP = absolute_structure.bijvoet_differences_probability_plot(
          t_analysis, use_students_t_distribution=True, students_t_nu=nu_calc)
        if tPP.distribution.degrees_of_freedom() < 100:
          assert tPP.correlation.coefficient() > NPP.correlation.coefficient()
      else:
        assert approx_equal(NPP.correlation.coefficient(), 1, 0.005)
      for analysis in analyses:
        assert approx_equal(analysis.hooft_y, 0, 1e-2)
        assert approx_equal(analysis.p2_true, 1)
        assert approx_equal(analysis.p2_false, 0)
        assert approx_equal(analysis.p3_true, 1)
        assert approx_equal(analysis.p3_false, 0)
        assert approx_equal(analysis.p3_racemic_twin, 0)
      if debug:
        csv_utils.writer(open('npp.csv', 'wb'), (NPP.x,NPP.y))
        if self.use_students_t_errors:
          csv_utils.writer(open('tpp.csv', 'wb'), (tPP.x,tPP.y))
      assert approx_equal(NPP.fit.y_intercept(), 0)
      # and now with the wrong absolute structure
      gaussian = absolute_structure.hooft_analysis(self.fo2, self.fc_i)
      analyses = [gaussian]
      NPP = absolute_structure.bijvoet_differences_probability_plot(gaussian)
      if self.use_students_t_errors:
        nu_calc = absolute_structure.maximise_students_t_correlation_coefficient(
          NPP.y, min_nu=1, max_nu=200)
        t_analysis = absolute_structure.students_t_hooft_analysis(
          self.fo2, self.fc_i, nu_calc, probability_plot_slope=NPP.fit.slope())
        analyses.append(gaussian)
        tPP = absolute_structure.bijvoet_differences_probability_plot(
          t_analysis, use_students_t_distribution=True)
        if tPP.distribution.degrees_of_freedom() < 100:
          assert tPP.correlation.coefficient() > NPP.correlation.coefficient()
      else:
        assert approx_equal(NPP.correlation.coefficient(), 1, 0.002)
        assert approx_equal(NPP.fit.y_intercept(), 0)
      for analysis in analyses:
        assert approx_equal(analysis.hooft_y, 1, 1e-2)
        assert approx_equal(analysis.p2_true, 0)
        assert approx_equal(analysis.p2_false, 1)
        assert approx_equal(analysis.p3_true, 0)
        assert approx_equal(analysis.p3_false, 1)
        assert approx_equal(analysis.p3_racemic_twin, 0)
      # test for the case of a racemic twin
      gaussian = absolute_structure.hooft_analysis(self.fo2_twin, self.fc)
      analyses = [gaussian]
      NPP = absolute_structure.bijvoet_differences_probability_plot(gaussian)
      if self.use_students_t_errors:
        nu_calc = absolute_structure.maximise_students_t_correlation_coefficient(
          NPP.y, min_nu=1, max_nu=200)
        t_analysis = absolute_structure.students_t_hooft_analysis(
          self.fo2_twin, self.fc, nu_calc, probability_plot_slope=NPP.fit.slope())
        tPP = absolute_structure.bijvoet_differences_probability_plot(
          t_analysis, use_students_t_distribution=True)
        if tPP.distribution.degrees_of_freedom() < 100:
          assert tPP.correlation.coefficient() > NPP.correlation.coefficient()
      else:
        assert approx_equal(NPP.correlation.coefficient(), 1, 0.002)
        assert approx_equal(NPP.fit.y_intercept(), 0)
      for analysis in analyses:
        assert approx_equal(analysis.hooft_y, 0.5, 1e-2)
        assert approx_equal(analysis.p3_true, 0)
        assert approx_equal(analysis.p3_false, 0)
        assert approx_equal(analysis.p3_racemic_twin, 1)
  
  
  def run_call_back(flags, space_group_info):
    if not space_group_info.group().is_centric():
      if flags.fix_random_seeds:
        random.seed(1)
        flex.set_random_seed(1)
        scitbx.random.set_random_seed(1)
      hooft_analysis_test_case(
        space_group_info,use_students_t_errors=False).exercise(debug=flags.Debug)
      if students_t_available:
        hooft_analysis_test_case(
          space_group_info,use_students_t_errors=True).exercise(debug=flags.Debug)
  
  def run():
    debug_utils.parse_options_loop_space_groups(
      sys.argv[1:],
      run_call_back,
      keywords=("fix_random_seeds",),
    )
  
  if __name__ == '__main__':
    libtbx.utils.show_times_at_exit()
    run()
  