from __future__ import division
from cctbx.array_family import flex
import sys
from libtbx import adopt_init_args
from cctbx import maptbx
import libtbx.phil
from libtbx.str_utils import format_value
from mmtbx.refinement import print_statistics

peak_search_params_str = """
peak_search_level = 1
  .type=int
max_peaks = 0
  .type=int
  .short_caption=Maximum peaks
interpolate = True
  .type=bool
min_distance_sym_equiv = None
  .type=float
  .short_caption=Minimum distance between symmetry-equivalent atoms
general_positions_only = False
  .type=bool
min_cross_distance = 1.8
  .type=float
  .short_caption=Minimum cross distance
min_cubicle_edge = 5
  .type=float
  .short_caption=Minimum edge length of cubicles used for \
    fast neighbor search
  .expert_level=2
"""

master_params = libtbx.phil.parse("""\
  use_sigma_scaled_maps = True
    .type=bool
    .help = Default is sigma scaled map, map in absolute scale is used \
            otherwise.
  resolution_factor = 1./4.
    .type=float
  map_next_to_model
    .expert_level=2
    .style = noauto
  {
    min_model_peak_dist = 1.8
      .type=float
      .short_caption = Minimum distance from model
    max_model_peak_dist = 6.0
      .type=float
      .short_caption = Maximum distance from model
    min_peak_peak_dist = 1.8
      .type=float
      .short_caption = Minimum distance between peaks
    use_hydrogens = False
      .type = bool
      .short_caption = Use hydrogens
  }
  max_number_of_peaks = None
    .type=int
    .expert_level=1
  peak_search
    .expert_level=1
    .short_caption = Search settings
    .style = box auto_align
  {
    %s
  }
""" % peak_search_params_str)

class peaks_holder(object):
  def __init__(self, heights, sites, iseqs_of_closest_atoms = None):
    assert (len(heights) == len(sites))
    self.heights = heights
    self.sites = sites
    self.iseqs_of_closest_atoms = iseqs_of_closest_atoms

  def filter_by_secondary_map(self, map, min_value):
    n_deleted = 0
    k = 0
    while (k < len(self.sites)):
      map_value = map.tricubic_interpolation(self.sites[k])
      if (map_value < min_value):
        del self.sites[k]
        del self.heights[k]
        n_deleted += 1
      else :
        k += 1
    return n_deleted

  def sort(self, reverse=False):
    from scitbx.array_family import flex
    selection = flex.sort_permutation(self.heights, reverse=reverse)
    heights_sorted = self.heights.select(selection)
    sites_sorted = self.sites.select(selection)
    iseqs_sorted = None
    if (self.iseqs_of_closest_atoms is not None):
      iseqs_sorted = self.iseqs_of_closest_atoms.select(selection)
    self.heights = heights_sorted
    self.sites = sites_sorted
    self.iseqs_of_closest_atoms = iseqs_sorted

