# TODO: more helpful output?

from __future__ import print_function, division # Python 3 stuff
import phaser
from phaser import pickle_support
from phaser.phenix_interface import master_params, phaser_style, match_modes
from phaser.phenix_interface import phaser_result
from phaser.phenix_interface import callbacks
from iotbx import file_reader
import iotbx.phil
from cctbx import crystal, sgtbx
from scitbx.array_family import flex # also required for unpickling?
from libtbx import object_oriented_patterns as oop
from libtbx.utils import Sorry, import_python_object, multi_out, null_out
from libtbx.phil import interface as phil_interface
from libtbx.str_utils import StringIO
from libtbx import easy_pickle
import os
import re
import sys

pickle_support.enable()

mr_modes = ["MR_AUTO","MR_FRF","MR_FTF","MR_LLG","MR_PAK","MR_RNP", "MR_BRF","MR_BTF"]
ep_modes = ["EP_AUTO","EP_SAD"]

def translate_keyword (name, style, prefix=None, multiple=False) :
  method_name = "set"
  if (multiple) :
    method_name = "add"
  if prefix is not None :
    method_name = "%s_" % prefix
  fields = [ s[:4] for s in name.split("_") ]
  method_name += "_".join([ s.upper() for s in fields ])
  if (style.as_multiple) and (method_name.startswith("set")) :
    method_name = method_name.replace("set", "add")
  return method_name

def convert_keywords (input,
                      keyword_phil,
                      mode,
                      script_out,
                      inherit_mode=False,
                      default_method_base=None,
                      offset_start=0) :
  s = script_out
  for keyword in keyword_phil.objects[offset_start:] :
    if keyword.is_template != 0 :
      continue
    style = phaser_style(keyword)
    if style.ignore :
      continue
    elif match_modes(mode, style.mode, inherit_mode) :
      if keyword.is_scope :
        method_base = default_method_base
        if method_base is None :
          method_base = translate_keyword(name=keyword.name,
                                              style=style,
                                              multiple=keyword.multiple)
        elif (keyword.multiple) :
          method_base = "add" + method_base[3:]
        if style.listargs or style.array :
          # setXXXX_YYYY(param1, ...)
          args = []
          for subkeywd in keyword.objects :
            assert not subkeywd.multiple
            value = subkeywd.extract()
            if value is None :
              keyword.show()
              raise Sorry("%s must not be None." % subkeywd.full_path())
            args.append(value)
          method = getattr(input, method_base)
          if style.listargs :
            print("input.%s(%s)" % (method_base,
              ", ".join([str(a) for a in args])),
              file=s)
            method(*args)
          else :
            print("input.%s([%s])" % (method_base,
              ", ".join([str(a) for a in args])),
              file=s)
            method(args)
        elif style.kwargs :
          kwargs = {}
          kwargs_out = []
          for subkeywd in keyword.objects :
            assert not subkeywd.multiple
            value = subkeywd.extract()
            if value is None :
              keyword.show()
              raise Sorry("%s must not be None." % subkeywd.full_path())
            kwargs[subkeywd.name] = value
            value_out = str(value)
            if isinstance(value, str) :
              value_out = "'%s'" % value
            kwargs_out.append("%s=%s" % (subkeywd.name, value_out))
          method = getattr(input, method_base)
          print("input.%s(%s)" % (method_base, ", ".join(kwargs_out)), file=s)
          try :
            method(**kwargs)
          except Exception, e :
            if (type(e).__name__ == "ArgumentError") :
              raise RuntimeError(
                "invalid boost.python arguments:\n%s\nKeywords:\n%s" %
                  (str(e), str(kwargs)))
            else :
              raise
        elif style.silent :
          # XXX: this should work recursively...
          convert_keywords(input=input,
                           keyword_phil=keyword,
                           mode=mode,
                           script_out=script_out,
                           inherit_mode=True)
        else :
          # setXXXX_YYYY(param1)
          # setXXXX_ZZZZ(param2)
          # ...
          current_offset = 0
          for subkeywd in keyword.objects :
            if subkeywd.is_template != 0 :
              continue
            style2 = phaser_style(subkeywd)
            if style2.ignore :
              current_offset += 1
              continue
            if (not match_modes(mode, style2.mode, True)) :
              continue
            if subkeywd.is_scope :
              subkeywd_base = method_base
              if style2.parent_silent :
                subkeywd_base = None
              convert_keywords(input=input,
                               keyword_phil=keyword,
                               mode=mode,
                               script_out=script_out,
                               inherit_mode=True,
                               default_method_base=subkeywd_base,
                               offset_start=current_offset)
              current_offset += 1
              break
              #continue
            value = subkeywd.extract()
            if value is None :
              current_offset += 1
              continue
            method_base = translate_keyword(name=keyword.name,
                                            style=style,
                                            prefix=default_method_base)
            if subkeywd.multiple or style2.as_multiple :
              method_base = "add" + method_base[3:]
            if style2.silent :
              method_name = method_base
            else :
              method_name = translate_keyword(name=subkeywd.name,
                                              style=style2,
                                              prefix=method_base)
            method = getattr(input, method_name)
            if value is not None :
              if isinstance(value, str) :
                print("input.%s('%s')" % (method_name, value), file=s)
                method(value)
              else :
                if isinstance(value, sgtbx.space_group_info) :
                  print("input.%s(%s)" % (method_name, str(value)), file=s)
                  method(str(value))
                elif (isinstance(value, list)) and (style2.as_multiple) :
                  for val_ in value :
                    print("input.%s(%s)" % (method_name, str(val_)), file=s)
                    method(val_)
                else :
                  print("input.%s(%s)" % (method_name, str(value)), file=s)
                  method(value)
            else :
              print("input.%s()" % method_name, file=s)
              method()
            current_offset += 1
      elif keyword.is_definition :
        value = keyword.extract()
        if value is None :
          continue
        method_name = translate_keyword(name=keyword.name,
                                        style=style,
                                        prefix=default_method_base,
                                        multiple=keyword.multiple)
        #print(keyword.full_path())
        method = getattr(input, method_name)
        if isinstance(value, str) :
          print("input.%s('%s')" % (method_name, value), file=s)
        else :
          print("input.%s(%s)" % (method_name, str(value)), file=s)
        method(value)
    else :
      continue

#-----------------------------------------------------------------------
# REFLECTION DATA HANDLING
def get_miller_array (miller_arrays, labin, file_name) :
  for miller_array in miller_arrays :
    array_labels = miller_array.info().label_string()
    if array_labels == labin :
      return miller_array
  raise Sorry("Couldn't find array %s in file %s." % (labin, file_name))

def crystal_sanity_check (space_group, unit_cell, master_space_group,
    master_unit_cell, array_name) :
  if type(space_group).__name__ == "space_group_info" :
    space_group = space_group.group()
  if type(master_space_group).__name__ == "space_group_info" :
    master_space_group = master_space_group.group()
  if not None in [space_group, master_space_group] :
    derived_point_group = space_group.build_derived_point_group()
    if master_space_group is not None :
      master_point_group=master_space_group.build_derived_point_group()
      if derived_point_group != master_point_group :
        raise Sorry(("The space group for the data in %s (%s) does "+
          "not match the point group of the overall space group (%s).") %
          (array_name, str(space_group), str(master_space_group)))
  if not None in [unit_cell, master_unit_cell] :
    if not master_unit_cell.is_similar_to(unit_cell) :
      raise Sorry(("The unit cell for the data in %s (%s) is "+
        "significantly different than the overall unit cell (%s).") %
        (array_name, str(unit_cell), str(master_unit_cell)))

class process_mr_hklin (object) :
  def __init__ (self, params, symmetry, out=sys.stdout) :
    array_name="%s:%s" % (params.phaser.hklin,params.phaser.labin)
    hkl_file = file_reader.any_file(params.phaser.hklin)
    hkl_file.check_file_type("hkl")
    input_array = get_miller_array(
      miller_arrays=hkl_file.file_object.as_miller_arrays(),
      labin=params.phaser.labin,
      file_name=params.phaser.hklin)
    if (input_array.sigmas() is None) :
      raise Sorry("Experimental sigmas are required.")
    is_changed = False
    if (symmetry is not None) :
      input_array = input_array.customized_copy(
        crystal_symmetry=symmetry)
      is_changed = True
    elif (params.phaser.crystal_symmetry.space_group is not None) :
      input_array = input_array.customized_copy(
        space_group_info=params.phaser.crystal_symmetry.space_group)
      is_changed = True
    elif (params.phaser.crystal_symmetry.unit_cell is not None) :
      input_array = input_array.customized_copy(
        unit_cell=params.phaser.crystal_symmetry.unit_cell)
      is_changed = True
    self.space_group = input_array.space_group_info()
    self.unit_cell = input_array.unit_cell()
    crystal_sanity_check(
      space_group=input_array.space_group(),
      unit_cell=self.unit_cell,
      master_space_group=params.phaser.crystal_symmetry.space_group,
      master_unit_cell=params.phaser.crystal_symmetry.unit_cell,
      array_name=array_name)
    if (self.space_group is None) or (self.unit_cell is None) :
      raise Sorry(("No build-in crystal information for the data in %s:%s; "+
        "please define the crystal symmetry parameters first.") %
        (params.phaser.hklin, params.phaser.labin))
    if input_array.is_xray_intensity_array() :
      from cctbx import french_wilson
      input_array = french_wilson.french_wilson_scale(
        miller_array=input_array,
        merge=True,
        log=out)
      if (input_array is None) :
        raise Sorry("French-Wilson treatment of intensities could not be "+
          "performed because suitable resolution binning could not be "+
          "determined.  This may indicate a severely incomplete dataset or "+
          "corrupted input file; please inspect your data and convert to "+
          "amplitudes manually if necessary.")
      is_changed = True
    if input_array.anomalous_flag() :
      input_array = input_array.average_bijvoet_mates()
      is_changed = True
    elif (not input_array.is_xray_amplitude_array()) :
      raise Sorry("%s is not X-ray data." % array_labels)
    if (not input_array.is_unique_set_under_symmetry()) :
      input_array = input_array.merge_equivalents().array()
      is_changed = True
    bad_sigmas = (input_array.sigmas() <= 0)
    if (bad_sigmas.count(True) > 0) :
      input_array = input_array.select(~bad_sigmas)
      is_changed = True
    #if input_array.info() is None : # XXX: will this work?
    if is_changed :
      hklin = "%s_data.mtz" % params.phaser.keywords.general.root
      if not params.test_mode :
        hklin = os.path.abspath(hklin)
      input_array = input_array.map_to_asu()
      # Save final data as MTZ file
      mtz_dataset = input_array.as_mtz_dataset(column_root_label="F")
      mtz_object = mtz_dataset.mtz_object()
      mtz_object.write(file_name=hklin)
      del mtz_object
      params.phaser.hklin = hklin
      labout = "F,SIGF"
      params.phaser.labin = labout
    (self.d_max, self.d_min) = input_array.d_max_min()
    self.hklin = params.phaser.hklin

