<gh_stars>100-1000
"Reading and writing of scalepack merge reflection files."
from __future__ import absolute_import, division, print_function

# Sample scalepack OUTPUT FILE
#    1
# -987
#    34.698    71.491    54.740    90.000   106.549    90.000 p21
#   0   0   4  3617.6   287.2
#   0   1   6 12951.7  1583.6 12039.2  1665.8
#
# Format: (3I4, 4F8.1)

from cctbx import uctbx
from cctbx import sgtbx
from cctbx import crystal
from cctbx import miller
from cctbx.array_family import flex
from libtbx.math_utils import iround
from libtbx import easy_pickle
import os
import sys
from six.moves import range
from six.moves import zip

class FormatError(Exception): pass

class reader(object):

  def __init__(self, file_handle, header_only=False):
    line = file_handle.readline()
    if (line.rstrip() != "    1"):
      raise FormatError("line 1: expecting '    1'")
    file_handle.readline() # ignore line 2
    line_error = "line 3: expecting unit cell parameters and space group label"
    line = file_handle.readline()
    if (len(line) < 63 or line[60] != ' '):
      raise FormatError(line_error)
    try:
      uc_params = [float(line[i * 10 : (i + 1) * 10]) for i in range(6)]
    except KeyboardInterrupt: raise
    except Exception:
      raise FormatError(line_error)
    self.unit_cell = uctbx.unit_cell(uc_params)
    self.space_group_symbol = line[61:].strip()
    if (len(self.space_group_symbol) == 0):
      raise FormatError(line_error)
    try:
      self.space_group_info = sgtbx.space_group_info(self.space_group_symbol)
    except KeyboardInterrupt: raise
    except Exception:
      self.space_group_info = None
    if (header_only): return
    self.miller_indices = flex.miller_index()
    self.i_obs = flex.double()
    self.sigmas = flex.double()
    self.anomalous = False
    line_count = 3
    while 1:
      line = file_handle.readline()
      line_count += 1
      line_error = "line %d: expecting (3I4, 4F8.1)" % line_count
      if (line == ""): break
      line = line.rstrip()
      if (len(line) == 0): continue
      line += (" " * 44)
      flds = []
      used = 0
      for width in (4,4,4,8,8,8,8):
        next_used = used + width
        flds.append(line[used:next_used].strip())
        used = next_used
      try:
        h = [int(flds[i]) for i in range(3)]
      except KeyboardInterrupt: raise
      except Exception:
        raise FormatError(line_error)
      for i in (0,1):
        j = 3+2*i
        if (len(flds[j])):
          try:
            i_obs, sigma = (float(flds[j]), float(flds[j+1]))
          except KeyboardInterrupt: raise
          except Exception:
            if (flds[j] == "*"*len(flds[j])) : # XXX overload
              continue
            raise FormatError(line_error)
          # XXX scalepack uses I=0, sigmaI=-1 to denote a missing Friedel
          # mate
          if (i_obs == 0) and (sigma == -1):
            continue
          if (i):
            h = [-e for e in h]
            self.anomalous = True
          if (sigma != 0 or (i_obs != 0 and i_obs != 1)):
            self.miller_indices.append(h)
            self.i_obs.append(i_obs)
            self.sigmas.append(sigma)

  def show_summary(self, out=None, prefix=""):
    if (out is None): out = sys.stdout
    self.unit_cell.show_parameters(f=out, prefix=prefix+"Unit cell: ")
    print(prefix + "Space group symbol:", self.space_group_symbol, file=out)
    print(prefix + "Anomalous flag:", self.anomalous, file=out)
    print(prefix + "Number of reflections:", self.miller_indices.size(), file=out)

  def crystal_symmetry(self):
    def have_r_with_axes_info():
      s = self.space_group_symbol.upper()
      return (s.startswith("R") and (s.endswith("R") or s.endswith("H")))
    return crystal.symmetry(
      unit_cell=self.unit_cell,
      space_group_info=self.space_group_info,
      correct_rhombohedral_setting_if_necessary=not have_r_with_axes_info())

  def as_miller_array(self,
        crystal_symmetry=None,
        force_symmetry=False,
        merge_equivalents=True,
        base_array_info=None,
        anomalous=None):
    if base_array_info is None:
      base_array_info = miller.array_info(source_type="scalepack_merge")
    crystal_symmetry_from_file = self.crystal_symmetry()
    if anomalous or self.anomalous:
      labels = ["I(+)", "SIGI(+)", "I(-)", "SIGI(-)"]
    else:
      labels = ["I", "SIGI"]
    return (miller.array(
      miller_set=miller.set(
        crystal_symmetry=crystal_symmetry_from_file.join_symmetry(
          other_symmetry=crystal_symmetry,
          force=force_symmetry),
        indices=self.miller_indices,
        anomalous_flag=anomalous or self.anomalous),
      data=self.i_obs,
      sigmas=self.sigmas)
      .set_info(base_array_info.customized_copy(
        labels=labels,
        crystal_symmetry_from_file=crystal_symmetry_from_file))
      .set_observation_type_xray_intensity())

  def as_miller_arrays(self,
        crystal_symmetry=None,
        force_symmetry=False,
        merge_equivalents=True,
        base_array_info=None,
        anomalous=None):
    return [self.as_miller_array(
      crystal_symmetry=crystal_symmetry,
      force_symmetry=force_symmetry,
      merge_equivalents=merge_equivalents,
      base_array_info=base_array_info,
      anomalous=anomalous,
    )]

def format_f8_1_or_i8(h, label, value):
  if (value < 1.e6): return "%8.1f" % value
  result = "%8d" % iround(value)
  if (len(result) > 8):
    raise ValueError(
      "Value is too large for scalepack merge format: hkl=%s, %s=%.6g" % (
        str(h).replace(" ",""), label, value))
  return result

def scale_intensities_if_necessary(miller_array,out=sys.stdout):
    value=miller_array.data().min_max_mean().max
    max_value=0.9999e+08
    if value > max_value: # scale it
      miller_array=miller_array.deep_copy().apply_scaling(target_max=max_value)
      print("NOTE: Scaled array (max=%s) to fit in Scalepack format" %(
         str(max_value)), file=out)
    return miller_array

def write(file_name=None, file_object=None, miller_array=None,
          space_group_symbol=None,
          line_1="    1",
          line_2=" -987",
          scale_intensities_for_scalepack_merge=False,
          out=sys.stdout):

  assert [file_name, file_object].count(None) == 1
  assert miller_array.is_xray_intensity_array() or \
    miller_array.is_xray_amplitude_array(), \
    "Input data must be amplitudes or intensities."
  assert miller_array.sigmas() is not None
  if (file_object is None):
    file_object = open(file_name, "w")
  if (not miller_array.is_xray_intensity_array()):
    miller_array = miller_array.f_as_f_sq()
  if (space_group_symbol is None):
    space_group_symbol = str(miller_array.space_group_info())
    if (not space_group_symbol.endswith(")")): # not universal Hermann-Mauguin
      space_group_symbol = space_group_symbol.replace(" ", "").lower()

  if scale_intensities_for_scalepack_merge: # 2014-01-07 TT
    miller_array=scale_intensities_if_necessary(miller_array,out=out)

  print(line_1, file=file_object)
  print(line_2, file=file_object)
  print(("%10.3f"*6) % miller_array.unit_cell().parameters(), end=' ', file=file_object)
  print(space_group_symbol, file=file_object)
  if (not miller_array.anomalous_flag()):
    for h,f,s in zip(miller_array.indices(),
                     miller_array.data(),
                     miller_array.sigmas()):
      print(((("%4d"*3) % h)
        + format_f8_1_or_i8(h, "intensity", f)
        + format_f8_1_or_i8(h, "sigma", s)), file=file_object)
  else:
    asu, matches = miller_array.match_bijvoet_mates()
    sel_pairs_plus = matches.pairs_hemisphere_selection("+")
    sel_pairs_minus = matches.pairs_hemisphere_selection("-")
    indices = asu.indices().select(sel_pairs_plus)
    data_plus = asu.data().select(sel_pairs_plus)
    data_minus = asu.data().select(sel_pairs_minus)
    sigmas_plus = asu.sigmas().select(sel_pairs_plus)
    sigmas_minus = asu.sigmas().select(sel_pairs_minus)
    for h,fp,sp,fm,sm in zip(indices, data_plus, sigmas_plus,
                                      data_minus, sigmas_minus):
      print(((("%4d"*3) % h)
        + format_f8_1_or_i8(h, "intensity", fp)
        + format_f8_1_or_i8(h, "sigma", sp)
        + format_f8_1_or_i8(h, "intensity", fm)
        + format_f8_1_or_i8(h, "sigma", sm)), file=file_object)
    sel_singles = matches.singles_hemisphere_selection("+")
    indices = asu.indices().select(sel_singles)
    data = asu.data().select(sel_singles)
    sigmas = asu.sigmas().select(sel_singles)
    for h,f,s in zip(indices, data, sigmas):
      print(((("%4d"*3) % h)
        + format_f8_1_or_i8(h, "intensity", f)
        + format_f8_1_or_i8(h, "sigma", s)), file=file_object)
    sel_singles = matches.singles_hemisphere_selection("-")
    indices = -asu.indices().select(sel_singles)
    data = asu.data().select(sel_singles)
    sigmas = asu.sigmas().select(sel_singles)
    for h,f,s in zip(indices, data, sigmas):
      print(((("%4d"*3) % h) + (" "*16)
        + format_f8_1_or_i8(h, "intensity", f)
        + format_f8_1_or_i8(h, "sigma", s)), file=file_object)
  file_object.close()

def run(args):
  to_pickle = "--pickle" in args
  for file_name in args:
    if (file_name.startswith("--")): continue
    s = reader(open(file_name, "r"))
    miller_array = s.as_miller_array(info="From file: "+file_name)
    miller_array.show_summary()
    if (to_pickle):
      pickle_file_name = os.path.split(file_name)[1] + ".pickle"
      print("Writing:", pickle_file_name)
      easy_pickle.dump(pickle_file_name, miller_array)
    print()

"Reading and writing of scalepack merge reflection files."
from __future__ import absolute_import, division, print_function

# Sample scalepack OUTPUT FILE
#    1
# -987
#    34.698    71.491    54.740    90.000   106.549    90.000 p21
#   0   0   4  3617.6   287.2
#   0   1   6 12951.7  1583.6 12039.2  1665.8
#
# Format: (3I4, 4F8.1)

from cctbx import uctbx
from cctbx import sgtbx
from cctbx import crystal
from cctbx import miller
from cctbx.array_family import flex
from libtbx.math_utils import iround
from libtbx import easy_pickle
import os
import sys
from six.moves import range
from six.moves import zip

class FormatError(Exception): pass

class reader(object):

  def __init__(self, file_handle, header_only=False):
    line = file_handle.readline()
    if (line.rstrip() != "    1"):
      raise FormatError("line 1: expecting '    1'")
    file_handle.readline() # ignore line 2
    line_error = "line 3: expecting unit cell parameters and space group label"
    line = file_handle.readline()
    if (len(line) < 63 or line[60] != ' '):
      raise FormatError(line_error)
    try:
      uc_params = [float(line[i * 10 : (i + 1) * 10]) for i in range(6)]
    except KeyboardInterrupt: raise
    except Exception:
      raise FormatError(line_error)
    self.unit_cell = uctbx.unit_cell(uc_params)
    self.space_group_symbol = line[61:].strip()
    if (len(self.space_group_symbol) == 0):
      raise FormatError(line_error)
    try:
      self.space_group_info = sgtbx.space_group_info(self.space_group_symbol)
    except KeyboardInterrupt: raise
    except Exception:
      self.space_group_info = None
    if (header_only): return
    self.miller_indices = flex.miller_index()
    self.i_obs = flex.double()
    self.sigmas = flex.double()
    self.anomalous = False
    line_count = 3
    while 1:
      line = file_handle.readline()
      line_count += 1
      line_error = "line %d: expecting (3I4, 4F8.1)" % line_count
      if (line == ""): break
      line = line.rstrip()
      if (len(line) == 0): continue
      line += (" " * 44)
      flds = []
      used = 0
      for width in (4,4,4,8,8,8,8):
        next_used = used + width
        flds.append(line[used:next_used].strip())
        used = next_used
      try:
        h = [int(flds[i]) for i in range(3)]
      except KeyboardInterrupt: raise
      except Exception:
        raise FormatError(line_error)
      for i in (0,1):
        j = 3+2*i
        if (len(flds[j])):
          try:
            i_obs, sigma = (float(flds[j]), float(flds[j+1]))
          except KeyboardInterrupt: raise
          except Exception:
            if (flds[j] == "*"*len(flds[j])) : # XXX overload
              continue
            raise FormatError(line_error)
          # XXX scalepack uses I=0, sigmaI=-1 to denote a missing Friedel
          # mate
          if (i_obs == 0) and (sigma == -1):
            continue
          if (i):
            h = [-e for e in h]
            self.anomalous = True
          if (sigma != 0 or (i_obs != 0 and i_obs != 1)):
            self.miller_indices.append(h)
            self.i_obs.append(i_obs)
            self.sigmas.append(sigma)

  def show_summary(self, out=None, prefix=""):
    if (out is None): out = sys.stdout
    self.unit_cell.show_parameters(f=out, prefix=prefix+"Unit cell: ")
    print(prefix + "Space group symbol:", self.space_group_symbol, file=out)
    print(prefix + "Anomalous flag:", self.anomalous, file=out)
    print(prefix + "Number of reflections:", self.miller_indices.size(), file=out)

  def crystal_symmetry(self):
    def have_r_with_axes_info():
      s = self.space_group_symbol.upper()
      return (s.startswith("R") and (s.endswith("R") or s.endswith("H")))
    return crystal.symmetry(
      unit_cell=self.unit_cell,
      space_group_info=self.space_group_info,
      correct_rhombohedral_setting_if_necessary=not have_r_with_axes_info())

  def as_miller_array(self,
        crystal_symmetry=None,
        force_symmetry=False,
        merge_equivalents=True,
        base_array_info=None,
        anomalous=None):
    if base_array_info is None:
      base_array_info = miller.array_info(source_type="scalepack_merge")
    crystal_symmetry_from_file = self.crystal_symmetry()
    if anomalous or self.anomalous:
      labels = ["I(+)", "SIGI(+)", "I(-)", "SIGI(-)"]
    else:
      labels = ["I", "SIGI"]
    return (miller.array(
      miller_set=miller.set(
        crystal_symmetry=crystal_symmetry_from_file.join_symmetry(
          other_symmetry=crystal_symmetry,
          force=force_symmetry),
        indices=self.miller_indices,
        anomalous_flag=anomalous or self.anomalous),
      data=self.i_obs,
      sigmas=self.sigmas)
      .set_info(base_array_info.customized_copy(
        labels=labels,
        crystal_symmetry_from_file=crystal_symmetry_from_file))
      .set_observation_type_xray_intensity())

  def as_miller_arrays(self,
        crystal_symmetry=None,
        force_symmetry=False,
        merge_equivalents=True,
        base_array_info=None,
        anomalous=None):
    return [self.as_miller_array(
      crystal_symmetry=crystal_symmetry,
      force_symmetry=force_symmetry,
      merge_equivalents=merge_equivalents,
      base_array_info=base_array_info,
      anomalous=anomalous,
    )]

