text: from __future__ import absolute_import, division, print_function
  # LIBTBX_SET_DISPATCHER_NAME cxi.pixel_histograms
  # LIBTBX_PRE_DISPATCHER_INCLUDE_SH export PHENIX_GUI_ENVIRONMENT=1
  # LIBTBX_PRE_DISPATCHER_INCLUDE_SH export BOOST_ADAPTBX_FPE_DEFAULT=1
  
  import sys
  
  from libtbx import easy_pickle
  from libtbx.option_parser import option_parser
  from scitbx.array_family import flex
  from scitbx.math import curve_fitting
  
  from xfel.cxi.cspad_ana import cspad_tbx
  from six.moves import range
  
  def run(args):
    assert len(args) > 0
    command_line = (option_parser()
                    .option("--roi",
                            type="string",
                            help="Region of interest for summing up histograms"
                            "from neighbouring pixels.")
                    .option("--log_scale",
                            action="store_true",
                            default=False,
                            help="Draw y-axis on a log scale.")
                    .option("--normalise",
                            action="store_true",
                            default=False,
                            help="Normalise by number of member images.")
                    .option("--save",
                            action="store_true",
                            default=False,
                            help="Save each plot as a png.")
                    .option("--start",
                            type="string",
                            help="Starting pixel coordinates")
                    .option("--fit_gaussians",
                            action="store_true",
                            default=False,
                            help="Fit gaussians to the peaks.")
                    .option("--n_gaussians",
                            type="int",
                            default=2,
                            help="Number of gaussians to fit.")
                    .option("--estimated_gain",
                            type="float",
                            default=30,
                            help="The approximate position of the one photon peak.")
                    ).process(args=args)
    log_scale = command_line.options.log_scale
    fit_gaussians = command_line.options.fit_gaussians
    roi = cspad_tbx.getOptROI(command_line.options.roi)
    normalise = command_line.options.normalise
    save_image = command_line.options.save
    starting_pixel = command_line.options.start
    n_gaussians = command_line.options.n_gaussians
    estimated_gain = command_line.options.estimated_gain
    if starting_pixel is not None:
      starting_pixel = eval(starting_pixel)
      assert isinstance(starting_pixel, tuple)
    args = command_line.args
  
    path = args[0]
    window_title = path
    d = easy_pickle.load(path)
    args = args[1:]
    pixels = None
    if len(args) > 0:
      pixels = [eval(arg) for arg in args]
      for pixel in pixels:
        assert isinstance(pixel, tuple)
        assert len(pixel) == 2
    if roi is not None:
  
      summed_hist = {}
      for i in range(roi[2], roi[3]):
        for j in range(roi[0], roi[1]):
          if (i,j) not in summed_hist:
            summed_hist[(0,0)] = d[(i,j)]
          else:
            summed_hist[(0,0)].update(d[(i,j)])
      d = summed_hist
  
    #if roi is not None:
  
      #summed_hist = None
      #for i in range(roi[2], roi[3]):
        #for j in range(roi[0], roi[1]):
          #if summed_hist is None:
            #summed_hist = d[(i,j)]
          #else:
            #summed_hist.update(d[(i,j)])
  
      #title = str(roi)
      #plot(hist, window_title=window_title, title=title,log_scale=log_scale,
           #normalise=normalise, save_image=save_image, fit_gaussians=fit_gaussians)
      #return
  
  
    histograms = pixel_histograms(d, estimated_gain=estimated_gain)
    histograms.plot(
       pixels=pixels, starting_pixel=starting_pixel, fit_gaussians=fit_gaussians,
      n_gaussians=n_gaussians, window_title=window_title, log_scale=log_scale,
      save_image=save_image)
  
  
  class pixel_histograms(object):
  
    def __init__(self, histograms, estimated_gain=30):
      self.histograms = histograms
      self.estimated_gain = estimated_gain
  
    def plot(self, pixels=None, starting_pixel=None, fit_gaussians=True,
             n_gaussians=2, window_title=None, log_scale=False, save_image=False):
      from matplotlib import pyplot
      normalise=False # XXX
      assert [pixels, starting_pixel].count(None) > 0
      print("n_images:", flex.sum(self.histograms.values()[0].slots()))
      print("n_pixels:", len(self.histograms))
      if pixels is None:
        pixels = sorted(self.histograms.keys())
        if starting_pixel is not None:
          pixels = pixels[pixels.index(starting_pixel):]
      for pixel in pixels:
        hist = self.histograms[pixel]
        print(pixel)
        title = str(pixel)
        if fit_gaussians:
          pyplot.subplot(211)
        self.plot_one_histogram(
          hist, window_title=window_title, title=title,log_scale=log_scale,
          normalise=normalise, save_image=save_image)
        fontsize = 18
        pyplot.ylabel("Counts", fontsize=fontsize)
        pyplot.xlabel("ADUs", fontsize=fontsize)
        if fit_gaussians:
          self.plot_gaussians(pixel, n_gaussians=n_gaussians, log_scale=log_scale)
        pyplot.ylabel("Counts", fontsize=fontsize)
        pyplot.xlabel("ADUs", fontsize=fontsize)
        #axes = pyplot.axes()
        #for tick in axes.xaxis.get_ticklabels():
          #tick.set_fontsize(fontsize)
        #for tick in axes.yaxis.get_ticklabels():
          #tick.set_fontsize(fontsize)
  
        if save_image:
          pyplot.savefig("%s.png" %title)
        else:
          pyplot.show()
  
    def plot_gaussians(self, pixel, n_gaussians=2, log_scale=False):
      from matplotlib import pyplot
      if log_scale:
        pyplot.ylim(ymin=0.1)
      hist = self.histograms[pixel]
      #pyplot.ylim(ymax=flex.max(hist.slots()))
      gaussians = self.fit_one_histogram(pixel, n_gaussians=n_gaussians)
      x = hist.slot_centers()
      y_calc = flex.double(x.size(), 0)
      for g in gaussians:
        print(g.params)
        y = g(x)
        y_calc += y
        pyplot.plot(x, y, linewidth=2)
  
      #print "Peak height ratio: %.2f" %(gaussians[0].params[0]/gaussians[1].params[0])
      pyplot.plot(x, y_calc)
      # Plot the fit residuals
      xlim = pyplot.xlim() # store for reuse below
      pyplot.subplot(212)
      residual = hist.slots().as_double() - y_calc
      pyplot.plot(x, residual, linewidth=2)
      pyplot.xlim(xlim)
  
    def plot_one_histogram(self, histogram,
                           window_title=None, title=None,
                           log_scale=False, normalise=False, save_image=False):
      from matplotlib import pyplot
      slots = histogram.slots().as_double()
      if normalise:
        normalisation = (flex.sum(slots) + histogram.n_out_of_slot_range()) / 1e5
        print("normalising by factor: ", normalisation)
        slots /= normalisation
      bins, data = hist_outline(histogram)
      if log_scale:
        data.set_selected(data == 0, 0.1) # otherwise lines don't get drawn when we have some empty bins
        pyplot.yscale("log")
      pyplot.plot(bins, data, '-k', linewidth=2)
      pyplot.suptitle(title)
      data_min = min([slot.low_cutoff for slot in histogram.slot_infos() if slot.n > 0])
      data_max = max([slot.low_cutoff for slot in histogram.slot_infos() if slot.n > 0])
      pyplot.xlim(data_min, data_max+histogram.slot_width())
  
    def fit_one_histogram(self, pixel, n_gaussians=2):
      histogram = self.histograms[pixel]
      fitted_gaussians = []
  
      slot_centers = histogram.slot_centers()
      slots = histogram.slots().as_double()
  
      zero_peak_gaussian = None
      for i in range(n_gaussians):
        if i == 0:
          lower_threshold = -1000
          upper_threshold = 0.4 * self.estimated_gain
          mean = 0
          fit = self.single_peak_fit(histogram, lower_threshold, upper_threshold, mean,
                                     zero_peak_gaussian=zero_peak_gaussian)
          hist_max = flex.max(histogram.slots())
          if abs(fit.functions[0].params[0] - hist_max)/hist_max > 0.1:
            upper_threshold = 0.3 * self.estimated_gain
            fit = self.single_peak_fit(histogram, lower_threshold, upper_threshold, mean,
                                       zero_peak_gaussian=zero_peak_gaussian)
        else:
          y_obs = histogram.slots().as_double()
          x = histogram.slot_centers()
          y_calc = flex.double(y_obs.size(), 0)
          for g in fitted_gaussians:
            y_calc += g(x)
          residual = y_obs - y_calc
          # triangular smoothing of residual to find peak position
          residual = sliding_average(residual)
          residual = sliding_average(residual)
          for n in (4, 5, 6, 7, 8):
          #for n in (5, 6, 7, 8):
            # we assume that the peaks are separated by at least n sigma
            n_sigma = abs(n * fitted_gaussians[0].params[2])
            slot_i = histogram.get_i_slot(fitted_gaussians[i-1].params[1]+n_sigma)
            max_slot_i = flex.max_index(residual[slot_i:]) + slot_i
            mean = slot_centers[max_slot_i]
            lower_threshold = mean - 0.3 * (mean - fitted_gaussians[0].params[1])
            upper_threshold = mean + 0.4 * (mean - fitted_gaussians[0].params[1])
            #print lower_threshold, mean, upper_threshold
            #zero_peak_gaussian = None
            fit = self.single_peak_fit(histogram, lower_threshold, upper_threshold, mean,
                                       zero_peak_gaussian=zero_peak_gaussian)
            if (fit.functions[0].params[1] > fitted_gaussians[-1].params[1]
                and fit.functions[0].sigma > 0.5 * fitted_gaussians[-1].sigma
                and (fit.functions[0].params[1] - fitted_gaussians[-1].params[1]) > n_sigma):
              break
        fitted_gaussians += fit.functions
        if i == 0: zero_peak_gaussian = fit.functions[0]
  
  
      if len(fitted_gaussians) > 1:
        try:
          check_pixel_histogram_fit(histogram, fitted_gaussians)
        except PixelFitError as e:
          print("PixelFitError:", str(e))
        gain = fitted_gaussians[1].params[1] - fitted_gaussians[0].params[1]
        print("gain: %s" %gain)
        zero_peak = fitted_gaussians[0].params[1]
        photon_threshold = 2/3
        n_single_photons = flex.sum(
          histogram.slots()[histogram.get_i_slot(photon_threshold * gain + zero_peak):])
        n_double_photons = flex.sum(
          histogram.slots()[histogram.get_i_slot((1+photon_threshold) * gain + zero_peak):])
        n_single_photons -= n_double_photons
        print("n_single_photons: %i" %n_single_photons)
        print("n_double_photons: %i" %n_double_photons)
        #print n_double_photons/(n_single_photons+n_double_photons)
      return fitted_gaussians
  
    def single_peak_fit(self, hist, lower_threshold, upper_threshold, mean,
                        zero_peak_gaussian=None):
      lower_slot = 0
      for slot in hist.slot_centers():
        lower_slot += 1
        if slot > lower_threshold: break
      upper_slot = 0
      for slot in hist.slot_centers():
        upper_slot += 1
        if slot > upper_threshold: break
  
      x = hist.slot_centers()
      y = hist.slots().as_double()
      starting_gaussians = [curve_fitting.gaussian(
        a=flex.max(y[lower_slot:upper_slot]), b=mean, c=3)]
     # print starting_gaussians
      #mamin: fit gaussian will take the maximum between starting point (lower_slot) and ending (upper_slot) as a
      if zero_peak_gaussian is not None:
        y -= zero_peak_gaussian(x)
      if 1:
        fit = curve_fitting.lbfgs_minimiser(
          starting_gaussians, x[lower_slot:upper_slot], y[lower_slot:upper_slot])
        sigma = abs(fit.functions[0].params[2])
        if sigma < 1 or sigma > 10:
          if flex.sum(y[lower_slot:upper_slot]) < 15: #mamin I changed 15 to 5
            # No point wasting time attempting to fit a gaussian if there aren't any counts
            #raise PixelFitError("Not enough counts to fit gaussian")
            return fit
          print("using cma_es:", sigma)
          fit = curve_fitting.cma_es_minimiser(
            starting_gaussians, x[lower_slot:upper_slot], y[lower_slot:upper_slot])
      else:
        fit = curve_fitting.cma_es_minimiser(
          starting_gaussians, x[lower_slot:upper_slot], y[lower_slot:upper_slot])
      return fit
  
    def pixels(self):
      for pixel in self.histograms.keys():
        yield pixel
  
  def sliding_average(y, n=3):
    assert n % 2 == 1 # n must be odd
    n_either_side = n//2
    summed = y[n_either_side:-n_either_side]
    for i in range(n_either_side):
      summed += y[n_either_side-i-1:-(n_either_side+i+1)]
      end = -(n_either_side-i-1)
      if end == 0: end = None
      summed += y[(n_either_side+i+1):end]
    averaged = summed / n
    for i in range(n_either_side):
      averaged.insert(i, y[i])
    for i in range(n_either_side, 0, -1):
      averaged.append(y[-i])
    return averaged
  
  def hist_outline(hist):
  
    step_size = hist.slot_width()
    half_step_size = 0.5 * step_size
    n_slots = len(hist.slots())
  
    bins = flex.double(n_slots * 2 + 2, 0)
    data = flex.double(n_slots * 2 + 2, 0)
    for i in range(n_slots):
      bins[2 * i + 1] = hist.slot_centers()[i] - half_step_size
      bins[2 * i + 2] = hist.slot_centers()[i] + half_step_size
      data[2 * i + 1] = hist.slots()[i]
      data[2 * i + 2] = hist.slots()[i]
  
    bins[0] = bins[1] - step_size
    bins[-1] = bins[-2] + step_size
    data[0] = 0
    data[-1] = 0
  
    return (bins, data)
  
  
  class PixelFitError(RuntimeError):
    pass
  
  def check_pixel_histogram_fit(hist, gaussians):
    assert gaussians is not None
    #if gaussians is None:
      ## Presumably the peak fitting failed in some way
      #print "Skipping pixel %s" %str(pixel)
      #continue
    zero_peak_diff = gaussians[0].params[1]
    if len(gaussians) < 2:
      raise PixelFitError("Only one gaussian!")
    y_obs = hist.slots().as_double()
    x = hist.slot_centers()
    y_calc = flex.double(y_obs.size(), 0)
    for g in gaussians:
      y_calc += g(x)
    residual = y_obs - y_calc
  
    # check the overall residual
    if flex.max(residual)/flex.sum(hist.slots()) > 0.015:
      raise PixelFitError("Bad fit residual: %f" %(flex.max(residual)/flex.sum(hist.slots())))
  
    # check the residual around the zero photon peak
    zero_gaussian = gaussians[0]
    selection = ((x < zero_gaussian.params[1] + 1 * zero_gaussian.sigma))
    #if ((flex.max(residual.select(selection))/flex.sum(hist.slots()) > 0.008)
        #or (flex.min(residual.select(selection))/flex.sum(hist.slots()) < -0.0067)):
      #raise PixelFitError("Bad fit residual around zero photon peak")
  
    # check the residual around the one photon peak
    one_gaussian = gaussians[1]
    selection = ((x > one_gaussian.params[1] - 1.4 * one_gaussian.sigma)) # &
    if selection.count(True) == 0:
      raise PixelFitError("Bad fit residual around one photon peak")
    max_residual_sel = flex.max(residual.select(selection))
    if max_residual_sel > 20 and max_residual_sel > 1.2 * one_gaussian.params[0]:
      raise PixelFitError("Bad fit residual: %f" %max_residual_sel)
  
    gain = gaussians[1].params[1] - gaussians[0].params[1]
    if 0 and estimated_gain is not None and abs(gain - estimated_gain) > 0.5 * estimated_gain:
      print("bad gain!!!!!", pixel, gain)
    #elif (one_gaussian.sigma / zero_gaussian.sigma) > 1.9:
      #raise PixelFitError("Bad sigma ratio: %.1f, %.1f" %(one_gaussian.sigma, zero_gaussian.sigma))
    elif gaussians[1].sigma < (0.5 * gaussians[0].sigma):
      raise PixelFitError("Bad sigma: %f" %gaussians[1].sigma)
    elif gain < (4 * gaussians[0].sigma):
      raise PixelFitError("Bad gain: %f" %gain)
    elif gain > (20 * gaussians[0].sigma): # XXX is 20 to low?
      raise PixelFitError("Bad gain: %f" %gain)
  
  if __name__ == '__main__':
    run(sys.argv[1:])
  