class process_ep_crystal (object) :
  def __init__ (self, params, symmetry, out=sys.stdout) :
    self.d_min = params.phaser.keywords.resolution.high
    self.d_max = params.phaser.keywords.resolution.low
    self.space_group = params.phaser.crystal_symmetry.space_group
    self.unit_cell = params.phaser.crystal_symmetry.unit_cell
    if self.space_group is not None :
      derived_point_group = None
    self.miller_arrays = {}
    input_files = {}
    # Pass 1: collect arrays symmetry information
    for xtal in params.phaser.crystal :
      xtal_id = xtal.xtal_id
      for dataset in xtal.dataset :
        if None in [dataset.wave_id, dataset.hklin, dataset.labin] :
          continue
        wave_id = dataset.wave_id
        if not dataset.hklin in input_files :
          hkl_file = file_reader.any_file(dataset.hklin)
          hkl_file.check_file_type("hkl")
          input_files[dataset.hklin] = hkl_file
        hkl_file = input_files[dataset.hklin]
        array_name = "%s:%s" % (dataset.hklin, dataset.labin)
        input_array = get_miller_array(
          miller_arrays=hkl_file.file_object.as_miller_arrays(),
          labin=dataset.labin,
          file_name=dataset.hklin)
        if not input_array.anomalous_flag() :
          raise Sorry(("Anomalous data required in mode '%s' (offending "+
            "data: %s).") % (params.phaser.mode, array_name))
        if (input_array.sigmas() is None) :
          raise Sorry("Experimental sigmas are required.")
        array_sg = input_array.space_group_info()
        array_uc = input_array.unit_cell()
        crystal_sanity_check(
          space_group=array_sg,
          unit_cell=array_uc,
          master_space_group=self.space_group,
          master_unit_cell=self.unit_cell,
          array_name=array_name)
        if (self.space_group is None) :
          self.space_group = array_sg
        if (self.unit_cell is None) :
          self.unit_cell = array_uc
        self.miller_arrays[(xtal_id, wave_id)] = input_array
    if symmetry is None :
      if (self.space_group is None) or (self.unit_cell is None) :
        raise Sorry("Couldn't determine automatic crystal symmetry "+
          "information from the provided reflection file(s).  Please "+
          "specify the space group and unit cell.")
      symmetry = crystal.symmetry(
        unit_cell=self.unit_cell,
        space_group_info=self.space_group,
        raise_sorry_if_incompatible_unit_cell=True)
    not_alnum = re.compile("[^a-zA-Z0-9_\-]{1,}")
    mtz_dataset = None
    self.data = None
    hklin = "%s_data.mtz" % params.phaser.keywords.general.root
    if not params.test_mode :
      hklin = os.path.abspath(hklin)
    # Pass 2: standardize data and create MTZ file
    for xtal in params.phaser.crystal :
      xtal_id = xtal.xtal_id
      for dataset in xtal.dataset :
        wave_id = dataset.wave_id
        column_root_label = "F_%s_%s" % (xtal_id, wave_id)
        if not_alnum.search(column_root_label) is not None :
          raise Sorry(("The crystal ID and dataset ID must only contain "+
            "alphanumeric characters or underscores (current: '%s','%s').") %
            (xtal_id, wave_id))
        elif len(column_root_label) > 24 :
          raise Sorry(("The crystal ID and dataset ID combined must not be "+
            "more than 21 characters (current: '%s'+'%s' = %dc.).") %
            (xtal_id, wave_id, len(column_root_label) - 3))
        input_arry = self.miller_arrays[(xtal_id, wave_id)]
        input_array = input_array.customized_copy(crystal_symmetry=symmetry)
        (d_max, d_min) = input_array.d_max_min()
        if self.d_max is None or self.d_max < d_max :
          self.d_max = d_max
        if self.d_min is None or self.d_min > d_min :
          self.d_min = d_min
        if input_array.is_xray_intensity_array() :
          from cctbx import french_wilson
          input_array = french_wilson.french_wilson_scale(
            miller_array=input_array,
            merge=True,
            log=out)
        elif (not input_array.is_xray_amplitude_array()) :
          raise Sorry("%s is not X-ray data." % array_name)
        input_array = input_array.map_to_asu()
        if (not input_array.is_unique_set_under_symmetry()) :
          input_array = input_array.merge_equivalents().array()
        bad_sigmas = (input_array.sigmas() <= 0)
        input_array = input_array.select(~bad_sigmas)
        if (self.data is None) :
          self.data = input_array
        if mtz_dataset is None :
          mtz_dataset = input_array.as_mtz_dataset(
            column_root_label=column_root_label)
        else :
          mtz_dataset.add_miller_array(
            miller_array=input_array,
            column_root_label=column_root_label)
        labels = "%s(+),SIG%s(+),%s(-),SIG%s(-)" % tuple([column_root_label]*4)
        dataset.hklin = hklin
        dataset.labin = labels
    mtz_object = mtz_dataset.mtz_object()
    mtz_object.write(file_name=hklin)
    del mtz_object
    self.hklin = hklin
