Python Debug.banner方法代码示例

# 需要导入模块: from xia2.Handlers.Streams import Debug [as 别名]
# 或者: from xia2.Handlers.Streams.Debug import banner [as 别名]
def get_unit_cell_errors(stop_after=None):
  '''Actually process something...'''
  wd = os.getcwd()

  all_miller_indices, all_two_thetas_obs, reference_cell, reference_lattice, reference_wavelength = load_sweeps_with_common_indexing()

  Chatter.banner('Unit cell sampling')
  Debug.banner('Unit cell sampling')
  span = miller.index_span(all_miller_indices)
  Chatter.write("Found %d reflections in 2theta range %.3f - %.3f deg" % (len(all_miller_indices), min(all_two_thetas_obs), max(all_two_thetas_obs)))
  Chatter.write("Miller index range: %s - %s" % (str(span.min()), str(span.max())))
  unit_cell_info = { 'reflections':
                     { 'count': len(all_miller_indices),
                       'min_2theta': min(all_two_thetas_obs),
                       'max_2theta': max(all_two_thetas_obs),
                       'min_miller': list(span.min()),
                       'max_miller': list(span.max())
                     } }

  # Exclude 1% of reflections to remove potential outliers
  # eg. indexed/integrated high angle noise
  two_theta_cutoff = sorted(all_two_thetas_obs)[-int(len(all_two_thetas_obs) * 0.01)-1]
  Chatter.write("Excluding outermost 1%% of reflections (2theta >= %.3f)" % two_theta_cutoff)
  two_thetas_select = all_two_thetas_obs < two_theta_cutoff
  all_two_thetas_obs = all_two_thetas_obs.select(two_thetas_select)
  all_miller_indices = all_miller_indices.select(two_thetas_select)
  Chatter.write("Kept %d reflections in 2theta range %.3f - %.3f deg" % (len(all_miller_indices), min(all_two_thetas_obs), max(all_two_thetas_obs)))
  span = miller.index_span(all_miller_indices)
  unit_cell_info['reflections_filtered'] = \
                     { 'count': len(all_miller_indices),
                       'min_2theta': min(all_two_thetas_obs),
                       'max_2theta': max(all_two_thetas_obs),
                       'min_miller': list(span.min()),
                       'max_miller': list(span.max())
                     }

  # prepare MonteCarlo sampling
  mc_runs = 50
  sample_size = min(len(all_miller_indices) // 2, 100)
  unit_cell_info['sampling'] = { 'method': 'montecarlo', 'runs': mc_runs, 'used_per_run': sample_size }
  unit_cell_info['reference'] = { 'cell': reference_cell.parameters(), 'cell_volume': reference_cell.volume(),
                                  'lattice': reference_lattice, 'wavelength': reference_wavelength }

  Chatter.write("\nRandomly sampling %d x %d reflections for Monte Carlo iterations" % (mc_runs, sample_size))
  Debug.write("Refinements start with reference unit cell: %s" % reference_cell)

  MC = []
  MCconstrained = []
  used_index_range = flex.miller_index()
  used_two_theta_range_min = 1e300
  used_two_theta_range_max = 0
  used_reflections = set()

  for n in range(mc_runs): # MC sampling
    # Select sample_size reflections
    sample = flex.size_t(random.sample(range(len(all_miller_indices)), sample_size))
    used_reflections = used_reflections.union(set(sample))
    miller_indices = all_miller_indices.select(sample)
    two_thetas_obs = all_two_thetas_obs.select(sample)

    # Record
    span = miller.index_span(miller_indices)
    used_index_range.append(span.min())
    used_index_range.append(span.max())
    used_two_theta_range_min = min(used_two_theta_range_min, min(two_thetas_obs))
    used_two_theta_range_max = max(used_two_theta_range_max, max(two_thetas_obs))

    refined = _refinery(two_thetas_obs, miller_indices, reference_wavelength, reference_cell)
    MC.append(refined.unit_cell().parameters() + (refined.unit_cell().volume(),))
    Debug.write('Run %d refined to: %s' % (n, str(refined.unit_cell())))
    if reference_lattice is not None and reference_lattice is not 'aP':
      refined = _refinery(two_thetas_obs, miller_indices, reference_wavelength, reference_cell, reference_lattice[0])
      MCconstrained.append(refined.unit_cell().parameters() + (refined.unit_cell().volume(),))
      Debug.write('Run %d (constrained %s) refined to: %s' % (n, reference_lattice[0], str(refined.unit_cell())))

    if (n % 50) == 0:
      sys.stdout.write("\n%5s ." % (str(n) if n > 0 else ''))
    else:
      sys.stdout.write(".")
    sys.stdout.flush()

  assert used_two_theta_range_min < used_two_theta_range_max

  def stats_summary(l):
    mean = sum(l) / len(l)
    var = 0
    for y in l:
      var = var + ((y - mean) ** 2)
    popvar = var / (len(l)-1)
    popstddev = math.sqrt(popvar)
    stderr = popstddev / math.sqrt(len(l))
    return { 'mean': mean, 'variance': var, 'population_variance': popvar,
        'population_standard_deviation': popstddev, 'standard_error': stderr }

  print
  Chatter.write("")
  Chatter.write("Unit cell estimation based on %d Monte Carlo runs," % len(MC))
  span = miller.index_span(used_index_range)
  Chatter.write("drawn from miller indices between %s and %s" % (str(span.min()), str(span.max())))
  Chatter.write("with associated 2theta angles between %.3f and %.3f deg" % (used_two_theta_range_min, used_two_theta_range_max))
#.........这里部分代码省略.........