<reponame>dperl-sol/cctbx_project<gh_stars>100-1000
from __future__ import absolute_import, division, print_function
from six.moves import range
from rstbx.apps.stills.simple_integration import IntegrationMetaProcedure
from rstbx.apps import simple_integration
from scitbx import matrix
import math,copy
from dials.array_family import flex
from six.moves import zip

class integrate_one_frame(IntegrationMetaProcedure):
  def __init__(self, triclinic):
    simple_integration.__init__(self)
    IntegrationMetaProcedure.__init__(self)
    self.triclinic_pairs = triclinic.triclinic_pairs
    self.triclinic = triclinic

  def integration_concept(self, image_number, cb_op_to_primitive, verbose=False, **kwargs):
    if kwargs.get("user-reentrant") != None:
      return self.integration_concept_detail(experiments=kwargs.get("reentrant_experiments"),
                                    reflections=kwargs.get("reentrant_reflections"),
                                    spots=self.triclinic.get_observations_with_outlier_removal(),
                                    image_number=image_number, cb_op_to_primitive=cb_op_to_primitive, **kwargs)

    experiments = self.prepare_dxtbx_models(setting_specific_ai = self.triclinic.inputai, sg="P1")
    reflections = self.prepare_reflection_list(experiments[0].detector)
    self.refine(experiments=experiments, reflections=reflections)

    if self.horizons_phil.isoforms == []: # normal behavior; no isoform constraints
      setting_experiments = self.prepare_dxtbx_models(setting_specific_ai = self.inputai,
      sg = self.inputpd["symmetry"].space_group_info().type().lookup_symbol())
      setting_reflections = copy.deepcopy(reflections)
      setting_reflections["miller_index"] = cb_op_to_primitive.inverse().apply(reflections["miller_index"])
      R = self.refine(experiments=setting_experiments, reflections=setting_reflections)
    else:
      isoform_refineries = []
      setting_reflections = copy.deepcopy(reflections)
      setting_reflections["miller_index"] = cb_op_to_primitive.inverse().apply(reflections["miller_index"])
      look_symbol = self.inputpd["symmetry"].space_group_info().type().lookup_symbol()

      for isoform in self.horizons_phil.isoforms:
        print("Testing isoform %s"%isoform.name,isoform.cell.parameters())
        print("asserting", look_symbol ,"==", isoform.lookup_symbol)

        assert look_symbol == isoform.lookup_symbol
        setting_experiments = self.prepare_dxtbx_models(setting_specific_ai = self.inputai,
          sg = look_symbol, isoform = isoform.cell )
        P = self.refine(experiments=setting_experiments, reflections=setting_reflections, isoform=True)
        print(P.rmsds())
        isoform_refineries.append(P)
      positional_rmsds = [math.sqrt(P.rmsds()[0]**2 + P.rmsds()[1]**2) for P in isoform_refineries]
      print("Positional rmsds for all isoforms:", positional_rmsds)
      minrmsd_mm = min(positional_rmsds)
      minindex = positional_rmsds.index(minrmsd_mm)
      print("The smallest rmsd is %5.1f um from isoform %s"%(1000.*minrmsd_mm,self.horizons_phil.isoforms[minindex].name))
      if self.horizons_phil.isoforms[minindex].rmsd_target_mm is not None:
        print("asserting",minrmsd_mm ,"<", self.horizons_phil.isoforms[minindex].rmsd_target_mm)
        assert minrmsd_mm < self.horizons_phil.isoforms[minindex].rmsd_target_mm
      print("Acceptable rmsd for isoform %s."%(self.horizons_phil.isoforms[minindex].name), end=' ')
      if len (self.horizons_phil.isoforms)==2:
        print("Rmsd gain over the other isoform %5.1f um."%(1000.*abs(positional_rmsds[0] - positional_rmsds[1])))
      else:
        print()
      R = isoform_refineries[minindex]
      # Now one last check to see if direct beam is out of bounds
      if self.horizons_phil.isoforms[minindex].beam_restraint is not None:
        refined_beam = matrix.col(R.get_experiments()[0].detector[0].get_beam_centre(experiments[0].beam.get_s0()))
        known_beam = matrix.col(self.horizons_phil.isoforms[minindex].beam_restraint)
        print("asserting",(refined_beam-known_beam).length(),"<",self.horizons_phil.isoforms[minindex].rmsd_target_mm)
        assert (refined_beam-known_beam).length() < self.horizons_phil.isoforms[minindex].rmsd_target_mm
        # future--circle of confusion could be given as a separate length in mm instead of reusing rmsd_target
      self.identified_isoform = self.horizons_phil.isoforms[minindex].name

    self.integration_concept_detail(experiments=R.get_experiments(), reflections=setting_reflections,
                                    spots=self.triclinic.get_observations_with_outlier_removal(),
                                    image_number=image_number, cb_op_to_primitive=cb_op_to_primitive, **kwargs)

  def integration_concept_detail(self, experiments, reflections, spots,image_number,cb_op_to_primitive,**kwargs):
    detector = experiments[0].detector
    crystal = experiments[0].crystal
    from cctbx.crystal import symmetry
    c_symmetry = symmetry(space_group = crystal.get_space_group(), unit_cell = crystal.get_unit_cell())

    self.image_number = image_number
    NEAR = 10
    pxlsz = detector[0].get_pixel_size()

    Predicted = self.get_predictions_accounting_for_centering(experiments,reflections,cb_op_to_primitive,**kwargs)

    FWMOSAICITY = self.inputai.getMosaicity()
    self.DOMAIN_SZ_ANG = kwargs.get("domain_size_ang",  self.__dict__.get("actual",0)  )
    refineflag = {True:0,False:1}[kwargs.get("domain_size_ang",0)==0]
    c_symmetry.show_summary(prefix="EXCURSION%1d REPORT FWMOS= %6.4f DOMAIN= %6.1f "%(refineflag,FWMOSAICITY,self.DOMAIN_SZ_ANG))
    from annlib_ext import AnnAdaptor
    self.cell = c_symmetry.unit_cell()

    query = flex.double()
    print(len(self.predicted))

    for pred in self.predicted: # predicted spot coord in pixels
      query.append(pred[0]/pxlsz[0])
      query.append(pred[1]/pxlsz[1])

    self.reserve_hkllist_for_signal_search = self.hkllist

    reference = flex.double()

    assert self.length>NEAR# Can't do spot/pred matching with too few spots
    for spot in spots:
      reference.append(spot.ctr_mass_x())
      reference.append(spot.ctr_mass_y())

    IS_adapt = AnnAdaptor(data=reference,dim=2,k=NEAR)
    IS_adapt.query(query)
    idx_cutoff = float(min(self.mask_focus[image_number]))

    from rstbx.apps.slip_helpers import slip_callbacks
    cache_refinement_spots = getattr(slip_callbacks.slip_callback,"requires_refinement_spots",False)

    indexed_pairs_provisional = []
    correction_vectors_provisional = []
    c_v_p_flex = flex.vec3_double()
    this_setting_matched_indices = reflections["miller_index"]
    for j,item in enumerate(this_setting_matched_indices):
      this_setting_index = self.hkllist.first_index(item)
      if this_setting_index:
        Match = dict(spot=j,pred=this_setting_index)
        indexed_pairs_provisional.append(Match)
        vector = matrix.col(
            [reflections["xyzobs.px.value"][j][0] - self.predicted[Match["pred"]][0]/pxlsz[0],
             reflections["xyzobs.px.value"][j][1] - self.predicted[Match["pred"]][1]/pxlsz[1]])
        correction_vectors_provisional.append(vector)
        c_v_p_flex.append((vector[0],vector[1],0.))
    self.N_correction_vectors = len(correction_vectors_provisional)
    self.rmsd_px = math.sqrt(flex.mean(c_v_p_flex.dot(c_v_p_flex)))
    print("... %d provisional matches"%self.N_correction_vectors, end=' ')
    print("r.m.s.d. in pixels: %6.3f"%(self.rmsd_px))

    if self.horizons_phil.integration.enable_residual_scatter:
      from matplotlib import pyplot as plt
      fig = plt.figure()
      for cv in correction_vectors_provisional:
        plt.plot([cv[1]],[-cv[0]],"r.")
      plt.title(" %d matches, r.m.s.d. %5.2f pixels"%(len(correction_vectors_provisional),math.sqrt(flex.mean(c_v_p_flex.dot(c_v_p_flex)))))
      plt.axes().set_aspect("equal")
      self.show_figure(plt,fig,"res")
      plt.close()

    if self.horizons_phil.integration.enable_residual_map:
      from matplotlib import pyplot as plt
      PX = reflections["xyzobs.px.value"]
      fig = plt.figure()
      for match,cv in zip(indexed_pairs_provisional,correction_vectors_provisional):
        plt.plot([PX[match["spot"]][1]],[-PX[match["spot"]][0]],"r.")
        plt.plot([self.predicted[match["pred"]][1]/pxlsz[1]],[-self.predicted[match["pred"]][0]/pxlsz[0]],"g.")
        plt.plot([PX[match["spot"]][1], PX[match["spot"]][1] + 10.*cv[1]],
                 [-PX[match["spot"]][0], -PX[match["spot"]][0] - 10.*cv[0]],'r-')
      if kwargs.get("user-reentrant") != None and self.horizons_phil.integration.spot_prediction == "dials" \
             and self.horizons_phil.integration.enable_residual_map_deltapsi:
        from rstbx.apps.stills.util import residual_map_special_deltapsi_add_on
        residual_map_special_deltapsi_add_on(
          reflections = self.dials_spot_prediction,
          matches = indexed_pairs_provisional, experiments=experiments,
          hkllist = self.hkllist,
          predicted = self.predicted, plot=plt, eta_deg=FWMOSAICITY, deff=self.DOMAIN_SZ_ANG
          )
      plt.xlim([0,detector[0].get_image_size()[1]])
      plt.ylim([-detector[0].get_image_size()[0],0])
      plt.title(" %d matches, r.m.s.d. %5.2f pixels"%(len(correction_vectors_provisional),math.sqrt(flex.mean(c_v_p_flex.dot(c_v_p_flex)))))
      plt.axes().set_aspect("equal")
      self.show_figure(plt,fig,"map")
      plt.close()

    indexed_pairs = indexed_pairs_provisional
    correction_vectors = correction_vectors_provisional
    ########### skip outlier rejection for this derived class

    ### However must retain the ability to write out correction vectiors.
    if True: # at Aaron's request; test later
      correction_lengths = flex.double([v.length() for v in correction_vectors_provisional])
      clorder = flex.sort_permutation(correction_lengths)
      sorted_cl = correction_lengths.select(clorder)
      indexed_pairs = []
      correction_vectors = []
      self.correction_vectors = []
      for icand in range(len(sorted_cl)):
        # somewhat arbitrary sigma = 1.0 cutoff for outliers
        indexed_pairs.append(indexed_pairs_provisional[clorder[icand]])
        correction_vectors.append(correction_vectors_provisional[clorder[icand]])
        if cache_refinement_spots:
          self.spotfinder.images[self.frame_numbers[self.image_number]]["refinement_spots"].append(
          spots[reflections[indexed_pairs[-1]["spot"]]['spotfinder_lookup']])
        if kwargs.get("verbose_cv")==True:
            print("CV OBSCENTER %7.2f %7.2f REFINEDCENTER %7.2f %7.2f"%(
              float(self.inputpd["size1"])/2.,float(self.inputpd["size2"])/2.,
              self.inputai.xbeam()/pxlsz[0], self.inputai.ybeam()/pxlsz[1]), end=' ')
            print("OBSSPOT %7.2f %7.2f PREDSPOT %7.2f %7.2f"%(
              reflections[indexed_pairs[-1]["spot"]]['xyzobs.px.value'][0],
              reflections[indexed_pairs[-1]["spot"]]['xyzobs.px.value'][1],
              self.predicted[indexed_pairs[-1]["pred"]][0]/pxlsz[0],
              self.predicted[indexed_pairs[-1]["pred"]][1]/pxlsz[1]), end=' ')
            the_hkl = self.hkllist[indexed_pairs[-1]["pred"]]
            print("HKL %4d %4d %4d"%the_hkl,"%2d"%self.setting_id, end=' ')
            radial, azimuthal = spots[indexed_pairs[-1]["spot"]].get_radial_and_azimuthal_size(
              self.inputai.xbeam()/pxlsz[0], self.inputai.ybeam()/pxlsz[1])
            print("RADIALpx %5.3f AZIMUTpx %5.3f"%(radial,azimuthal))

        # Store a list of correction vectors in self.
        radial, azimuthal = spots[indexed_pairs[-1]['spot']].get_radial_and_azimuthal_size(
          self.inputai.xbeam()/pxlsz[0], self.inputai.ybeam()/pxlsz[1])
        self.correction_vectors.append(
          dict(obscenter=(float(self.inputpd['size1']) / 2,
                          float(self.inputpd['size2']) / 2),
               refinedcenter=(self.inputai.xbeam() / pxlsz[0],
                              self.inputai.ybeam() / pxlsz[1]),
               obsspot=(reflections[indexed_pairs[-1]['spot']]['xyzobs.px.value'][0],
                        reflections[indexed_pairs[-1]['spot']]['xyzobs.px.value'][1]),
               predspot=(self.predicted[indexed_pairs[-1]['pred']][0] / pxlsz[0],
                         self.predicted[indexed_pairs[-1]['pred']][1] / pxlsz[1]),
               hkl=(self.hkllist[indexed_pairs[-1]['pred']][0],
                    self.hkllist[indexed_pairs[-1]['pred']][1],
                    self.hkllist[indexed_pairs[-1]['pred']][2]),
               setting_id=self.setting_id,
               radial=radial,
               azimuthal=azimuthal))


    self.inputpd["symmetry"] = c_symmetry
    self.inputpd["symmetry"].show_summary(prefix="SETTING ")


    if self.horizons_phil.integration.model == "user_supplied":
      # Not certain of whether the reentrant_* dictionary keys create a memory leak
      if kwargs.get("user-reentrant",None)==None:
        kwargs["reentrant_experiments"] = experiments
        kwargs["reentrant_reflections"] = reflections
        from cxi_user import post_outlier_rejection
        self.indexed_pairs = indexed_pairs
        self.spots = spots
        post_outlier_rejection(self,image_number,cb_op_to_primitive,self.horizons_phil,kwargs)
        return
    ########### finished with user-supplied code


    correction_lengths=flex.double([v.length() for v in correction_vectors])

    self.r_residual = pxlsz[0]*flex.mean(correction_lengths)

    #assert len(indexed_pairs)>NEAR # must have enough indexed spots
    if (len(indexed_pairs) <= NEAR):
      raise Sorry("Not enough indexed spots, only found %d, need %d" % (len(indexed_pairs), NEAR))

    reference = flex.double()
    for item in indexed_pairs:
      reference.append(spots[item["spot"]].ctr_mass_x())
      reference.append(spots[item["spot"]].ctr_mass_y())

    PS_adapt = AnnAdaptor(data=reference,dim=2,k=NEAR)
    PS_adapt.query(query)

    self.BSmasks = []
    # do not use null: self.null_correction_mapping( predicted=self.predicted,
    self.positional_correction_mapping( predicted=self.predicted,
                                        correction_vectors = correction_vectors,
                                        PS_adapt = PS_adapt,
                                        IS_adapt = IS_adapt,
                                        spots = spots)

    # which spots are close enough to interfere with background?
    MAXOVER=6
    OS_adapt = AnnAdaptor(data=query,dim=2,k=MAXOVER) #six near nbrs
    OS_adapt.query(query)
    if self.mask_focus[image_number] is None:
      raise Sorry("No observed/predicted spot agreement; no Spotfinder masks; skip integration")
    nbr_cutoff = 2.0* max(self.mask_focus[image_number])
    FRAME = int(nbr_cutoff/2)
    #print "The overlap cutoff is %d pixels"%nbr_cutoff
    nbr_cutoff_sq = nbr_cutoff * nbr_cutoff

    #print "Optimized C++ section...",
    self.set_frame(FRAME)
    self.set_background_factor(kwargs["background_factor"])
    self.set_nbr_cutoff_sq(nbr_cutoff_sq)
    self.set_guard_width_sq(self.horizons_phil.integration.guard_width_sq)
    self.set_detector_gain(self.horizons_phil.integration.detector_gain)
    flex_sorted = flex.int()
    for item in self.sorted:
      flex_sorted.append(item[0]);flex_sorted.append(item[1]);

    if self.horizons_phil.integration.mask_pixel_value is not None:
      self.set_mask_pixel_val(self.horizons_phil.integration.mask_pixel_value)

    image_obj = self.imagefiles.imageindex(self.frame_numbers[self.image_number])
    image_obj.read()
    rawdata = image_obj.linearintdata # assume image #1

    if self.inputai.active_areas != None:
      self.detector_xy_draft = self.safe_background( rawdata=rawdata,
                          predicted=self.predicted,
                          OS_adapt=OS_adapt,
                          sorted=flex_sorted,
                          tiles=self.inputai.active_areas.IT,
                          tile_id=self.inputai.active_areas.tile_id);
    else:
      self.detector_xy_draft = self.safe_background( rawdata=rawdata,
                          predicted=self.predicted,
                          OS_adapt=OS_adapt,
                          sorted=flex_sorted);
    for i in range(len(self.predicted)): # loop over predicteds
      B_S_mask = {}
      keys = self.get_bsmask(i)
      for k in range(0,len(keys),2):
        B_S_mask[(keys[k],keys[k+1])]=True
      self.BSmasks.append(B_S_mask)
    #print "Done"
    return

  def get_predictions_accounting_for_centering(self,experiments,reflections,cb_op_to_primitive,**kwargs):
    # interface requires this function to set current_orientation
    # in the actual setting used for Miller index calculation
    detector = experiments[0].detector
    crystal = experiments[0].crystal

    if self.horizons_phil.integration.model == "user_supplied":

      lower_limit_domain_size = math.pow(
       crystal.get_unit_cell().volume(),
       1./3.)*self.horizons_phil.integration.mosaic.domain_size_lower_limit # default 10-unit cell block size minimum reasonable domain
      actual_used_domain_size = kwargs.get("domain_size_ang",lower_limit_domain_size)

      self.block_counter+=1
      rot_mat = matrix.sqr(cb_op_to_primitive.c().r().as_double()).transpose()

      from cctbx.crystal_orientation import crystal_orientation, basis_type
      centered_orientation = crystal_orientation(crystal.get_A(),basis_type.reciprocal)
      self.current_orientation = centered_orientation
      self.current_cb_op_to_primitive = cb_op_to_primitive
      primitive_orientation = centered_orientation.change_basis(rot_mat)

      self.inputai.setOrientation(primitive_orientation)
      from cxi_user import pre_get_predictions
      if self.block_counter < 2:
        KLUDGE = self.horizons_phil.integration.mosaic.kludge1 # bugfix 1 of 2 for protocol 6, equation 2
        self.inputai.setMosaicity(KLUDGE*self.inputai.getMosaicity())

      oldbase = self.inputai.getBase()

      #print oldbase.xbeam, oldbase.ybeam
      newbeam = detector[0].get_beam_centre(experiments[0].beam.get_s0())
      newdistance = -detector[0].get_beam_centre_lab(experiments[0].beam.get_s0())[2]

      from labelit.dptbx import Parameters
      base = Parameters(xbeam = newbeam[0], ybeam = newbeam[1], #oldbase.xbeam, ybeam = oldbase.ybeam,
                        distance = newdistance, twotheta = 0.0)
      self.inputai.setBase(base)
      self.inputai.setWavelength(experiments[0].beam.get_wavelength())

      self.bp3_wrapper = pre_get_predictions(self.inputai, self.horizons_phil,
        raw_image = self.imagefiles.images[self.image_number],
        imageindex = self.frame_numbers[self.image_number],
        spotfinder = self.spotfinder,
        limiting_resolution = self.limiting_resolution,
        domain_size_ang = actual_used_domain_size,
        )

      BPpredicted = self.bp3_wrapper.ucbp3.selected_predictions_labelit_format()
      BPhkllist = self.bp3_wrapper.ucbp3.selected_hkls()

      self.actual = actual_used_domain_size
      primitive_hkllist = BPhkllist
      #not sure if matrix needs to be transposed first for outputting HKL's???:
      self.hkllist = cb_op_to_primitive.inverse().apply(primitive_hkllist)

      if self.horizons_phil.integration.spot_prediction == "dials":
        from dials.algorithms.spot_prediction import StillsReflectionPredictor
        predictor = StillsReflectionPredictor(experiments[0])
        Rcalc = flex.reflection_table.empty_standard(len(self.hkllist))
        Rcalc['miller_index'] = self.hkllist
        predictor.for_reflection_table(Rcalc, crystal.get_A())
        self.predicted = Rcalc['xyzcal.mm']
        self.dials_spot_prediction = Rcalc
        self.dials_model = experiments

      elif self.horizons_phil.integration.spot_prediction == "ucbp3":
        self.predicted = BPpredicted

      self.inputai.setOrientation(centered_orientation)
      if self.inputai.active_areas != None:
        self.predicted,self.hkllist = self.inputai.active_areas(
                                      self.predicted,self.hkllist,self.pixel_size)
      if self.block_counter < 2:
         down = self.inputai.getMosaicity()/KLUDGE
         print("Readjusting mosaicity back down to ",down)
         self.inputai.setMosaicity(down)
      return

  def refine(self, experiments, reflections, isoform=None):
    from dials.algorithms.refinement.refiner import phil_scope
    from libtbx.phil import parse

    params = phil_scope.fetch(source=parse('')).extract()
    params.refinement.reflections.weighting_strategy.delpsi_constant=100000.
    params.refinement.reflections.weighting_strategy.override="stills"
    params.refinement.parameterisation.auto_reduction.action="fix"
    #params.refinement.reflections.do_outlier_rejection=True
    #params.refinement.reflections.iqr_multiplier=0.5
    #params.refinement.reflections.minimum_sample_size=50
    #params.refinement.reflections.maximum_sample_size=50
    #params.refinement.reflections.random_seed=1
    if self.horizons_phil.integration.dials_refinement.strategy=="distance":
      params.refinement.parameterisation.beam.fix="all"
      params.refinement.parameterisation.detector.fix_list=["Tau1"] # fix detector rotz, allow distance to refine
    elif self.horizons_phil.integration.dials_refinement.strategy=="wavelength":
      params.refinement.parameterisation.beam.fix="in_spindle_plane,out_spindle_plane"
      params.refinement.parameterisation.detector.fix_list=["Dist","Tau1"] # fix detector rotz and distance
    elif self.horizons_phil.integration.dials_refinement.strategy=="fix":
      params.refinement.parameterisation.beam.fix="all"
      params.refinement.parameterisation.detector.fix_list=["Dist","Tau1"] # fix detector rotz and distance
    if isoform is not None:
      params.refinement.reflections.outlier.algorithm="null"
      params.refinement.parameterisation.crystal.fix="cell"
    from dials.algorithms.refinement.refiner import RefinerFactory
    refiner = RefinerFactory.from_parameters_data_experiments(params,
      reflections, experiments)

    history = refiner.run()
    print(history.keys())
    for item in history["rmsd"]:
      print("%5.2f %5.2f %8.5f"%(item[0],item[1],180.*item[2]/math.pi))

    #for item in history["parameter_vector"]:
    #  print ["%8.5f"%a for a in item]
    print(refiner.selection_used_for_refinement().count(True),"spots used for refinement")
    print(refiner.get_experiments()[0].beam)
    print(refiner.get_experiments()[0].detector)
    print("Distance:", -refiner.get_experiments()[0].detector[0].get_beam_centre_lab(refiner.get_experiments()[0].beam.get_s0())[2])
    print(refiner.get_experiments()[0].crystal)
    return refiner

  def prepare_reflection_list(self,detector):

    spots = self.triclinic.get_observations_with_outlier_removal()
    ordinary_python_list_of_indexed_observations = [
      {
        "id":0,
        "panel":0,
        "miller_index":item["pred"],
        "xyzobs.px.value":(spots[item["spot"]].ctr_mass_x(),spots[item["spot"]].ctr_mass_y(),0.0),
        "xyzobs.px.variance":(0.25,0.25,0.25),
        "spotfinder_lookup":item["spot"]
      }
      for item in self.triclinic_pairs
    ]

    self.length = len(ordinary_python_list_of_indexed_observations)
    R= flex.reflection_table.empty_standard(self.length)

    R['miller_index'] = flex.miller_index([item["miller_index"] for item in ordinary_python_list_of_indexed_observations])
    R['xyzobs.px.value'] = flex.vec3_double([item["xyzobs.px.value"] for item in ordinary_python_list_of_indexed_observations])
    R['xyzobs.px.variance'] = flex.vec3_double([item["xyzobs.px.variance"] for item in ordinary_python_list_of_indexed_observations])
    R['spotfinder_lookup'] = flex.int([item["spotfinder_lookup"] for item in ordinary_python_list_of_indexed_observations])

    R['xyzobs.mm.value'] = flex.vec3_double(self.length)
    R['xyzobs.mm.variance'] = flex.vec3_double(self.length)

    pxlsz = detector[0].get_pixel_size()

    for idx in range(self.length):
      R['xyzobs.mm.value'][idx] = (R['xyzobs.px.value'][idx][0]*pxlsz[0], R['xyzobs.px.value'][idx][1]*pxlsz[1], R['xyzobs.px.value'][idx][2])
      R['xyzobs.mm.variance'][idx] = (R['xyzobs.px.variance'][idx][0]*pxlsz[0], R['xyzobs.px.variance'][idx][1]*pxlsz[1], R['xyzobs.px.variance'][idx][2])

    return R

  def prepare_dxtbx_models(self,setting_specific_ai,sg,isoform=None):

    from dxtbx.model import BeamFactory
    beam = BeamFactory.simple(wavelength = self.inputai.wavelength)

    from dxtbx.model import DetectorFactory
    detector = DetectorFactory.simple(
      sensor = DetectorFactory.sensor("PAD"),
      distance = setting_specific_ai.distance(),
      beam_centre = [setting_specific_ai.xbeam(), setting_specific_ai.ybeam()],
      fast_direction = "+x",
      slow_direction = "+y",
      pixel_size = [self.pixel_size,self.pixel_size],
      image_size = [self.inputpd['size1'],self.inputpd['size1']],
      )

    direct = matrix.sqr(setting_specific_ai.getOrientation().direct_matrix())
    from dxtbx.model import MosaicCrystalKabsch2010
    crystal = MosaicCrystalKabsch2010(
      real_space_a = matrix.row(direct[0:3]),
      real_space_b = matrix.row(direct[3:6]),
      real_space_c = matrix.row(direct[6:9]),
      space_group_symbol = sg,
    )
    crystal.set_mosaicity(setting_specific_ai.getMosaicity())
    if isoform is not None:
      newB = matrix.sqr(isoform.fractionalization_matrix()).transpose()
      crystal.set_B(newB)

    from dxtbx.model import Experiment, ExperimentList
    experiments = ExperimentList()
    experiments.append(Experiment(beam=beam,
                                  detector=detector,
                                  crystal=crystal))

    print(beam)
    print(detector)
    print(crystal)
    return experiments

