FOX/ObjCryst++  1.10.X (development)
1 /* ObjCryst++ Object-Oriented Crystallographic Library
2  (c) 2006- Vincent Favre-Nicolin
4  This program is free software; you can redistribute it and/or modify
5  it under the terms of the GNU General Public License as published by
6  the Free Software Foundation; either version 2 of the License, or
7  (at your option) any later version.
9  This program is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  GNU General Public License for more details.
14  You should have received a copy of the GNU General Public License
15  along with this program; if not, write to the Free Software
16  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 */
18 /*
19 * source file for Indexing classes & functions
20 *
21 */
25 #include <math.h>
26 #include <fstream>
27 #include <iostream>
28 #include <vector>
29 #include <list>
30 #include <time.h>
31 #include "ObjCryst/RefinableObj/RefinableObj.h"
32 #include "ObjCryst/RefinableObj/LSQNumObj.h"
33 namespace ObjCryst
34 {
41 enum CrystalCentering
51 float EstimateCellVolume(const float dmin, const float dmax, const float nbrefl,
52  const CrystalSystem system,const CrystalCentering centering,const float kappa=1);
59 {
60  public:
61  RecUnitCell(const float zero=0,const float par0=0,const float par1=0,const float par2=0,
62  const float par3=0,const float par4=0,const float par5=0,CrystalSystem lattice=CUBIC,
63  const CrystalCentering cent=LATTICE_P);
64  RecUnitCell(const RecUnitCell &old);
65  void operator=(const RecUnitCell &rhs);
66  // access to ith parameter
67  //float operator[](const unsigned int i);
72  float hkl2d(const float h,const float k,const float l,REAL *derivpar=NULL,const unsigned int derivhkl=0) const;
78  void hkl2d_delta(const float h,const float k,const float l,const RecUnitCell &delta, float & dmin, float &dmax) const;
85  std::vector<float> DirectUnitCell(const bool equiv=false)const;
104  REAL par[7];
105  float zero;
106  CrystalSystem mlattice;
107  CrystalCentering mCentering;
108 };
114 class PeakList
115 {
116  public:
117  PeakList();
118  PeakList(const PeakList &old);
119  void operator=(const PeakList &rhs);
120  ~PeakList();
121  void ImportDhklDSigmaIntensity(std::istream &is,float defaultsigma=.001);
122  void ImportDhklIntensity(std::istream &is);
123  void ImportDhkl(std::istream &is);
124  void Import2ThetaIntensity(std::istream &is, const float wavelength=1.5418);
137  float Simulate(float zero, float a, float b, float c,
138  float alpha, float beta, float gamma,
139  bool deg, unsigned int nb=20, unsigned int nbspurious=0,
140  float sigma=0, float percentMissing=0, const bool verbose=false);
141  void ExportDhklDSigmaIntensity(std::ostream &out)const;
144  void AddPeak(const float d, const float iobs=1.0,const float dobssigma=0.0,const float iobssigma=0.0,
145  const int h=0,const int k=0, const int l=0,const float d2calc=0);
146  void RemovePeak(unsigned int i);
147  void Print(std::ostream &os) const;
149  struct hkl0
150  {
151  hkl0(const int h=0,const int k=0, const int l=0);
153  int h,k,l;
154  };
158  struct hkl
159  {
160  hkl(const float dobs=1.0,const float iobs=0.0,const float dobssigma=0.0,const float iobssigma=0.0,
161  const int h=0,const int k=0, const int l=0,const float d2calc=0);
163  float dobs;
165  float dobssigma;
167  float d2obs;
169  float d2obsmin;
171  float d2obsmax;
173  float iobs;
175  float iobssigma;
177  mutable int h,k,l;
179  mutable std::list<hkl0> vDicVolHKL;
181  mutable bool isIndexed;
183  mutable bool isSpurious;
185  mutable unsigned long stats;
187  mutable float d2calc;
189  mutable float d2diff;
190  };
192  vector<hkl> & GetPeakList();
194  const vector<hkl> & GetPeakList()const ;
198  mutable vector<hkl>mvHKL;
201  mutable list<hkl> mvPredictedHKL;
202 };
205 float Score(const PeakList &dhkl, const RecUnitCell &ruc, const unsigned int nbSpurious=0,
206  const bool verbose=false,const bool storehkl=false,
207  const bool storePredictedHKL=false);
213 {
214  public:
215  CellExplorer(const PeakList &dhkl, const CrystalSystem lattice, const unsigned int nbSpurious);
216  void Evolution(unsigned int ng,const bool randomize=true,const float f=0.7,const float cr=0.5,unsigned int np=100);
217  void SetLengthMinMax(const float min,const float max);
218  void SetAngleMinMax(const float min,const float max);
219  void SetVolumeMinMax(const float min,const float max);
220  void SetNbSpurious(const unsigned int nb);
222  void SetD2Error(const float err);
223  void SetMinMaxZeroShift(const float min,const float max);
224  void SetCrystalSystem(const CrystalSystem system);
225  void SetCrystalCentering(const CrystalCentering cent);
226  virtual const string& GetClassName() const;
227  virtual const string& GetName() const;
228  virtual void Print() const;
229  virtual unsigned int GetNbLSQFunction() const;
230  virtual const CrystVector_REAL& GetLSQCalc(const unsigned int) const;
231  virtual const CrystVector_REAL & GetLSQObs(const unsigned int) const;
232  virtual const CrystVector_REAL & GetLSQWeight(const unsigned int) const;
233  virtual const CrystVector_REAL & GetLSQDeriv(const unsigned int, RefinablePar &);
234  virtual void BeginOptimization(const bool allowApproximations=false, const bool enableRestraints=false);
235  void LSQRefine(int nbCycle=1, bool useLevenbergMarquardt=true, const bool silent=false);
241  void DicVol(const float minScore=10,const unsigned int minDepth=3,const float stopOnScore=50.0,const unsigned int stopOnDepth=6);
247  void ReduceSolutions(const bool updateReportThreshold=false);
248  float GetBestScore()const;
249  const std::list<std::pair<RecUnitCell,float> >& GetSolutions()const;
250  std::list<std::pair<RecUnitCell,float> >& GetSolutions();
251  private:
252  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>());
253  void Init();
255  std::list<std::pair<RecUnitCell,float> > mvSolution;
256  unsigned int mnpar;
257  const PeakList *mpPeakList;
258  float mLengthMin,mLengthMax;
259  float mAngleMin,mAngleMax;
260  float mVolumeMin,mVolumeMax;
261  float mZeroShiftMin,mZeroShiftMax;
263  float mMin[7];
266  float mAmp[7];
270  CrystalCentering mCentering;
271  unsigned int mNbSpurious;
272  float mD2Error;
273  LSQNumObj mLSQObj;
274  mutable CrystVector_REAL mObs;
275  mutable CrystVector_REAL mCalc;
276  mutable CrystVector_REAL mWeight;
277  mutable CrystVector_REAL mDeriv;
281  float mBestScore;
283  std::vector<unsigned int> mvNbSolutionDepth;
284  float mMinScoreReport;
285  unsigned int mMaxDicVolDepth,mDicVolDepthReport;
288  mutable float mCosAngMax;
290  unsigned int mNbLSQExcept;
291 };
294 }//namespace
295 #endif
unsigned int mNbLSQExcept
Number of exceptions caught during LSQ, in a given search - above 20 LSQ is disabled.
Definition: Indexing.h:290
vector< hkl > & GetPeakList()
Get peak list.
Definition: Indexing.cpp:930
virtual const string & GetName() const
Name of the object.
Definition: Indexing.cpp:1415
void SetD2Error(const float err)
Allowed error on 1/d (squared!), used for dicvol.
Definition: Indexing.cpp:1408
RecUnitCell mRecUnitCell
Reciprocal unit cell used for least squares refinement.
Definition: Indexing.h:279
float mCosAngMax
Stored value of cos(max ang) for tricilinic search - we do not want to recompute the cos at every dic...
Definition: Indexing.h:288
list< hkl > mvPredictedHKL
Full list of calculated HKL positions for a given solution, up to a given resolution After finding a ...
Definition: Indexing.h:201
std::list< std::pair< RecUnitCell, float > > mvSolution
Max number of obs reflections to use.
Definition: Indexing.h:255
Algorithm class to find the correct indexing from observed peak positions.
Definition: Indexing.h:212
unsigned long stats
Indexing statistics.
Definition: Indexing.h:185
float EstimateCellVolume(const float dmin, const float dmax, const float nbrefl, const CrystalSystem system, const CrystalCentering centering, const float kappa)
Estimate volume from number of peaks at a given dmin See J.
Definition: Indexing.cpp:45
Different lattice types.
Definition: Indexing.h:38
CrystalCentering mCentering
Centering type.
Definition: Indexing.h:270
float iobs
Observed peak intensity.
Definition: Indexing.h:173
float hkl2d(const float h, const float k, const float l, REAL *derivpar=NULL, const unsigned int derivhkl=0) const
Compute d*^2 for hkl reflection if deriv != -1, compute derivate versus the corresponding parameter...
Definition: Indexing.cpp:106
Generic Refinable Object.
Definition: RefinableObj.h:752
Class to store positions of observed reflections.
Definition: Indexing.h:114
bool isIndexed
Is this line indexed ?
Definition: Indexing.h:181
float Score(const PeakList &dhkl, const RecUnitCell &rpar, const unsigned int nbSpurious, const bool verbose, const bool storehkl, const bool storePredictedHKL)
Compute score for a candidate RecUnitCell and a PeakList.
Definition: Indexing.cpp:936
float d2obs
Observed peak position 1/d^2.
Definition: Indexing.h:167
float d2obsmin
Min value for observed peak position 1/(d+disgma/2)^2.
Definition: Indexing.h:169
int h
Miller indices.
Definition: Indexing.h:153
vector< hkl > mvHKL
Predict peak positions Best h,k,l for each observed peak (for least-squares refinement) This is store...
Definition: Indexing.h:198
virtual const CrystVector_REAL & GetLSQDeriv(const unsigned int, RefinablePar &)
Get the first derivative values for the LSQ function, for a given parameter.
Definition: Indexing.cpp:1451
One set of Miller indices, a possible indexation for a reflection.
Definition: Indexing.h:149
float mAmp[7]
Max amplitude (max=min+amplitude) for all parameters All parameters are treated as periodic for DE (...
Definition: Indexing.h:266
float d2calc
Calculated position, 1/d^2.
Definition: Indexing.h:187
std::list< hkl0 > vDicVolHKL
Possible Miller indices, stored during a dichotomy search.
Definition: Indexing.h:179
float d2obsmax
Min value for observed peak position 1/(d-disgma/2)^2.
Definition: Indexing.h:171
CrystalSystem mlattice
Lattice type for which we search.
Definition: Indexing.h:268
void ReduceSolutions(const bool updateReportThreshold=false)
Sort all solutions by score, remove duplicates.
Definition: Indexing.cpp:2628
virtual const CrystVector_REAL & GetLSQCalc(const unsigned int) const
Get the current calculated value for the LSQ function.
Definition: Indexing.cpp:1427
void DicVol(const float minScore=10, const unsigned int minDepth=3, const float stopOnScore=50.0, const unsigned int stopOnDepth=6)
Run DicVOl algorithm, store only solutions with score >minScore or depth>=minDepth, stop at the end of one volume interval (~400 A^3) if best score>stopOnScore, or if one solution was found at depth>=stopOnDepth.
Definition: Indexing.cpp:2126
virtual const CrystVector_REAL & GetLSQObs(const unsigned int) const
Get the observed values for the LSQ function.
Definition: Indexing.cpp:1440
REAL par[7]
The 6 parameters defining 1/d_hkl^2 = d*_hkl^2, for different crystal classes, from: d*_hkl^2 = zero ...
Definition: Indexing.h:104
float d2diff
1/d^2 difference, obs-calc
Definition: Indexing.h:189
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)
Add one peak.
Definition: Indexing.cpp:886
virtual void BeginOptimization(const bool allowApproximations=false, const bool enableRestraints=false)
This should be called by any optimization class at the begining of an optimization.
Definition: Indexing.cpp:1481
float iobssigma
Error on observed peak intensity.
Definition: Indexing.h:175
virtual unsigned int GetNbLSQFunction() const
Number of LSQ functions.
Definition: Indexing.cpp:1424
Lightweight class describing the reciprocal unit cell, for the fast computation of d*_hkl^2...
Definition: Indexing.h:58
float mBestScore
Current best score.
Definition: Indexing.h:281
float dobssigma
Error on peak position.
Definition: Indexing.h:165
One observed diffraction line, to be indexed.
Definition: Indexing.h:158
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
Definition: Indexing.cpp:1410
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)
Generate a list of simulated peak positions, from given lattice parameters.
Definition: Indexing.cpp:804
The namespace which includes all objects (crystallographic and algorithmic) in ObjCryst++.
Definition: Atom.cpp:47
float dobs
Observed peak position 1/d.
Definition: Indexing.h:163
bool isSpurious
Is this an impurity line ?
Definition: Indexing.h:183
Generic class for parameters of refinable objects.
Definition: RefinableObj.h:223
int h
Miller indices, after line is indexed.
Definition: Indexing.h:177
float mMin[7]
Min values for all parameters (7=unit cell +zero)
Definition: Indexing.h:263
std::vector< unsigned int > mvNbSolutionDepth
Number of solutions found during dicvol search, at each depth.
Definition: Indexing.h:283
void hkl2d_delta(const float h, const float k, const float l, const RecUnitCell &delta, float &dmin, float &dmax) const
Compute d*^2 for one hkl reflection: this functions computes a d*^2 range (min,max) for a given range...
Definition: Indexing.cpp:454
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)
light-weight class storing the reciprocal space unitcell
Definition: Indexing.cpp:80
(Quick & dirty) Least-Squares Refinement Object with Numerical derivatives
Definition: LSQNumObj.h:37
std::vector< float > DirectUnitCell(const bool equiv=false) const
Compute real space unit cell from reciprocal one.
Definition: Indexing.cpp:562
virtual const CrystVector_REAL & GetLSQWeight(const unsigned int) const
Get the weight values for the LSQ function.
Definition: Indexing.cpp:1445