hkl_unmerge.hh

// hkl_unmerge.hh
// Phil Evans August 2003-2004 etc
//
// Classes for unmerged hkl lists
//  at present just for the items used by Mosflm & Scala

// These data structures are primarily directed at Scala so are
// in its namespace

#ifndef HKL_UNMERGE_HEADER
#define HKL_UNMERGE_HEADER

#include <map>
#include "hkl_datatypes.hh"
#include "hkl_controls.hh"
#include "controls.hh"
#include "hkl_symmetry.hh"
#include "icering.hh"
#include "observationflags.hh"
#include "range.hh"
#include "hash.hh"
#include "runthings.hh"

namespace scala
{
  //==============================================================
  class data_flags
  // True if column is present
  {
  public:
    data_flags();  
    void print() const; // for debugging
    
    bool is_h, is_k, is_l, is_misym, is_batch,
      is_I, is_sigI, is_Ipr, is_sigIpr, is_fractioncalc,
      is_Xdet, is_Ydet, is_Rot, is_Width, is_LP, is_Mpart,
      is_ObsFlag, is_BgPkRatio, is_scale, is_sigscale, is_time;
  };  // data_flags
  //===================================================================
  // *******************   observation, reflections etc
  class SelectI 
  {
    //

  public:
    SelectI(){}
    //  Algorithm for selection as in Scala
    //   INTEGRATED   integrated intensity I                  (SelectIcolFlag=0)
    //   PROFILE      profile-fitted intensity IPR            (SelectIcolFlag=-1)
    //   COMBINE      weighted mean after combining partials  (SelectIcolFlag=+1)
    //       combine integrated intensity Iint and
    //       profile-fitted intensity Ipr as
    //         I = w * Ipr  +  (1-w) * Iint
    //       where w = 1/(1+(Iint/Imid)**IPOWER)
    //       Imid (= IcolFlag)  is the point of equal weight, and
    //       should be about the mean intensity.
    //  Ipower default = 3

    // return true if COMBINE
    static bool Combine() {return (selecticolflag > 0);}
    // If IcolFlag > 0, set selecticolflag and set imid = Imid, ie COMBINE
    // If IcolFlag = 0, select Iint, < 0 select Ipr
    static void SetIcolFlag(const int& IcolFlag, const double& Imid, const int Ipower=3);
    // Set SelectIcolFlag
    static int& SelectIcolFlag() {return selecticolflag;}
    // store imid = overall <I> if needed and not set
    static void SetAverageIntensity(const double& meanI);

    static IsigI GetCombinedI(const Rtype& Ic, const Rtype& varIc,
                              const Rtype& Ipr, const Rtype& varIpr);

    static std::string format();

  private:
    static int selecticolflag;
    static int ipowercomb;
    static double imid;
  //===================================================================
  };  // SelectI

// PartFlagSwitch
//   EMPTY  nothing stored
//   FULL   fully recorded reflection (1 part)
//   COMPLETE_CHECKED  complete, all Mpart flags consistent
//   COMPLETE passed total fraction checks
//   SCALE passed checks to be scaled
//   INCOMPLETE
//   EXCLUDED
 enum PartFlagSwitch
   {EMPTY, FULL, COMPLETE_CHECKED, COMPLETE, SCALE, INCOMPLETE, EXCLUDED};
  //===================================================================
  class observation_part
  // hkl is the reduced indices in the "current" spacegroup
  // run is set by call from hkl_unmerge_list::organise
  {
  public:

    // constructors
    observation_part();  // don't use 
    observation_part(const Hkl& hkl_in,
                     const int& isym_in, const int& batch_in,
                     const Rtype& I_in, const Rtype& sigI_in,
                     const Rtype& I_pr, const Rtype& sigIpr_in,
                     const Rtype& Xdet_in, const Rtype& Ydet_in,
                     const Rtype& phi_in, const Rtype& time_in,
                     const Rtype& fraction_calc_in, const Rtype& width_in,
                     const Rtype& LP_in,
                     const int& Npart_in, const int& Ipart_in,
                     const ObservationFlag& ObsFlag_in);
    // Npart is number of parts, from input (MPART)
    //       = 1 for full, = -1 unknown but partial
    // Ipart serial number in parts, from input (MPART)
    //       = 1 for full or unknown

    void set_run(const int& run) {run_ = run;}  


    inline  IsigI I_sigI() const {return IsigI(I_, sigI_);}
    inline  Rtype Ic() const {return I_;}       // actual I column
    inline  Rtype sigIc() const {return sigI_;} 
    inline  IsigI I_sigIpr() const {return IsigI(Ipr_, sigIpr_);}
    inline  Rtype Ipr() const {return Ipr_;}       //  Ipr column
    inline  Rtype sigIpr() const {return sigIpr_;} 

    void StoreIsigI(const IsigI& Is) {I_ = Is.I(); sigI_ = Is.sigI();}
    void StoreIsigIpr(const IsigI& Ipr) {Ipr_ = Ipr.I(); sigIpr_ = Ipr.sigI();}

    // selected I column
    // selIcolFlag = 0 Ic, = -1 Ipr
    inline  Rtype Ic(const int& selIcolFlag) const
    {return (selIcolFlag == 0) ? I_ : Ipr_;}
    inline  Rtype sigIc(const int& selIcolFlag) const
    {return (selIcolFlag == 0) ? sigI_ : sigIpr_;}


    inline  Rtype Xdet() const {return Xdet_;}
    inline  Rtype Ydet() const {return Ydet_;}
    inline  Rtype phi()  const {return phi_;}
    inline  Rtype time() const {return time_;}
    inline  Rtype fraction_calc() const {return fraction_calc_;}
    inline  Rtype width() const {return width_;}
    inline  Rtype LP() const {return LP_;}   // divide by this to get raw intensity

    inline  int isym() const {return isym_;}
    inline  int batch() const {return batch_;}
    inline  int Npart() const {return Npart_;}
    inline  int Ipart() const {return Ipart_;}
    inline  ObservationFlag ObsFlag() const {return ObsFlag_;}
    inline  int run() const {return run_;}
    inline  Hkl hkl()  const {return hkl_;};

    void set_isym(const int& isym) {isym_ = isym;}  // set isym
    void set_hkl(const Hkl& hkl) {hkl_ = hkl;} // set hkl
    void set_batch(const int& Batch) {batch_ = Batch;} // set batch number
    void set_phi(const Rtype& Phi) {phi_ = Phi;} // set Phi
    void offset_phi(const Rtype& Offset) {phi_ += Offset;} // offset Phi
    void offset_time(const Rtype& Offset) {time_ += Offset;} // offset Time
    void negate_time() {time_ = -time_;} // negate Time

  private:
    Hkl hkl_;
    int isym_, batch_;
    Rtype I_, sigI_;           // from column I
    Rtype Ipr_, sigIpr_;           // from column Ipr (if present)
    Rtype Xdet_, Ydet_, phi_, time_;
    Rtype fraction_calc_, width_, LP_;
    int  Npart_, Ipart_;
    ObservationFlag ObsFlag_;
    int run_;
  };  // class observation_part

  //===================================================================
  class observation
  {
    // 
  public:
    observation();
    observation(const Hkl hkl_in,
                const int& isym_in,
                const int& run_in,
                const int& datasetIndex_in,
                const int& Npart_in,
                observation_part ** const part1_in, 
                const Rtype& TotFrac,
                const PartFlagSwitch& partialstatus_in,
                const ObservationFlag& obsflag_in);

    // Accessors
    Rtype I() const {return I_;}    
    Rtype sigI() const {return sigI_;} 
    IsigI I_sigI() const {return IsigI(I_,sigI_);}    

    Rtype kI() const {return I_/gscale;}    
    Rtype ksigI() const {return sigI_/gscale;} 
    IsigI kI_sigI() const {return IsigI(I_/gscale,sigI_/gscale);}    


    Rtype phi() const {return phi_;}  
    Rtype time() const {return time_;} 
    int Isym() const {return isym_;}  
    Rtype LP() const {return LP_;}

    int run() const {return run_;};

    int datasetIndex() const {return datasetIndex_;}    
    bool IsAccepted() const {return obs_status.IsAccepted();} 

    int num_parts() const;  
    observation_part get_part(const int& kpart) const; 
    void replace_part(const int& kpart, const observation_part& obs_part); 

    Hkl hkl_original() const {return hkl_original_;} 
    bool IsFull() const {return part_flag == FULL;} 
    Rtype Gscale() const {return gscale;}  
    std::pair<float,float> XYdet() const; 
    float TotalFraction() const {return totalfraction;} 
    int Batch() const;  

    // Partial status flag
    PartFlagSwitch PartFlag() const {return part_flag;} 

    bool set_IsigI_phi_time(const Rtype& I, const Rtype& sigI,
                            const Rtype& phi, const Rtype& time,
                            const Rtype& LP);
    void set_batch(const int& batch) {batch_ = batch;}
    void set_sigI(Rtype& sigI) {sigI_ = sigI;}
    void SetGscale(const Rtype& g) {gscale = g;}
    void StoreS2(const float& Thetap, const float& Phip)
    {thetap=Thetap; phip=Phip;}
    void GetS2(double& Thetap, double& Phip) const 
    {Thetap=thetap; Phip=phip;}
    void StoreS(const FVect3& dStarvec) {s_dif = dStarvec;}
    FVect3 GetS() const {return s_dif;}

    void ResetObsAccept(ObservationFlagControl& ObsFlagControl);        
    ObservationFlag Observationflag() const {return obs_flag;}

    ObservationStatus ObsStatus() const {return obs_status;}
    void UpdateStatus(const ObservationStatus& Status) {obs_status = Status;}

  private:
    Hkl hkl_original_;
    int isym_;
    int run_;
    int datasetIndex_;
    int Npart_;
    observation_part ** part1;
    int batch_;  // batch number of "central" batch
    Rtype totalfraction;
    PartFlagSwitch part_flag;
    ObservationFlag obs_flag;
    Rtype gscale;         // inverse scale g
    Rtype I_, sigI_;
    Rtype phi_;
    Rtype time_;
    Rtype LP_;
    ObservationStatus obs_status;
    Rtype thetap, phip;   // secondary beam direction polar angles
    FVect3 s_dif;             // diffraction vector at Phi setting, 1/A units
  }; // class observation

  //===================================================================
  class reflection
  {
  public:
    reflection();  


    reflection(const Hkl& hkl, const int& index_obs1, const Dtype& s);
    reflection(const reflection& refl); 

    reflection& operator= (const reflection& refl); 

    // Storage
    void store_last_index(const int& index_obs2);
    inline int first_index() const {return index_first_obs_;} 
    inline int last_index() const {return index_last_obs_;} 
    inline Hkl hkl() const {return hkl_reduced_;} 

    void add_observation_list(const std::vector<observation>& obs_in);


    Rtype sum_partials(int& Nfull, int& Npart, int& Nscaled);

    // Retrieval
    int num_observations() const; 
    int NvalidObservations() const; 
    observation get_observation(const int& lobs) const; 

    int next_observation(observation& obs) const; 
    void reset() const {NextObs = -1;} // (NextObs is mutable)
    void replace_observation(const observation& obs);
    void replace_observation(const observation& obs, const int& lobs);
    // Count valid observations
    void CountNValid();

    Rtype invresolsq() const {return invresolsq_;} 

    // Counts observations reclassified in ObsFlagControl
    void ResetObsAccept(ObservationFlagControl& ObsFlagControl);

    void SetStatus(const int& istat) {statusflag = istat;} 
    int Status() const {return statusflag;} 

  private:
    Hkl hkl_reduced_; // reduced indices
    std::vector<observation> observations;
    int index_first_obs_, index_last_obs_;
    Rtype invresolsq_;
    mutable int NextObs;
    int NvalidObs;
    int statusflag;  // 0 OK accept, !=0  reject
  }; //class reflection

  //===================================================================
  class hkl_unmerge_list

  {
  public:
    hkl_unmerge_list(): status(EMPTY) {};
    void init(const std::string& Title,
              const int& NreflReserve, 
              const hkl_symmetry& symmetry,
              const all_controls& controls,
              const std::vector<Xdataset>& DataSets,
              const std::vector<Batch>& Batches);
    void init(const std::string& Title,
              const int& NreflReserve, 
              const hkl_symmetry& symmetry,
              const all_controls& controls);
    void clear();

    bool IsEmpty() const {return (status == EMPTY);}

    bool IsReady() const {return (status == SUMMED);}

    void StoreDatasetBatch(const std::vector<Xdataset>& DataSets,
                           const std::vector<Batch>& Batches);


    void AddDatasetBatch(const std::vector<Xdataset>& DataSets,
                                           const std::vector<Batch>& Batches);


    int append(const hkl_unmerge_list& OtherList);

    //----
    hkl_unmerge_list(const hkl_unmerge_list& List);
    hkl_unmerge_list& operator= (const hkl_unmerge_list& List);


    // Set controls, limits, etc ----------
    void SetResoLimits(const float& LowReso, const float& HighReso); 

    void ResetResoLimits();
    void SetIceRings(const Rings& rings);

    // list and flag observations which are outside these limits
    void ImposeResoByRunLimits();

    // Retrieve information  -------------------------------
    int num_reflections() const;
    int num_reflections_valid() const {return Nref_valid;}
    int num_reflections_icering() const {return Nref_icering;}

    int num_observations() const;
    int num_observations_full() const {return Nobs_full;}
    int num_observations_part() const {return Nobs_partial;}
    int num_observations_scaled() const {return Nobs_scaled;}
    int num_observations_rejected_FracTooSmall() const
    {return partial_flags.NrejFractionTooSmall();}
    int num_observations_rejected_FracTooLarge() const
    {return partial_flags.NrejFractionTooLarge();}
    int num_observations_rejected_Gap() const {return partial_flags.NrejGap();}
    Rtype MinSigma() const {return sigmamin;}
    Scell cell(const PxdName& PXDsetName = PxdName()) const;
    Scell Cell() const {return cell(PxdName());}  

    hkl_symmetry symmetry() const {return refl_symm;}

    void SetMtzSym(const CMtz::SYMGRP& Mtzsym) {mtzsym = Mtzsym;}
    CMtz::SYMGRP MtzSym() const {return mtzsym;}

    ReindexOp TotalReindex() const {return totalreindex;}
    // Return resolution
    ResoRange ResRange() const {return ResolutionRange;}      
    Range Srange() const {return ResolutionRange.SRange();}   
    Range RRange() const {return ResolutionRange.RRange();}   
    ResoRange ResLimRange() const {return ResoLimRange;}      
    Rtype DstarMax() const;                                   
    // Return dataset stuff
    int num_datasets() const {return ndatasets;} 
    Xdataset xdataset(const int& jset) const {return datasets.at(jset);}  
    std::vector<Xdataset> AllXdatasets() const {return datasets;} 
    // Return batch stuff
    int num_batches() const {return nbatches;} 
    std::vector<Batch> Batches() const {return batches;}  

    Batch batch(const int& jbat) const {return batches.at(jbat);}
    int batch_serial(const int& batch) const {return batch_lookup.lookup(batch);}

    // Return batch serial number for batch batchnum, or if this one is
    // not present, search upwards until one is found or maxbatchnum is reached.
    // Return -1 if nothing found
    int NextBatchSerial(const int& batchnum, const int& maxbatchnum) const;
    // Return batch serial number for batch batchnum, or if this one is
    // not present, search backwards until one is found or 0 is reached.
    // Return -1 if nothing found
    int LastBatchSerial(const int& batchnum) const;

    // Run stuff
    int num_runs() const {return runlist.size();} 
    std::vector<Run> RunList() const {return runlist;}  

    void OffsetBatchNumbers(const std::vector<int>& runOffsets);

    void RejectBatch(const int& batch);

    void RejectBatchSerial(const int& jbat);

    int PurgeRejectedBatches();

    int EliminateBatches(const std::vector<int> RejectedBatches);

    std::string Title() const {return FileTitle;}

    void AppendFileName(const std::string& Name);
    std::string Filename() const {return filename;} 

    reflection get_reflection(const int& lref) const;

    int get_reflection(reflection& refl, const Hkl& hkl) const;
    int next_reflection(reflection& refl) const;
    void replace_reflection(const reflection& refl);

    void replace_observation(const observation& obs, const int& lobs);


    void rewind() const;

    observation_part& find_part(const int& i) const;
    int num_parts() const {return N_part_list;} // length of part list

    // Do things             -------------------------------


    int change_symmetry(const hkl_symmetry& new_symm,
                         const ReindexOp& reindex_op,
                         const bool& AllowFractIndex = false);


    int prepare();

    int sum_partials();


    void CalcSecondaryBeams(const int& pole);

    void SetPoles(const int& pole);


    std::pair<float, float> CalcSecondaryBeamPolar 
    (const int& batchNum, const Hkl& hkl_original, const float& phi,
     DVect3& sPhi) const;


    DVect3 CalcSecondaryBeam
    (const int& batchNum, const Hkl& hkl_original, const float& phi,
     DVect3& sPhi) const;
      

    void ResetObsAccept(ObservationFlagControl& ObsFlagControl);
    void ResetReflAccept();

    void store_part(const Hkl& hkl_in,
                    const int& isym_in, const int& batch_in,
                    const Rtype& I_in, const Rtype& sigI_in,
                    const Rtype& Ipr, const Rtype& sigIpr,
                    const Rtype& Xdet_in, const Rtype& Ydet_in,
                    const Rtype& phi_in, const Rtype& time_in,
                    const Rtype& fraction_calc_in, const Rtype& width_in,
                    const Rtype& LP_in, 
                    const int& Npart_in, const int& Ipart_in,
                    const ObservationFlag& ObsFlag_in);
    int close_part_list(const ResoRange& RRange,
                        const bool& Sorted);  // returns number of observations