#    print(("Converted reflection data for use in Phaser; new file is "+
#          "%s with labels %s.") % (params.phaser.hklin, labout), file=out)

#-----------------------------------------------------------------------
# MAIN DRIVER
class phaser_parameter_interpreter (object) :
  def __init__ (self, index, working_dir, out=sys.stdout) :
    params = index.working_phil.extract()
    self.params = params
    mode = params.phaser.mode
    s = StringIO()
    self._script_out = s
    self.input_data = None # XXX hack for density modification
    print("import phaser", file=s)
    if mode.startswith("MR_") :
      print("data = phaser.InputMR_DAT()", file=s)
      data = phaser.InputMR_DAT()
      run_data_class = phaser.runMR_DAT
    elif mode.startswith("EP_") :
      print("data = phaser.InputEP_DAT()", file=s)
      data = phaser.InputEP_DAT()
      run_data_class = phaser.runEP_DAT
    else :
      data = None
      run_data_class = None
    phaser_out = phaser.Output()
    phaser_out.setPackagePhenix(out)
    if data is not None :
      if params.dry_run or params.show_script :
        data.setMUTE(True)
    not_alnum = re.compile("[^a-zA-Z0-9_\-\.]{1,}")
    if params.phaser.keywords.general.root is None :
      raise Sorry("Please define a root filename for this run.")
    elif not_alnum.search(params.phaser.keywords.general.root) is not None :
      raise Sorry("The root filename must be alphanumeric (underscores OK).")
    symmetry = None
    if (params.phaser.crystal_symmetry.unit_cell is not None and
        params.phaser.crystal_symmetry.space_group is not None) :
      symmetry = crystal.symmetry(
        unit_cell=params.phaser.crystal_symmetry.unit_cell,
        space_group_info=params.phaser.crystal_symmetry.space_group,
        raise_sorry_if_incompatible_unit_cell=True)
    #-------------------------------------------------------------------
    # MR DATA (HKLIN)
    process_data = None
    if match_modes(mode, index.get_phaser_modes("phaser.hklin")) : # XXX: MR_*
      if params.phaser.hklin is None :
        raise Sorry("You must specify a reflections file in mode %s." % mode)
      elif params.phaser.labin is None :
        # MR modes
        raise Sorry("You must specify reflection data labels in mode %s."%mode)
      process_data = process_mr_hklin(params=params,
                                      symmetry=symmetry,
                                      out=out)
      print("data.setHKLI('%s')" % params.phaser.hklin, file=s)
      data.setHKLI(params.phaser.hklin)
      labels = params.phaser.labin.split(",")
      assert len(labels) == 2 # should have been taken care of already
      print("data.setLABI('%s','%s')" % (labels[0], labels[1]), file=s)
      data.setLABI(labels[0], labels[1])
    #-------------------------------------------------------------------
    # EP DATA
    elif match_modes(mode, index.get_phaser_modes("phaser.crystal")) :
      # TODO: eventually, allow multiples of everything
      if len(params.phaser.crystal) == 0 :
        raise Sorry("Crystal data required for mode %s." % mode)
      elif len(params.phaser.crystal) > 1 :
        raise Sorry("Only one crystal supported at this time.")
      xtal = params.phaser.crystal[0]
      if len(xtal.dataset) > 1 :
        raise Sorry("Only one wavelength supported at this time.")
      process_data = process_ep_crystal(params=params,
                                        symmetry=symmetry,
                                        out=out)
      # XXX: for now, still have one global hklin
      print("data.setHKLI('%s')" % process_data.hklin, file=s)
      data.setHKLI(process_data.hklin)
      self.input_data = process_data.data
      for xtal in params.phaser.crystal :
        xtalid = xtal.xtal_id
        if not None in [xtal.pdb_file, xtal.ha_file] :
          raise Sorry("For each crystal, you may supply heavy-atom sites as "+
            "EITHER a PDB file OR 'ha' file, but not both.")
        for dataset in xtal.dataset :
          waveid = dataset.wave_id
          labels = dataset.labin.split(",")
          if len(labels) != 4 :
            raise Sorry(("The data labels %s do not appear to correspond to "+
             "the appropriate data type.  For experimental phasing, data must "+
              "be merged, anomalous amplitudes in MTZ format.  The graphical "+
              "reflection file editor (phenix.reflection_file_editor) can be "+
              "used to convert the data as neeeded.") % dataset.labin)
          print("data.addCRYS_ANOM_LABI('%s', '%s', '%s', '%s', '%s', '%s')" %
            (xtalid, waveid, labels[0], labels[1], labels[2],labels[3]), file=s)
          data.addCRYS_ANOM_LABI(xtalid, waveid, labels[0], labels[1],
                                 labels[2], labels[3])
    #-------------------------------------------------------------------
    # SYMMETRY
    if match_modes(mode, index.get_phaser_modes("phaser.crystal_symmetry")) :
      if params.phaser.crystal_symmetry.space_group is None :
        if process_data.space_group is not None :
          params.phaser.crystal_symmetry.space_group = process_data.space_group
        else :
          raise Sorry("Couldn't determine space group information - please "+
                      "supply crystal symmetry parameters.")
      space_group = str(params.phaser.crystal_symmetry.space_group)
      print("data.setSPAC_NAME('%s')" % space_group, file=s)
      data.setSPAC_NAME(space_group)
      if params.phaser.crystal_symmetry.unit_cell is None :
        if process_data.unit_cell is not None :
          params.phaser.crystal_symmetry.unit_cell = process_data.unit_cell
        else :
          raise Sorry("Couldn't determine unit cell information - please "+
                      "supply crystal symmetry parameters.")
      unit_cell = params.phaser.crystal_symmetry.unit_cell.parameters()
      print("data.setCELL6([%s])" %
            ", ".join([ "%.3f" % x for x in unit_cell ]), file=s)
      data.setCELL6(unit_cell)
    if process_data is not None :
      if params.phaser.keywords.resolution.high is None :
        params.phaser.keywords.resolution.high = process_data.d_min
      if params.phaser.keywords.resolution.low is None :
        params.phaser.keywords.resolution.low = process_data.d_max
    #-------------------------------------------------------------------
    # MAIN INPUT OBJECT
    input_class = import_python_object(import_path="phaser.Input%s" % mode,
      error_prefix="",
      target_must_be="",
      where_str="").object
    print("input = phaser.Input%s()" % mode, file=s)
    input = input_class()
    self.input = input
    if data is not None :
      print("r = phaser.%s(data)" % run_data_class.__name__, file=s)
      r = run_data_class(data, phaser_out)
      handle_runtime_error(r, "setting up data inputs")
      print("input.setCELL6(r.getUnitCell())", file=s)
      print("input.setSPAC_HALL(r.getSpaceGroupHall())", file=s)
      input.setCELL6(r.getUnitCell())
      input.setSPAC_HALL(r.getSpaceGroupHall())
      if (params.density_modification.dm) :
        if (params.density_modification.solvent_fraction is None) :
          params.density_modification.solvent_fraction = \
            self.get_solvent_content(r)
      if mode.startswith("MR") and not mode in ["MR_PAK"] :
        print("input.setREFL(r.getMiller(),r.getFobs(),r.getSigFobs())",file=s)
        input.setREFL(r.getMiller(),r.getFobs(),r.getSigFobs())
      elif mode.startswith("EP") :
        print("input.setCRYS_MILLER(r.getMiller())", file=s)
        input.setCRYS_MILLER(r.getMiller())
        for xtal in params.phaser.crystal :
          xtalid = xtal.xtal_id
          for wave in xtal.dataset :
            waveid = wave.wave_id
            print("Fpos = r.getFpos('%s', '%s')" % (xtalid,waveid), file=s)
            print("Spos = r.getSIGFpos('%s', '%s')" % (xtalid,waveid), file=s)
            print("Ppos = r.getPpos('%s', '%s')" % (xtalid,waveid), file=s)
            print("Fneg = r.getFneg('%s', '%s')" % (xtalid,waveid), file=s)
            print("Sneg = r.getSIGFneg('%s', '%s')" % (xtalid,waveid), file=s)
            print("Pneg = r.getPneg('%s', '%s')" % (xtalid,waveid), file=s)
            Fpos = r.getFpos(xtalid,waveid)
            Spos = r.getSIGFpos(xtalid,waveid)
            Ppos = r.getPpos(xtalid,waveid)
            Fneg = r.getFneg(xtalid,waveid)
            Sneg = r.getSIGFneg(xtalid,waveid)
            Pneg = r.getPneg(xtalid,waveid)
            print(("input.addCRYS_ANOM_DATA('%s','%s',Fpos,Spos,Ppos,Fneg,"+
                   "Sneg,Pneg)") % (xtalid, waveid), file=s)
            input.addCRYS_ANOM_DATA(xtalid,waveid,Fpos,Spos,Ppos,Fneg,Sneg,Pneg)
            if [dataset.wavelength, dataset.energy].count(None) != 1 :
              raise Sorry(("Please enter EITHER a wavelength (Angstroms) "+
                  "OR energy (eV), but not both."))
            elif dataset.wavelength is not None : # TODO: more than 1?
              print("input.setWAVE(%.5f)" % dataset.wavelength, file=s)
              input.setWAVE(dataset.wavelength)
            elif dataset.energy is not None :
              wavelength = dataset.energy / 12398.0
              print("input.setWAVE(%.5f)" % wavelength, file=s)
              input.setWAVE(wavelength)
          if xtal.pdb_file is not None :
            print("input.setATOM_PDB('%s','%s')"%(xtalid,xtal.pdb_file),file=s)
            input.setATOM_PDB(xtalid, xtal.pdb_file)
          elif xtal.ha_file is not None :
            print("input.setATOM_HA('%s','%s')"%(xtalid,xtal.ha_file),file=s)
            input.setATOM_HA(xtalid, xtal.ha_file)
      del r
      del data
    #-------------------------------------------------------------------
    # special case: streamlined single-model search (like phenix.automr)
    have_composition = have_search = False
    if (match_modes(mode, index.get_phaser_modes("phaser.model"))) :
      if (len(params.phaser.model) > 0) :
        assert (len(params.phaser.ensemble) == 0)
        assert (len(params.phaser.composition.chain) == 0)
        assert (len(params.phaser.search) == 0)
        assert ((params.phaser.component_copies is not None) and
                (params.phaser.component_copies > 0))
        assert ((params.phaser.search_copies is not None) and
                (params.phaser.search_copies > 0))
        if (params.phaser.model_rmsd is not None) :
          assert (len(params.phaser.model) == len(params.phaser.model_rmsd))
          for pdb_file, rmsd in zip(params.phaser.model,
                                    params.phaser.model_rmsd) :
            input.addENSE_PDB_RMS("ense_1", pdb_file, rmsd)
        elif (params.phaser.model_identity is not None) :
          for pdb_file, identity in zip(params.phaser.model,
                                        params.phaser.model_identity) :
            input.addENSE_PDB_ID("ense_1", pdb_file, identity*0.01)
        if (params.phaser.seq_file is not None) :
          if (params.phaser.chain_type == "protein") :
            input.addCOMP_PROT_SEQ_NUM(params.phaser.seq_file,
              params.phaser.component_copies)
          else :
            input.addCOMP_NUCL_SEQ_NUM(params.phaser.seq_file,
              params.phaser.component_copies)
        elif (params.phaser.mol_weight is not None) :
          if (params.phaser.chain_type == "protein") :
            input.addCOMP_PROT_MW_NUM(params.phaser.mol_weight,
              params.phaser.component_copies)
          else :
            input.addCOMP_NUCL_MW_NUM(params.phaser.mol_weight,
              params.phaser.component_copies)
        input.addSEAR_ENSE_NUM("ense_1", params.phaser.search_copies)
        have_composition = have_search = True
    #-------------------------------------------------------------------
    # ENSEMBLES
    if match_modes(mode, index.get_phaser_modes("phaser.ensemble")) :
      for i, ensemble in enumerate(params.phaser.ensemble) :
        model_id = ensemble.model_id
        if model_id is None :
          raise Sorry("Missing model ID for ensemble #%d." % (i+1))
        if len(ensemble.coordinates) > 0 :
          if ensemble.map_hklin is not None :
            # XXX: does Phaser already check this?
            raise Sorry(("Please specify either PDB files or an MTZ file "+
              "for ensemble %s, but not both.") % model_id)
          for j, coords in enumerate(ensemble.coordinates) :
            if coords.pdb is None :
              raise Sorry(("Missing PDB file for ensemble %s coordinate set "+
                "%d.") % (model_id, j+1))
            if (not coords.read_variance_from_pdb_remarks) :
              if [coords.rmsd, coords.identity].count(None) != 1 :
                raise Sorry(("Please specify EITHER an RMSD OR fractional "+
                  "identity for PDB file %s.") % coords.pdb)
              elif coords.rmsd is not None :
                print("input.addENSE_PDB_RMS('%s','%s',%s)" % (model_id,
                  coords.pdb, str(coords.rmsd)), file=s)
                input.addENSE_PDB_RMS(model_id, coords.pdb, coords.rmsd)
              else :
                print("input.addENSE_PDB_ID('%s','%s',%s)" % (model_id,
                  coords.pdb, str(coords.identity)), file=s)
                input.addENSE_PDB_ID(model_id, coords.pdb, coords.identity)
            else :
              print("input.addENSE_PDB_CARD('%s', '%s', True)" % (model_id,
                coords.pdb), file=s)
              input.addENSE_PDB_CARD(model_id, coords.pdb, True)
            # TODO: addENSE_CARD
          if ensemble.disable_check :
            input.setENSE_DISA_CHEC(model_id, True)
          if ensemble.use_hetatm :
            input.setENSE_HETA(model_id, True)
        elif ensemble.map_hklin is not None :
          ens = ensemble
          if ens.map_labels is None :
            raise Sorry("You must specify labels for map data for %s." %
                        model_id)
          if ((ens.map_extent is None or len(ens.map_extent) != 3) or
                (ens.map_centre is None or len(ens.map_centre) != 3)) :
            # TODO: confirm this
            raise Sorry(("You must specify the extent and center of the map "+
              "used for ensemble %s; this must both take the form of three "+
              "decimal numbers between zero and one.") % model_id)
          if (None in [ens.map_protein_mw, ens.map_nucleic_mw] or
              (ens.map_protein_mw == 0 and ens.map_nucleic_mw == 0)) :
            raise Sorry(("Please specify protein and/or nucleic acid masses "+
              "for ensemble %s.") % model_id)
          (f_label, phi_label) = ens.map_labels.split(",")
          print(("input.addENSE_MAP('%s', '%s', '%s', '%s', %s, %.3f, %s, "+
                 "%.1f, %.1f") % (model_id, ens.map_hklin, f_label, phi_label,
                  str(ens.map_extent), ens.map_rms, str(ens.map_centre),
                  ens.map_protein_mw, ens.map_nucleic_mw), file=s)
          input.addENSE_MAP(model_id, ens.map_hklin, f_label, phi_label,
            ens.map_extent, ens.map_rms, ens.map_centre, ens.map_protein_mw,
            ens.map_nucleic_mw)
        if ensemble.fractional_comp is not None :
          # XXX: should I check if every ensemble has this?
          print("input.setENSE_FRAC('%s',%s)" % (ensemble.model_id,
            ensemble.fractional_comp), file=s)
          input.addCOMP_ENSE_FRAC(ensemble.model_id, ensemble.fractional_comp)
          have_composition = True
        if ensemble.solution_at_origin :
          print("input.addSOLU_ORIG_ENSE('%s')" % ensemble.model_id, file=s)
          input.addSOLU_ORIG_ENSE(ensemble.model_id)
    #-------------------------------------------------------------------
    # SEARCH PROCEDURE
    if match_modes(mode, index.get_phaser_modes("phaser.search")) :
      if (not have_search) and (len(params.phaser.search) == 0) :
        raise Sorry("No search procedure defined.")
      for i, component in enumerate(params.phaser.search) :
        if component.ensembles is None :
          raise Sorry("No model IDs supplied for search component %d." % (i+1))
        elif len(component.ensembles) == 1:
          print("input.addSEAR_ENSE_NUM('%s',%s)" % (component.ensembles[0],
              str(component.copies)), file=s)
          input.addSEAR_ENSE_NUM(component.ensembles[0], component.copies)
        else :
          print("input.addSEAR_ENSE_OR_ENSE_NUM(%s,%s)" %
            (str(component.ensembles), str(component.copies)), file=s)
          input.addSEAR_ENSE_OR_ENSE_NUM(component.ensembles, component.copies)
    if match_modes(mode, index.get_phaser_modes("phaser.composition")) :
      if (not have_composition) :
        self.set_composition(input)
    #-------------------------------------------------------------------
    # OLD SOLUTION
    if params.phaser.solution is not None :
      assert os.path.isfile(params.phaser.solution)
      old_result = easy_pickle.load(params.phaser.solution)
      # XXX preserve backwards compatibility
      if hasattr(old_result, "get_phaser_object") :
        old_result = old_result.get_phaser_object()
      if type(old_result).__name__ == "ResultMR_RF" :
        solution = old_result.getDotRlist()
      else :
        solution = old_result.getDotSol()
      if len(solution) == 0 :
        raise Sorry("The specified previous solutions file is empty.")
      if mode in ["MR_FTF", "MR_BTF"] :
        if not solution.is_rlist() :
          raise Sorry("Input solutions for the translation functions must be "+
            "the output from a rotation function.")
      else :
        if solution.is_rlist() :
          raise Sorry("The specified solutions file is from a rotation "+
            "function; both rotation and translation are required when "+
            "supplying an existing solution in this mode.")
      input.setSOLU(solution)
    #-------------------------------------------------------------------
    # KEYWORD CONVERSION
    keywords = params.phaser.keywords
    keyword_phil = index.get_scope_by_name("phaser.keywords")
    convert_keywords(input, keyword_phil, mode, s)
    index.update_from_python(params)
    self.params = params
    print("result = phaser.run%s(input)" % mode, file=s)
    self.script = s.getvalue()
    self.input = input
    self.mode = mode

  def get_solvent_content (self, data_result) :
    params = self.params
    if (params.phaser.composition.solvent is not None) :
      return params.phaser.composition.solvent
    else :
      input = phaser.InputCCA()
      input.setSPAC_HALL(data_result.getSpaceGroupHall())
      input.setCELL6(data_result.getUnitCell())
      self.set_composition(input, script_out=null_out())
      input.setMUTE(True)
      r = phaser.runCCA(input)
      if (not r.Success()) :
        raise Sorry("Could not determine solvent fraction automatically.  If "+
          "you want to use density modification, we recommend you explicitly "+
          "set the solvent fraction (available under \"Density modification\""+
          "in the Settings menu).")
      vm = r.getBestVM()
      return (1.0 - (1.23 / vm))

  #-------------------------------------------------------------------
  # COMPOSITION
  def set_composition (self, input=None, script_out=None) :
    params = self.params
    if (script_out is None) :
      script_out = self._script_out
    s = script_out
    if (input is None) :
      input = self.input
    if params.phaser.composition.solvent is not None :
      if len(params.phaser.composition.chain) > 0 :
        raise Sorry("Please specify EITHER fractional solvent content OR "+
          "individual components, but not both.")
      print("input.setCOMP_SOLV()", file=s)
      input.setCOMP_SOLV()
      # XXX: is there currently no way to set this?
    elif len(params.phaser.composition.chain) > 0 :
      print("input.setCOMP_ASU()", file=s)
      input.setCOMP_ASU()
      for i, chain in enumerate(params.phaser.composition.chain) :
        mass_attr = ["sequence_file","sequence","nres","mw"]
        masses = [ getattr(chain, name) for name in mass_attr]
        if masses.count(None) != 3 :
          raise Sorry(("Please specify ONE source of mass information for "+
            "chain #%d: sequence file, sequence string, # of residues, or "+
            "molecular weight (Da).") % (i+1))
        if chain.num is None :
          raise Sorry("Please specify a number of copies for component #%d."%
            (i+1))
        if chain.chain_type == "protein" :
          method_name = "addCOMP_PROT_"
        else :
          method_name = "addCOMP_NUCL_"
        if chain.sequence_file is not None :
          method_name += "SEQ_NUM"
          mass_value = chain.sequence_file
        elif chain.sequence is not None :
          method_name += "STR_NUM"
          mass_value = chain.sequence
        elif chain.nres is not None :
          method_name += "NRES_NUM"
          mass_value = chain.nres
        else : # chain.mw is not None
          method_name += "MW_NUM"
          mass_value = chain.mw
        method = getattr(input, method_name)
        if isinstance(mass_value, str) :
          out_value = "'%s'" % mass_value
        else :
          out_value = str(mass_value)
        print("input.%s(%s,%s)"%(method_name,out_value,str(chain.num)),
              file=s)
        method(mass_value, chain.num)
    else :
      print("input.setCOMP_AVER()", file=s)
      input.setCOMP_AVER()
    for i, atom in enumerate(params.phaser.composition.atom) :
      if atom.element is None or atom.num is None :
        raise Sorry(("Please specify an element and number of atoms "+
          "for heavy atom %d.") % (i+1))
      print("input.addCOMP_ATOM_NUM('%s', %d)" % (atom.element,atom.num),
        file=s)
      try :
        input.addCOMP_ATOM_NUM(atom.element, atom.num)
      except ValueError, e :
        raise Sorry(str(e))

  def run (self, out=sys.stdout, use_callback=True) :
    run_class = import_python_object(import_path="phaser.run%s" % self.mode,
      error_prefix="",
      target_must_be="",
      where_str="").object
    phaser_out = phaser.Output()
    if use_callback :
      callback = callbacks.manager()
    else :
      callback = oop.null()
    out_wrapper = multi_out()
    log_file = open(self.params.phaser.keywords.general.root + ".log", "w")
    out_wrapper.register("Primary output stream", out)
    out_wrapper.register("Phaser logfile", log_file)
    phaser_out.setPackagePhenix(out_wrapper)
    if use_callback :
      phaser_out.setPhenixCallback(callback)
    try :
      r = run_class(self.input, phaser_out)
      callback.flush_summary()
      handle_runtime_error(r, "running main function")
      self.result = r
    finally :
      log_file.close()
    if (self.params.density_modification.dm) :
      from phaser.phenix_interface import dm
      mtz_files = []
      if (self.params.phaser.mode.startswith("MR")) :
        mtz_files = [self.params.phaser.keywords.general.root + ".1.mtz"]
      elif (self.params.phaser.mode.startswith("EP")) :
        mtz_files = [self.params.phaser.keywords.general.root + ".1.mtz",
                     self.params.phaser.keywords.general.root + ".1.hand.mtz" ]
      callback.call_back(message="phaser_start", data=None)
      for file_name in mtz_files :
        if os.path.isfile(file_name) :
          callback.call_back(message="MODE",
            data="density modification with %s" % file_name)
          dm.run_density_modification(
            phases_file=file_name,
            f_obs=self.input_data,
            params=self.params.density_modification,
            log=out)
      callback.call_back(message="phaser_end", data=None)