from __future__ import absolute_import, division, print_function
from six.moves import range
from rstbx.apps.stills.simple_integration import IntegrationMetaProcedure
from rstbx.apps import simple_integration
from scitbx import matrix
import math,copy
from dials.array_family import flex
from six.moves import zip

class integrate_one_frame(IntegrationMetaProcedure):
  def __init__(self, triclinic):
    simple_integration.__init__(self)
    IntegrationMetaProcedure.__init__(self)
    self.triclinic_pairs = triclinic.triclinic_pairs
    self.triclinic = triclinic

  def integration_concept(self, image_number, cb_op_to_primitive, verbose=False, **kwargs):
    if kwargs.get("user-reentrant") != None:
      return self.integration_concept_detail(experiments=kwargs.get("reentrant_experiments"),
                                    reflections=kwargs.get("reentrant_reflections"),
                                    spots=self.triclinic.get_observations_with_outlier_removal(),
                                    image_number=image_number, cb_op_to_primitive=cb_op_to_primitive, **kwargs)

    experiments = self.prepare_dxtbx_models(setting_specific_ai = self.triclinic.inputai, sg="P1")
    reflections = self.prepare_reflection_list(experiments[0].detector)
    self.refine(experiments=experiments, reflections=reflections)

    if self.horizons_phil.isoforms == []: # normal behavior; no isoform constraints
      setting_experiments = self.prepare_dxtbx_models(setting_specific_ai = self.inputai,
      sg = self.inputpd["symmetry"].space_group_info().type().lookup_symbol())
      setting_reflections = copy.deepcopy(reflections)
      setting_reflections["miller_index"] = cb_op_to_primitive.inverse().apply(reflections["miller_index"])
      R = self.refine(experiments=setting_experiments, reflections=setting_reflections)
    else:
      isoform_refineries = []
      setting_reflections = copy.deepcopy(reflections)
      setting_reflections["miller_index"] = cb_op_to_primitive.inverse().apply(reflections["miller_index"])
      look_symbol = self.inputpd["symmetry"].space_group_info().type().lookup_symbol()

      for isoform in self.horizons_phil.isoforms:
        print("Testing isoform %s"%isoform.name,isoform.cell.parameters())
        print("asserting", look_symbol ,"==", isoform.lookup_symbol)

        assert look_symbol == isoform.lookup_symbol
        setting_experiments = self.prepare_dxtbx_models(setting_specific_ai = self.inputai,
          sg = look_symbol, isoform = isoform.cell )
        P = self.refine(experiments=setting_experiments, reflections=setting_reflections, isoform=True)
        print(P.rmsds())
        isoform_refineries.append(P)
      positional_rmsds = [math.sqrt(P.rmsds()[0]**2 + P.rmsds()[1]**2) for P in isoform_refineries]
      print("Positional rmsds for all isoforms:", positional_rmsds)
      minrmsd_mm = min(positional_rmsds)
      minindex = positional_rmsds.index(minrmsd_mm)
      print("The smallest rmsd is %5.1f um from isoform %s"%(1000.*minrmsd_mm,self.horizons_phil.isoforms[minindex].name))
      if self.horizons_phil.isoforms[minindex].rmsd_target_mm is not None:
        print("asserting",minrmsd_mm ,"<", self.horizons_phil.isoforms[minindex].rmsd_target_mm)
        assert minrmsd_mm < self.horizons_phil.isoforms[minindex].rmsd_target_mm
      print("Acceptable rmsd for isoform %s."%(self.horizons_phil.isoforms[minindex].name), end=' ')
      if len (self.horizons_phil.isoforms)==2:
        print("Rmsd gain over the other isoform %5.1f um."%(1000.*abs(positional_rmsds[0] - positional_rmsds[1])))
      else:
        print()
      R = isoform_refineries[minindex]
      # Now one last check to see if direct beam is out of bounds
      if self.horizons_phil.isoforms[minindex].beam_restraint is not None:
        refined_beam = matrix.col(R.get_experiments()[0].detector[0].get_beam_centre(experiments[0].beam.get_s0()))
        known_beam = matrix.col(self.horizons_phil.isoforms[minindex].beam_restraint)
        print("asserting",(refined_beam-known_beam).length(),"<",self.horizons_phil.isoforms[minindex].rmsd_target_mm)
        assert (refined_beam-known_beam).length() < self.horizons_phil.isoforms[minindex].rmsd_target_mm
        # future--circle of confusion could be given as a separate length in mm instead of reusing rmsd_target
      self.identified_isoform = self.horizons_phil.isoforms[minindex].name

    self.integration_concept_detail(experiments=R.get_experiments(), reflections=setting_reflections,
                                    spots=self.triclinic.get_observations_with_outlier_removal(),
                                    image_number=image_number, cb_op_to_primitive=cb_op_to_primitive, **kwargs)

  def integration_concept_detail(self, experiments, reflections, spots,image_number,cb_op_to_primitive,**kwargs):
    detector = experiments[0].detector
    crystal = experiments[0].crystal
    from cctbx.crystal import symmetry
    c_symmetry = symmetry(space_group = crystal.get_space_group(), unit_cell = crystal.get_unit_cell())

    self.image_number = image_number
    NEAR = 10
    pxlsz = detector[0].get_pixel_size()

    Predicted = self.get_predictions_accounting_for_centering(experiments,reflections,cb_op_to_primitive,**kwargs)

    FWMOSAICITY = self.inputai.getMosaicity()
    self.DOMAIN_SZ_ANG = kwargs.get("domain_size_ang",  self.__dict__.get("actual",0)  )
    refineflag = {True:0,False:1}[kwargs.get("domain_size_ang",0)==0]
    c_symmetry.show_summary(prefix="EXCURSION%1d REPORT FWMOS= %6.4f DOMAIN= %6.1f "%(refineflag,FWMOSAICITY,self.DOMAIN_SZ_ANG))
    from annlib_ext import AnnAdaptor
    self.cell = c_symmetry.unit_cell()

    query = flex.double()
    print(len(self.predicted))

    for pred in self.predicted: # predicted spot coord in pixels
      query.append(pred[0]/pxlsz[0])
      query.append(pred[1]/pxlsz[1])

    self.reserve_hkllist_for_signal_search = self.hkllist

    reference = flex.double()

    assert self.length>NEAR# Can't do spot/pred matching with too few spots
    for spot in spots:
      reference.append(spot.ctr_mass_x())
      reference.append(spot.ctr_mass_y())

    IS_adapt = AnnAdaptor(data=reference,dim=2,k=NEAR)
    IS_adapt.query(query)
    idx_cutoff = float(min(self.mask_focus[image_number]))

    from rstbx.apps.slip_helpers import slip_callbacks
    cache_refinement_spots = getattr(slip_callbacks.slip_callback,"requires_refinement_spots",False)

    indexed_pairs_provisional = []
    correction_vectors_provisional = []
    c_v_p_flex = flex.vec3_double()
    this_setting_matched_indices = reflections["miller_index"]
    for j,item in enumerate(this_setting_matched_indices):
      this_setting_index = self.hkllist.first_index(item)
      if this_setting_index:
        Match = dict(spot=j,pred=this_setting_index)
        indexed_pairs_provisional.append(Match)
        vector = matrix.col(
            [reflections["xyzobs.px.value"][j][0] - self.predicted[Match["pred"]][0]/pxlsz[0],
             reflections["xyzobs.px.value"][j][1] - self.predicted[Match["pred"]][1]/pxlsz[1]])
        correction_vectors_provisional.append(vector)
        c_v_p_flex.append((vector[0],vector[1],0.))
    self.N_correction_vectors = len(correction_vectors_provisional)
    self.rmsd_px = math.sqrt(flex.mean(c_v_p_flex.dot(c_v_p_flex)))
    print("... %d provisional matches"%self.N_correction_vectors, end=' ')
    print("r.m.s.d. in pixels: %6.3f"%(self.rmsd_px))

    if self.horizons_phil.integration.enable_residual_scatter:
      from matplotlib import pyplot as plt
      fig = plt.figure()
      for cv in correction_vectors_provisional:
        plt.plot([cv[1]],[-cv[0]],"r.")
      plt.title(" %d matches, r.m.s.d. %5.2f pixels"%(len(correction_vectors_provisional),math.sqrt(flex.mean(c_v_p_flex.dot(c_v_p_flex)))))
      plt.axes().set_aspect("equal")
      self.show_figure(plt,fig,"res")
      plt.close()

    if self.horizons_phil.integration.enable_residual_map:
      from matplotlib import pyplot as plt
      PX = reflections["xyzobs.px.value"]
      fig = plt.figure()
      for match,cv in zip(indexed_pairs_provisional,correction_vectors_provisional):
        plt.plot([PX[match["spot"]][1]],[-PX[match["spot"]][0]],"r.")
        plt.plot([self.predicted[match["pred"]][1]/pxlsz[1]],[-self.predicted[match["pred"]][0]/pxlsz[0]],"g.")
        plt.plot([PX[match["spot"]][1], PX[match["spot"]][1] + 10.*cv[1]],
                 [-PX[match["spot"]][0], -PX[match["spot"]][0] - 10.*cv[0]],'r-')
      if kwargs.get("user-reentrant") != None and self.horizons_phil.integration.spot_prediction == "dials" \
             and self.horizons_phil.integration.enable_residual_map_deltapsi:
        from rstbx.apps.stills.util import residual_map_special_deltapsi_add_on
        residual_map_special_deltapsi_add_on(
          reflections = self.dials_spot_prediction,
          matches = indexed_pairs_provisional, experiments=experiments,
          hkllist = self.hkllist,
          predicted = self.predicted, plot=plt, eta_deg=FWMOSAICITY, deff=self.DOMAIN_SZ_ANG
          )
      plt.xlim([0,detector[0].get_image_size()[1]])
      plt.ylim([-detector[0].get_image_size()[0],0])
      plt.title(" %d matches, r.m.s.d. %5.2f pixels"%(len(correction_vectors_provisional),math.sqrt(flex.mean(c_v_p_flex.dot(c_v_p_flex)))))
      plt.axes().set_aspect("equal")
      self.show_figure(plt,fig,"map")
      plt.close()

    indexed_pairs = indexed_pairs_provisional
    correction_vectors = correction_vectors_provisional
    ########### skip outlier rejection for this derived class

    ### However must retain the ability to write out correction vectiors.
    if True: # at Aaron's request; test later
      correction_lengths = flex.double([v.length() for v in correction_vectors_provisional])
      clorder = flex.sort_permutation(correction_lengths)
      sorted_cl = correction_lengths.select(clorder)
      indexed_pairs = []
      correction_vectors = []
      self.correction_vectors = []
      for icand in range(len(sorted_cl)):
        # somewhat arbitrary sigma = 1.0 cutoff for outliers
        indexed_pairs.append(indexed_pairs_provisional[clorder[icand]])
        correction_vectors.append(correction_vectors_provisional[clorder[icand]])
        if cache_refinement_spots:
          self.spotfinder.images[self.frame_numbers[self.image_number]]["refinement_spots"].append(
          spots[reflections[indexed_pairs[-1]["spot"]]['spotfinder_lookup']])
        if kwargs.get("verbose_cv")==True:
            print("CV OBSCENTER %7.2f %7.2f REFINEDCENTER %7.2f %7.2f"%(
              float(self.inputpd["size1"])/2.,float(self.inputpd["size2"])/2.,
              self.inputai.xbeam()/pxlsz[0], self.inputai.ybeam()/pxlsz[1]), end=' ')
            print("OBSSPOT %7.2f %7.2f PREDSPOT %7.2f %7.2f"%(
              reflections[indexed_pairs[-1]["spot"]]['xyzobs.px.value'][0],
              reflections[indexed_pairs[-1]["spot"]]['xyzobs.px.value'][1],
              self.predicted[indexed_pairs[-1]["pred"]][0]/pxlsz[0],
              self.predicted[indexed_pairs[-1]["pred"]][1]/pxlsz[1]), end=' ')
            the_hkl = self.hkllist[indexed_pairs[-1]["pred"]]
            print("HKL %4d %4d %4d"%the_hkl,"%2d"%self.setting_id, end=' ')
            radial, azimuthal = spots[indexed_pairs[-1]["spot"]].get_radial_and_azimuthal_size(
              self.inputai.xbeam()/pxlsz[0], self.inputai.ybeam()/pxlsz[1])
            print("RADIALpx %5.3f AZIMUTpx %5.3f"%(radial,azimuthal))

        # Store a list of correction vectors in self.
        radial, azimuthal = spots[indexed_pairs[-1]['spot']].get_radial_and_azimuthal_size(
          self.inputai.xbeam()/pxlsz[0], self.inputai.ybeam()/pxlsz[1])
        self.correction_vectors.append(
          dict(obscenter=(float(self.inputpd['size1']) / 2,
                          float(self.inputpd['size2']) / 2),
               refinedcenter=(self.inputai.xbeam() / pxlsz[0],
                              self.inputai.ybeam() / pxlsz[1]),
               obsspot=(reflections[indexed_pairs[-1]['spot']]['xyzobs.px.value'][0],
                        reflections[indexed_pairs[-1]['spot']]['xyzobs.px.value'][1]),
               predspot=(self.predicted[indexed_pairs[-1]['pred']][0] / pxlsz[0],
                         self.predicted[indexed_pairs[-1]['pred']][1] / pxlsz[1]),
               hkl=(self.hkllist[indexed_pairs[-1]['pred']][0],
                    self.hkllist[indexed_pairs[-1]['pred']][1],
                    self.hkllist[indexed_pairs[-1]['pred']][2]),
               setting_id=self.setting_id,
               radial=radial,
               azimuthal=azimuthal))


    self.inputpd["symmetry"] = c_symmetry
    self.inputpd["symmetry"].show_summary(prefix="SETTING ")


    if self.horizons_phil.integration.model == "user_supplied":
      # Not certain of whether the reentrant_* dictionary keys create a memory leak
      if kwargs.get("user-reentrant",None)==None:
        kwargs["reentrant_experiments"] = experiments
        kwargs["reentrant_reflections"] = reflections
        from cxi_user import post_outlier_rejection
        self.indexed_pairs = indexed_pairs
        self.spots = spots
        post_outlier_rejection(self,image_number,cb_op_to_primitive,self.horizons_phil,kwargs)
        return
    ########### finished with user-supplied code


    correction_lengths=flex.double([v.length() for v in correction_vectors])

    self.r_residual = pxlsz[0]*flex.mean(correction_lengths)

    #assert len(indexed_pairs)>NEAR # must have enough indexed spots
    if (len(indexed_pairs) <= NEAR):
      raise Sorry("Not enough indexed spots, only found %d, need %d" % (len(indexed_pairs), NEAR))

    reference = flex.double()
    for item in indexed_pairs:
      reference.append(spots[item["spot"]].ctr_mass_x())
      reference.append(spots[item["spot"]].ctr_mass_y())

    PS_adapt = AnnAdaptor(data=reference,dim=2,k=NEAR)
    PS_adapt.query(query)

    self.BSmasks = []
    # do not use null: self.null_correction_mapping( predicted=self.predicted,
    self.positional_correction_mapping( predicted=self.predicted,
                                        correction_vectors = correction_vectors,
                                        PS_adapt = PS_adapt,
                                        IS_adapt = IS_adapt,
                                        spots = spots)

    # which spots are close enough to interfere with background?
    MAXOVER=6
    OS_adapt = AnnAdaptor(data=query,dim=2,k=MAXOVER) #six near nbrs
    OS_adapt.query(query)
    if self.mask_focus[image_number] is None:
      raise Sorry("No observed/predicted spot agreement; no Spotfinder masks; skip integration")
    nbr_cutoff = 2.0* max(self.mask_focus[image_number])
    FRAME = int(nbr_cutoff/2)
    #print "The overlap cutoff is %d pixels"%nbr_cutoff
    nbr_cutoff_sq = nbr_cutoff * nbr_cutoff

    #print "Optimized C++ section...",
    self.set_frame(FRAME)
    self.set_background_factor(kwargs["background_factor"])
    self.set_nbr_cutoff_sq(nbr_cutoff_sq)
    self.set_guard_width_sq(self.horizons_phil.integration.guard_width_sq)
    self.set_detector_gain(self.horizons_phil.integration.detector_gain)
    flex_sorted = flex.int()
    for item in self.sorted:
      flex_sorted.append(item[0]);flex_sorted.append(item[1]);

    if self.horizons_phil.integration.mask_pixel_value is not None:
      self.set_mask_pixel_val(self.horizons_phil.integration.mask_pixel_value)

    image_obj = self.imagefiles.imageindex(self.frame_numbers[self.image_number])
    image_obj.read()
    rawdata = image_obj.linearintdata # assume image #1

    if self.inputai.active_areas != None:
      self.detector_xy_draft = self.safe_background( rawdata=rawdata,
                          predicted=self.predicted,
                          OS_adapt=OS_adapt,
                          sorted=flex_sorted,
                          tiles=self.inputai.active_areas.IT,
                          tile_id=self.inputai.active_areas.tile_id);
    else:
      self.detector_xy_draft = self.safe_background( rawdata=rawdata,
                          predicted=self.predicted,
                          OS_adapt=OS_adapt,
                          sorted=flex_sorted);
    for i in range(len(self.predicted)): # loop over predicteds
      B_S_mask = {}
      keys = self.get_bsmask(i)
      for k in range(0,len(keys),2):
        B_S_mask[(keys[k],keys[k+1])]=True
      self.BSmasks.append(B_S_mask)
    #print "Done"
    return

  def get_predictions_accounting_for_centering(self,experiments,reflections,cb_op_to_primitive,**kwargs):
    # interface requires this function to set current_orientation
    # in the actual setting used for Miller index calculation
    detector = experiments[0].detector
    crystal = experiments[0].crystal

    if self.horizons_phil.integration.model == "user_supplied":

      lower_limit_domain_size = math.pow(
       crystal.get_unit_cell().volume(),
       1./3.)*self.horizons_phil.integration.mosaic.domain_size_lower_limit # default 10-unit cell block size minimum reasonable domain
      actual_used_domain_size = kwargs.get("domain_size_ang",lower_limit_domain_size)

      self.block_counter+=1
      rot_mat = matrix.sqr(cb_op_to_primitive.c().r().as_double()).transpose()

      from cctbx.crystal_orientation import crystal_orientation, basis_type
      centered_orientation = crystal_orientation(crystal.get_A(),basis_type.reciprocal)
      self.current_orientation = centered_orientation
      self.current_cb_op_to_primitive = cb_op_to_primitive
      primitive_orientation = centered_orientation.change_basis(rot_mat)

      self.inputai.setOrientation(primitive_orientation)
      from cxi_user import pre_get_predictions
      if self.block_counter < 2:
        KLUDGE = self.horizons_phil.integration.mosaic.kludge1 # bugfix 1 of 2 for protocol 6, equation 2
        self.inputai.setMosaicity(KLUDGE*self.inputai.getMosaicity())

      oldbase = self.inputai.getBase()

      #print oldbase.xbeam, oldbase.ybeam
      newbeam = detector[0].get_beam_centre(experiments[0].beam.get_s0())
      newdistance = -detector[0].get_beam_centre_lab(experiments[0].beam.get_s0())[2]

      from labelit.dptbx import Parameters
      base = Parameters(xbeam = newbeam[0], ybeam = newbeam[1], #oldbase.xbeam, ybeam = oldbase.ybeam,
                        distance = newdistance, twotheta = 0.0)
      self.inputai.setBase(base)
      self.inputai.setWavelength(experiments[0].beam.get_wavelength())

      self.bp3_wrapper = pre_get_predictions(self.inputai, self.horizons_phil,
        raw_image = self.imagefiles.images[self.image_number],
        imageindex = self.frame_numbers[self.image_number],
        spotfinder = self.spotfinder,
        limiting_resolution = self.limiting_resolution,
        domain_size_ang = actual_used_domain_size,
        )

      BPpredicted = self.bp3_wrapper.ucbp3.selected_predictions_labelit_format()
      BPhkllist = self.bp3_wrapper.ucbp3.selected_hkls()

      self.actual = actual_used_domain_size
      primitive_hkllist = BPhkllist
      #not sure if matrix needs to be transposed first for outputting HKL's???:
      self.hkllist = cb_op_to_primitive.inverse().apply(primitive_hkllist)

      if self.horizons_phil.integration.spot_prediction == "dials":
        from dials.algorithms.spot_prediction import StillsReflectionPredictor
        predictor = StillsReflectionPredictor(experiments[0])
        Rcalc = flex.reflection_table.empty_standard(len(self.hkllist))
        Rcalc['miller_index'] = self.hkllist
        predictor.for_reflection_table(Rcalc, crystal.get_A())
        self.predicted = Rcalc['xyzcal.mm']
        self.dials_spot_prediction = Rcalc
        self.dials_model = experiments

      elif self.horizons_phil.integration.spot_prediction == "ucbp3":
        self.predicted = BPpredicted

      self.inputai.setOrientation(centered_orientation)
      if self.inputai.active_areas != None:
        self.predicted,self.hkllist = self.inputai.active_areas(
                                      self.predicted,self.hkllist,self.pixel_size)
      if self.block_counter < 2:
         down = self.inputai.getMosaicity()/KLUDGE
         print("Readjusting mosaicity back down to ",down)
         self.inputai.setMosaicity(down)
      return

  def refine(self, experiments, reflections, isoform=None):
    from dials.algorithms.refinement.refiner import phil_scope
    from libtbx.phil import parse

    params = phil_scope.fetch(source=parse('')).extract()
    params.refinement.reflections.weighting_strategy.delpsi_constant=100000.
    params.refinement.reflections.weighting_strategy.override="stills"
    params.refinement.parameterisation.auto_reduction.action="fix"
    #params.refinement.reflections.do_outlier_rejection=True
    #params.refinement.reflections.iqr_multiplier=0.5
    #params.refinement.reflections.minimum_sample_size=50
    #params.refinement.reflections.maximum_sample_size=50
    #params.refinement.reflections.random_seed=1
    if self.horizons_phil.integration.dials_refinement.strategy=="distance":
      params.refinement.parameterisation.beam.fix="all"
      params.refinement.parameterisation.detector.fix_list=["Tau1"] # fix detector rotz, allow distance to refine
    elif self.horizons_phil.integration.dials_refinement.strategy=="wavelength":
      params.refinement.parameterisation.beam.fix="in_spindle_plane,out_spindle_plane"
      params.refinement.parameterisation.detector.fix_list=["Dist","Tau1"] # fix detector rotz and distance
    elif self.horizons_phil.integration.dials_refinement.strategy=="fix":
      params.refinement.parameterisation.beam.fix="all"
      params.refinement.parameterisation.detector.fix_list=["Dist","Tau1"] # fix detector rotz and distance
    if isoform is not None:
      params.refinement.reflections.outlier.algorithm="null"
      params.refinement.parameterisation.crystal.fix="cell"
    from dials.algorithms.refinement.refiner import RefinerFactory
    refiner = RefinerFactory.from_parameters_data_experiments(params,
      reflections, experiments)

    history = refiner.run()
    print(history.keys())
    for item in history["rmsd"]:
      print("%5.2f %5.2f %8.5f"%(item[0],item[1],180.*item[2]/math.pi))

    #for item in history["parameter_vector"]:
    #  print ["%8.5f"%a for a in item]
    print(refiner.selection_used_for_refinement().count(True),"spots used for refinement")
    print(refiner.get_experiments()[0].beam)
    print(refiner.get_experiments()[0].detector)
    print("Distance:", -refiner.get_experiments()[0].detector[0].get_beam_centre_lab(refiner.get_experiments()[0].beam.get_s0())[2])
    print(refiner.get_experiments()[0].crystal)
    return refiner

  def prepare_reflection_list(self,detector):

    spots = self.triclinic.get_observations_with_outlier_removal()
    ordinary_python_list_of_indexed_observations = [
      {
        "id":0,
        "panel":0,
        "miller_index":item["pred"],
        "xyzobs.px.value":(spots[item["spot"]].ctr_mass_x(),spots[item["spot"]].ctr_mass_y(),0.0),
        "xyzobs.px.variance":(0.25,0.25,0.25),
        "spotfinder_lookup":item["spot"]
      }
      for item in self.triclinic_pairs
    ]

    self.length = len(ordinary_python_list_of_indexed_observations)
    R= flex.reflection_table.empty_standard(self.length)

    R['miller_index'] = flex.miller_index([item["miller_index"] for item in ordinary_python_list_of_indexed_observations])
    R['xyzobs.px.value'] = flex.vec3_double([item["xyzobs.px.value"] for item in ordinary_python_list_of_indexed_observations])
    R['xyzobs.px.variance'] = flex.vec3_double([item["xyzobs.px.variance"] for item in ordinary_python_list_of_indexed_observations])
    R['spotfinder_lookup'] = flex.int([item["spotfinder_lookup"] for item in ordinary_python_list_of_indexed_observations])

    R['xyzobs.mm.value'] = flex.vec3_double(self.length)
    R['xyzobs.mm.variance'] = flex.vec3_double(self.length)

    pxlsz = detector[0].get_pixel_size()

    for idx in range(self.length):
      R['xyzobs.mm.value'][idx] = (R['xyzobs.px.value'][idx][0]*pxlsz[0], R['xyzobs.px.value'][idx][1]*pxlsz[1], R['xyzobs.px.value'][idx][2])
      R['xyzobs.mm.variance'][idx] = (R['xyzobs.px.variance'][idx][0]*pxlsz[0], R['xyzobs.px.variance'][idx][1]*pxlsz[1], R['xyzobs.px.variance'][idx][2])

    return R

  def prepare_dxtbx_models(self,setting_specific_ai,sg,isoform=None):

    from dxtbx.model import BeamFactory
    beam = BeamFactory.simple(wavelength = self.inputai.wavelength)

    from dxtbx.model import DetectorFactory
    detector = DetectorFactory.simple(
      sensor = DetectorFactory.sensor("PAD"),
      distance = setting_specific_ai.distance(),
      beam_centre = [setting_specific_ai.xbeam(), setting_specific_ai.ybeam()],
      fast_direction = "+x",
      slow_direction = "+y",
      pixel_size = [self.pixel_size,self.pixel_size],
      image_size = [self.inputpd['size1'],self.inputpd['size1']],
      )

    direct = matrix.sqr(setting_specific_ai.getOrientation().direct_matrix())
    from dxtbx.model import MosaicCrystalKabsch2010
    crystal = MosaicCrystalKabsch2010(
      real_space_a = matrix.row(direct[0:3]),
      real_space_b = matrix.row(direct[3:6]),
      real_space_c = matrix.row(direct[6:9]),
      space_group_symbol = sg,
    )
    crystal.set_mosaicity(setting_specific_ai.getMosaicity())
    if isoform is not None:
      newB = matrix.sqr(isoform.fractionalization_matrix()).transpose()
      crystal.set_B(newB)

    from dxtbx.model import Experiment, ExperimentList
    experiments = ExperimentList()
    experiments.append(Experiment(beam=beam,
                                  detector=detector,
                                  crystal=crystal))

    print(beam)
    print(detector)
    print(crystal)
    return experiments