    void StoreDataFlags(const data_flags& DataFlags) {dataflags = DataFlags;}
    data_flags DataFlags() const {return dataflags;}

    // dump given reflection, for debugging
    void dump_reflection(const Hkl& hkl) const;

    // Set runs
    void AutoSetRun();

    // ......................................................
  private:
    void initialise(const int NreflReserve, 
                    const hkl_symmetry& symmetry);

    // organise observations into reflections
    int organise();   // returns number of reflections
    // assemble partials into observations, returns number of observations
    int partials();
    void MakeHklLookup() const;

    void SetBatchList();
    void AverageBatchData();

    // Status:-
    //   EMPTY          initial state
    //   RAWLIST        raw observation list read in
    //   SORTED         sorted
    //   ORGANISED      organised into reflections
    //   PREPARED       reflections split into observations
    //                   (ie partials checked and assembled,
    //                   but not summed)
    //   SUMMED         partials summed, reflection.observations valid
    enum rfl_status
      {EMPTY, RAWLIST, SORTED, ORGANISED, PREPARED, SUMMED};
    rfl_status status;
    bool ChangeIndex;
    ReindexOp totalreindex;   // cumulative from original HKLIN setting

    std::vector <observation_part> obs_part_list;
    std::vector <observation_part *> obs_part_pointer; 
    std::vector <reflection> refl_list;
    Rtype sigmamin;
    int N_part_list;
    int Nref;
    int Nref_valid;

    int Nobservations;
    int Nobs_full;
    int Nobs_partial;
    int Nobs_scaled;
    // hkl index list: only generated if needed
    mutable clipper::HKL_lookup hkl_lookup;
    mutable bool is_hkl_lookup;   // initially false, true if hkl_lookup has been generated

    mutable int Nref_icering;

    mutable int NextRefNum;  // for next_reflection

    hkl_symmetry refl_symm;
    CMtz::SYMGRP mtzsym;     // Mtz-style symmetry, may be null
                             //  null is mtzsym.spcgrp = -1

    // Average cell (over all datasets)
    Scell averagecell;
    ResoRange ResolutionRange; // low, high, in stored array
    ResoRange ResoLimRange;    // resolution range limits
                               // reflections outside this range will
                               // not be "accept"ed

    // Controls
    int run_set;   // has run be set in part list?
                   // = -1 no controls, = 0 not set, = +1 set
    run_controls run_flags;

    bool partial_set;  // Has partial_flags been set?
    partial_controls partial_flags;

    // Ice rings
    Rings Icerings;

    // Datasets & batches
    std::vector<Xdataset> datasets;
    int ndatasets;
    std::vector<Batch> batches;      // list of batches
    int nbatches;
    hash_table batch_lookup;

    // IsPhiOffset initially false, reset in AutoSetRun
    // true if batch list set up & offset needed
    // Used in store_part to apply offset to data as it is stored,
    // or if StoreDatasetBatch (which calls AutoSetRun) is called
    // after storing the data, then offsets are applied in AutoSetRun
    // Reset to false in close_part_list
    bool IsPhiOffset;

    // Runs
    std::vector<Run> runlist;

    std::string FileTitle;
    std::string filename;

    // Flags for which data items are actually present
    data_flags dataflags;

    // ****  Private member functions
    // sort part list
    void sort();
    // for the organise method
    void add_refl(const Hkl& hkl_red, const int& index, const Dtype& s);
    void end_refl(const int& index);


    void set_run();   // set all runs in part list

    bool accept() const;    // true if reflection(NextRefNum) is acceptable
                       // also counts ice ring

    void SetUpRuns();
    void CheckAllRuns();
    void SetRunInput();

    // Check that all batches in each run have the same or similar orientations
    // Return false if different
    // Store orientation in all runs
    bool StoreOrientationRun();

    // Are new datasets the same as any old ones?
    // For each dataset from other list, store equivalent dataset
    // index in present list, if it is the same dataset
    void MergeDatasetLists(const std::vector<Xdataset>& otherDatasets,
                           const std::vector<Batch>& otherBatches);

    struct ComparePartOrder
    // This construct seems to be a way of getting pointer-to-function
    // into the argument for sort. Copied from the Web.
    {
        bool operator()(const observation_part * part1,
                        const observation_part * part2);
    };
  };  // class hkl_unmerge_list
} // end namespace scala

#endif
// end hkl_unmerge

---