class manager(object):
  def __init__(self,
               fmodel,
               map_type,
               map_cutoff,
               params = None,
               log = None,
               use_all_data = True,
               silent = False,
               map_coeffs=None):
    adopt_init_args(self, locals())
    assert (map_type is not None) or (map_coeffs is not None)
    self.mapped = False
    self.peaks_ = None
    if(self.log is None): self.log = sys.stdout
    if(self.params is None): self.params = master_params.extract()
    if (map_coeffs is not None):
      fft_map = map_coeffs.fft_map(
        resolution_factor=self.params.resolution_factor,
        symmetry_flags=maptbx.use_space_group_symmetry)
      if(self.params.use_sigma_scaled_maps):
        fft_map.apply_sigma_scaling()
        map_units = "sigma"
      else:
        fft_map.apply_volume_scaling()
        map_units = "e/A**3"
      fft_map_data = fft_map.real_map_unpadded()
    else:
      fft_map = self.fmodel.electron_density_map().\
          fft_map(resolution_factor = self.params.resolution_factor,
                  symmetry_flags    = maptbx.use_space_group_symmetry,
                  map_type          = self.map_type,
                  use_all_data      = use_all_data)
      if(self.params.use_sigma_scaled_maps):
        fft_map.apply_sigma_scaling()
        map_units = "sigma"
      else:
        fft_map.apply_volume_scaling()
        map_units = "e/A**3"
      fft_map_data = fft_map.real_map_unpadded()
    crystal_gridding_tags = fft_map.tags()
    max_number_of_peaks = self.params.max_number_of_peaks
    if(self.params.max_number_of_peaks is None):
      max_number_of_peaks = self.fmodel.xray_structure.scatterers().size() * 5
    negative = False
    if(self.map_cutoff < 0):
      self.map_cutoff *= -1
      negative = True
      fft_map_data *= -1
    min_distance_sym_equiv = self.params.peak_search.min_distance_sym_equiv
    if(min_distance_sym_equiv is None):
      min_distance_sym_equiv = \
        self.fmodel.xray_structure.min_distance_sym_equiv()
    peak_search_parameters = maptbx.peak_search_parameters(
      peak_search_level      = self.params.peak_search.peak_search_level,
      max_peaks              = self.params.peak_search.max_peaks,
      peak_cutoff            = self.map_cutoff,
      interpolate            = self.params.peak_search.interpolate,
      min_distance_sym_equiv = min_distance_sym_equiv,
      general_positions_only = self.params.peak_search.general_positions_only,
      min_cross_distance     = self.params.peak_search.min_cross_distance,
      min_cubicle_edge       = self.params.peak_search.min_cubicle_edge)
    if(self.fmodel.r_work() > 0.00001 and self.fmodel.r_free() > 0.00001):
      cluster_analysis = crystal_gridding_tags.peak_search(
        parameters = peak_search_parameters,
        map = fft_map_data).all(max_clusters = max_number_of_peaks)
      heights = cluster_analysis.heights()
      if(negative):
        heights *= -1.
      self.peaks_ = peaks_holder(heights = heights,
                                 sites   = cluster_analysis.sites())
      if(not self.silent):
        print >>self.log,"Number of peaks found at %s map (map cutoff=%s %s)= %s"%(
          self.map_type, format_value("%-5.2f", self.map_cutoff).strip(),
          map_units, format_value("%-12d", self.peaks_.sites.size()))

  def peaks(self):
    return self.peaks_

  def peaks_mapped(self):
    if(self.peaks_ is None): return None
    assert self.mapped == False
    max_dist = self.params.map_next_to_model.max_model_peak_dist
    min_dist = self.params.map_next_to_model.min_model_peak_dist
    if (min_dist is None):
      min_dist = 0.
    if (max_dist is None):
      max_dist = float(sys.maxint)
    xray_structure = self.fmodel.xray_structure.deep_copy_scatterers()
    use_selection = None
    if(not self.params.map_next_to_model.use_hydrogens):
      use_selection = ~xray_structure.hd_selection()
    initial_number_of_sites = self.peaks_.sites.size()
    if(not self.silent):
      print >> self.log, "Filter by distance & map next to the model:"
    result = xray_structure.closest_distances(sites_frac = self.peaks_.sites,
      distance_cutoff = max_dist, use_selection = use_selection)
    smallest_distances_sq = result.smallest_distances_sq
    smallest_distances = result.smallest_distances
    in_box = smallest_distances_sq > 0
    not_too_far = smallest_distances_sq <= max_dist**2
    not_too_close = smallest_distances_sq >= min_dist**2
    selection = (not_too_far & not_too_close & in_box)
    iseqs_of_closest_atoms = result.i_seqs.select(selection)
    peaks = peaks_holder(
      heights                = self.peaks_.heights.select(selection),
      sites                  = result.sites_frac.select(selection),
      iseqs_of_closest_atoms = iseqs_of_closest_atoms)
    sd = flex.sqrt(smallest_distances_sq.select(in_box))
    d_min = flex.min_default(sd, 0)
    d_max = flex.max_default(sd, 0)
    if(not self.silent):
      print >> self.log,"   mapped sites are within: %5.3f - %5.3f"%(d_min,d_max)
      print >> self.log, "   number of sites selected in [dist_min=%5.2f, " \
        "dist_max=%5.2f]: %d from: %d" % (min_dist, max_dist, peaks.sites.size(),
        initial_number_of_sites)
    smallest_distances = flex.sqrt(smallest_distances_sq.select(selection))
    d_min = flex.min_default(smallest_distances, 0)
    d_max = flex.max_default(smallest_distances, 0)
    if(not self.silent):
      print >> self.log,"   mapped sites are within: %5.3f - %5.3f"%(d_min,d_max)
    self.mapped = True
    self.peaks_ = peaks
    return peaks

  def show_mapped(self, pdb_atoms):
    if(self.peaks_ is None): return None
    peaks = self.peaks()
    if(peaks.iseqs_of_closest_atoms is None):
      raise RuntimeError("iseqs_of_closest_atoms is None")
    scatterers = self.fmodel.xray_structure.scatterers()
    assert scatterers.size() == pdb_atoms.size()
    assert peaks.sites.size() == peaks.heights.size()
    assert peaks.heights.size() == peaks.iseqs_of_closest_atoms.size()
    print >> self.log
    dist = self.fmodel.xray_structure.unit_cell().distance
    for i in flex.sort_permutation(data=peaks.iseqs_of_closest_atoms):
      s = peaks.sites[i]
      h = peaks.heights[i]
      i_seq = peaks.iseqs_of_closest_atoms[i]
      sc = scatterers[i_seq]
      d = dist(s, sc.site)
      element = sc.element_symbol()
      print >> self.log, "peak= %8.3f closest distance to %s = %8.3f" % (
        h, pdb_atoms[i_seq].id_str(), d)
      assert d <= self.params.map_next_to_model.max_model_peak_dist
      assert d >= self.params.map_next_to_model.min_model_peak_dist

def show_highest_peaks_and_deepest_holes(fmodel,
                                         pdb_atoms,
                                         map_type,
                                         map_cutoff_plus,
                                         map_cutoff_minus,
                                         log = None):
  if(log is None): log = sys.stdout
  print_statistics.make_header(
    "residual map %s: highest peaks and deepst holes"%map_type, out = log)
  fp_params = master_params.extract()
  fp_params.map_next_to_model.min_model_peak_dist = 0.7
  fp_params.map_next_to_model.min_peak_peak_dist = 0.7
  fp_params.peak_search.min_cross_distance = 0.7
  fp_params.map_next_to_model.use_hydrogens = True
  for par in [(map_cutoff_plus,"peaks"), (map_cutoff_minus,"holes")]:
    print_statistics.make_sub_header(par[1], out = log)
    result = manager(fmodel     = fmodel,
                     map_type   = "mFobs-DFmodel",
                     map_cutoff = par[0],
                     params     = fp_params,
                     log        = log)
    result.peaks_mapped()
    result.show_mapped(pdb_atoms = pdb_atoms)