def format_f8_1_or_i8(h, label, value):
  if (value < 1.e6): return "%8.1f" % value
  result = "%8d" % iround(value)
  if (len(result) > 8):
    raise ValueError(
      "Value is too large for scalepack merge format: hkl=%s, %s=%.6g" % (
        str(h).replace(" ",""), label, value))
  return result

def scale_intensities_if_necessary(miller_array,out=sys.stdout):
    value=miller_array.data().min_max_mean().max
    max_value=0.9999e+08
    if value > max_value: # scale it
      miller_array=miller_array.deep_copy().apply_scaling(target_max=max_value)
      print("NOTE: Scaled array (max=%s) to fit in Scalepack format" %(
         str(max_value)), file=out)
    return miller_array

def write(file_name=None, file_object=None, miller_array=None,
          space_group_symbol=None,
          line_1="    1",
          line_2=" -987",
          scale_intensities_for_scalepack_merge=False,
          out=sys.stdout):

  assert [file_name, file_object].count(None) == 1
  assert miller_array.is_xray_intensity_array() or \
    miller_array.is_xray_amplitude_array(), \
    "Input data must be amplitudes or intensities."
  assert miller_array.sigmas() is not None
  if (file_object is None):
    file_object = open(file_name, "w")
  if (not miller_array.is_xray_intensity_array()):
    miller_array = miller_array.f_as_f_sq()
  if (space_group_symbol is None):
    space_group_symbol = str(miller_array.space_group_info())
    if (not space_group_symbol.endswith(")")): # not universal Hermann-Mauguin
      space_group_symbol = space_group_symbol.replace(" ", "").lower()

  if scale_intensities_for_scalepack_merge: # 2014-01-07 TT
    miller_array=scale_intensities_if_necessary(miller_array,out=out)

  print(line_1, file=file_object)
  print(line_2, file=file_object)
  print(("%10.3f"*6) % miller_array.unit_cell().parameters(), end=' ', file=file_object)
  print(space_group_symbol, file=file_object)
  if (not miller_array.anomalous_flag()):
    for h,f,s in zip(miller_array.indices(),
                     miller_array.data(),
                     miller_array.sigmas()):
      print(((("%4d"*3) % h)
        + format_f8_1_or_i8(h, "intensity", f)
        + format_f8_1_or_i8(h, "sigma", s)), file=file_object)
  else:
    asu, matches = miller_array.match_bijvoet_mates()
    sel_pairs_plus = matches.pairs_hemisphere_selection("+")
    sel_pairs_minus = matches.pairs_hemisphere_selection("-")
    indices = asu.indices().select(sel_pairs_plus)
    data_plus = asu.data().select(sel_pairs_plus)
    data_minus = asu.data().select(sel_pairs_minus)
    sigmas_plus = asu.sigmas().select(sel_pairs_plus)
    sigmas_minus = asu.sigmas().select(sel_pairs_minus)
    for h,fp,sp,fm,sm in zip(indices, data_plus, sigmas_plus,
                                      data_minus, sigmas_minus):
      print(((("%4d"*3) % h)
        + format_f8_1_or_i8(h, "intensity", fp)
        + format_f8_1_or_i8(h, "sigma", sp)
        + format_f8_1_or_i8(h, "intensity", fm)
        + format_f8_1_or_i8(h, "sigma", sm)), file=file_object)
    sel_singles = matches.singles_hemisphere_selection("+")
    indices = asu.indices().select(sel_singles)
    data = asu.data().select(sel_singles)
    sigmas = asu.sigmas().select(sel_singles)
    for h,f,s in zip(indices, data, sigmas):
      print(((("%4d"*3) % h)
        + format_f8_1_or_i8(h, "intensity", f)
        + format_f8_1_or_i8(h, "sigma", s)), file=file_object)
    sel_singles = matches.singles_hemisphere_selection("-")
    indices = -asu.indices().select(sel_singles)
    data = asu.data().select(sel_singles)
    sigmas = asu.sigmas().select(sel_singles)
    for h,f,s in zip(indices, data, sigmas):
      print(((("%4d"*3) % h) + (" "*16)
        + format_f8_1_or_i8(h, "intensity", f)
        + format_f8_1_or_i8(h, "sigma", s)), file=file_object)
  file_object.close()

def run(args):
  to_pickle = "--pickle" in args
  for file_name in args:
    if (file_name.startswith("--")): continue
    s = reader(open(file_name, "r"))
    miller_array = s.as_miller_array(info="From file: "+file_name)
    miller_array.show_summary()
    if (to_pickle):
      pickle_file_name = os.path.split(file_name)[1] + ".pickle"
      print("Writing:", pickle_file_name)
      easy_pickle.dump(pickle_file_name, miller_array)
    print()
# Sample scalepack OUTPUT FILE #    1 # -987 #    34.698    71.491    54.740    90.000   106.549    90.000 p21 #   0   0   4  3617.6   287.2 #   0   1   6 12951.7  1583.6 12039.2  1665.8 # # Format: (3I4, 4F8.1) # ignore line 2 # XXX overload # XXX scalepack uses I=0, sigmaI=-1 to denote a missing Friedel # mate # scale it # not universal Hermann-Mauguin # 2014-01-07 TT