#!/usr/bin/env python
#
# dials.command_line.index_symops.py
#
#  Copyright (C) 2015 Diamond Light Source
#
#  Author: Graeme Winter
#
#  This code is distributed under the BSD license, a copy of which is
#  included in the root directory of this package.

from __future__ import division

from libtbx.phil import command_line
import iotbx.phil
from cctbx import sgtbx
from dials.util.options import OptionParser
from dials.util.options import flatten_reflections, flatten_experiments

help_message = '''

This program can be used to (i) generate the symmetry from the experiment
and apply to the input reflections, and (ii) for each symop for that symmetry
attempt to calculate the CC on that operation within the strong spot list.

  dials.index_symops experiment.json indexed.pickle [d_min=3.0] [d_max=10.0]

'''

phil_scope = iotbx.phil.parse("""
d_min = 0
  .type = float
  .help = "High resolution limit to use for analysis"
d_max = 0
  .type = float
  .help = "Low resolution limit to use for analysis"
""", process_includes=True)


def run(args):
  import libtbx.load_env
  usage = "%s [options] experiment.json indexed.pickle" % \
    libtbx.env.dispatcher_name

  parser = OptionParser(
    usage=usage,
    phil=phil_scope,
    read_reflections=True,
    read_experiments=True,
    check_format=False,
    epilog=help_message)

  params, options = parser.parse_args(show_diff_phil=True)

  reflections = flatten_reflections(params.input.reflections)
  experiments = flatten_experiments(params.input.experiments)
  if len(reflections) == 0 or len(experiments) == 0:
    parser.print_help()
    return
  assert(len(reflections) == 1)
  assert(len(experiments) == 1)
  experiment = experiments[0]
  crystal = experiment.crystal
  reflections = reflections[0]

  from copy import deepcopy

  original_miller_indices = reflections['miller_index']

  space_group = crystal.get_space_group()
  unit_cell = crystal.get_unit_cell()
  from cctbx.crystal import symmetry as crystal_symmetry
  cs = crystal_symmetry(unit_cell, space_group.type().lookup_symbol())

  from cctbx.miller import set as miller_set

  ms = miller_set(cs, original_miller_indices)
  ms = ms.array(reflections['intensity.sum.value'])

  if params.d_min or params.d_max:
    ms = ms.resolution_filter(d_min=params.d_min, d_max=params.d_max)

  print '%d reflections analysed' % ms.size()
  print '%10s %6s %5s' % ('Symop', 'Nref', 'CC')

  for smx in space_group.smx():
    reindexed = deepcopy(reflections)
    miller_indices = reflections['miller_index']
    reindexed_miller_indices = sgtbx.change_of_basis_op(smx).apply(
      miller_indices)
    rms = miller_set(cs, reindexed_miller_indices)
    rms = rms.array(reflections['intensity.sum.value'])
    if params.d_min or params.d_max:
      rms = rms.resolution_filter(d_min=params.d_min, d_max=params.d_max)
    intensity, intensity_rdx = rms.common_sets(ms)
    cc = intensity.correlation(intensity_rdx).coefficient()

    print '%10s %6d %.3f' % (smx, intensity.size(), cc)

if __name__ == '__main__':
  import sys
  from libtbx.utils import show_times_at_exit
  show_times_at_exit()
  run(sys.argv[1:])

@@ -1003,16 +1003,125 @@
 alysis%22%0A
+symop_threshold = 0%0A  .type = float%0A  .help = %22Thereshold above which we consider a symmetry operator true.%22%0A
 %22%22%22, pro
@@ -1897,16 +1897,54 @@
 deepcopy
+%0A  from dials.array_family import flex
 %0A%0A  orig
@@ -2346,16 +2346,83 @@
 .value'%5D
+ /%0A                flex.sqrt(reflections%5B'intensity.sum.variance'%5D)
 )%0A%0A  if 
@@ -2618,16 +2618,36 @@
  'CC')%0A%0A
+  true_symops = %5B%5D%0A%0A
   for sm
@@ -2952,16 +2952,87 @@
 .value'%5D
+ /%0A                    flex.sqrt(reflections%5B'intensity.sum.variance'%5D)
 )%0A    if
@@ -3248,16 +3248,158 @@
 ient()%0A%0A
+    accept = ''%0A%0A    if params.symop_threshold:%0A      if cc %3E params.symop_threshold:%0A        true_symops.append(smx)%0A        accept = '***'%0A%0A
     prin
@@ -3414,16 +3414,19 @@
 %256d %25.3f
+ %25s
 ' %25 (smx
@@ -3447,16 +3447,326 @@
 ze(), cc
+, accept)%0A%0A  if params.symop_threshold:%0A    from cctbx.sgtbx import space_group%0A    sg = space_group()%0A    for symop in true_symops:%0A      sg = sg.expand_smx(symop)%0A    sg_symbols = sg.match_tabulated_settings()%0A    print 'Derived point group from symmetry operations: %25s' %25 %5C%0A      sg_symbols.hermann_mauguin(
 )%0A%0Aif __

Added symop threshold for analysis (user set at the moment) and derivation of the apparent rotational symmetry group / point group: The following parameters have been modified: