from __future__ import division
from libtbx.utils import Sorry
import os.path
import sys

def run_density_modification (phases_file, f_obs, params, log=None) :
  """
  Run simple solvent-flipping in mmtbx.density_modification on Phaser output.
  f_obs will be None if molecular replacement was run, since the output MTZ
  file will include anisotropy-corrected amplitudes.  Unfortunatley the EP
  mode does not currently output these data.
  """
  if (log is None) :
    log = sys.stdout
  from mmtbx import density_modification
  from iotbx import file_reader
  import iotbx.mtz
  phaser_out = file_reader.any_file(phases_file, force_type="hkl")
  phaser_out.assert_file_type("hkl")
  miller_arrays = phaser_out.file_server.miller_arrays
  hl_coeffs = None
  for array in miller_arrays :
    labels = array.info().label_string()
    if (labels == "F,SIGF") :
      f_obs = array
    elif (labels == "HLA,HLB,HLC,HLD") :
      hl_coeffs = array
  if (f_obs is None) or (hl_coeffs is None) :
    raise Sorry("Missing F-obs or HL coefficients array.")
  if (f_obs.anomalous_flag()) :
    f_obs = f_obs.average_bijvoet_mates()
  if (hl_coeffs.anomalous_flag()) :
    hl_coeffs = hl_coeffs.average_bijvoet_mates()
  f_obs = f_obs.common_set(hl_coeffs)
  hl_coeffs = hl_coeffs.common_set(f_obs)
  assert (f_obs.indices() == hl_coeffs.indices())
  run_dm = density_modification.density_modification(
    f_obs=f_obs,
    hl_coeffs_start=hl_coeffs,
    params=params,
    log=log)
  map_coeffs = run_dm.map_coeffs_in_original_setting
  mtz_dataset = map_coeffs.as_mtz_dataset(
    column_root_label="FWT",
    label_decorator=iotbx.mtz.ccp4_label_decorator())
  mtz_object = mtz_dataset.mtz_object()
  mtz_out = os.path.splitext(phases_file)[0] + "_dm.mtz"
  mtz_object.write(file_name=mtz_out)
  return mtz_out