def handle_runtime_error (result, procedure) :
  if not result.Success() :
    raise Sorry("ERROR %s (%s): %s" % (procedure, result.ErrorName(),
                result.ErrorMessage()))

########################################################################
def run (args, out=sys.stdout, phil_str=None, working_phil=None) :
  try :
    from phenix.utilities import citations
  except ImportError :
    citations = oop.null()
  import iotbx.phil
  citations.clear_citations()
  sources = []
  if phil_str is not None :
    input_phil = iotbx.phil.parse(phil_str)
    sources.append(phil_str)
  if working_phil is not None :
    sources.append(working_phil)
  master_phil = master_params()
  cmdline = iotbx.phil.process_command_line_with_files(
    args=args,
    master_phil=master_phil,
    pdb_file_def="phaser.model",
    reflection_file_def="phaser.hklin",
    seq_file_def="phaser.seq_file")
  working_phil = cmdline.work
  index = phil_interface.index(master_phil, working_phil,
    parse=iotbx.phil.parse)
  params = index.get_python_object()
  if params.show_defaults :
    master_phil.show(out=out)
    return True
  # if labin is None and this was started from the command line, figure out
  # default Miller array to use
  if ((params.phaser.hklin is not None) and (params.phaser.labin is None) and
      (not params.phaser.run_control.gui_run)) :
    from iotbx import file_reader
    hklin = file_reader.any_file(params.phaser.hklin, force_type="hkl")
    hklin.check_file_type("hkl")
    f_label = i_label = None
    for array in hklin.file_server.miller_arrays :
      label = array.info().label_string()
      if (label in ["F,SIGF", "FOBS,SIGFOBS", "F(+),SIGF(+),F(-),SIGF(-)"]) :
        f_label = label
        break
      elif ((array.is_xray_amplitude_array()) and
          (not label.startswith("FC")) and (f_label is None)) :
        f_label = label
      elif (array.is_xray_intensity_array()) and (i_label is None) :
        i_label = label
    if (f_label is not None) :
      params.phaser.labin = f_label
    elif (i_label is not None) :
      params.phaser.labin = i_label
    index.update_from_python(params)
  #if params.dry_run :
  #  index.update("phaser.keywords.mute = True")
  tmp_out = out
  if not params.verbose :
    tmp_out = StringIO()
  run_phaser = phaser_parameter_interpreter(index, os.getcwd(), out=tmp_out)
  if params.show_script :
    out.write(run_phaser.script)
    return True
  elif params.dry_run :
    run_phaser.params.dry_run = False
    index.update_from_python(run_phaser.params)
    index.get_diff().show(out=out)
    return True
  if params.verbose :
    index.get_diff().show(out=out)
  run_phaser.run(out=out)
  params = run_phaser.params # XXX this has updated solvent_fraction
  citations.add_citation("phaser", source="Phaser")
  citations.show(out=out, source="Phaser")
  result = phaser_result(
    phaser_object=run_phaser.result,
    params=params,
    run_dir=os.getcwd())
  return result