# normal behavior; no isoform constraints # Now one last check to see if direct beam is out of bounds # future--circle of confusion could be given as a separate length in mm instead of reusing rmsd_target # predicted spot coord in pixels # Can't do spot/pred matching with too few spots ########### skip outlier rejection for this derived class ### However must retain the ability to write out correction vectiors. # at Aaron's request; test later # somewhat arbitrary sigma = 1.0 cutoff for outliers # Store a list of correction vectors in self. # Not certain of whether the reentrant_* dictionary keys create a memory leak ########### finished with user-supplied code #assert len(indexed_pairs)>NEAR # must have enough indexed spots # do not use null: self.null_correction_mapping( predicted=self.predicted, # which spots are close enough to interfere with background? #six near nbrs #print "The overlap cutoff is %d pixels"%nbr_cutoff #print "Optimized C++ section...", # assume image #1 # loop over predicteds #print "Done" # interface requires this function to set current_orientation # in the actual setting used for Miller index calculation # default 10-unit cell block size minimum reasonable domain # bugfix 1 of 2 for protocol 6, equation 2 #print oldbase.xbeam, oldbase.ybeam #oldbase.xbeam, ybeam = oldbase.ybeam, #not sure if matrix needs to be transposed first for outputting HKL's???: #params.refinement.reflections.do_outlier_rejection=True #params.refinement.reflections.iqr_multiplier=0.5 #params.refinement.reflections.minimum_sample_size=50 #params.refinement.reflections.maximum_sample_size=50 #params.refinement.reflections.random_seed=1 # fix detector rotz, allow distance to refine # fix detector rotz and distance # fix detector rotz and distance #for item in history["parameter_vector"]: #  print ["%8.5f"%a for a in item]