from __future__ import division
from libtbx.utils import Sorry, null_out
import sys

def find_anomalous_scatterers (
    fmodel,
    pdb_hierarchy,
    wavelength,
    prefix = "phaser_AX",
    verbose = True,
    log = None,
    n_cycles=None) :
  """
  Runs LLG completion (searching for 'AX') to collect a list of anomalous
  scatterers around the unit cell.
  """
  if (wavelength is None) :
    raise Sorry("Please specify the X-ray wavelength.")
  if (log is None) : log = sys.stdout
  from phaser.phenix_adaptors import sad_target
  import phaser
  assert (fmodel.f_obs().anomalous_flag())
  f_obs = fmodel.f_obs()
  r_free_flags = fmodel.r_free_flags()
  f_obs = f_obs.select(f_obs.data() > 0)
  f_obs, r_free_flags = f_obs.common_sets(other=r_free_flags)
  adaptor = sad_target.data_adaptor(
    f_obs = f_obs,
    r_free_flags = r_free_flags,
    verbose = True)
  phaser_input = adaptor.set_ep_auto_data(
    xray_structure = fmodel.xray_structure,
    pdb_hierarchy = pdb_hierarchy)
  phaser_input.setLLGC_COMP(True)
  phaser_input.addLLGC_SCAT("AX")
  phaser_input.setWAVE(wavelength)
  phaser_input.setROOT(prefix)
  phaser_input.setMUTE(not verbose)
  phaser_input.setHKLO(False)
  phaser_input.setXYZO(False)
  if (n_cycles is not None) :
    assert (n_cycles >= 1)
    phaser_input.setLLGC_NCYC(n_cycles)
  out = phaser.Output()
  if (not verbose) :
    out.setPackagePhenix(null_out())
  else :
    out.setPackagePhenix(log)
  result = phaser.runEP_AUTO(phaser_input, out)
  # FIXME attempting to call result.getIotbx() will cause a segfault if the
  # job was unsuccessful
  if (result.Failure()) :
    return None
  # XXX will these already be arranged to be closest to the partial model?
  return result.getIotbx() # pdb hierarchy object