Indexing.h

/*  ObjCryst++ Object-Oriented Crystallographic Library
    (c) 2006- Vincent Favre-Nicolin vincefn@users.sourceforge.net

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/
/*
*  source file for Indexing classes & functions
*
*/
#ifndef _OBJCRYST_INDEXING_H_
#define _OBJCRYST_INDEXING_H_

#include <math.h>
#include <fstream>
#include <iostream>
#include <vector>
#include <list>
#include <time.h>
#include "ObjCryst/RefinableObj/RefinableObj.h"
#include "ObjCryst/RefinableObj/LSQNumObj.h"
namespace ObjCryst
{
enum CrystalSystem
{ TRICLINIC, MONOCLINIC, ORTHOROMBIC, HEXAGONAL, RHOMBOEDRAL, TETRAGONAL, CUBIC};

enum CrystalCentering
{ LATTICE_P,LATTICE_I,LATTICE_A,LATTICE_B,LATTICE_C,LATTICE_F};

float EstimateCellVolume(const float dmin, const float dmax, const float nbrefl,
                         const CrystalSystem system,const CrystalCentering centering,const float kappa=1);

class RecUnitCell
{
   public:
      RecUnitCell(const float zero=0,const float par0=0,const float par1=0,const float par2=0,
                  const float par3=0,const float par4=0,const float par5=0,CrystalSystem lattice=CUBIC,
                  const CrystalCentering cent=LATTICE_P);
      RecUnitCell(const RecUnitCell &old);
      void operator=(const RecUnitCell &rhs);
      // access to ith parameter
      //float operator[](const unsigned int i);
      float hkl2d(const float h,const float k,const float l,REAL *derivpar=NULL,const unsigned int derivhkl=0) const;
      void hkl2d_delta(const float h,const float k,const float l,const RecUnitCell &delta, float & dmin, float &dmax) const;
      std::vector<float> DirectUnitCell(const bool equiv=false)const;
      REAL par[7];
      float zero;
      CrystalSystem mlattice;
      CrystalCentering mCentering;
};

class PeakList
{
   public:
      PeakList();
      PeakList(const PeakList &old);
      void operator=(const PeakList &rhs);
      ~PeakList();
      void ImportDhklDSigmaIntensity(std::istream &is,float defaultsigma=.001);
      void ImportDhklIntensity(std::istream &is);
      void ImportDhkl(std::istream &is);
      void Import2ThetaIntensity(std::istream &is, const float wavelength=1.5418);
      float Simulate(float zero, float a, float b, float c,
                    float alpha, float beta, float gamma,
                    bool deg, unsigned int nb=20, unsigned int nbspurious=0,
                    float sigma=0, float percentMissing=0, const bool verbose=false);
      void ExportDhklDSigmaIntensity(std::ostream &out)const;
      void AddPeak(const float d, const float iobs=1.0,const float dobssigma=0.0,const float iobssigma=0.0,
                   const int h=0,const int k=0, const int l=0,const float d2calc=0);
      void RemovePeak(unsigned int i);
      void Print(std::ostream &os) const;
      struct hkl0
      {
         hkl0(const int h=0,const int k=0, const int l=0);
         int h,k,l;
      };
      struct hkl
      {
         hkl(const float dobs=1.0,const float iobs=0.0,const float dobssigma=0.0,const float iobssigma=0.0,
             const int h=0,const int k=0, const int l=0,const float d2calc=0);
         float dobs;
         float dobssigma;
         float d2obs;
         float d2obsmin;
         float d2obsmax;
         float iobs;
         float iobssigma;
         mutable int h,k,l;
         mutable std::list<hkl0> vDicVolHKL;
         mutable bool isIndexed;
         mutable bool isSpurious;
         mutable unsigned long stats;
         mutable float d2calc;
         mutable float d2diff;
      };
      vector<hkl> & GetPeakList();
      const vector<hkl> & GetPeakList()const ;
      mutable vector<hkl>mvHKL;
      mutable list<hkl> mvPredictedHKL;
};

float Score(const PeakList &dhkl, const RecUnitCell &ruc, const unsigned int nbSpurious=0,
            const bool verbose=false,const bool storehkl=false,
            const bool storePredictedHKL=false);

class CellExplorer:public RefinableObj
{
   public:
      CellExplorer(const PeakList &dhkl, const CrystalSystem lattice, const unsigned int nbSpurious);
      void Evolution(unsigned int ng,const bool randomize=true,const float f=0.7,const float cr=0.5,unsigned int np=100);
      void SetLengthMinMax(const float min,const float max);
      void SetAngleMinMax(const float min,const float max);
      void SetVolumeMinMax(const float min,const float max);
      void SetNbSpurious(const unsigned int nb);
      void SetD2Error(const float err);
      void SetMinMaxZeroShift(const float min,const float max);
      void SetCrystalSystem(const CrystalSystem system);
      void SetCrystalCentering(const CrystalCentering cent);
      virtual const string& GetClassName() const;
      virtual const string& GetName() const;
      virtual void Print() const;
      virtual unsigned int GetNbLSQFunction() const;
      virtual const CrystVector_REAL& GetLSQCalc(const unsigned int) const;
      virtual const CrystVector_REAL & GetLSQObs(const unsigned int) const;
      virtual const CrystVector_REAL & GetLSQWeight(const unsigned int) const;
      virtual const CrystVector_REAL & GetLSQDeriv(const unsigned int, RefinablePar &);
      virtual void BeginOptimization(const bool allowApproximations=false, const bool enableRestraints=false);
      void LSQRefine(int nbCycle=1, bool useLevenbergMarquardt=true, const bool silent=false);
      void DicVol(const float minScore=10,const unsigned int minDepth=3,const float stopOnScore=50.0,const unsigned int stopOnDepth=6);
      void ReduceSolutions(const bool updateReportThreshold=false);
      float GetBestScore()const;
      const std::list<std::pair<RecUnitCell,float> >& GetSolutions()const;
      std::list<std::pair<RecUnitCell,float> >& GetSolutions();
   private:
      unsigned int RDicVol(RecUnitCell uc0, RecUnitCell uc1, unsigned int depth,unsigned long &nbCalc,const float minV,const float maxV,vector<unsigned int> vdepth=vector<unsigned int>());
      void Init();
      std::list<std::pair<RecUnitCell,float> > mvSolution;
      unsigned int mnpar;
      const PeakList *mpPeakList;
      float mLengthMin,mLengthMax;
      float mAngleMin,mAngleMax;
      float mVolumeMin,mVolumeMax;
      float mZeroShiftMin,mZeroShiftMax;
      float mMin[7];
      float mAmp[7];
      CrystalSystem mlattice;
      CrystalCentering mCentering;
      unsigned int mNbSpurious;
      float mD2Error;
      LSQNumObj mLSQObj;
      mutable CrystVector_REAL mObs;
      mutable CrystVector_REAL mCalc;
      mutable CrystVector_REAL mWeight;
      mutable CrystVector_REAL mDeriv;
      RecUnitCell mRecUnitCell;
      float mBestScore;
      std::vector<unsigned int> mvNbSolutionDepth;
      float mMinScoreReport;
      unsigned int mMaxDicVolDepth,mDicVolDepthReport;
      mutable float mCosAngMax;
      unsigned int mNbLSQExcept;
};


}//namespace
#endif