FOX/ObjCryst++  2022
Indexing.cpp
1 /* ObjCryst++ Object-Oriented Crystallographic Library
2  (c) 2006- Vincent Favre-Nicolin vincefn@users.sourceforge.net
3 
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.
8 
9  This program is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  GNU General Public License for more details.
13 
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 */
22 #include <algorithm>
23 #include <iomanip>
24 
25 #include "ObjCryst/ObjCryst/Indexing.h"
26 #include "ObjCryst/Quirks/VFNDebug.h"
27 #include "ObjCryst/Quirks/VFNStreamFormat.h"
28 #include "ObjCryst/Quirks/Chronometer.h"
29 
30 using namespace std;
31 
32 #ifndef M_PI
33 #define M_PI 3.14159265358979323846264338327950288
34 #endif
35 
36 #ifndef DEG2RAD
37 #define DEG2RAD (M_PI/180.)
38 #endif
39 #ifndef RAD2DEG
40 #define RAD2DEG (180./M_PI)
41 #endif
42 
43 namespace ObjCryst
44 {
45 
46 float EstimateCellVolume(const float dmin, const float dmax, const float nbrefl,
47  const CrystalSystem system,const CrystalCentering centering,const float kappa)
48 {
49  const float q1=dmin*dmin*dmin-dmax*dmax*dmax;
50  const float q2=dmin*dmin -dmax*dmax;
51  float D0,C0;
52  if(system==TRICLINIC)
53  {
54  C0=2.095;
55  return nbrefl/(C0*kappa*q1);
56  }
57  if(system==CUBIC)
58  {
59  if(centering==LATTICE_P) D0=0.862;
60  if(centering==LATTICE_I) D0=0.475;
61  if(centering==LATTICE_F) D0=0.354;
62  return pow(nbrefl/(D0*kappa*q2),(float)1.5);
63  }
64  // "*.85" means using D0_min rather than D0
65  if(system==MONOCLINIC) {C0=1.047;D0=0.786*.85;}
66  if(system==ORTHORHOMBIC){C0=0.524;D0=1.36 *.85;}
67  if(system==HEXAGONAL) {C0=0.150;D0=1.04 *.85;}
68  if(system==RHOMBOEDRAL){C0=0.230;D0=1.04 *.85;}
69  if(system==TETRAGONAL) {C0=0.214;D0=1.25 *.85;}
70  if((centering==LATTICE_I)||(centering==LATTICE_A)||(centering==LATTICE_B)||(centering==LATTICE_C)) {C0/=2;D0/=2;}
71  if(centering==LATTICE_F){C0/=4;D0/=4;}
72  const float alpha=D0*q2/(3*C0*q1), beta=nbrefl/(2*kappa*C0*q1);
73  const float eta=beta-alpha*alpha*alpha,gamma=sqrt(beta*beta-2*beta*alpha*alpha*alpha);
74  const float v=pow(pow(eta+gamma,(float)(1/3.))+pow(fabs(eta-gamma),(float)(1/3.))-alpha,(float)3);
75  return v;
76 }
77 
81 RecUnitCell::RecUnitCell(const float zero,const float p0,const float p1,const float p2,
82  const float p3,const float p4,const float p5,CrystalSystem lattice,
83  const CrystalCentering cent, const unsigned int nbspurious):
84 mlattice(lattice),mCentering(cent),mNbSpurious(nbspurious)
85 {
86  this->par[0]=zero;
87  this->par[1]=p0;
88  this->par[2]=p1;
89  this->par[3]=p2;
90  this->par[4]=p3;
91  this->par[5]=p4;
92  this->par[6]=p5;
93 }
94 
95 RecUnitCell::RecUnitCell(const RecUnitCell &old)
96 {
97  *this=old;
98 }
99 
100 void RecUnitCell::operator=(const RecUnitCell &rhs)
101 {
102  for(unsigned int i=0;i<7;++i) par[i]=rhs.par[i];
103  mlattice=rhs.mlattice;
104  mCentering=rhs.mCentering;
105  mNbSpurious=rhs.mNbSpurious;
106 }
107 
108 float RecUnitCell::hkl2d(const float h,const float k,const float l,REAL *derivpar,const unsigned int derivhkl) const
109 {
110  if((derivpar==NULL)&&(derivhkl==0))
111  {
112  switch(mlattice)
113  {
114  case TRICLINIC:
115  {
116  return par[0]+par[1]*h*h + par[2]*k*k + par[3]*l*l + par[4]*h*k + par[5]*k*l + par[6]*h*l;
117  break;
118  }
119  case MONOCLINIC:
120  {
121  return par[0]+par[1]*par[1]*h*h + par[2]*par[2]*k*k + par[3]*par[3]*l*l + 2*par[1]*par[3]*par[4]*h*l;
122  break;
123  }
124  case ORTHORHOMBIC:
125  {
126  return par[0]+par[1]*par[1]*h*h + par[2]*par[2]*k*k + par[3]*par[3]*l*l;
127  break;
128  }
129  case HEXAGONAL:
130  {
131  return par[0]+par[1]*par[1]*(h*h + k*k + h*k)+ par[2]*par[2]*l*l ;
132  break;
133  }
134  case RHOMBOEDRAL:
135  {
136  return par[0]+par[1]*par[1]*(h*h + k*k + l*l + 2*par[2]*(h*k + k*l + h*l));
137  break;
138  }
139  case TETRAGONAL:
140  {
141  return par[0]+par[1]*par[1]*(h*h + k*k) + par[2]*par[2]*l*l;
142  break;
143  }
144  case CUBIC:
145  {
146  return par[0]+par[1]*par[1]*(h*h+k*k+l*l);
147  break;
148  }
149  }
150  }
151  if(derivhkl==1)
152  {
153  switch(mlattice)
154  {
155  case TRICLINIC:
156  {
157  return 2*par[1]*h + par[4]*k + par[6]*l;
158  break;
159  }
160  case MONOCLINIC:
161  {
162  return 2*par[1]*par[1]*h + 2*par[1]*par[3]*par[4]*l;
163  break;
164  }
165  case ORTHORHOMBIC:
166  {
167  return 2*par[1]*par[1]*h;
168  break;
169  }
170  case HEXAGONAL:
171  {
172  return par[1]*par[1]*(2*h + k);
173  break;
174  }
175  case RHOMBOEDRAL:
176  {
177  return par[1]*par[1]*(2*h + 2*par[2]*(k + l));
178  break;
179  }
180  case TETRAGONAL:
181  {
182  return 2*par[1]*par[1]*h;
183  break;
184  }
185  case CUBIC:
186  {
187  return 2*par[1]*par[1]*h;
188  break;
189  }
190  }
191  }
192  if(derivhkl==2)
193  {
194  switch(mlattice)
195  {
196  case TRICLINIC:
197  {
198  return 2*par[2]*k + par[4]*h + par[5]*l;
199  break;
200  }
201  case MONOCLINIC:
202  {
203  return 2*par[2]*par[2]*k;
204  break;
205  }
206  case ORTHORHOMBIC:
207  {
208  return 2*par[2]*par[2]*k;
209  break;
210  }
211  case HEXAGONAL:
212  {
213  return par[1]*par[1]*(2*k + h);
214  break;
215  }
216  case RHOMBOEDRAL:
217  {
218  return par[1]*par[1]*(2*k + l*l + 2*par[2]*(h + l));
219  break;
220  }
221  case TETRAGONAL:
222  {
223  return 2*par[1]*par[1]*k;
224  break;
225  }
226  case CUBIC:
227  {
228  return 2*par[1]*par[1]*k;
229  break;
230  }
231  }
232  }
233  if(derivhkl==3)
234  {
235  switch(mlattice)
236  {
237  case TRICLINIC:
238  {
239  return 2*par[3]*l + par[5]*k + par[6]*h;
240  break;
241  }
242  case MONOCLINIC:
243  {
244  return 2*par[3]*par[3]*l + 2*par[1]*par[3]*par[4]*h;
245  break;
246  }
247  case ORTHORHOMBIC:
248  {
249  return 2*par[3]*par[3]*l;
250  break;
251  }
252  case HEXAGONAL:
253  {
254  return 2*par[2]*par[2]*l;
255  break;
256  }
257  case RHOMBOEDRAL:
258  {
259  return par[1]*par[1]*(2*l + 2*par[2]*(k + h));
260  break;
261  }
262  case TETRAGONAL:
263  {
264  return 2*par[2]*par[2]*l;
265  break;
266  }
267  case CUBIC:
268  {
269  return 2*par[1]*par[1]*l;
270  break;
271  }
272  }
273  }
274 
275  if(derivpar==&par[0]) return 1.0;
276  if(derivpar==(par+1))
277  {
278  switch(mlattice)
279  {
280  case TRICLINIC:
281  {
282  return h*h;
283  break;
284  }
285  case MONOCLINIC:
286  {
287  return 2*par[1]*h*h + 2*par[3]*par[4]*h*l;
288  break;
289  }
290  case ORTHORHOMBIC:
291  {
292  return 2*par[1]*h*h;
293  break;
294  }
295  case HEXAGONAL:
296  {
297  return 2*par[1]*(h*h + k*k + h*k);
298  break;
299  }
300  case RHOMBOEDRAL:
301  {
302  return 2*par[1]*(h*h + k*k + l*l + 2*par[2]*(h*k + k*l + h*l));
303  break;
304  }
305  case TETRAGONAL:
306  {
307  return 2*par[1]*(h*h + k*k);
308  break;
309  }
310  case CUBIC:
311  {
312  return 2*par[1]*(h*h+k*k+l*l);
313  break;
314  }
315  }
316  }
317  if(derivpar==(par+2))
318  {
319  switch(mlattice)
320  {
321  case TRICLINIC:
322  {
323  return k*k;
324  break;
325  }
326  case MONOCLINIC:
327  {
328  return 2*par[2]*k*k;
329  break;
330  }
331  case ORTHORHOMBIC:
332  {
333  return 2*par[2]*k*k;
334  break;
335  }
336  case HEXAGONAL:
337  {
338  return 2*par[2]*l*l ;
339  break;
340  }
341  case RHOMBOEDRAL:
342  {
343  return par[1]*par[1]*(h*h + k*k + l*l + 2*(h*k + k*l + h*l));
344  break;
345  }
346  case TETRAGONAL:
347  {
348  return 2*par[2]*l*l;
349  break;
350  }
351  case CUBIC:
352  {
353  throw 0;
354  break;
355  }
356  }
357  }
358  if(derivpar==(par+3))
359  {
360  switch(mlattice)
361  {
362  case TRICLINIC:
363  {
364  return l*l;
365  break;
366  }
367  case MONOCLINIC:
368  {
369  return 2*par[3]*l*l + 2*par[1]*par[4]*h*l;
370  break;
371  }
372  case ORTHORHOMBIC:
373  {
374  return 2*par[3]*l*l;
375  break;
376  }
377  case HEXAGONAL:
378  {
379  throw 0;
380  break;
381  }
382  case RHOMBOEDRAL:
383  {
384  throw 0;
385  break;
386  }
387  case TETRAGONAL:
388  {
389  throw 0;
390  break;
391  }
392  case CUBIC:
393  {
394  throw 0;
395  break;
396  }
397  }
398  }
399  if(derivpar==(par+4))
400  {
401  switch(mlattice)
402  {
403  case TRICLINIC:
404  {
405  return h*k;
406  break;
407  }
408  case MONOCLINIC:
409  {
410  return 2*par[1]*par[3]*h*l;
411  break;
412  }
413  default:
414  {
415  throw 0;
416  break;
417  }
418  }
419  }
420  if(derivpar==(par+5))
421  {
422  switch(mlattice)
423  {
424  case TRICLINIC:
425  {
426  return k*l;
427  break;
428  }
429  default:
430  {
431  throw 0;
432  break;
433  }
434  }
435  }
436  if(derivpar==(par+6))
437  {
438  switch(mlattice)
439  {
440  case TRICLINIC:
441  {
442  return h*l;
443  break;
444  }
445  default:
446  {
447  throw 0;
448  break;
449  }
450  }
451  }
452  throw 0;
453  return 0.0;
454 }
455 
456 void RecUnitCell::hkl2d_delta(const float h,const float k,const float l,
457  const RecUnitCell &delta, float & dmin, float &dmax) const
458 {
459  const float p0m=par[0]-delta.par[0] , p0p=par[0]+delta.par[0],
460  p1m=par[1]-delta.par[1] , p1p=par[1]+delta.par[1],
461  p2m=par[2]-delta.par[2] , p2p=par[2]+delta.par[2],
462  p3m=par[3]-delta.par[3] , p3p=par[3]+delta.par[3],
463  p4m=par[4]-delta.par[4] , p4p=par[4]+delta.par[4],
464  p5m=par[5]-delta.par[5] , p5p=par[5]+delta.par[5],
465  p6m=par[6]-delta.par[6] , p6p=par[6]+delta.par[6];
466  switch(mlattice)
467  {
468  case TRICLINIC:
469  {//par[0]+par[1]*h*h + par[2]*k*k + par[3]*l*l + par[4]*h*k + par[5]*k*l + par[6]*h*l;
470  float p4mm,p5mm,p6mm,p4pp,p5pp,p6pp;
471  if((h*k)>0){p4mm=p4m;p4pp=p4p;}else{p4mm=p4p;p4pp=p4m;}
472  if((k*l)>0){p5mm=p5m;p5pp=p5p;}else{p5mm=p5p;p5pp=p5m;}
473  if((h*l)>0){p6mm=p6m;p6pp=p6p;}else{p6mm=p6p;p6pp=p6m;}
474  dmin=p0m+p1m*h*h+p2m*k*k+p3m*l*l+p4mm*h*k+p5mm*k*l+p6mm*h*l;
475  dmax=p0p+p1p*h*h+p2p*k*k+p3p*l*l+p4pp*h*k+p5pp*k*l+p6pp*h*l;
476  /*
477  if(dmin<0)
478  {
479  cout<<"hkl2d_delta: dmin<0 ! "<<int(h)<<","<<int(k)<<","<<int(l)<<endl;
480  for(unsigned int i=0;i<7;++i) cout<<par[i]<<" +/-"<<delta.par[i]<<endl;
481  exit(0);
482  }*/
483  return;
484  }
485  case MONOCLINIC: //OK
486  {
487  if((h*l)>0)
488  {
489  dmin = p0m + p1m*p1m*h*h + p2m*p2m*k*k + p3m*p3m*l*l + 2*p1m*p3m*p4m*h*l;
490  dmax = p0p + p1p*p1p*h*h + p2p*p2p*k*k + p3p*p3p*l*l + 2*p1p*p3p*p4p*h*l;
491  return;
492  }
493  else
494  {
495  const bool b1=(h*(par[1]*h+par[3]*par[4]*l))>0;// d(d*^2)/dp1
496  const bool b3=(l*(par[3]*l+par[1]*par[4]*h))>0;// d(d*^2)/dp2
497  if(b1 && b3)
498  {
499  dmin = p0m + p1m*p1m*h*h + p2m*p2m*k*k + p3m*p3m*l*l + 2*p1m*p3m*p4p*h*l;
500  dmax = p0p + p1p*p1p*h*h + p2p*p2p*k*k + p3p*p3p*l*l + 2*p1p*p3p*p4m*h*l;
501  return;
502  }
503  else if(b1 && (!b3))
504  {
505  dmin = p0m + p1m*p1m*h*h + p2m*p2m*k*k + p3p*p3p*l*l + 2*p1m*p3p*p4p*h*l;
506  dmax = p0p + p1p*p1p*h*h + p2p*p2p*k*k + p3m*p3m*l*l + 2*p1p*p3m*p4m*h*l;
507  return;
508  }
509  else if((!b1) && b3)
510  {
511  dmin = p0m + p1p*p1p*h*h + p2m*p2m*k*k + p3m*p3m*l*l + 2*p1p*p3m*p4p*h*l;
512  dmax = p0p + p1m*p1m*h*h + p2p*p2p*k*k + p3p*p3p*l*l + 2*p1m*p3p*p4m*h*l;
513  return;
514  }
515  else
516  {
517  dmin = p0m + p1p*p1p*h*h + p2m*p2m*k*k + p3p*p3p*l*l + 2*p1p*p3p*p4p*h*l;
518  dmax = p0p + p1m*p1m*h*h + p2p*p2p*k*k + p3m*p3m*l*l + 2*p1m*p3m*p4m*h*l;
519  return;
520  }
521  }
522  }
523  case ORTHORHOMBIC: //OK
524  {
525  dmin= p0m + p1m*p1m*h*h + p2m*p2m*k*k + p3m*p3m*l*l;
526  dmax= p0p + p1p*p1p*h*h + p2p*p2p*k*k + p3p*p3p*l*l;
527  return;
528  }
529  case HEXAGONAL: //OK
530  {
531  dmin=p0m+p1m*p1m*(h*h + k*k + h*k)+ p2m*p2m*l*l ;
532  dmax=p0p+p1p*p1p*(h*h + k*k + h*k)+ p2p*p2p*l*l ;
533  return;
534  }
535  case RHOMBOEDRAL:
536  {
537  if((h*k + k*l + h*l)>=0)
538  {
539  dmin= p0m+p1m*p1m*(h*h + k*k + l*l + 2*p2m*(h*k + k*l + h*l));
540  dmax= p0p+p1p*p1p*(h*h + k*k + l*l + 2*p2p*(h*k + k*l + h*l));
541  }
542  else
543  {
544  dmin= p0m+p1m*p1m*(h*h + k*k + l*l + 2*p2p*(h*k + k*l + h*l));
545  dmax= p0p+p1p*p1p*(h*h + k*k + l*l + 2*p2m*(h*k + k*l + h*l));
546  }
547  return;
548  }
549  case TETRAGONAL: //OK
550  {
551  dmin= p0m + p1m*p1m*(h*h + k*k) + p2m*p2m*l*l;
552  dmax= p0p + p1p*p1p*(h*h + k*k) + p2p*p2p*l*l;
553  return;
554  }
555  case CUBIC: //OK
556  {
557  dmin=p0m + p1m*p1m*(h*h+k*k+l*l);
558  dmax=p0p + p1p*p1p*(h*h+k*k+l*l);
559  return;
560  }
561  }
562 }
563 
564 vector<float> RecUnitCell::DirectUnitCell(const bool equiv)const
565 {
566  // reciprocal unit cell parameters
567  float aa,bb,cc,calphaa,cbetaa,cgammaa,salphaa,sbetaa,sgammaa;
568  switch(mlattice)
569  {
570  case TRICLINIC:
571  {
572  aa=sqrt(par[1]);
573  bb=sqrt(par[2]);
574  cc=sqrt(par[3]);
575  calphaa=par[5]/(2*bb*cc);
576  cbetaa =par[6]/(2*aa*cc);
577  cgammaa=par[4]/(2*aa*bb);
578  salphaa=sqrt(abs(1-calphaa*calphaa));
579  sbetaa =sqrt(abs(1-cbetaa *cbetaa));
580  sgammaa=sqrt(abs(1-cgammaa*cgammaa));
581  break;
582  }
583  case MONOCLINIC:
584  {
585  aa=par[1];
586  bb=par[2];
587  cc=par[3];
588  calphaa=0;
589  cbetaa=par[4];
590  cgammaa=0;
591  salphaa=1;
592  sbetaa =sqrt(abs(1-cbetaa *cbetaa));
593  sgammaa=1;
594  break;
595  }
596  case ORTHORHOMBIC:
597  {
598  aa=par[1];
599  bb=par[2];
600  cc=par[3];
601  calphaa=0;
602  cbetaa =0;
603  cgammaa=0;
604  salphaa=1;
605  sbetaa =1;
606  sgammaa=1;
607  break;
608  }
609  case HEXAGONAL:
610  {
611  aa=par[1];
612  bb=par[1];
613  cc=par[2];
614  calphaa=0;
615  cbetaa =0;
616  cgammaa=0.5;
617  salphaa=1;
618  sbetaa =1;
619  sgammaa=0.8660254037844386;
620  break;
621  }
622  case RHOMBOEDRAL:
623  {
624  aa=par[1];
625  bb=par[1];
626  cc=par[1];
627  calphaa=par[4];
628  cbetaa =par[4];
629  cgammaa=par[4];
630  salphaa=sqrt(abs(1-calphaa *calphaa));
631  sbetaa =salphaa;
632  sgammaa=salphaa;
633  break;
634  }
635  case TETRAGONAL:
636  {
637  aa=par[1];
638  bb=par[1];
639  cc=par[2];
640  calphaa=0;
641  cbetaa =0;
642  cgammaa=0;
643  salphaa=1;
644  sbetaa =1;
645  sgammaa=1;
646  break;
647  }
648  case CUBIC:
649  {
650  aa=par[1];
651  bb=par[1];
652  cc=par[1];
653  calphaa=0;
654  cbetaa =0;
655  cgammaa=0;
656  salphaa=1;
657  sbetaa =1;
658  sgammaa=1;
659  break;
660  }
661  // This should never happen. Avoid using unitialized cell parameters.
662  default:
663  throw 0;
664  }
665  // Volume of reciprocal unit cell
666  const float vv=sqrt(abs(1-calphaa*calphaa-cbetaa*cbetaa-cgammaa*cgammaa+2*calphaa*cbetaa*cgammaa));
667 
668  const float a=salphaa/(aa*vv);
669  const float b=sbetaa /(bb*vv);
670  const float c=sgammaa/(cc*vv);
671 
672  const float calpha=(cbetaa *cgammaa-calphaa)/(sbetaa *sgammaa);
673  const float cbeta =(calphaa*cgammaa-cbetaa )/(salphaa*sgammaa);
674  const float cgamma=(calphaa*cbetaa -cgammaa)/(salphaa*sbetaa );
675 
676  const float alpha=acos( calpha );
677  const float beta =acos( cbeta );
678  const float gamma=acos( cgamma );
679 
680  const float v=a*b*c*sqrt(1-calpha*calpha-cbeta*cbeta-cgamma*cgamma+2*calpha*cbeta*cgamma);
681 
682  vector<float> par(7);
683  par[0]=a;
684  par[1]=b;
685  par[2]=c;
686  par[3]=alpha;
687  par[4]=beta;
688  par[5]=gamma;
689  par[6]=v;
690  return par;
691 }
693 PeakList::hkl0::hkl0(const int h0,const int k0, const int l0):
694 h(h0),k(k0),l(l0)
695 {}
696 
698 PeakList::hkl::hkl(const float d,const float i,const float ds,const float is,
699  const int h0,const int k0, const int l0,const float dc0):
700 dobs(d),dobssigma(ds),d2obs(d*d),d2obsmin((d-ds/2)*(d-ds/2)),d2obsmax((d+ds/2)*(d+ds/2)),iobs(i),iobssigma(is),
701 h(h0),k(k0),l(l0),isIndexed(false),isSpurious(false),stats(0),
702 d2calc(dc0),d2diff(0)
703 {}
704 
705 bool compareHKL_d(const PeakList::hkl &d1, const PeakList::hkl &d2)
706 {
707  return d1.dobs < d2.dobs;
708 }
709 
710 
712 
713 PeakList::PeakList()
714 {}
715 
716 PeakList::PeakList(const PeakList &old)
717 {
718  *this=old;
719 }
720 
721 void PeakList::operator=(const PeakList &rhs)
722 {
723  VFN_DEBUG_ENTRY("PeakList::operator=(PeakList &old)",10);
724  mvHKL=rhs.mvHKL;
725  VFN_DEBUG_EXIT("PeakList::operator=(PeakList &old)",10);
726 }
727 
728 PeakList::~PeakList()
729 {}
730 
731 void PeakList::ImportDhklDSigmaIntensity(istream &is,float defaultsigma)
732 {
733  float d,sigma,iobs;
734  while(true)
735  {// :TODO: use readline to make sure when the end is reached
736  is >>d;
737  cout<<__FILE__<<":"<<__LINE__<<" "<<mvHKL.size()<<":d="<<d;
738  if(is.good()==false) break;
739  is>>sigma;
740  if(is.good()==false) break;
741  is>>iobs;
742  if(sigma<=0) sigma=d*defaultsigma;
743  if(iobs<=0) iobs=1.0;
744  mvHKL.push_back(hkl(1/d,iobs,1/(d-sigma/2)-1/(d+sigma/2)));
745  cout<<"+/-"<<sigma<<", I="<<iobs<<" 1/d="<<1/d<<endl;
746  if(is.good()==false) break;
747  }
748  sort(mvHKL.begin(),mvHKL.end(),compareHKL_d);
749  cout<<"Imported "<<mvHKL.size()<<" observed reflection positions."<<endl;
750 }
751 
752 void PeakList::ImportDhklIntensity(istream &is)
753 {
754  float d,iobs;
755  while(true)
756  {// :TODO: use readline to make sure when the end is reached
757  is >>d;
758  if(is.eof()) break;
759  is>>iobs;
760  mvHKL.push_back(hkl(1/d,iobs));
761  cout<<__FILE__<<":"<<__LINE__<<" "<<mvHKL.size()<<":d="<<d<<", I="<<iobs<<" 1/d="<<1/d<<endl;
762  if(is.eof()) break;
763  }
764  sort(mvHKL.begin(),mvHKL.end(),compareHKL_d);
765  cout<<"Imported "<<mvHKL.size()<<" observed reflection positions."<<endl;
766 }
767 
768 void PeakList::ImportDhkl(istream &is)
769 {
770  std::vector<std::pair<float,float> > v;
771  float d;
772  while(true)
773  {// :TODO: use readline to make sure when the end is reached
774  is >>d;
775  if(is.eof()) break;
776  mvHKL.push_back(hkl(1/d));
777  cout<<__FILE__<<":"<<__LINE__<<" "<<mvHKL.size()<<":d="<<d<<" 1/d="<<1/d<<endl;
778  if(is.eof()) break;
779  }
780  sort(mvHKL.begin(),mvHKL.end(),compareHKL_d);
781  cout<<"Imported "<<v.size()<<" observed reflection positions."<<endl;
782 }
783 
784 
785 template<class T,class U> bool comparePairFirst(std::pair<T,U> &p1, std::pair<T,U> &p2)
786 {
787  return p1.first < p2.first;
788 }
789 
790 void PeakList::Import2ThetaIntensity(istream &is, const float wavelength)
791 {
792  std::list<std::pair<float,float> > v;
793  float d,iobs;
794  while(true)
795  {// :TODO: use readline to make sure when the end is reached
796  is >>d;
797  if(is.eof()) break;
798  is>>iobs;
799  d=2*sin(d/2*DEG2RAD)/wavelength;
800  mvHKL.push_back(hkl(1/d,iobs));
801  cout<<__FILE__<<":"<<__LINE__<<" "<<mvHKL.size()<<":d="<<1/d<<", I="<<iobs<<" 1/d="<<d<<endl;
802  if((is.eof())||v.size()>=20) break;
803  }
804  sort(mvHKL.begin(),mvHKL.end(),compareHKL_d);
805  cout<<"Imported "<<v.size()<<" observed reflection positions."<<endl;
806 }
807 float PeakList::Simulate(float zero, float a, float b, float c,
808  float alpha, float beta, float gamma,
809  bool deg, unsigned int nb, unsigned int nbspurious,
810  float sigma,float percentMissing, const bool verbose,
811  const bool merge)
812 {
813  if(deg){alpha*=DEG2RAD;beta*=DEG2RAD;gamma*=DEG2RAD;}
814  const float v=sqrt(1-cos(alpha)*cos(alpha)-cos(beta)*cos(beta)-cos(gamma)*cos(gamma)
815  +2*cos(alpha)*cos(beta)*cos(gamma));
816 
817  const float aa=sin(alpha)/a/v;
818  const float bb=sin(beta )/b/v;
819  const float cc=sin(gamma)/c/v;
820 
821  const float alphaa=acos( (cos(beta )*cos(gamma)-cos(alpha))/sin(beta )/sin(gamma) );
822  const float betaa =acos( (cos(alpha)*cos(gamma)-cos(beta ))/sin(alpha)/sin(gamma) );
823  const float gammaa=acos( (cos(alpha)*cos(beta )-cos(gamma))/sin(alpha)/sin(beta ) );
824 
825  RecUnitCell ruc(zero,aa*aa,bb*bb,cc*cc,2*aa*bb*cos(gammaa),2*bb*cc*cos(alphaa),2*aa*cc*cos(betaa),TRICLINIC,LATTICE_P);
826  std::list<float> vd2;
827 
828  for(int h=0;h<=20;h++)
829  for(int k=-20;k<=20;k++)
830  {
831  if((h==0) && (k<0)) k=0;
832  for(int l=-20;l<=20;l++)
833  {
834  if((h==0) && (k==0) && (l<=0)) l=1;
835  vd2.push_back(sqrt(ruc.hkl2d(h,k,l)));
836  }
837  }
838  //
839  std::list<float>::iterator pos=vd2.begin();
840  if(percentMissing>0.90) percentMissing=0.90;
841  for(;pos!=vd2.end();++pos)
842  {
843  if((rand()/float(RAND_MAX))<percentMissing) *pos=1e10;
844  }
845  vd2.sort();
846  if(merge)
847  {// Remove duplicate reflections (e.g. symmetry-equivalent) with (near-)identical d-spacing,
848  // so that they are counted only once amongst the 'nb' observed peaks.
849  const float atol=1e-5;
850  for(std::list<float>::iterator p1=vd2.begin(), p2=++vd2.begin(); p2!=vd2.end();)
851  {
852  if(fabs(*p2-*p1)<atol) p2=vd2.erase(p2);
853  else {++p1;++p2;}
854  }
855  }
856  pos=vd2.begin();
857  const float dmin=*pos/2;
858  for(unsigned int i=0;i<nb;++i)pos++;
859  const float dmax=*pos;
860 
861  for(unsigned int i=0;i<nbspurious;++i)
862  {
863  const unsigned int idx=1+i*nb/nbspurious+(rand()%nbspurious);
864  pos=vd2.begin();
865  for(unsigned int j=0;j<idx;++j) pos++;
866  *pos=dmin+rand()/float(RAND_MAX)*(dmax-dmin);
867  }
868 
869  pos=vd2.begin();
870  for(unsigned int i=0;i<nb;++i)
871  {
872  float d=*pos++;
873  const float ds=d*sigma;
874  float d1=d+ds*(rand()/float(RAND_MAX)*2-1);
875  //cout<<d<<" "<<ds<<" "<<d1<<" "<<sigma<<endl;
876  mvHKL.push_back(hkl(d1,1.0,ds));
877  }
878 
879  if(verbose)
880  {
881  char buf[200];
882  sprintf(buf,"a=%6.3f b=%6.3f c=%6.3f alpha=%6.2f beta=%6.2f gamma=%6.2f V=%8.2f",a,b,c,alpha*RAD2DEG,beta*RAD2DEG,gamma*RAD2DEG,v*a*b*c);
883  cout<<"PeakList::Simulate():"<<buf<<endl;
884  }
885  sort(mvHKL.begin(),mvHKL.end(),compareHKL_d);
886  return v*a*b*c;
887 }
888 
889 void PeakList::ExportDhklDSigmaIntensity(std::ostream &os)const
890 {
891  for(vector<PeakList::hkl>::const_iterator pos=mvHKL.begin();pos!=mvHKL.end();++pos)
892  {
893  const float sigma=1/(pos->dobs-pos->dobssigma/2)-1/(pos->dobs+pos->dobssigma/2);
894  os<< std::fixed << setw(6) << setprecision(3) << 1/pos->dobs <<" "<< sigma <<" "<< std::scientific << pos->iobs <<endl;
895  }
896 }
897 
898 void PeakList::AddPeak(const float d, const float iobs,const float dobssigma,const float iobssigma,
899  const int h,const int k, const int l,const float d2calc)
900 {
901  if(dobssigma<=0)
902  {// Manually added peak ? Use other reflection's sigmas to evaluate sigma for this reflection
903  float s=0;
904  for(vector<hkl>::const_iterator pos=mvHKL.begin();pos!=mvHKL.end();++pos)
905  s+= pos->dobssigma;
906  s/=mvHKL.size();
907  if(s>0) mvHKL.push_back(hkl(d,iobs,s,iobssigma,h,k,l,d2calc));
908  else mvHKL.push_back(hkl(d,iobs,1e-4,iobssigma,h,k,l,d2calc));
909  }
910  else mvHKL.push_back(hkl(d,iobs,dobssigma,iobssigma,h,k,l,d2calc));
911  sort(mvHKL.begin(),mvHKL.end(),compareHKL_d);
912  //this->Print(cout);
913 }
914 
915 void PeakList::RemovePeak(unsigned int idx)
916 {
917  for(unsigned int i=idx;i<(mvHKL.size()-1);++i) mvHKL[i]=mvHKL[i+1];
918  mvHKL.pop_back();
919 }
920 
921 void PeakList::Print(std::ostream &os) const
922 {
923  unsigned int i=0;
924  char buf[200];
925  os<<"PeakList, with "<<mvHKL.size()<<" peaks"<<endl;
926  for(vector<PeakList::hkl>::const_iterator pos=mvHKL.begin();pos!=mvHKL.end();++pos)
927  {
928  const float sigma=1/(pos->dobs-pos->dobssigma/2)-1/(pos->dobs+pos->dobssigma/2);
929  if(pos->isIndexed)
930  sprintf(buf,"#%3d d=%6.3f+/-%7.4f dcalc=%6.3f, diff=%7.4f, iobs=%6.3f HKL=%2d %2d %2d Spurious=%1d stats=%6lu",
931  i++,1/pos->dobs,sigma,
932  1/sqrt(abs(pos->d2calc)),1/sqrt(abs(pos->d2calc))-1/pos->dobs,
933  pos->iobs,pos->h,pos->k,pos->l,pos->isSpurious,pos->stats);
934  else
935  sprintf(buf,"#%3d d=%6.3f+/-%6.3f iobs=%6.3f UNINDEXED Spurious=%1d stats=%6lu",
936  i++,1/pos->dobs,1/(pos->dobs-pos->dobssigma/2)-1/(pos->dobs+pos->dobssigma/2),
937  pos->iobs,pos->isSpurious,pos->stats);
938  os<<buf<<endl;
939  }
940 }
941 
942 vector<PeakList::hkl> & PeakList::GetPeakList(){return mvHKL;}
943 
944 const vector<PeakList::hkl> & PeakList::GetPeakList()const {return mvHKL;}
945 
947 
948 float Score(const PeakList &dhkl, const RecUnitCell &rpar, const unsigned int nbSpurious,
949  const bool verbose,const bool storehkl,const bool storePredictedHKL)
950 {
951  const bool autozero=false;
952  vector<PeakList::hkl>::const_iterator pos,first,last;
953  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
954  {
955  if(storehkl) pos->isIndexed=false;
956  pos->d2calc=0;
957  pos->d2diff=1000;
958  }
959  const unsigned long nb=dhkl.GetPeakList().size();
960  if(storePredictedHKL) dhkl.mvPredictedHKL.clear();
961 
962  unsigned long nbCalc=0;
963  int h,k,l;
964  float predict_coeff=1;
965  if(storePredictedHKL)predict_coeff=2;
966  const float dmax=dhkl.mvHKL[nb-1].d2obs*predict_coeff*1.05;
967  int sk0,sl0;// do we need >0 *and* <0 indices for k,l ?
968  switch(rpar.mlattice)
969  {
970  case TRICLINIC:
971  sk0=-1;sl0=-1;
972  break;
973  case MONOCLINIC:
974  sk0=1;sl0=-1;
975  break;
976  case ORTHORHOMBIC:
977  sk0=1;sl0=1;
978  break;
979  case HEXAGONAL:
980  sk0=-1;sl0=1;
981  break;
982  case RHOMBOEDRAL:
983  sk0=-1;sl0=-1;
984  break;
985  case TETRAGONAL:
986  sk0=1;sl0=1;
987  break;
988  case CUBIC:
989  sk0=1;sl0=1;
990  break;
991  // This should never happen. Avoid using unitialized values.
992  default:
993  throw 0;
994  }
995  int stepk,stepl;// steps in k,l to use for centered lattices
996  switch(rpar.mCentering)
997  {
998  case LATTICE_P:stepk=1;stepl=1;break;
999  case LATTICE_I:stepk=1;stepl=2;break;
1000  case LATTICE_A:stepk=1;stepl=2;break;
1001  case LATTICE_B:stepk=1;stepl=2;break;
1002  case LATTICE_C:stepk=2;stepl=1;break;
1003  case LATTICE_F:stepk=2;stepl=2;break;
1004  // This should never happen. Avoid using unitialized values.
1005  default: throw 0;
1006  }
1007  first=dhkl.GetPeakList().begin();last=dhkl.GetPeakList().end();
1008  unsigned long nbCalcH,nbCalcK;// Number of calculated lines below dmax for each h,k
1009  for(h=0;;++h)
1010  {
1011  nbCalcH=0;
1012  for(int sk=sk0;sk<=1;sk+=2)
1013  {
1014  if(h==0) sk=1;// no need to explore 0kl with both sk -1 and 1
1015  if(stepk==2) k=(h%2);// For LATTICE_C,LATTICE_F: h odd => k odd
1016  else k=0;
1017  for(;;k+=stepk)
1018  {
1019  nbCalcK=0;
1020  for(int sl=sl0;sl<=1;sl+=2)
1021  {
1022  if((h+k)==0)
1023  {
1024  sl=1;// No need to list 0 0 l with l<0
1025  l=1;
1026  }
1027  else
1028  {
1029  if(h==0)
1030  {
1031  if(rpar.mlattice==MONOCLINIC) sl=1;// 0 k l and 0 k -l are equivalent
1032  if((sk<0)||(sl<0)) l=1;// Do not list 0 k 0 with k<0
1033  else l=0;// h==k==0 already covered
1034  }
1035  else
1036  {
1037  if(sl<0) l=1;// Do not list h k 0 twice
1038  else l=0;
1039  }
1040  }
1041  if(stepl==2)
1042  {
1043  if(rpar.mCentering==LATTICE_I) l+=(h+k+l)%2;
1044  if(rpar.mCentering==LATTICE_A) l+=(k+l)%2;// Start at hk1 if k odd
1045  if( (rpar.mCentering==LATTICE_B)
1046  ||(rpar.mCentering==LATTICE_F)) l+=(h+l)%2;// Start at hk1 if h odd
1047  }
1048  for(;;l+=stepl)
1049  {
1050  const float d2=rpar.hkl2d(h,sk*k,sl*l);
1051  if(d2>dmax)
1052  {
1053  //cout<<__FILE__<<":"<<__LINE__<<" hkl: "<<h<<" "<<sk*k<<" "<<sl*l<<":"<<sqrt(d2)<<" deriv="<<sl*rpar.hkl2d(h,sk*k,sl*l,NULL,3)<<"/"<<sqrt(dmax)<<endl;
1054  // Only break if d is increasing with l
1055  if((sl*rpar.hkl2d(h,sk*k,sl*l,NULL,3))>=0) break;
1056  else continue;
1057  }
1058  nbCalc++;nbCalcK++;nbCalcH++;
1059  if(storePredictedHKL)
1060  {
1061  dhkl.mvPredictedHKL.push_back(PeakList::hkl(0,0,0,0,h,sk*k,sl*l,d2));
1062  //continue;
1063  }
1064  for(pos=first;pos!=last;++pos)
1065  {
1066  const float tmp=d2-pos->d2obs;
1067  if(tmp<.1)
1068  {
1069  if(tmp<-.1) break;
1070  if(fabs(tmp)<fabs(pos->d2diff))
1071  {
1072  pos->d2diff=tmp;
1073  if(storehkl)
1074  {
1075  pos->h=h;
1076  pos->k=sk*k;
1077  pos->l=sl*l;
1078  pos->isIndexed=true;
1079  pos->d2calc=d2;
1080  }
1081  }
1082  /*
1083  if((verbose)&&(fabs(tmp)<.01))
1084  cout<<__FILE__<<":"<<__LINE__<<" hkl: "<<h<<" "<<k<<" "<<l
1085  <<"#"<<i<<": calc="<<sqrt(d2)<<", obs="<<sqrt(*pd2)<<", min_epsilon="<<*pdq2<<", tmp="<<tmp<<endl;
1086  */
1087  }
1088  }
1089  }
1090  }
1091  if(nbCalcK==0) //d(hk0)>dmax
1092  {
1093  //cout<<__FILE__<<":"<<__LINE__<<" hkl: "<<h<<" "<<sk*k<<" "<<0<<" deriv="<<sk*rpar.hkl2d(h,sk*k,0,NULL,2)<<endl;
1094  if((sk*rpar.hkl2d(h,sk*k,0,NULL,2))>=0) break;
1095  }
1096  }
1097  }
1098  if(nbCalcH==0) break;//h00 beyond limit
1099  }
1100  float epsilon=0.0,zero=0.0;
1101  if(autozero)
1102  {
1103  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos) zero+=pos->d2diff;
1104  zero/=nb;
1105  }
1106  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
1107  {
1108  epsilon +=fabs(pos->d2diff-zero);
1109  }
1110  if(nbSpurious>0)
1111  {// find worst fitting lines and remove them from epsilon calculation
1112  list<pair<float,unsigned int> > vdiff_idx;
1113  unsigned int i=0;
1114  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
1115  vdiff_idx.push_back(make_pair(fabs(pos->d2diff),i++));
1116  vdiff_idx.sort(comparePairFirst<float,unsigned int>);
1117  i=0;
1118  for(list<pair<float,unsigned int> >::reverse_iterator rpos=vdiff_idx.rbegin();rpos!=vdiff_idx.rend();++rpos)
1119  {// :TODO: correct zero after removing spurious lines
1120  epsilon -= fabs(rpos->first-zero);
1121  if(storehkl) dhkl.GetPeakList()[rpos->second].isIndexed=false;
1122  dhkl.GetPeakList()[rpos->second].stats++;
1123  if(++i==nbSpurious) break;
1124  }
1125  }
1126  if(verbose)
1127  {
1128  float epstmp=0;
1129  //unsigned long i=0;
1130  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
1131  {
1132  epstmp+=fabs(pos->d2diff-zero);
1133  //cout<<"Line #"<<i++<<": obs="<<pos->d2obs<<", diff="<<pos->d2diff<<" -> epsilon="<<epstmp<<endl;
1134  }
1135  cout<<"epsilon="<<epstmp<<", dmax="<<dmax<<" ,nb="<<nb<<" ,nbcalc="<<nbCalc<<endl;
1136  }
1137  /*
1138  else
1139  {//Only stat+ the worst
1140  float max=-1;
1141  unsigned int worst=0;
1142  unsigned int i=0;
1143  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
1144  if(abs(pos->d2diff)>max) {worst=i++;max=abs(pos->d2diff);}
1145  else i++;
1146  dhkl.GetPeakList()[worst].stats++;
1147  }
1148  */
1149  if(nbCalc==0) return 0;
1150  const float score=(dmax/predict_coeff)*nb/(2*epsilon*nbCalc);
1151  if(verbose)
1152  {
1153  dhkl.Print(cout);
1154  cout<<"Final score:"<<score<<", nbCalc="<<nbCalc<<" ,<epsilon>="<<epsilon<<" nb="<<nb<<" Qn="<<sqrt(dmax)<<endl;
1155  }
1156  return score;
1157 }
1158 
1160 
1161 CellExplorer::CellExplorer(const PeakList &dhkl, const CrystalSystem lattice, const unsigned int nbSpurious):
1162 mnpar(3),mpPeakList(&dhkl),
1163 mLengthMin(4),mLengthMax(25),
1164 mAngleMin(M_PI),mAngleMax(2*M_PI/3),
1165 mVolumeMin(0),mVolumeMax(1600),
1166 mZeroShiftMin(0),mZeroShiftMax(0),
1167 mlattice(lattice),mCentering(LATTICE_P),mNbSpurious(nbSpurious),
1168 mObs(0),mCalc(0),mWeight(0),mDeriv(0),mBestScore(0.0),
1169 mMinScoreReport(10),mMaxDicVolDepth(6),mDicVolDepthReport(6),
1170 mNbLSQExcept(0)
1171 {
1172  this->Init();
1173 }
1174 
1175 void CellExplorer::Evolution(unsigned int ng,const bool randomize,const float f,const float cr,unsigned int np)
1176 {
1177  this->Init();
1178  const bool autozero=true;
1179  //cout<<__FILE__<<":"<<__LINE__<<"<CellExplorer::Evolution(...): randomizing,ng="<<ng
1180  // <<"random="<<randomize<<"f="<<f<<"cr="<<cr<<"np="<<np<<endl;
1181  vector<pair<RecUnitCell,float> > vRUC(np);
1182  vector<pair<RecUnitCell,float> > vTrial(np);
1183  float bestScore=-1e20;
1184  vector<pair<RecUnitCell,float> >::iterator bestpos=vRUC.begin();
1185 
1186  const clock_t mTime0=clock();
1187 
1188  if(randomize)
1189  {
1190  for(unsigned int i=0;i<vRUC.size();++i)
1191  {
1192  vRUC[i].first.mlattice=mlattice;
1193  vTrial[i].first.mlattice=mlattice;
1194  for(unsigned int k=0;k<mnpar;++k) vRUC[i].first.par[k]=mMin[k]+mAmp[k]*rand()/(float)RAND_MAX;
1195  vRUC[i].second=Score(*mpPeakList,vRUC[i].first,mNbSpurious);
1196  }
1197  }
1198 
1199  for(unsigned long i=ng;i>0;--i)
1200  {
1201  for(unsigned j=0;j<np;j++)
1202  {
1203  if(true)
1204  {// DE/rand/1/exp
1205  unsigned int r1=j,r2=j,r3=j;
1206  while(r1==j)r1=rand()%np;
1207  while((r2==j)||(r1==r2))r2=rand()%np;
1208  while((r3==j)||(r3==r1)||(r3==r2))r3=rand()%np;
1209  unsigned int ncr=1+(int)(cr*mnpar*rand()/(float)RAND_MAX);
1210  unsigned int ncr0=rand()%mnpar;
1211  RecUnitCell *t0=&(vTrial[j].first);
1212  RecUnitCell *c0=&(vRUC[j].first);
1213  RecUnitCell *c1=&(vRUC[r1].first);
1214  RecUnitCell *c2=&(vRUC[r2].first);
1215  RecUnitCell *c3=&(vRUC[r3].first);
1216  for(unsigned int k=0;k<mnpar;++k)t0->par[k] = c0->par[k];
1217  for(unsigned int k=0;k<ncr;++k)
1218  {
1219  const unsigned l=(ncr0+k)%mnpar;
1220  const float v1=c1->par[l]-mMin[l];
1221  const float v2=c2->par[l]-mMin[l];
1222  const float v3=c3->par[l]-mMin[l];
1223  t0->par[l]=mMin[l]+fmod(v1+f*(v2-v3)+3*mAmp[l],mAmp[l]);
1224  }
1225  }
1226  if(false)
1227  {// DE/rand-to-best/1/exp
1228  unsigned int r1=j,r2=j,r3=j;
1229  while(r1==j)r1=rand()%np;
1230  while((r2==j)||(r1==r2))r2=rand()%np;
1231  while((r3==j)||(r3==r1)||(r3==r2))r3=rand()%np;
1232  unsigned int ncr=1+(int)(cr*(mnpar-1)*rand()/(float)RAND_MAX);
1233  unsigned int ncr0=rand()%mnpar;
1234  RecUnitCell *t0=&(vTrial[j].first);
1235  RecUnitCell *c0=&(vRUC[j].first);
1236  //RecUnitCell *c1=&(vRUC[r1].first);
1237  RecUnitCell *c2=&(vRUC[r2].first);
1238  RecUnitCell *c3=&(vRUC[r3].first);
1239  RecUnitCell *best=&(bestpos->first);
1240  for(unsigned int k=0;k<6;++k)t0->par[k] = c0->par[k];//mMin[k]+mAmp[k]*rand()/(float)RAND_MAX;
1241  for(unsigned int k=0;k<ncr;++k)
1242  {
1243  const unsigned l=(ncr0+k)%mnpar;
1244  const float v0=c0->par[l]-mMin[l];
1245  //const float v1=c1->par[l]-mMin[l];
1246  const float v2=c2->par[l]-mMin[l];
1247  const float v3=c3->par[l]-mMin[l];
1248  const float vb=best->par[l]-mMin[l];
1249  t0->par[l]=mMin[l]+fmod(vb+f*(vb-v0)+f*(v2-v3)+5*mAmp[l],mAmp[l]);
1250  }
1251  }
1252  if(false)
1253  {// MC
1254  const float amp=.05/(1+i*.01);
1255  RecUnitCell *t0=&(vTrial[j].first);
1256  for(unsigned int k=0;k<6;++k)
1257  {
1258 
1259  t0->par[k] = mMin[k]+ fmod((float)(amp*mAmp[k]*(rand()/(float)RAND_MAX-0.5)+5*mAmp[k]),(float)mAmp[k]);
1260  }
1261  }
1262  }
1263  // Compute cost for all trials and select best
1264  vector<pair<RecUnitCell,float> >::iterator posTrial,pos;
1265  posTrial=vTrial.begin();
1266  pos=vRUC.begin();
1267  for(;posTrial!=vTrial.end();)
1268  {
1269  // If using auto-zero, fix zero parameter
1270  if(autozero) posTrial->first.par[0]=0;
1271  // Did we go beyond allowed volume ?
1272  switch(mlattice)
1273  {
1274  case TRICLINIC:
1275  break;
1276  case MONOCLINIC:
1277  {
1278  float v0=posTrial->first.par[1]*posTrial->first.par[2]*posTrial->first.par[3];
1279  while(v0<1/mVolumeMax)
1280  {
1281  const unsigned int i=rand()%3+1;
1282  posTrial->first.par[i]*=1/(mVolumeMax*v0)+1e-4;
1283  if(posTrial->first.par[i]>(mMin[i]+mAmp[i])) posTrial->first.par[i]=mMin[i]+mAmp[i];
1284  v0=posTrial->first.par[1]*posTrial->first.par[2]*posTrial->first.par[3];
1285  }
1286  break;
1287  }
1288  case ORTHORHOMBIC:
1289  break;
1290  case HEXAGONAL:
1291  break;
1292  case RHOMBOEDRAL:
1293  break;
1294  case TETRAGONAL:
1295  break;
1296  case CUBIC:
1297  break;
1298  }
1299 
1300  const float score=Score(*mpPeakList,posTrial->first,mNbSpurious);
1301  if(score > pos->second)
1302  {
1303  pos->second=score;
1304  const REAL *p0=posTrial->first.par;
1305  REAL *p1=pos->first.par;
1306  for(unsigned int k=0;k<mnpar;++k) *p1++ = *p0++;
1307  if(score>bestScore)
1308  {
1309  bestScore=score;
1310  bestpos=pos;
1311  }
1312  }
1313  /*
1314  else
1315  {
1316  if(log(rand()/(float)RAND_MAX)>(-(score-pos->second)))
1317  {
1318  pos->second=score;
1319  const float *p0=posTrial->first.par;
1320  float *p1=pos->first.par;
1321  for(unsigned int k=0;k<mnpar;++k) *p1++ = *p0++;
1322  }
1323  }
1324  */
1325  ++pos;++posTrial;
1326  }
1327  if((i%100000)==0)
1328  {
1329  vector<float> par=bestpos->first.DirectUnitCell();
1330  cout<<"Generation #"<<ng-i<<", Best score="<<bestScore
1331  <<" Trial: a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]<<", alpha="
1332  <<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG<<", V="<<par[6]
1333  <<" "<<(ng-i)*np/((clock()-mTime0)/(float)CLOCKS_PER_SEC)<<" trials/s"
1334  <<endl;
1335  }
1336  if(false)//((i%10000)==0)
1337  {// Randomize periodically
1338  for(vector<pair<RecUnitCell,float> >::iterator pos=vRUC.begin();pos!=vRUC.end();++pos)
1339  {
1340  if(pos==bestpos) continue;
1341  for(unsigned int k=0;k<mnpar;++k) pos->first.par[k]=mMin[k]+mAmp[k]*rand()/(float)RAND_MAX;
1342  }
1343  }
1344  }
1345  /*
1346  for(vector<pair<RecUnitCell,float> >::iterator pos=vRUC.begin();pos!=vRUC.end();++pos)
1347  {
1348  // Final cost
1349  vector<float> par=pos->first.DirectUnitCell();
1350  cout<<__FILE__<<":"<<__LINE__<<" Trial: a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]<<", alpha="
1351  <<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG<<", V="<<par[6]
1352  <<", score="<<pos->second<<endl;
1353  }
1354  Score(*mpPeakList,bestpos->first,mNbSpurious,true);
1355  */
1356 
1357  //this->ReduceSolutions(true);
1358 
1359  mRecUnitCell=bestpos->first;
1360  float score=Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true);
1361  vector<float> par=mRecUnitCell.DirectUnitCell();
1362  cout<<__FILE__<<":"<<__LINE__<<" Best-DE : a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]<<", alpha="
1363  <<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG<<", V="<<par[6]
1364  <<", score="<<bestpos->second
1365  <<" ("<<ng*np/((clock()-mTime0)/(float)CLOCKS_PER_SEC)<<" trials/s)"<<endl;
1366  if(score>mMinScoreReport*.5)
1367  {
1368  // Now, do a least-squares refinement on best
1369  mRecUnitCell=bestpos->first;
1370  this->LSQRefine(10,true,true);
1371  par=mRecUnitCell.DirectUnitCell();
1372  score=Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true);
1373  cout<<__FILE__<<":"<<__LINE__<<" Best-LSQ: a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]<<", alpha="
1374  <<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG<<", V="<<par[6]
1375  <<", score="<<score<<endl;
1376  if((score>mMinScoreReport)&&(score>(mBestScore/3)))
1377  {
1378  if(score>mBestScore) mBestScore=score;
1379  mvSolution.push_back(make_pair(mRecUnitCell,score));
1380  mvSolution.back().first.mNbSpurious = mNbSpurious;
1381  this->ReduceSolutions(true);// We may have solutions from previous runs
1382  }
1383  }
1384 }
1385 
1386 void CellExplorer::SetLengthMinMax(const float min,const float max)
1387 {
1388  mLengthMin=min;
1389  mLengthMax=max;
1390 }
1391 void CellExplorer::SetAngleMinMax(const float min,const float max)
1392 {
1393  mAngleMin=min;
1394  mAngleMax=max;
1395 }
1396 void CellExplorer::SetVolumeMinMax(const float min,const float max)
1397 {
1398  mVolumeMin=min;
1399  mVolumeMax=max;
1400 }
1401 void CellExplorer::SetNbSpurious(const unsigned int nb)
1402 {
1403  mNbSpurious=nb;
1404 }
1405 void CellExplorer::SetMinMaxZeroShift(const float min,const float max)
1406 {
1407  mZeroShiftMin=min;
1408  mZeroShiftMax=max;
1409 }
1410 
1411 void CellExplorer::SetCrystalSystem(const CrystalSystem system)
1412 {
1413  mlattice=system;
1414 }
1415 
1416 void CellExplorer::SetCrystalCentering(const CrystalCentering cent)
1417 {
1418  mCentering=cent;
1419 }
1420 
1421 void CellExplorer::SetD2Error(const float err){mD2Error=err;}
1422 
1423 const string& CellExplorer::GetClassName() const
1424 {
1425  const static string className="CellExplorer";
1426  return className;
1427 }
1428 const string& CellExplorer::GetName() const
1429 {
1430  const static string name="Some CellExplorer Object";
1431  return name;
1432 }
1433 void CellExplorer::Print() const
1434 {
1435  this->RefinableObj::Print();
1436 }
1437 unsigned int CellExplorer::GetNbLSQFunction() const
1438 {return 1;}
1439 
1440 const CrystVector_REAL& CellExplorer::GetLSQCalc(const unsigned int) const
1441 {
1442  VFN_DEBUG_ENTRY("CellExplorer::GetLSQCalc()",2)
1443  unsigned int j=0;
1444  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();++pos)
1445  {
1446  if(pos->isIndexed)
1447  mCalc(j++)=mRecUnitCell.hkl2d(pos->h,pos->k,pos->l);
1448  }
1449  //cout<<__FILE__<<":"<<__LINE__<<"LSQCalc : Score:"<<Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true,false)<<endl;
1450  VFN_DEBUG_EXIT("CellExplorer::GetLSQCalc()",2)
1451  return mCalc;
1452 }
1453 const CrystVector_REAL& CellExplorer::GetLSQObs(const unsigned int) const
1454 {
1455  VFN_DEBUG_MESSAGE("CellExplorer::GetLSQObs()",2)
1456  return mObs;
1457 }
1458 const CrystVector_REAL& CellExplorer::GetLSQWeight(const unsigned int) const
1459 {
1460  VFN_DEBUG_MESSAGE("CellExplorer::GetLSQWeight()",2)
1461  //:TODO: exclude the worst points (user-chosen number)
1462  return mWeight;
1463 }
1464 const CrystVector_REAL& CellExplorer::GetLSQDeriv(const unsigned int, RefinablePar &refpar)
1465 {
1466  VFN_DEBUG_ENTRY("CellExplorer::GetLSQDeriv()",2)
1467  REAL *par=NULL;
1468  if(refpar.GetName()=="Reciprocal unit cell par #0") par=mRecUnitCell.par+1;
1469  else
1470  if(refpar.GetName()=="Reciprocal unit cell par #1") par=mRecUnitCell.par+2;
1471  else
1472  if(refpar.GetName()=="Reciprocal unit cell par #2") par=mRecUnitCell.par+3;
1473  else
1474  if(refpar.GetName()=="Reciprocal unit cell par #3") par=mRecUnitCell.par+4;
1475  else
1476  if(refpar.GetName()=="Reciprocal unit cell par #4") par=mRecUnitCell.par+5;
1477  else
1478  if(refpar.GetName()=="Reciprocal unit cell par #5") par=mRecUnitCell.par+6;
1479  else
1480  if(refpar.GetName()=="Zero") par=mRecUnitCell.par+0;
1481  else cout<<__FILE__<<":"<<__LINE__<<":Parameter not found:"<<refpar.GetName()<<endl;
1482  unsigned int j=0;
1483  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();++pos)
1484  {
1485  VFN_DEBUG_MESSAGE("CellExplorer::GetLSQDeriv():"<<j<<"/"<<mpPeakList->GetPeakList().size(),2)
1486  VFN_DEBUG_MESSAGE("CellExplorer::GetLSQDeriv():"<<pos->h<<","<<pos->k<<","<<pos->l,2)
1487  if(pos->isIndexed)
1488  mDeriv(j++)=mRecUnitCell.hkl2d(pos->h,pos->k,pos->l,par);
1489  }
1490  VFN_DEBUG_EXIT("CellExplorer::GetLSQDeriv()",2)
1491  return mDeriv;
1492 }
1493 
1494 void CellExplorer::BeginOptimization(const bool allowApproximations, const bool enableRestraints)
1495 {
1496  VFN_DEBUG_ENTRY("CellExplorer::BeginOptimization()",5)
1497  Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true,false);
1498  const unsigned int nb=mpPeakList->GetPeakList().size();
1499  mCalc.resize(nb-mNbSpurious);
1500  mObs.resize(nb-mNbSpurious);
1501  mWeight.resize(nb-mNbSpurious);
1502  mDeriv.resize(nb-mNbSpurious);
1503  int j=0;
1504  float thres=0.0;
1505  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();++pos)
1506  if(thres<pos->iobs) thres=pos->iobs;
1507  thres/=10;// weight=1 for intensities up to Imax/10
1508 
1509  //cout <<"Beginning optimization with reflexions:"<<endl;
1510  //char buf[100];
1511  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();++pos)
1512  {
1513  if(pos->isIndexed)
1514  {
1515  mObs(j)=pos->d2obs;
1516  if(mObs(j)>thres) mWeight(j)=1;
1517  else mWeight(j)=mObs(j)/thres;
1518  /*
1519  sprintf(buf,"#%2d (%3d %3d %3d) dobs=%6.3f dcalc=%6.3f iobs=%6.3f weight=%6.4f",
1520  i,mpPeakList->mvHKL[i].h,mpPeakList->mvHKL[i].k,mpPeakList->mvHKL[i].l,
1521  1/mpPeakList->mvdobs[i],1/sqrt(mRecUnitCell.hkl2d(mpPeakList->mvHKL[i].h,mpPeakList->mvHKL[i].k,mpPeakList->mvHKL[i].l)),
1522  mpPeakList->mviobs[i],mWeight(j));
1523  */
1524  j++;
1525  }
1526  /*
1527  else
1528  {
1529  sprintf(buf,"#%2d (%3d %3d %3d) dobs=%6.3f dcalc=%6.3f iobs=%6.3f SPURIOUS",
1530  i,mpPeakList->mvHKL[i].h,mpPeakList->mvHKL[i].k,mpPeakList->mvHKL[i].l,
1531  1/mpPeakList->mvdobs[i],1/sqrt(mRecUnitCell.hkl2d(mpPeakList->mvHKL[i].h,mpPeakList->mvHKL[i].k,mpPeakList->mvHKL[i].l)),
1532  mpPeakList->mviobs[i]);
1533  }
1534  cout<<buf<<endl;
1535  */
1536  }
1537  this->RefinableObj::BeginOptimization(allowApproximations,enableRestraints);
1538  VFN_DEBUG_EXIT("CellExplorer::BeginOptimization()",5)
1539 }
1540 
1541 void CellExplorer::LSQRefine(int nbCycle, bool useLevenbergMarquardt, const bool silent)
1542 {
1543  if(mNbLSQExcept>100)
1544  {
1545  if(!silent) cout<<"CellExplorer::LSQRefine(): LSQ was disabled, too many (>100) exceptions caught. Weird unit cell parameters ?";
1546  return;
1547  }
1548  VFN_DEBUG_ENTRY("CellExplorer::LSQRefine()",5)
1549  mLSQObj.SetRefinedObj(*this);
1550  mLSQObj.PrepareRefParList(true);
1551  //this->BeginOptimization();
1552  //cout<<FormatVertVector<REAL>(this->GetLSQObs(0),this->GetLSQCalc(0),this->GetLSQWeight(0),this->GetLSQDeriv(0,this->GetPar((long)0)))<<endl;
1553  try {mLSQObj.Refine(nbCycle,useLevenbergMarquardt,silent);}
1554  catch(const ObjCrystException &except)
1555  {
1556  if(++mNbLSQExcept>100) cout<<"WARNING CellExplorer::LSQRefine(): LSQ was disabled, too many (>100) exceptions caught. Weird unit cell parameters ?"<<endl ;
1557  }
1558  if(!silent) mpPeakList->Print(cout);
1559  VFN_DEBUG_EXIT("CellExplorer::LSQRefine()",5)
1560 }
1561 
1575 bool DichoIndexed(const PeakList &dhkl, const RecUnitCell &par,const RecUnitCell &dpar,
1576  const unsigned int nbUnindexed=0,const bool verbose=false,unsigned int useStoredHKL=0,
1577  const unsigned int maxNbMissingBelow5=0)
1578 {
1579  const unsigned int nb=dhkl.GetPeakList().size();
1580  int nbIndexed=nb-nbUnindexed;// Number of reflections we require to be indexed
1581  float d5=0;
1582  if(maxNbMissingBelow5>0) d5=dhkl.GetPeakList()[4].d2obs;
1583  // number of missing reflections calculated below 5th observed line
1584  unsigned int nbMissingBelow5=0;
1585  // List of indexed reflections
1586  vector<PeakList::hkl>::const_iterator pos,first,last,end;
1587  if(useStoredHKL==1)
1588  {// We already now possible Miller indices for all reflections
1589  unsigned int nbUnIx = 0;
1590  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
1591  {
1592  pos->isIndexed=false;
1593  for(list<PeakList::hkl0>::const_iterator phkl0=pos->vDicVolHKL.begin();phkl0!=pos->vDicVolHKL.end();++phkl0)
1594  {
1595  float d0,d1;
1596  par.hkl2d_delta(phkl0->h,phkl0->k,phkl0->l,dpar,d0,d1);
1597  if((pos->d2obsmax>=d0) && (d1>=pos->d2obsmin))
1598  {
1599  pos->d2calc=(d0+d1)/2;
1600  pos->isIndexed=true;
1601  if(--nbIndexed==0) return true;
1602  break;
1603  }
1604  }
1605  if(!(pos->isIndexed)) if(++nbUnIx>nbUnindexed) return false;
1606  }
1607  return false;
1608  }
1609  const bool storePossibleHKL=(useStoredHKL==2);
1610 
1611  if(storePossibleHKL)
1612  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
1613  {
1614  pos->isIndexed=false;
1615  pos->vDicVolHKL.clear();
1616  }
1617  else
1618  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos) pos->isIndexed=false;
1619 
1620  int h,k,l;
1621  float dmax=dhkl.GetPeakList()[nb-1].d2obs;
1622  float dmin=dhkl.GetPeakList()[0 ].d2obs;
1623 
1624 
1625  int sk0,sl0;// do we need >0 *and* <0 indices for k,l ?
1626  switch(par.mlattice)
1627  {
1628  case TRICLINIC:
1629  sk0=-1;sl0=-1;
1630  break;
1631  case MONOCLINIC:
1632  {
1633  sk0=1;sl0=-1;
1634  break;
1635  }
1636  case ORTHORHOMBIC:
1637  sk0=1;sl0=1;
1638  break;
1639  case HEXAGONAL:
1640  sk0=-1;sl0=1;
1641  break;
1642  case RHOMBOEDRAL:
1643  sk0=-1;sl0=-1;
1644  break;
1645  case TETRAGONAL:
1646  sk0=1;sl0=1;
1647  break;
1648  case CUBIC:
1649  sk0=1;sl0=1;
1650  break;
1651  // This should never happen. Avoid using unitialized values.
1652  default:
1653  throw 0;
1654  }
1655  int stepk,stepl;// steps in k,l to use for centered lattices
1656  switch(par.mCentering)
1657  {
1658  case LATTICE_P:stepk=1;stepl=1;break;
1659  case LATTICE_I:stepk=1;stepl=2;break;
1660  case LATTICE_A:stepk=1;stepl=2;break;
1661  case LATTICE_B:stepk=1;stepl=2;break;
1662  case LATTICE_C:stepk=2;stepl=1;break;
1663  case LATTICE_F:stepk=2;stepl=2;break;
1664  // This should never happen. Avoid using unitialized values.
1665  default: throw 0;
1666  }
1667  //RecUnitCell par0(par),par1(par);
1668  //for(unsigned int i=0;i<7;++i) {par0.par[i]-=dpar.par[i];par1.par[i]+=dpar.par[i];}
1669 
1670  //currently first & last unindexed dhkl
1671  first=dhkl.GetPeakList().begin(),last=dhkl.GetPeakList().end(),end=dhkl.GetPeakList().end();
1672 
1673  unsigned long nbCalcH,nbCalcK;// Number of calculated lines below dmax for each h,k
1674  for(h=0;;++h)
1675  {
1676  if(verbose) cout<<"H="<<h<<endl;
1677  nbCalcH=0;
1678  for(int sk=sk0;sk<=1;sk+=2)
1679  {
1680  if(h==0) sk=1;
1681  if(stepk==2) k=(h%2);// For LATTICE_C,LATTICE_F: h odd => k odd
1682  else k=0;
1683  for(;;k+=stepk)
1684  {
1685  if(verbose) cout<<"K="<<k*sk<<endl;
1686  nbCalcK=0;
1687  for(int sl=sl0;sl<=1;sl+=2)
1688  {
1689  int l0=0;
1690  if((h+k)==0)
1691  {
1692  sl=1;// No need to list 0 0 l with l<0
1693  l0=1;
1694  }
1695  else
1696  {
1697  if(h==0)
1698  {
1699  if(par.mlattice==MONOCLINIC) sl=1;// 0 k l and 0 k -l are equivalent
1700  if((sk<0)||(sl<0)) l0=1;// Do not list 0 k 0 with k<0
1701  else l0=0;// h==k==0 already covered
1702  }
1703  else
1704  {
1705  if(sl<0) l0=1;// Do not list h k 0 twice
1706  else l0=0;
1707  }
1708  }
1709  if(stepl==2)
1710  {
1711  if(par.mCentering==LATTICE_I) l0+=(h+k+l0)%2;
1712  if(par.mCentering==LATTICE_A) l0+=(k+l0)%2;// Start at k+l even
1713  if( (par.mCentering==LATTICE_B)
1714  ||(par.mCentering==LATTICE_F)) l0+=(h+l0)%2;// Start at h+l even
1715  }
1716  if(verbose) cout<<"SL="<<sl<<", L0="<<l0<<", STEPL="<<stepl<<", Centering="<<par.mCentering<<endl;
1717  for(l=l0;;l+=stepl)
1718  {
1719  if(verbose) cout<<"L="<<l<<","<<sl<<endl;
1720  float d0,d1;
1721  par.hkl2d_delta(h,sk*k,sl*l,dpar,d0,d1);
1722  if(d0<dmax) {nbCalcH++;nbCalcK++;}
1723  if((d1<dmin)&&(maxNbMissingBelow5==0)) continue;
1724  if(d0>dmax)
1725  {
1726  if(par.mlattice==TRICLINIC)
1727  {
1728  // Must check that d is increasing with l, otherwise we still need to increase it
1729  if(verbose) cout<<"L?="<< par.hkl2d(h,sk*k,sl*l,NULL,3)*sl <<", dmax="<<dmax<<endl;
1730  if((par.hkl2d(h,sk*k,sl*l,NULL,3)*sl)>0) break;
1731  }
1732  else break;
1733  }
1734  bool missing=(d0<d5)&&(maxNbMissingBelow5>0);
1735  for(pos=first;pos!=end;++pos)
1736  {
1737  if(pos==last) break;
1738  if((!storePossibleHKL)&&(pos->isIndexed)&&missing) continue;
1739  const float d2obs=pos->d2obs,d2obsmin=pos->d2obsmin, d2obsmax=pos->d2obsmax;
1740  if((d2obsmax>=d0) && (d1>=d2obsmin))
1741  {
1742  missing=false;
1743  if(!(pos->isIndexed))
1744  {
1745  pos->d2calc=(d0+d1)/2;
1746  --nbIndexed;
1747  pos->isIndexed=true;
1748  }
1749  if(verbose) cout<<d1<<" < ? <"<<d0<<"("<<h<<","<<sk*k<<","<<sl*l<<"): "<<d2obs<<" (remaining to index:"<<nbIndexed<<")"<<endl;
1750  if(storePossibleHKL)
1751  pos->vDicVolHKL.push_back(PeakList::hkl0(h,sk*k,sl*l));
1752  else
1753  {
1754  if((!storePossibleHKL)&&(nbIndexed==0)) return true;
1755  if(pos==first){first++;dmin=first->d2obsmin;}
1756  if(pos==last){last--;dmax=last->d2obsmax;}
1757  }
1758  }
1759  }
1760  if(missing) if(++nbMissingBelow5>=maxNbMissingBelow5)return false;
1761  }
1762  }
1763  if(nbCalcK==0) break;// && ((par.hkl2d(h,sk*k,0,NULL,2)*sk)>0)) break; // k too large
1764  }
1765  }
1766  if(nbCalcH==0) break;//h too large
1767  }
1768  if(verbose)
1769  {
1770  dhkl.Print(cout);
1771  }
1772  return nbIndexed<=0;
1773 }
1774 
1775 float CellExplorer::GetBestScore()const{return mBestScore;}
1776 const std::list<std::pair<RecUnitCell,float> >& CellExplorer::GetSolutions()const {return mvSolution;}
1777 std::list<std::pair<RecUnitCell,float> >& CellExplorer::GetSolutions() {return mvSolution;}
1778 
1779 unsigned int CellExplorer::RDicVol(RecUnitCell par0,RecUnitCell dpar, unsigned int depth,unsigned long &nbCalc,const float minV,const float maxV,vector<unsigned int> vdepth)
1780 {
1781  static bool localverbose=false;
1782  if(mlattice==TRICLINIC)
1783  {
1784  const float p1=par0.par[1] , p2=par0.par[2] , p3=par0.par[3] , p4=par0.par[4] , p5=par0.par[5] , p6=par0.par[6];
1785  const float p1m=p1-dpar.par[1], p2m=p2-dpar.par[2], p3m=p3-dpar.par[3], p4m=p4-dpar.par[4], p5m=p5-dpar.par[5], p6m=p6-dpar.par[6];
1786  const float p1p=p1+dpar.par[1], p2p=p2+dpar.par[2], p3p=p3+dpar.par[3], p4p=p4+dpar.par[4], p5p=p5+dpar.par[5], p6p=p6+dpar.par[6];
1787 
1788  // a*<b*<c*
1789  if((p1m>p2p)||(p2m>p3p)) return 0;
1790 
1791  // max/min absolute values for p4,p5,p6
1792  if((p4m>p1p)||(-p4p>p1p)) return 0;//abs(p4)<p1 <=> b* < b*+/-a*
1793  if((p5m>p2p)||(-p5p>p2p)) return 0;//abs(p5)<p2 <=> c* < c*+/-b*
1794  if((p6m>p1p)||(-p6p>p1p)) return 0;//abs(p6)<p1 <=> c* < c*+/-a*
1795 
1796  const float max6=p1p+p2p-p4m-p5m;
1797  if((p6m>max6)||(-p6p>max6)) return 0;//abs(p6)<p1+p2-p4-p5 <=> c* < c*+/-a*+/-b*
1798 
1799  float p6mm,p6pp,p5mm,p5pp,p4mm,p4pp; // p6pp: smaller V*, larger V, etc..
1800  // derivative of V*^2 with p6
1801  if((p4*p5-2*p2*p6)>0) {p6pp=p6p;p6mm=p6m;}
1802  else{p6pp=p6m;p6mm=p6p;}
1803  // derivative of V*^2 with p5
1804  if((p4*p6-2*p1*p5)>0) {p5pp=p5p;p5mm=p5m;}
1805  else{p5pp=p5m;p5mm=p5p;}
1806  // derivative of V*^2 with p5
1807  if((p5*p6-2*p3*p4)>0) {p4pp=p4p;p4mm=p4m;}
1808  else{p4pp=p4m;p4mm=p4p;}
1809 
1810  //const float vmin0=1/sqrt(abs(p1p*p2p*p3p*(1-p5mm*p5mm/(4*p2p*p3p)-p6mm*p6mm/(4*p1p*p3p)-p4mm*p4mm/(4*p1p*p2p)+p4mm*p5mm*p6mm/(4*p1m*p2m*p3m))));
1811  //const float vmax0=1/sqrt(abs(p1m*p2m*p3m*(1-p5pp*p5pp/(4*p2m*p3m)-p6pp*p6pp/(4*p1m*p3m)-p4pp*p4pp/(4*p1m*p2m)+p4pp*p5pp*p6pp/(4*p1m*p2m*p3m))));
1812  const float vmin0=1/sqrt(abs(p1p*p2p*p3p*(1-p5mm*p5mm/(4*p2p*p3p)-p6mm*p6mm/(4*p1p*p3p)-p4mm*p4mm/(4*p1p*p2p)+p4mm*p5mm*p6mm/(4*p1m*p2m*p3m))));
1813  const float vmax0=1/sqrt(abs(p1m*p2m*p3m*(1-p5pp*p5pp/(4*p2m*p3m)-p6pp*p6pp/(4*p1m*p3m)-p4pp*p4pp/(4*p1m*p2m)+p4pp*p5pp*p6pp/(4*p1m*p2m*p3m))));
1814  if((vmin0>maxV)||(vmax0<minV)) return 0;
1815  }
1816  else
1817  if((depth>0)&&(depth<=2))// test if volume is within range
1818  {
1819  RecUnitCell parm=par0,parp=par0;
1820  for(unsigned int i=0;i<6;++i) {parm.par[i]-=dpar.par[i];parp.par[i]+=dpar.par[i];}
1821  vector<float> parmd=parm.DirectUnitCell();
1822  vector<float> parpd=parp.DirectUnitCell();
1823  if((parpd[6]>maxV)||(parmd[6]<minV))return 0;
1824  }
1825  unsigned int useStoredHKL=1;//Use already stored hkl
1826  if(depth==0) useStoredHKL=2; //Store possible hkl for all observed lines
1827 
1828  unsigned int maxMissingBelow5=0;
1829  // In the triclinic case, accept a maximum of 5 missing reflections below the 5th observed line
1830  if(mlattice==TRICLINIC) maxMissingBelow5=5;
1831 
1832  bool indexed=DichoIndexed(*mpPeakList,par0,dpar,mNbSpurious,localverbose,useStoredHKL,maxMissingBelow5);
1833 
1834  #if 0
1835  // If indexation failed but depth>=4, try adding a zero ?
1836  if( (!indexed) && (depth>=4))
1837  {//:TODO: Check if this is OK ! Vary value with depth
1838  dpar.par[0]=.0001;
1839  indexed=DichoIndexed(*mpPeakList,par0,dpar,mNbSpurious,false,useStoredHKL,maxMissingBelow5);
1840  //if(indexed) cout<<"Added zero - SUCCESS !"<<endl;
1841  }
1842  #endif
1843 
1844  if((indexed)&&(useStoredHKL==2))
1845  {
1846  // Test if two successive lines have been indexed exclusively with the same hkl
1847  unsigned int nbident=0;
1848  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();)
1849  {
1850  if(pos->vDicVolHKL.size()==1)
1851  {
1852  const PeakList::hkl0 h0=pos->vDicVolHKL.front();
1853  if(++pos==mpPeakList->GetPeakList().end()) break;
1854  if(pos->vDicVolHKL.size()==1)
1855  {
1856  const PeakList::hkl0 h1=pos->vDicVolHKL.front();
1857  if((h0.h==h1.h)&&(h0.k==h1.k)&&(h0.l==h1.l)) nbident++;
1858  if(nbident>mNbSpurious) {indexed=false;break;}
1859  }
1860  }
1861  else ++pos;
1862  }
1863  }
1864 
1865  // if we can zoom in for one parameter directly, we need per-parameter depth
1866  if(vdepth.size()==0)
1867  {
1868  vdepth.resize(mnpar-1);
1869  for(vector<unsigned int>::iterator pos=vdepth.begin();pos!=vdepth.end();) *pos++=depth;
1870  }
1871  else
1872  for(vector<unsigned int>::iterator pos=vdepth.begin();pos!=vdepth.end();++pos) if(*pos<depth)*pos=depth;
1873  #if 1
1874  if(false)//((useStoredHKL==2)&&(mNbSpurious==0)&&indexed)
1875  { // If high-d lines have been associated to a single reflection which is either h00, 0k0 or 00l,
1876  // jump the corresponding parameter to higher depth (mDicVolDepthReport, lowest depth report) immediately
1877  vector<pair<unsigned int,float> > vq0(3);
1878  for(unsigned int i=0;i<3;++i) {vq0[i].first=0;vq0[i].second=0.0;}
1879  RecUnitCell par0orig=par0,dparorig=dpar;
1880  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();++pos)
1881  {
1882  if(pos->vDicVolHKL.size()==1)
1883  {
1884  const PeakList::hkl0 h0=pos->vDicVolHKL.front();
1885  if((h0.k==0)&&(h0.l==0)) {vq0[0].first+=1 ; vq0[0].second+=pos->dobs/h0.h;}
1886  else
1887  if((h0.h==0)&&(h0.l==0)) {vq0[1].first+=1 ; vq0[1].second+=pos->dobs/h0.k;}
1888  else
1889  if((h0.h==0)&&(h0.k==0)) {vq0[2].first+=1 ; vq0[2].second+=pos->dobs/h0.l;if(localverbose) cout<<h0.h<<" "<<h0.k<<" "<<h0.l<<": d="<<pos->dobs<<endl;}
1890  }
1891  }
1892  switch(mlattice)
1893  {
1894  case TRICLINIC:
1895  {// In the triclinic case we may start with p1 and p2 already at depth>0 (triclinic quick tests)
1896  if(vq0[0].first>0) {par0.par[1]=vq0[0].second/vq0[0].first ; dpar.par[1]*=pow((float)0.5,float(mDicVolDepthReport-vdepth[0]));vdepth[0]=mDicVolDepthReport;}
1897  if(vq0[1].first>0) {par0.par[2]=vq0[1].second/vq0[1].first ; dpar.par[2]*=pow((float)0.5,float(mDicVolDepthReport-vdepth[1]));vdepth[1]=mDicVolDepthReport;}
1898  if(vq0[2].first>0) {par0.par[3]=vq0[2].second/vq0[2].first ; dpar.par[3]*=pow((float)0.5,float(mDicVolDepthReport-vdepth[2]));vdepth[2]=mDicVolDepthReport;}
1899  break;
1900  }
1901  case MONOCLINIC:
1902  {
1903  if(vq0[0].first>0) {par0.par[1]=vq0[0].second/vq0[0].first ; vdepth[0]=mDicVolDepthReport;dpar.par[1]*=.0625;}
1904  if(vq0[1].first>0) {par0.par[2]=vq0[1].second/vq0[1].first ; vdepth[1]=mDicVolDepthReport;dpar.par[2]*=.0625;}
1905  if(vq0[2].first>0) {par0.par[3]=vq0[2].second/vq0[2].first ; vdepth[2]=mDicVolDepthReport;dpar.par[3]*=.0625;}
1906  break;
1907  }
1908  case ORTHORHOMBIC:
1909  {
1910  if(vq0[0].first>0) {par0.par[1]=vq0[0].second/vq0[0].first ; vdepth[0]=mDicVolDepthReport;dpar.par[1]*=.0625;}//pow((float)0.5,(int)(mDicVolDepthReport-depth));}
1911  if(vq0[1].first>0) {par0.par[2]=vq0[1].second/vq0[1].first ; vdepth[1]=mDicVolDepthReport;dpar.par[2]*=.0625;}//pow((float)0.5,(int)(mDicVolDepthReport-depth));}
1912  if(vq0[2].first>0) {par0.par[3]=vq0[2].second/vq0[2].first ; vdepth[2]=mDicVolDepthReport;dpar.par[3]*=.0625;}//pow((float)0.5,(int)(mDicVolDepthReport-depth));}
1913  break;
1914  }
1915  case HEXAGONAL:break;
1916  case RHOMBOEDRAL:break;
1917  case TETRAGONAL:break;
1918  case CUBIC:break;
1919  }
1920  // If all parameters are at a higher depth, jump the global depth immediately
1921  unsigned int newdepth=40;
1922  for(vector<unsigned int>::iterator pos=vdepth.begin();pos!=vdepth.end();++pos) if(*pos<newdepth) newdepth=*pos;
1923  if(newdepth>depth) depth=newdepth;
1924  if((vq0[0].first>0)||(vq0[1].first>0)||(vq0[2].first>0))
1925  {
1926  indexed=DichoIndexed(*mpPeakList,par0,dpar,mNbSpurious,false,1,maxMissingBelow5);
1927  if(false)
1928  {
1929  {
1930  RecUnitCell parm=par0orig,parp=par0;
1931  for(unsigned int i=0;i<6;++i) {parm.par[i]-=dparorig.par[i];parp.par[i]+=dparorig.par[i];}
1932  vector<float> parmd=parm.DirectUnitCell();
1933  vector<float> parpd=parp.DirectUnitCell();
1934  char buf[200];
1935  sprintf(buf,"orig: a=%5.2f-%5.2f b=%5.2f-%5.2f c=%5.2f-%5.2f alpha=%5.2f-%5.2f beta=%5.2f-%5.2f gamma=%5.2f-%5.2f V=%5.2f-%5.2f",
1936  parpd[0],parmd[0],parpd[1],parmd[1],parpd[2],parmd[2],parpd[3]*RAD2DEG,parmd[3]*RAD2DEG,
1937  parpd[4]*RAD2DEG,parmd[4]*RAD2DEG,parpd[5]*RAD2DEG,parmd[5]*RAD2DEG,parpd[6],parmd[6]);
1938  for(unsigned int i = 0; i < depth; ++i) cout << " ";
1939  cout<<buf<<"level="<<depth<<", indexed="<<indexed<<endl;
1940  }
1941  {
1942  RecUnitCell parm=par0,parp=par0;
1943  for(unsigned int i=0;i<6;++i) {parm.par[i]-=dpar.par[i];parp.par[i]+=dpar.par[i];}
1944  vector<float> parmd=parm.DirectUnitCell();
1945  vector<float> parpd=parp.DirectUnitCell();
1946  char buf[200];
1947  sprintf(buf,"bypass: a=%5.2f-%5.2f b=%5.2f-%5.2f c=%5.2f-%5.2f alpha=%5.2f-%5.2f beta=%5.2f-%5.2f gamma=%5.2f-%5.2f V=%5.2f-%5.2f",
1948  parpd[0],parmd[0],parpd[1],parmd[1],parpd[2],parmd[2],parpd[3]*RAD2DEG,parmd[3]*RAD2DEG,
1949  parpd[4]*RAD2DEG,parmd[4]*RAD2DEG,parpd[5]*RAD2DEG,parmd[5]*RAD2DEG,parpd[6],parmd[6]);
1950  for(unsigned int i = 0; i < depth; ++i) cout << " ";
1951  cout<<buf<<"level="<<depth<<", indexed="<<indexed<<endl;
1952  }
1953  }
1954  }
1955  }
1956  #endif
1957  if(false)//(depth==1)&&(rand()%10==0))
1958  {
1959  RecUnitCell parm=par0,parp=par0;
1960  for(unsigned int i=0;i<4;++i) {parm.par[i]-=dpar.par[i];parp.par[i]+=dpar.par[i];}
1961  for(unsigned int i=4;i<7;++i) {parm.par[i]+=dpar.par[i];parp.par[i]-=dpar.par[i];}
1962  vector<float> parmd=parm.DirectUnitCell();
1963  vector<float> parpd=parp.DirectUnitCell();
1964  char buf[200];
1965  sprintf(buf,"a=%5.2f-%5.2f b=%5.2f-%5.2f c=%5.2f-%5.2f alpha=%6.2f-%6.2f beta=%6.2f-%6.2f gamma=%6.2f-%6.2f V=%6.2f-%6.2f",
1966  parpd[0],parmd[0],parpd[1],parmd[1],parpd[2],parmd[2],parpd[3]*RAD2DEG,parmd[3]*RAD2DEG,
1967  parpd[4]*RAD2DEG,parmd[4]*RAD2DEG,parpd[5]*RAD2DEG,parmd[5]*RAD2DEG,parpd[6],parmd[6]);
1968  for(unsigned int i = 0; i < depth; ++i) cout << " ";
1969  cout<<buf<<"level="<<depth<<", indexed="<<indexed<<"("<<mvSolution.size()<<" sol.)"<<endl;
1970  }
1971  nbCalc++;
1972  // :TODO: if we failed the dichotomy and reached some depth, try guessing a zero shift from the indexed reflections
1973  /*
1974  if((!indexed)&&(depth>=2))
1975  {
1976  vector<float> shifts(mpPeakList->GetPeakList().size());
1977  vector<PeakList::hkl>::const_iterator peakpos=mpPeakList->GetPeakList().begin();
1978  for(vector<float>::iterator spos=shifts.begin();spos!=shifts.end();)
1979  { *spos++ = peakpos->d2diff * (float)(peakpos->isIndexed&&(!peakpos->isSpurious));peakpos++;}
1980  sort(shifts.begin(),shifts.end());
1981  par0.par[0]=shifts[mpPeakList->GetPeakList().size()/2];//use median value
1982  indexed=DichoIndexed(*mpPeakList,par0,dpar,mNbSpurious);
1983  if(indexed) cout<<"Failed Dicho ? Trying auto-zero shifting :Worked !"<<endl;
1984  }
1985  */
1986  if(indexed)
1987  {
1988  unsigned int deeperSolutions=0;
1989  if(depth<mMaxDicVolDepth)
1990  {
1991  if(false)//depth>=5)
1992  {
1993  RecUnitCell parm=par0,parp=par0;
1994  for(unsigned int i=0;i<6;++i) {parm.par[i]-=dpar.par[i];parp.par[i]+=dpar.par[i];}
1995  vector<float> parmd=parm.DirectUnitCell();
1996  vector<float> parpd=parp.DirectUnitCell();
1997  char buf[200];
1998  sprintf(buf,"a=%5.2f-%5.2f b=%5.2f-%5.2f c=%5.2f-%5.2f alpha=%5.2f-%5.2f beta=%5.2f-%5.2f gamma=%5.2f-%5.2f V=%5.2f-%5.2f",
1999  parpd[0],parmd[0],parpd[1],parmd[1],parpd[2],parmd[2],parpd[3]*RAD2DEG,parmd[3]*RAD2DEG,
2000  parpd[4]*RAD2DEG,parmd[4]*RAD2DEG,parpd[5]*RAD2DEG,parmd[5]*RAD2DEG,parpd[6],parmd[6]);
2001  for(unsigned int i=0;i<depth;++i) cout<<" ";
2002  cout<<"Depth="<<depth<<" "<<buf<<endl;
2003  for(int i=0;i<=6;++i)cout<<par0.par[i]<<",";
2004  for(int i=0;i<=6;++i)cout<<dpar.par[i]<<",";
2005  cout<<endl;
2006  }
2007  RecUnitCell par=par0;
2008  // zero (if used...)
2009  dpar.par[0]=0.5*dpar.par[0];
2010  // Divide interval by 2, except if this parameter is already at a higher depth
2011  // because a main axis has been indexed already.
2012  for(unsigned int i=1;i<mnpar;++i) dpar.par[i]*=(0.5+0.5*(vdepth[i-1]>depth));
2013 
2014  for(int i0=-1;i0<=1;i0+=2)
2015  {
2016  //:TODO: dichotomy on zero shift ?
2017  if(localverbose) cout<<__FILE__<<":"<<__LINE__<<":"<<par.par[3]<<" +/- "<<dpar.par[3]<<" ("<<vdepth[2]<<")"<<endl;
2018  // Don't change parameter if it is already determined at a higher depth
2019  if(vdepth[0]==depth) {par.par[1]=par0.par[1]+i0*dpar.par[1];}
2020  else {i0=2;}// no need to dicho this parameter which is already at higher depth
2021  if(mnpar==2)
2022  deeperSolutions+=RDicVol(par,dpar, depth+1,nbCalc,minV,maxV,vdepth);
2023  else
2024  for(int i1=-1;i1<=1;i1+=2)
2025  {
2026  if(vdepth[1]==depth) {par.par[2]=par0.par[2]+i1*dpar.par[2];}
2027  else {i1=2;}// no need to dicho this parameter which is already at higher depth
2028  if(mnpar==3)
2029  deeperSolutions+=RDicVol(par,dpar, depth+1,nbCalc,minV,maxV,vdepth);
2030  else
2031  for(int i2=-1;i2<=1;i2+=2)
2032  {
2033  if(vdepth[2]==depth) {par.par[3]=par0.par[3]+i2*dpar.par[3];}
2034  else {i2=2;}// no need to dicho this parameter which is already at higher depth
2035  if(mnpar==4)
2036  deeperSolutions+=RDicVol(par,dpar, depth+1,nbCalc,minV,maxV,vdepth);
2037  else
2038  for(int i3=-1;i3<=1;i3+=2)
2039  {
2040  if(vdepth[3]==depth)par.par[4]=par0.par[4]+i3*dpar.par[4];
2041  else i3=2;
2042  if(mnpar==5)
2043  deeperSolutions+=RDicVol(par,dpar, depth+1,nbCalc,minV,maxV,vdepth);
2044  else
2045  for(int i4=-1;i4<=1;i4+=2)
2046  {
2047  par.par[5]=par0.par[5]+i4*dpar.par[5];
2048  //if(mnpar==7)
2049  // deeperSolutions+=RDicVol(par,dpar, depth+1,nbCalc,minV,maxV,vdepth);
2050  //else
2051  for(int i5=-1;i5<=1;i5+=2)
2052  {
2053  par.par[6]=par0.par[6]+i5*dpar.par[6];
2054  //if(localverbose) cout<<__FILE__<<":"<<__LINE__<<":"<<par.par[3]<<" +/- "<<dpar.par[3]<<" ("<<vdepth[2]<<")"<<endl;
2055  deeperSolutions+=RDicVol(par,dpar, depth+1,nbCalc,minV,maxV,vdepth);
2056  }
2057  }
2058  }
2059  }
2060  }
2061  }
2062  }
2063  if((deeperSolutions==0) &&(depth>=mDicVolDepthReport))
2064  {
2065  mRecUnitCell=par0;
2066  vector<float> par=mRecUnitCell.DirectUnitCell();
2067  float score=Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true,false);
2068  // If we already have enough reports at higher depths (depth+2), don't bother record this one
2069  bool report=true;
2070  if(depth<(mMaxDicVolDepth-1))
2071  if(mvNbSolutionDepth[depth+2]>100)report=false;
2072  if(report && (((score>(mMinScoreReport*.5))&&(depth>=mDicVolDepthReport)) || (depth>=mMaxDicVolDepth)))
2073  {
2074  if(false)//score>) mBestScore//((score>mMinScoreReport)||(depth>=mDicVolDepthReport))
2075  cout<<__FILE__<<":"<<__LINE__<<" Depth="<<depth<<" (DIC) ! a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]<<", alpha="
2076  <<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG<<", V="<<par[6]
2077  <<", score="<<score<<endl;
2078  this->LSQRefine(5,true,true);
2079 
2080  // Re-score (may change to a better hkl indexing), and refine again
2081  score=Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true,false);
2082  this->LSQRefine(5,true,true);
2083 
2084  par=mRecUnitCell.DirectUnitCell();
2085  score=Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true,false);
2086  if( ((score>mMinScoreReport)||(depth>=mDicVolDepthReport))
2087  &&((mvSolution.size()<50)||(score>(mBestScore/3)))
2088  &&((mvSolution.size()<50)||(score>mMinScoreReport)))
2089  {
2090  if((score>(mBestScore))||((score>(mBestScore*0.8))&&(mvSolution.size()<50)))//||(rand()%100==0))
2091  {
2092  char buf[200];
2093  {
2094  RecUnitCell parm=par0,parp=par0;
2095  for(unsigned int i=0;i<4;++i) {parm.par[i]-=dpar.par[i];parp.par[i]+=dpar.par[i];}
2096  for(unsigned int i=4;i<7;++i) {parm.par[i]+=dpar.par[i];parp.par[i]-=dpar.par[i];}
2097  vector<float> parmd=parm.DirectUnitCell();
2098  vector<float> parpd=parp.DirectUnitCell();
2099  sprintf(buf,"a=%5.2f-%5.2f b=%5.2f-%5.2f c=%5.2f-%5.2f alpha=%6.2f-%6.2f beta=%6.2f-%6.2f gamma=%6.2f-%6.2f V=%6.2f-%6.2f",
2100  parpd[0],parmd[0],parpd[1],parmd[1],parpd[2],parmd[2],parpd[3]*RAD2DEG,parmd[3]*RAD2DEG,
2101  parpd[4]*RAD2DEG,parmd[4]*RAD2DEG,parpd[5]*RAD2DEG,parmd[5]*RAD2DEG,parpd[6],parmd[6]);
2102  for(unsigned int i = 0; i < depth; ++i) cout << " ";
2103 
2104  cout<<buf<<"level="<<depth<<", indexed="<<indexed<<"("<<mvSolution.size()<<" sol.)"<<endl;
2105  sprintf(buf,"a=%7.5f-%7.5f b=%7.5f-%7.5f c=%7.5f-%7.5f alpha=%7.5f-%7.5f beta=%7.5f-%7.5f gamma=%7.5f-%7.5f",
2106  parp.par[1],parm.par[1],parp.par[2],parm.par[2],parp.par[3],parm.par[3],parp.par[4],parm.par[4],
2107  parp.par[5],parm.par[5],parp.par[6],parm.par[6]);
2108  for(unsigned int i = 0; i < depth; ++i) cout << " ";
2109  cout<<buf<<"level="<<depth<<", indexed="<<indexed<<"("<<mvSolution.size()<<" sol.)"<<endl;
2110  }
2111  sprintf(buf," Solution ? a=%7.3f b=%7.3f c=%7.3f alpha=%7.3f beta=%7.3f gamma=%7.3f V=%8.2f score=%6.2f #%4lu",
2112  par[0],par[1],par[2],par[3]*RAD2DEG,par[4]*RAD2DEG,par[5]*RAD2DEG,par[6],score,mvSolution.size());
2113  cout<<buf<<endl;
2114  mBestScore=score;
2115  }
2116  mvSolution.push_back(make_pair(mRecUnitCell,score));
2117  mvSolution.back().first.mNbSpurious = mNbSpurious;
2118  mvNbSolutionDepth[depth]+=1;
2119  if((mvSolution.size()>1100)&&(rand()%1000==0))
2120  {
2121  cout<<mvSolution.size()<<" solutions ! Redparing..."<<endl;
2122  this->ReduceSolutions(true);// This will update the min report score
2123  cout<<"-> "<<mvSolution.size()<<" remaining"<<endl;
2124  }
2125  }
2126  }
2127  }
2128  return deeperSolutions+1;
2129  }
2130  return 0;
2131 }
2132 
2133 vector<float> linspace(float min, float max,unsigned int nb)
2134 {
2135  vector<float> v(nb);
2136  for(unsigned int i=0;i<nb;++i) v[i]=min+(max-min)*i/(nb-1);
2137  return v;
2138 }
2139 
2140 void CellExplorer::DicVol(const float minScore,const unsigned int minDepth,const float stopOnScore,const unsigned int stopOnDepth)
2141 {
2142  mNbLSQExcept=0;
2143  mDicVolDepthReport=minDepth;
2144  mMinScoreReport=minScore;
2145  this->Init();
2146  if(minDepth>mMaxDicVolDepth) mMaxDicVolDepth=minDepth;
2147  mvNbSolutionDepth.resize(mMaxDicVolDepth+1);
2148  for(unsigned int i=0;i<=mMaxDicVolDepth;++i) mvNbSolutionDepth[i]=0;
2149 
2150  float latstep=0.5,
2151  vstep=(mVolumeMax-mVolumeMin)/(ceil((mVolumeMax-mVolumeMin)/500)-0.0001);
2152  mCosAngMax=abs(cos(mAngleMax));
2153  const float cosangstep=mCosAngMax/(ceil(mCosAngMax/.08)-.0001);
2154  if(((mVolumeMax-mVolumeMin)/vstep)>10) vstep=(mVolumeMax-mVolumeMin)/9.999;
2155  if(((mLengthMax-mLengthMin)/latstep)>25) latstep=(mLengthMax-mLengthMin)/24.9999;
2156 
2157  cout<<mLengthMin<<"->"<<mLengthMax<<","<<latstep<<","<<(mLengthMax-mLengthMin)/latstep<<endl;
2158  cout<<mAngleMin<<"->"<<mAngleMax<<","<<cosangstep<<","<<mCosAngMax<<","<<(mAngleMax-mAngleMin)/cosangstep<<endl;
2159  cout<<mVolumeMin<<"->"<<mVolumeMax<<","<<vstep<<","<<(mVolumeMax-mVolumeMin)/vstep<<endl;
2160  RecUnitCell par0,dpar;
2161  par0.mlattice=mlattice;
2162  dpar.mlattice=mlattice;
2163  par0.mCentering=mCentering;
2164  dpar.mCentering=mCentering;
2165  //Zero shift parameter - not used for dicvol right now ? :TODO:
2166  par0.par[0]=0.0;
2167  dpar.par[0]=0.0;
2168  unsigned long nbCalc=0;
2169  Chronometer chrono;
2170  float bestscore=0;
2171  list<pair<RecUnitCell,float> >::iterator bestpos;
2172  bool breakDepth=false;
2173  // In the triclinic case, first try assigning a* and b* from the first reflections
2174  if(false) //mlattice==TRICLINIC)
2175  for(float minv=mVolumeMin;minv<mVolumeMax;minv+=vstep)
2176  {
2177  float maxv=minv+vstep;
2178  if(maxv>mVolumeMax)maxv=mVolumeMax;
2179  cout<<"Starting: V="<<minv<<"->"<<maxv<<endl;
2180  const float minr=1/(mLengthMax*mLengthMax);
2181  const float maxr=1/(mLengthMin*mLengthMin);
2182  const float stepr=(maxr-minr)/24.999;
2183  float p1,p2;
2184  for(unsigned int i=0;i<=5;i++)
2185  {
2186  switch(i)
2187  {// Try to find a and b from the first observed reflections
2188  case 0: p1=mpPeakList->GetPeakList()[0].d2obs ;p2=mpPeakList->GetPeakList()[1].d2obs ; break;
2189  case 1: p1=mpPeakList->GetPeakList()[0].d2obs ;p2=mpPeakList->GetPeakList()[2].d2obs ; break;
2190  case 2: p1=mpPeakList->GetPeakList()[1].d2obs/2;p2=mpPeakList->GetPeakList()[0].d2obs ; break;
2191  case 3: p1=mpPeakList->GetPeakList()[1].d2obs/2;p2=mpPeakList->GetPeakList()[2].d2obs ; break;
2192  case 4: p1=mpPeakList->GetPeakList()[2].d2obs/2;p2=mpPeakList->GetPeakList()[0].d2obs ; break;
2193  case 5: p1=mpPeakList->GetPeakList()[2].d2obs/2;p2=mpPeakList->GetPeakList()[1].d2obs ; break;
2194  }
2195  //if(i>0) exit(0);
2196  if(p1>p2) continue;
2197  cout<<"Trying #"<<i<<": a*="<<p1<<", b*="<<p2<<endl;
2198  float min3r=p2,
2199  max3r=maxr;//:TODO: use larger value to take angles into account ?
2200  const float step3r=(max3r-min3r)/(ceil((max3r-min3r)/stepr)-.001);
2201  vector<unsigned int> vdepth(mnpar-1);
2202  for(vector<unsigned int>::iterator pos=vdepth.begin();pos!=vdepth.end();) *pos++=0;
2203  vdepth[0]=3;
2204  vdepth[1]=3;
2205  for(float p3=min3r;p3<max3r;p3+=step3r)
2206  {
2207  //cout<<" p3="<<p3<<endl;
2208  const float max4r=p3+step3r;
2209  const float step4r=max4r/(ceil(max4r/stepr)-.001);
2210  for(float p4=0;p4<max4r;p4+=step4r)
2211  {
2212  //cout<<" p4="<<p4<<endl;
2213  float max5r=(p2+stepr);
2214  const float step5r=max5r/(ceil(max5r/stepr)-.001);
2215  for(float p5=0;p5<max5r;p5+=step5r)
2216  {
2217  float max6r=(p1+stepr);
2218  const float step6r=max6r/(ceil(max6r/stepr)-.001);
2219  for(float p6=-max6r;p6<max6r;p6+=step6r)
2220  {
2221  //cout<<" p6="<<p6<<"/"<<p1<<"/"<<p3<<endl;
2222  dpar.par[1]=stepr*pow(float(0.51),int(vdepth[0]));
2223  dpar.par[2]=stepr*pow(float(0.51),int(vdepth[1]));
2224  dpar.par[3]=step3r*0.51;
2225  dpar.par[4]=step4r*0.51;
2226  dpar.par[5]=step5r*0.51;
2227  dpar.par[6]=step6r*0.51;
2228 
2229  par0.par[0]=0;
2230  par0.par[1]=p1;
2231  par0.par[2]=p2;
2232  par0.par[3]=p3+step3r/2;
2233  par0.par[4]=p4+step4r/2;
2234  par0.par[5]=p5+step5r/2;
2235  par0.par[6]=p6+step6r/2;
2236  //for(int i=0;i<=6;++i)cout<<par0.par[i]<<",";
2237  //cout<<endl;
2238  //for(int i=0;i<=6;++i)cout<<dpar.par[i]<<",";
2239  //cout<<endl;
2240  RDicVol(par0,dpar,0,nbCalc,minv,maxv,vdepth);
2241  }
2242  }
2243  }
2244  }
2245  cout<<"Finished trying: a*="<<p1<<" A, b*="<<p2<<" A, "<<nbCalc
2246  <<" unit cells tested, "<<nbCalc/chrono.seconds()<<" tests/s, Elapsed time="
2247  <<chrono.seconds()<<"s, Best score="<<mBestScore<<", "<<stopOnScore<<", "<<breakDepth<<endl;
2248  breakDepth=false;
2249  if(stopOnDepth>0)
2250  for(unsigned int i=stopOnDepth; i<mvNbSolutionDepth.size();++i)
2251  if(mvNbSolutionDepth[i]>1) {breakDepth=true;break;}
2252  if((mBestScore>stopOnScore)&&(breakDepth)) break;
2253  }//cases
2254  cout<<"Finished triclinic QUICK tests for: V="<<minv<<"->"<<maxv<<" A^3, "<<nbCalc
2255  <<" unit cells tested, "<<nbCalc/chrono.seconds()<<" tests/s, Elapsed time="
2256  <<chrono.seconds()<<"s, Best score="<<mBestScore<<endl;
2257  if((mBestScore>stopOnScore)&&(breakDepth)) break;
2258  }//volume
2259  if((mBestScore<stopOnScore)||(!breakDepth))
2260  for(float minv=mVolumeMin;minv<mVolumeMax;minv+=vstep)
2261  {
2262  float maxv=minv+vstep;
2263  if(maxv>mVolumeMax)maxv=mVolumeMax;
2264  cout<<"Starting: V="<<minv<<"->"<<maxv<<endl;
2265  switch(mlattice)
2266  {
2267  case TRICLINIC:
2268  {
2269  const unsigned int nbstep=25;
2270  vector<float> v1=linspace(mLengthMin,mLengthMax,nbstep);
2271  const float lstep=v1[1]-v1[0];
2272  for(unsigned int i1=0;i1<(nbstep-1);++i1)
2273  {
2274  const float p1 =(1/(v1[i1]*v1[i1])+1/(v1[i1+1]*v1[i1+1]))/2;
2275  const float dp1=(1/(v1[i1]*v1[i1])-1/(v1[i1+1]*v1[i1+1]))/2;
2276  //cout<<"p1="<<p1<<endl;
2277  //cout<<(v1[i1+1]-mLengthMin)/lstep+2.001<<endl;
2278  const unsigned int nb2=int((v1[i1+1]-mLengthMin)/lstep+2.001);
2279  vector<float> v2=linspace(mLengthMin,v1[i1+1],nb2);
2280  //for(unsigned int i2=0;i2<nb2;++i2) cout<<1/v2[i2]/v2[i2]<<" ";
2281  //cout<<endl;
2282  for(unsigned int i2=0;i2<(nb2-1);++i2)
2283  {
2284  const float p2 =(1/(v2[i2]*v2[i2])+1/(v2[i2+1]*v2[i2+1]))/2;
2285  const float dp2=(1/(v2[i2]*v2[i2])-1/(v2[i2+1]*v2[i2+1]))/2;
2286  //cout<<" p2="<<p2<<endl;
2287  const unsigned int nb3=int((v2[i2+1]-mLengthMin)/lstep+2.001);
2288  vector<float> v3=linspace(mLengthMin,v2[i2+1],nb3);
2289  for(unsigned int i3=0;i3<(nb3-1);++i3)
2290  {
2291  const float p3 =(1/(v3[i3]*v3[i3])+1/(v3[i3+1]*v3[i3+1]))/2;
2292  const float dp3=(1/(v3[i3]*v3[i3])-1/(v3[i3+1]*v3[i3+1]))/2;
2293 
2294  //const float vmax3=v1[i1+1]*v2[i2+1]*v3[i3+1];
2295  //const float vmin3=v1[i1]*v2[i2]*v3[i3]/5;
2296  const float vmin3=1/sqrt((p1+dp1)*(p2+dp2)*(p3+dp3));
2297  const float vmax3=1/sqrt((p1-dp1)*(p2-dp2)*(p3-dp3))*2; // *2 - sufficient margin to take into account angles
2298  if(vmax3<minv) continue;
2299  if(vmin3>maxv) continue;
2300 
2301  //cout<<" p3="<<p3<<endl;
2302  //char buf[200];
2303  //sprintf(buf,"a1=%6.4f +/-%6.4f (%6.3f-%6.3f) a2=%6.4f +/-%6.4f (%6.3f-%6.3f) a3=%6.4f +/-%6.4f (%6.3f-%6.3f), V=%6.2f-%6.2f",
2304  // p1,dp1,1/sqrt(p1+dp1),1/sqrt(p1-dp1),p2,dp2,1/sqrt(p2+dp2),1/sqrt(p2-dp2),
2305  // p3,dp3,1/sqrt(p3+dp3),1/sqrt(p3-dp3),vmin3,vmax3);
2306  //cout<<buf<<endl;
2307  const unsigned int nb4=int((p1+dp1)/(4*dp1)+2.001);
2308  vector<float> v4=linspace(0,p1+dp1,nb4);
2309  for(unsigned int i4=0;i4<(nb4-1);++i4)
2310  {
2311  const float p4 =(v4[i4+1]+v4[i4])/2;
2312  const float dp4=(v4[i4+1]-v4[i4])/2;
2313  //cout<<" p4="<<p4<<endl;
2314  const unsigned int nb5=int((p2+dp2)/(4*dp2)+2.001);
2315  vector<float> v5=linspace(0,p2+dp2,nb5);
2316  for(unsigned int i5=0;i5<(nb5-1);++i5)
2317  {
2318  const float p5 =(v5[i5+1]+v5[i5])/2;
2319  const float dp5=(v5[i5+1]-v5[i5])/2;
2320  //cout<<" p5="<<p5<<endl;
2321 
2322  float vmax6=(p1+dp1)+(p2+dp2)-(p4-dp4)-(p5-dp5);
2323  if(vmax6>(p1+dp1)) vmax6=p1+dp1;
2324  if(vmax6<0) continue;
2325  const unsigned int nb6=int((2*vmax6)/(4*dp1)+2.001);
2326  vector<float> v6=linspace(-vmax6,vmax6,nb6);
2327  //cout<<" p6="<<nb6<<","<<vmax6<<endl;
2328  for(unsigned int i6=0;i6<(nb6-1);++i6)
2329  {
2330  const float p6 =(v6[i6+1]+v6[i6])/2;
2331  const float dp6=(v6[i6+1]-v6[i6])/2;
2332  //cout<<" p6="<<p6<<endl;
2333 
2334  //char buf[200];
2335  //sprintf(buf,"%6.4f-%6.4f %6.4f-%6.4f %6.4f-%6.4f %6.4f-%6.4f %6.4f-%6.4f %6.4f-%6.4f ",
2336  // p1-dp1,p1+dp1,p2-dp2,p2+dp2,p3-dp3,p3+dp3,p4-dp4,p4+dp4,p5-dp5,p5+dp5,p6-dp6,p6+dp6);
2337  //cout<<"Testing: "<<buf<<endl;
2338  //cout<<i1<<","<<i2<<","<<i3<<","<<i4<<","<<i5<<","<<i6<<endl;
2339  dpar.par[1]=dp1;
2340  dpar.par[2]=dp2;
2341  dpar.par[3]=dp3;
2342  dpar.par[4]=dp4;
2343  dpar.par[5]=dp5;
2344  dpar.par[6]=dp6;
2345 
2346  par0.par[0]=0;
2347  par0.par[1]=p1;
2348  par0.par[2]=p2;
2349  par0.par[3]=p3;
2350  par0.par[4]=p4;
2351  par0.par[5]=p5;
2352  par0.par[6]=p6;
2353 
2354  RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2355  }
2356  }
2357  }
2358  }
2359  }
2360  }
2361  break;
2362  }
2363  case MONOCLINIC:
2364  {
2365  RecUnitCell parlarge,//parm: smallest reciprocal, largest direct cell
2366  parsmall;//parp: largest reciprocal, smallest direct cell
2367  vector<float> parlarged,parsmalld;
2368  latstep=(mLengthMax-mLengthMin)/24.999;
2369  for(float x4=0;x4<mCosAngMax+cosangstep;x4+=cosangstep)
2370  {
2371  const float sinbeta=sqrt(abs(1-x4*x4));
2372  float x1=mLengthMin;
2373  for(;x1<mLengthMax;x1+=latstep)
2374  {
2375  float x2=mLengthMin;
2376  for(;x2<mLengthMax;x2+=latstep)
2377  {
2378  float x3=x1;
2379  const float x3step=(mLengthMax-x1)/(ceil((mLengthMax-x1)/latstep)-0.001);
2380  for(;x3<mLengthMax;x3+=x3step) //x3+=(latstep+x3*sin4)
2381  {
2382  if((x3*x4)>x1) break;// | c * cos(beta) | <a
2383  dpar.par[1]=(1/(x1)-1/(x1+latstep))*0.5/sinbeta;
2384  dpar.par[2]=(1/(x2)-1/(x2+latstep))*0.5/sinbeta;
2385  dpar.par[3]=(1/(x3)-1/(x3+x3step ))*0.5/sinbeta;
2386  dpar.par[4]=cosangstep*0.5;
2387 
2388  par0.par[0]=0;
2389  par0.par[1]=(1/(x1)+1/(x1+latstep))*0.5/sinbeta;
2390  par0.par[2]=(1/(x2)+1/(x2+latstep))*0.5/sinbeta;
2391  par0.par[3]=(1/(x3)+1/(x3+x3step ))*0.5/sinbeta;
2392  par0.par[4]=x4+cosangstep*.5;
2393 
2394  const float smallv=x1*x2*x3*sinbeta;
2395  if(smallv>maxv) break;
2396  const float largev=(x1+latstep)*(x2+latstep)*(x3+latstep)*(sinbeta+cosangstep);
2397  if(largev<minv) continue;
2398  /*
2399  char buf[200];
2400  sprintf(buf,"a=%5.2f-%5.2f b=%5.2f-%5.2f c=%5.2f-%5.2f alpha=%5.2f-%5.2f beta=%5.2f-%5.2f gamma=%5.2f-%5.2f V=%5.2f-%5.2f",
2401  parsmalld[0],parlarged[0],parsmalld[1],parlarged[1],parsmalld[2],parlarged[2],parsmalld[3]*RAD2DEG,parlarged[3]*RAD2DEG,
2402  parsmalld[4]*RAD2DEG,parlarged[4]*RAD2DEG,parsmalld[5]*RAD2DEG,parlarged[5]*RAD2DEG,parsmalld[6],parlarged[6]);
2403  cout<<buf<<" VM="<<maxv<<", x3="<<x3<<endl;
2404  */
2405  RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2406  }//x3
2407  //if(((parsmalld[6]>maxv)&&(x3==x1))||(parlarged[1]>mLengthMax)) break;
2408  }//x2
2409  }//x1
2410  // Test if we have one solution before going to the next angle range
2411  for(list<pair<RecUnitCell,float> >::iterator pos=mvSolution.begin();pos!=mvSolution.end();++pos)
2412  {
2413  const float score=pos->second;//Score(*mpPeakList,pos->first,mNbSpurious);
2414  if(score>bestscore) {bestscore=score;bestpos=pos;}
2415  }
2416  bool breakDepth=false;
2417  if(stopOnDepth>0)
2418  for(unsigned int i=stopOnDepth; i<mvNbSolutionDepth.size();++i)
2419  if(mvNbSolutionDepth[i]>1) {breakDepth=true;break;}
2420  if((bestscore>stopOnScore)&&(breakDepth)) break;
2421  }//x4
2422  break;
2423  }
2424  case ORTHORHOMBIC:
2425  {
2426  if(false)
2427  {
2428  // Test 7.677350 5.803770 10.313160 V=480
2429  //const float a=7.75,b=5.75,c=10.25;
2430  // Test 6.062000 16.779400 16.8881 v=1750
2431  const float a=6.25,b=16.75,c=16.75;
2432  dpar.par[1]=(1/(a-.25)-1/(a+.25))*0.5;
2433  dpar.par[2]=(1/(b-.25)-1/(b+.25))*0.5;
2434  dpar.par[3]=(1/(c-.25)-1/(c+.25))*0.5;
2435  par0.par[0]=0;
2436  par0.par[1]=1/a;
2437  par0.par[2]=1/b;
2438  par0.par[3]=1/c;
2439  RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2440  break;
2441  }
2442  latstep=(mLengthMax-mLengthMin)/24.999;
2443  for(float x1=mLengthMin;x1<mLengthMax;x1+=latstep)
2444  {
2445  for(float x2=x1;x2<mLengthMax;x2+=latstep)
2446  {
2447  for(float x3=x2;x3<mLengthMax;x3+=latstep)
2448  {
2449  dpar.par[1]=(1/(x1)-1/(x1+latstep))*0.5;
2450  dpar.par[2]=(1/(x2)-1/(x2+latstep))*0.5;
2451  dpar.par[3]=(1/(x3)-1/(x3+latstep))*0.5;
2452 
2453  par0.par[0]=0;
2454  par0.par[1]=(1/(x1)+1/(x1+latstep))*0.5;
2455  par0.par[2]=(1/(x2)+1/(x2+latstep))*0.5;
2456  par0.par[3]=(1/(x3)+1/(x3+latstep))*0.5;
2457 
2458  const float vmin=x1*x2*x3,vmax=(x1+latstep)*(x2+latstep)*(x3+latstep);
2459  if(vmin>maxv) break;
2460  if(vmax>=minv) RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2461  }
2462  if((x1*x2*x2)>maxv) break;
2463  }
2464  if((x1*x1*x1)>maxv) break;
2465  }
2466  break;
2467  }
2468  case HEXAGONAL:
2469  {
2470  vector<float> parlarged,parsmalld;// Small & large UC in direct space
2471  latstep=(mLengthMax-mLengthMin)/24.999;
2472  for(float x1=mLengthMin;;x1+=latstep)
2473  {
2474  for(float x2=mLengthMin;x2<(mLengthMax+latstep);x2+=latstep)
2475  {
2476  dpar.par[1]=(1/(x1)-1/(x1+latstep))*0.5;
2477  dpar.par[2]=(1/(x2)-1/(x2+latstep))*0.5;
2478 
2479  par0.par[0]=0;
2480  par0.par[1]=(1/(x1)+1/(x1+latstep))*0.5;
2481  par0.par[2]=(1/(x2)+1/(x2+latstep))*0.5;
2482 
2483  RecUnitCell parlarge=par0,parsmall=par0;
2484  for(unsigned int i=0;i<6;++i) {parlarge.par[i]-=dpar.par[i];parsmall.par[i]+=dpar.par[i];}
2485  parlarged=parlarge.DirectUnitCell();
2486  parsmalld=parsmall.DirectUnitCell();
2487  /*
2488  char buf[200];
2489  sprintf(buf,"a=%5.2f-%5.2f b=%5.2f-%5.2f c=%5.2f-%5.2f alpha=%5.2f-%5.2f beta=%5.2f-%5.2f gamma=%5.2f-%5.2f V=%5.2f-%5.2f",
2490  parsmalld[0],parlarged[0],parsmalld[1],parlarged[1],parsmalld[2],parlarged[2],parsmalld[3]*RAD2DEG,parlarged[3]*RAD2DEG,
2491  parsmalld[4]*RAD2DEG,parlarged[4]*RAD2DEG,parsmalld[5]*RAD2DEG,parlarged[5]*RAD2DEG,parsmalld[6],parlarged[6]);
2492  */
2493  if((parsmalld[6]<maxv)&&(parlarged[6]>minv))
2494  {
2495  //cout<<buf<<endl;
2496  RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2497  }
2498  //else cout<<buf<<" BREAK"<<endl;
2499  }
2500  if(parlarged[0]>mLengthMax) break;
2501  }
2502  break;
2503  }
2504  case RHOMBOEDRAL: //:TODO:
2505  {
2506  latstep=(mLengthMax-mLengthMin)/24.999;
2507  for(float x1=mLengthMin;x1<(mLengthMax+latstep);x1+=latstep)
2508  {
2509  for(float x2=0;x2<mCosAngMax+cosangstep;x2+=cosangstep)
2510  {
2511  dpar.par[1]=latstep/2*1.1;
2512  dpar.par[2]=cosangstep/2*1.1;
2513 
2514  par0.par[0]=0;
2515  par0.par[1]=x1-latstep/2*1.1;
2516  par0.par[2]=x2-cosangstep/2*1.1;
2517  vector<float> par=par0.DirectUnitCell();
2518  if((par[6]<maxv)&&(par[6]>minv))
2519  {
2520  RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2521  }
2522  }
2523  }
2524  break;
2525  }
2526  case TETRAGONAL:
2527  {
2528  vector<float> parlarged,parsmalld;// Small & large UC in direct space
2529  latstep=(mLengthMax-mLengthMin)/24.999;
2530  for(float x1=mLengthMin;x1<mLengthMax;x1+=latstep)
2531  {
2532  for(float x2=mLengthMin;x2<mLengthMax;x2+=latstep)
2533  {
2534  dpar.par[1]=(1/(x1)-1/(x1+latstep))*0.5;
2535  dpar.par[2]=(1/(x2)-1/(x2+latstep))*0.5;
2536 
2537  par0.par[0]=0;
2538  par0.par[1]=(1/(x1)+1/(x1+latstep))*0.5;
2539  par0.par[2]=(1/(x2)+1/(x2+latstep))*0.5;
2540 
2541  RecUnitCell parlarge=par0,parsmall=par0;
2542  for(unsigned int i=0;i<6;++i) {parlarge.par[i]-=dpar.par[i];parsmall.par[i]+=dpar.par[i];}
2543  parlarged=parlarge.DirectUnitCell();
2544  parsmalld=parsmall.DirectUnitCell();
2545  /*
2546  char buf[200];
2547  sprintf(buf,"a=%5.2f-%5.2f b=%5.2f-%5.2f c=%5.2f-%5.2f alpha=%5.2f-%5.2f beta=%5.2f-%5.2f gamma=%5.2f-%5.2f V=%5.2f-%5.2f",
2548  parsmalld[0],parlarged[0],parsmalld[1],parlarged[1],parsmalld[2],parlarged[2],parsmalld[3]*RAD2DEG,parlarged[3]*RAD2DEG,
2549  parsmalld[4]*RAD2DEG,parlarged[4]*RAD2DEG,parsmalld[5]*RAD2DEG,parlarged[5]*RAD2DEG,parsmalld[6],parlarged[6]);
2550  */
2551  if((parsmalld[6]<maxv)&&(parlarged[6]>minv))
2552  {
2553  RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2554  }
2555  if(parsmalld[6]>maxv) break;
2556  }
2557  if((x1*mLengthMin*mLengthMin)>maxv) break;
2558  }
2559  break;
2560  }
2561  case CUBIC:
2562  {
2563  latstep=(mLengthMax-mLengthMin)/24.999;
2564  cout<<mLengthMax<<","<<mLengthMin<<","<<latstep<<endl;
2565  for(float x1=mLengthMin;x1<(mLengthMax+latstep);x1+=latstep)
2566  {
2567  dpar.par[1]=(1/(x1)-1/(x1+latstep))*0.5;
2568 
2569  par0.par[0]=0;
2570  par0.par[1]=(1/(x1)+1/(x1+latstep))*0.5;
2571 
2572  const float vmin=x1*x1*x1,vmax=(x1+latstep)*(x1+latstep)*(x1+latstep);
2573  if(vmin>maxv)break;
2574  if(vmax>minv) RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2575  }
2576  break;
2577  }
2578  }
2579  cout<<"Finished: V="<<minv<<"->"<<maxv<<" A^3, "<<nbCalc
2580  <<" unit cells tested, "<<nbCalc/chrono.seconds()<<" tests/s, Elapsed time="
2581  <<chrono.seconds()<<"s"<<endl;
2582  for(list<pair<RecUnitCell,float> >::iterator pos=mvSolution.begin();pos!=mvSolution.end();++pos)
2583  {
2584  const float score=pos->second;//Score(*mpPeakList,pos->first,mNbSpurious);
2585  if(score>bestscore) {bestscore=score;bestpos=pos;}
2586  }
2587  bool breakDepth=false;
2588  if(stopOnDepth>0)
2589  for(unsigned int i=stopOnDepth; i<mvNbSolutionDepth.size();++i)
2590  if(mvNbSolutionDepth[i]>1) {breakDepth=true;break;}
2591  if((bestscore>stopOnScore)&&(breakDepth)) break;
2592  }
2593  /*
2594  {// Tag spurious lines
2595  vector<int> vSpuriousScore;
2596  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();++pos)
2597  vSpuriousScore.push_back(pos->stats);
2598  sort(vSpuriousScore.begin(),vSpuriousScore.end());
2599  const int threshold=vSpuriousScore[vSpuriousScore.size()/2]*5;
2600  for(vector<PeakList::hkl>::iterator pos=mpPeakList->mvHKL.begin();pos!=mpPeakList->mvHKL.end();++pos)
2601  if(pos->stats > threshold) pos->isSpurious=true;
2602  else pos->isSpurious=false;
2603  mpPeakList->Print(cout);
2604  }
2605  */
2606  this->ReduceSolutions(true);
2607  bestscore=0;bestpos=mvSolution.end();
2608  for(list<pair<RecUnitCell,float> >::iterator pos=mvSolution.begin();pos!=mvSolution.end();++pos)
2609  {
2610  const float score=Score(*mpPeakList,pos->first,mNbSpurious);
2611  if(score>bestscore) {bestpos=pos;bestscore=score;}
2612  vector<float> par=pos->first.DirectUnitCell();
2613  cout<<__FILE__<<":"<<__LINE__<<" Solution ? a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]
2614  <<", alpha="<<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG
2615  <<", V="<<par[6]<<", score="<<score<<endl;
2616  }
2617  if(bestpos!=mvSolution.end())
2618  {
2619  vector<float> par=bestpos->first.DirectUnitCell();
2620  cout<<__FILE__<<":"<<__LINE__<<" BEST ? a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]
2621  <<", alpha="<<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG
2622  <<", V="<<par[6]<<", score="<<bestscore<<endl;
2623  cout<<nbCalc<<"unit cells tested, "<<nbCalc/chrono.seconds()<<" tests/s, Elapsed time="
2624  <<chrono.seconds()<<"s"<<endl;
2625  }
2626 }
2627 
2628 bool SimilarRUC(const RecUnitCell &c0,const RecUnitCell &c1, const float delta=0.005)
2629 {
2630  if(c0.mNbSpurious != c1.mNbSpurious) return false;
2631  vector<float> par0=c0.DirectUnitCell();
2632  vector<float> par1=c1.DirectUnitCell();
2633  float diff=0;
2634  for(unsigned int i=0;i<6;++i) diff += abs(par0[i]-par1[i]);
2635  return (diff/6)<delta;
2636 }
2637 
2638 bool compareRUCScore(std::pair<RecUnitCell,float> &p1, std::pair<RecUnitCell,float> &p2)
2639 {
2640  return p1.second > p2.second;
2641 }
2642 
2643 void CellExplorer::ReduceSolutions(const bool updateReportThreshold)
2644 {
2645  const bool verbose=false;
2646  std::list<std::pair<RecUnitCell,float> > vSolution2;
2647  // TODO: take into account number of spurious lines for cutoff value.
2648  // keep only solutions above mBestScore/5
2649  for(list<pair<RecUnitCell,float> >::iterator pos=mvSolution.begin();pos!=mvSolution.end();)
2650  {
2651  if(pos->second<(mBestScore/5)) pos=mvSolution.erase(pos);
2652  else ++pos;
2653  }
2654  if(updateReportThreshold&& ((mBestScore/5)>mMinScoreReport))
2655  {
2656  cout<<"CellExplorer::ReduceSolutions(): update threshold for report from "
2657  <<mMinScoreReport<<" to "<<mBestScore/5<<endl;
2658  mMinScoreReport=mBestScore/5;
2659  }
2660 
2661  while(mvSolution.size()>0)
2662  {
2663  vSolution2.push_back(mvSolution.front());
2664  mvSolution.pop_front();
2665  vector<float> par=vSolution2.back().first.DirectUnitCell();
2666  if(verbose)
2667  cout<<__FILE__<<":"<<__LINE__<<" SOLUTION: a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]
2668  <<", alpha="<<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG
2669  <<", V="<<par[6]<<", score="<<vSolution2.back().second<<", SIMILAR TO:"<<endl;
2670  for(list<pair<RecUnitCell,float> >::iterator pos=mvSolution.begin();pos!=mvSolution.end();)
2671  {
2672  if(SimilarRUC(pos->first,vSolution2.back().first))
2673  {
2674  par=pos->first.DirectUnitCell();
2675  if(verbose)
2676  cout<<__FILE__<<":"<<__LINE__<<" 1: a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]
2677  <<", alpha="<<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG
2678  <<", V="<<par[6]<<", score="<<pos->second<<" ("<<mvSolution.size()<<")"<<endl;
2679  if(vSolution2.back().first.mlattice==pos->first.mlattice)
2680  {
2681  if(pos->second>vSolution2.back().second) vSolution2.back()=*pos;
2682  }
2683  else if(vSolution2.back().first.mlattice<pos->first.mlattice) vSolution2.back()=*pos;
2684  pos=mvSolution.erase(pos);
2685  }
2686  else
2687  {
2688  par=pos->first.DirectUnitCell();
2689  if(verbose)
2690  cout<<__FILE__<<":"<<__LINE__<<" 0: a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]
2691  <<", alpha="<<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG
2692  <<", V="<<par[6]<<", score="<<pos->second<<" ("<<mvSolution.size()<<")"<<endl;
2693  ++pos;
2694  }
2695  }
2696  }
2697  mvSolution=vSolution2;
2698  mvSolution.sort(compareRUCScore);
2699 
2700  // keep at most 100 solutions, update mDicVolDepthReport and mMinScoreReport if necessary
2701  if(mvSolution.size()>100)
2702  {
2703  mvSolution.resize(100);
2704  if(updateReportThreshold && (mvSolution.back().second>mMinScoreReport))
2705  {
2706  cout<<"CellExplorer::ReduceSolutions(): update threshold for report from "
2707  <<mMinScoreReport<<" to "<<mvSolution.back().second<<endl;
2708  mMinScoreReport=mvSolution.back().second;
2709  }
2710  }
2711 }
2712 
2713 
2714 void CellExplorer::Init()
2715 {
2716  // Prepare global optimisation
2717  //for(unsigned int i=0;i<mpPeakList->nb;++i)
2718  // cout<<__FILE__<<":"<<__LINE__<<":d*="<<mpPeakList->mvdobs[i]<<", d*^2="<<mpPeakList->mvd2obs[i]<<endl;
2719  srand(time(NULL));
2720  vector<pair<RecUnitCell,float> >::iterator pos;
2721  const float min_latt=1./mLengthMax;
2722  const float max_latt=1./mLengthMin;
2723  const float amp_crossp=abs(cos(mAngleMax));
2724  //mMin[0]=-.002;mAmp[0]=.004;
2725  mMin[0]=.00;mAmp[0]=.00;
2726  switch(mlattice)
2727  {
2728  case TRICLINIC:
2729  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2730  mMin[2]=min_latt;mAmp[2]=max_latt-min_latt;
2731  mMin[3]=min_latt;mAmp[3]=max_latt-min_latt;
2732  mMin[4]=0;mAmp[4]=amp_crossp;
2733  mMin[5]=0;mAmp[5]=amp_crossp;
2734  mMin[6]=0;mAmp[6]=amp_crossp;
2735  mnpar=7;
2736  break;
2737  case MONOCLINIC:
2738  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2739  mMin[2]=min_latt;mAmp[2]=max_latt-min_latt;
2740  mMin[3]=min_latt;mAmp[3]=max_latt-min_latt;
2741  mMin[4]=0;mAmp[4]=amp_crossp;
2742  mnpar=5;
2743  break;
2744  case ORTHORHOMBIC:
2745  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2746  mMin[2]=min_latt;mAmp[2]=max_latt-min_latt;
2747  mMin[3]=min_latt;mAmp[3]=max_latt-min_latt;
2748  mnpar=4;
2749  break;
2750  case HEXAGONAL:
2751  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2752  mMin[2]=min_latt;mAmp[2]=max_latt-min_latt;
2753  mnpar=3;
2754  break;
2755  case RHOMBOEDRAL:
2756  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2757  mMin[2]=-amp_crossp;mAmp[2]=2*amp_crossp;
2758  mnpar=3;
2759  break;
2760  case TETRAGONAL:
2761  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2762  mMin[2]=min_latt;mAmp[2]=max_latt-min_latt;
2763  mnpar=3;
2764  break;
2765  case CUBIC:
2766  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2767  mnpar=2;
2768  break;
2769  }
2770  //for(unsigned int k=0;k<mnpar;++k) cout<<"par["<<k<<"]: "<<mMin[k]<<"->"<<mMin[k]+mAmp[k]<<endl;
2771 
2772  unsigned int nb1=0,nb2=0;
2773  switch(mlattice)
2774  {
2775  case TRICLINIC:
2776  {
2777  nb1=3;
2778  nb2=3;
2779  break;
2780  }
2781  case MONOCLINIC:
2782  {
2783  nb1=3;
2784  nb2=1;
2785  break;
2786  }
2787  case ORTHORHOMBIC:
2788  {
2789  nb1=3;
2790  break;
2791  }
2792  case HEXAGONAL:
2793  {
2794  nb1=2;
2795  break;
2796  }
2797  case RHOMBOEDRAL:
2798  {
2799  nb1=2;
2800  break;
2801  }
2802  case TETRAGONAL:
2803  {
2804  nb1=2;
2805  break;
2806  }
2807  case CUBIC:
2808  {
2809  nb1=1;
2810  break;
2811  }
2812  }
2813  this->ResetParList();
2814  {
2815  RefinablePar tmp("Zero",mRecUnitCell.par+0,-0.01,0.01,gpRefParTypeObjCryst,REFPAR_DERIV_STEP_ABSOLUTE,
2816  true,false,true,false);
2817  tmp.SetDerivStep(1e-4);
2818  this->AddPar(tmp);
2819  }
2820  char buf [50];
2821  string str="Reciprocal unit cell par #";
2822  for(unsigned int i=0;i<nb1;++i)
2823  {
2824  sprintf(buf,"%i",i);
2825  RefinablePar tmp(str+(string)buf,mRecUnitCell.par+i+1,
2826  0.01,1.,gpRefParTypeObjCryst,REFPAR_DERIV_STEP_ABSOLUTE,
2827  false,false,true,false);
2828  tmp.SetDerivStep(1e-4);
2829  this->AddPar(tmp);
2830  }
2831  for(unsigned int i=nb1;i<(nb1+nb2);++i)
2832  {
2833  sprintf(buf,"%i",i);
2834  RefinablePar tmp(str+(string)buf,mRecUnitCell.par+i+1,
2835  0.0,0.5,gpRefParTypeObjCryst,REFPAR_DERIV_STEP_ABSOLUTE,
2836  false,false,true,false);
2837  tmp.SetDerivStep(1e-4);
2838  this->AddPar(tmp);
2839  }
2840 }
2841 
2842 }//namespace
ObjCryst
The namespace which includes all objects (crystallographic and algorithmic) in ObjCryst++.
Definition: doc-main.h:25
ObjCryst::EstimateCellVolume
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:46
ObjCryst::DichoIndexed
bool DichoIndexed(const PeakList &dhkl, const RecUnitCell &par, const RecUnitCell &dpar, const unsigned int nbUnindexed=0, const bool verbose=false, unsigned int useStoredHKL=0, const unsigned int maxNbMissingBelow5=0)
Number of reflexions found in the intervals calculated between par+dpar and par-dpar.
Definition: Indexing.cpp:1575
ObjCryst::Score
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:948
ObjCryst::CrystalSystem
CrystalSystem
Different lattice types.
Definition: Indexing.h:40
ObjCryst::RecUnitCell
Lightweight class describing the reciprocal unit cell, for the fast computation of d*_hkl^2.
Definition: Indexing.h:62
ObjCryst::RecUnitCell::hkl2d_delta
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:456
ObjCryst::RecUnitCell::mNbSpurious
unsigned int mNbSpurious
The number of spurious lines used to match this RecUnitCell.
Definition: Indexing.h:111
ObjCryst::RecUnitCell::par
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:107
ObjCryst::RecUnitCell::hkl2d
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:108
ObjCryst::RecUnitCell::RecUnitCell
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, const unsigned int nbspurious=0)
light-weight class storing the reciprocal space unitcell
Definition: Indexing.cpp:81
ObjCryst::RecUnitCell::DirectUnitCell
std::vector< float > DirectUnitCell(const bool equiv=false) const
Compute real space unit cell from reciprocal one.
Definition: Indexing.cpp:564
ObjCryst::PeakList
Class to store positions of observed reflections.
Definition: Indexing.h:119
ObjCryst::PeakList::Simulate
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, const bool merge=false)
Generate a list of simulated peak positions, from given lattice parameters.
Definition: Indexing.cpp:807
ObjCryst::PeakList::GetPeakList
vector< hkl > & GetPeakList()
Get peak list.
Definition: Indexing.cpp:942
ObjCryst::PeakList::mvPredictedHKL
list< hkl > mvPredictedHKL
Full list of calculated HKL positions for a given solution, up to a given resolution After finding a ...
Definition: Indexing.h:211
ObjCryst::PeakList::mvHKL
vector< hkl > mvHKL
Predict peak positions Best h,k,l for each observed peak (for least-squares refinement) This is store...
Definition: Indexing.h:208
ObjCryst::PeakList::AddPeak
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:898
ObjCryst::PeakList::hkl0
One set of Miller indices, a possible indexation for a reflection.
Definition: Indexing.h:160
ObjCryst::PeakList::hkl
One observed diffraction line, to be indexed.
Definition: Indexing.h:169
ObjCryst::RefinablePar
Generic class for parameters of refinable objects.
Definition: RefinableObj.h:225
ObjCryst::RefinablePar::GetName
string GetName() const
Get the parameter's name.
Definition: RefinableObj.cpp:545
Chronometer
Simple chronometer class, with microsecond precision.
Definition: Chronometer.h:35