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==ORTHOROMBIC){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 ORTHOROMBIC:
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 ORTHOROMBIC:
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 ORTHOROMBIC:
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 ORTHOROMBIC:
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 ORTHOROMBIC:
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 ORTHOROMBIC:
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 ORTHOROMBIC:
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 ORTHOROMBIC: //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 ORTHOROMBIC:
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 {
812  if(deg){alpha*=DEG2RAD;beta*=DEG2RAD;gamma*=DEG2RAD;}
813  const float v=sqrt(1-cos(alpha)*cos(alpha)-cos(beta)*cos(beta)-cos(gamma)*cos(gamma)
814  +2*cos(alpha)*cos(beta)*cos(gamma));
815 
816  const float aa=sin(alpha)/a/v;
817  const float bb=sin(beta )/b/v;
818  const float cc=sin(gamma)/c/v;
819 
820  const float alphaa=acos( (cos(beta )*cos(gamma)-cos(alpha))/sin(beta )/sin(gamma) );
821  const float betaa =acos( (cos(alpha)*cos(gamma)-cos(beta ))/sin(alpha)/sin(gamma) );
822  const float gammaa=acos( (cos(alpha)*cos(beta )-cos(gamma))/sin(alpha)/sin(beta ) );
823 
824  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);
825  std::list<float> vd2;
826 
827  for(int h=0;h<=20;h++)
828  for(int k=-20;k<=20;k++)
829  {
830  if((h==0) && (k<0)) k=0;
831  for(int l=-20;l<=20;l++)
832  {
833  if((h==0) && (k==0) && (l<=0)) l=1;
834  vd2.push_back(sqrt(ruc.hkl2d(h,k,l)));
835  }
836  }
837  //
838  std::list<float>::iterator pos=vd2.begin();
839  if(percentMissing>0.90) percentMissing=0.90;
840  for(;pos!=vd2.end();++pos)
841  {
842  if((rand()/float(RAND_MAX))<percentMissing) *pos=1e10;
843  }
844  vd2.sort();
845  pos=vd2.begin();
846  const float dmin=*pos/2;
847  for(unsigned int i=0;i<nb;++i)pos++;
848  const float dmax=*pos;
849 
850  for(unsigned int i=0;i<nbspurious;++i)
851  {
852  const unsigned int idx=1+i*nb/nbspurious+(rand()%nbspurious);
853  pos=vd2.begin();
854  for(unsigned int j=0;j<idx;++j) pos++;
855  *pos=dmin+rand()/float(RAND_MAX)*(dmax-dmin);
856  }
857 
858  pos=vd2.begin();
859  for(unsigned int i=0;i<nb;++i)
860  {
861  float d=*pos++;
862  const float ds=d*sigma;
863  float d1=d+ds*(rand()/float(RAND_MAX)*2-1);
864  //cout<<d<<" "<<ds<<" "<<d1<<" "<<sigma<<endl;
865  mvHKL.push_back(hkl(d1,1.0,ds));
866  }
867 
868  if(verbose)
869  {
870  char buf[200];
871  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);
872  cout<<"PeakList::Simulate():"<<buf<<endl;
873  }
874  sort(mvHKL.begin(),mvHKL.end(),compareHKL_d);
875  return v*a*b*c;
876 }
877 
878 void PeakList::ExportDhklDSigmaIntensity(std::ostream &os)const
879 {
880  for(vector<PeakList::hkl>::const_iterator pos=mvHKL.begin();pos!=mvHKL.end();++pos)
881  {
882  const float sigma=1/(pos->dobs-pos->dobssigma/2)-1/(pos->dobs+pos->dobssigma/2);
883  os<< std::fixed << setw(6) << setprecision(3) << 1/pos->dobs <<" "<< sigma <<" "<< std::scientific << pos->iobs <<endl;
884  }
885 }
886 
887 void PeakList::AddPeak(const float d, const float iobs,const float dobssigma,const float iobssigma,
888  const int h,const int k, const int l,const float d2calc)
889 {
890  if(dobssigma<=0)
891  {// Manually added peak ? Use other reflection's sigmas to evaluate sigma for this reflection
892  float s=0;
893  for(vector<hkl>::const_iterator pos=mvHKL.begin();pos!=mvHKL.end();++pos)
894  s+= pos->dobssigma;
895  s/=mvHKL.size();
896  if(s>0) mvHKL.push_back(hkl(d,iobs,s,iobssigma,h,k,l,d2calc));
897  else mvHKL.push_back(hkl(d,iobs,1e-4,iobssigma,h,k,l,d2calc));
898  }
899  else mvHKL.push_back(hkl(d,iobs,dobssigma,iobssigma,h,k,l,d2calc));
900  sort(mvHKL.begin(),mvHKL.end(),compareHKL_d);
901  //this->Print(cout);
902 }
903 
904 void PeakList::RemovePeak(unsigned int idx)
905 {
906  for(unsigned int i=idx;i<(mvHKL.size()-1);++i) mvHKL[i]=mvHKL[i+1];
907  mvHKL.pop_back();
908 }
909 
910 void PeakList::Print(std::ostream &os) const
911 {
912  unsigned int i=0;
913  char buf[200];
914  os<<"PeakList, with "<<mvHKL.size()<<" peaks"<<endl;
915  for(vector<PeakList::hkl>::const_iterator pos=mvHKL.begin();pos!=mvHKL.end();++pos)
916  {
917  const float sigma=1/(pos->dobs-pos->dobssigma/2)-1/(pos->dobs+pos->dobssigma/2);
918  if(pos->isIndexed)
919  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",
920  i++,1/pos->dobs,sigma,
921  1/sqrt(abs(pos->d2calc)),1/sqrt(abs(pos->d2calc))-1/pos->dobs,
922  pos->iobs,pos->h,pos->k,pos->l,pos->isSpurious,pos->stats);
923  else
924  sprintf(buf,"#%3d d=%6.3f+/-%6.3f iobs=%6.3f UNINDEXED Spurious=%1d stats=%6lu",
925  i++,1/pos->dobs,1/(pos->dobs-pos->dobssigma/2)-1/(pos->dobs+pos->dobssigma/2),
926  pos->iobs,pos->isSpurious,pos->stats);
927  os<<buf<<endl;
928  }
929 }
930 
931 vector<PeakList::hkl> & PeakList::GetPeakList(){return mvHKL;}
932 
933 const vector<PeakList::hkl> & PeakList::GetPeakList()const {return mvHKL;}
934 
936 
937 float Score(const PeakList &dhkl, const RecUnitCell &rpar, const unsigned int nbSpurious,
938  const bool verbose,const bool storehkl,const bool storePredictedHKL)
939 {
940  const bool autozero=false;
941  vector<PeakList::hkl>::const_iterator pos,first,last;
942  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
943  {
944  if(storehkl) pos->isIndexed=false;
945  pos->d2calc=0;
946  pos->d2diff=1000;
947  }
948  const unsigned long nb=dhkl.GetPeakList().size();
949  if(storePredictedHKL) dhkl.mvPredictedHKL.clear();
950 
951  unsigned long nbCalc=0;
952  int h,k,l;
953  float predict_coeff=1;
954  if(storePredictedHKL)predict_coeff=2;
955  const float dmax=dhkl.mvHKL[nb-1].d2obs*predict_coeff*1.05;
956  int sk0,sl0;// do we need >0 *and* <0 indices for k,l ?
957  switch(rpar.mlattice)
958  {
959  case TRICLINIC:
960  sk0=-1;sl0=-1;
961  break;
962  case MONOCLINIC:
963  sk0=1;sl0=-1;
964  break;
965  case ORTHOROMBIC:
966  sk0=1;sl0=1;
967  break;
968  case HEXAGONAL:
969  sk0=-1;sl0=1;
970  break;
971  case RHOMBOEDRAL:
972  sk0=-1;sl0=-1;
973  break;
974  case TETRAGONAL:
975  sk0=1;sl0=1;
976  break;
977  case CUBIC:
978  sk0=1;sl0=1;
979  break;
980  // This should never happen. Avoid using unitialized values.
981  default:
982  throw 0;
983  }
984  int stepk,stepl;// steps in k,l to use for centered lattices
985  switch(rpar.mCentering)
986  {
987  case LATTICE_P:stepk=1;stepl=1;break;
988  case LATTICE_I:stepk=1;stepl=2;break;
989  case LATTICE_A:stepk=1;stepl=2;break;
990  case LATTICE_B:stepk=1;stepl=2;break;
991  case LATTICE_C:stepk=2;stepl=1;break;
992  case LATTICE_F:stepk=2;stepl=2;break;
993  // This should never happen. Avoid using unitialized values.
994  default: throw 0;
995  }
996  first=dhkl.GetPeakList().begin();last=dhkl.GetPeakList().end();
997  unsigned long nbCalcH,nbCalcK;// Number of calculated lines below dmax for each h,k
998  for(h=0;;++h)
999  {
1000  nbCalcH=0;
1001  for(int sk=sk0;sk<=1;sk+=2)
1002  {
1003  if(h==0) sk=1;// no need to explore 0kl with both sk -1 and 1
1004  if(stepk==2) k=(h%2);// For LATTICE_C,LATTICE_F: h odd => k odd
1005  else k=0;
1006  for(;;k+=stepk)
1007  {
1008  nbCalcK=0;
1009  for(int sl=sl0;sl<=1;sl+=2)
1010  {
1011  if((h+k)==0)
1012  {
1013  sl=1;// No need to list 0 0 l with l<0
1014  l=1;
1015  }
1016  else
1017  {
1018  if(h==0)
1019  {
1020  if(rpar.mlattice==MONOCLINIC) sl=1;// 0 k l and 0 k -l are equivalent
1021  if((sk<0)||(sl<0)) l=1;// Do not list 0 k 0 with k<0
1022  else l=0;// h==k==0 already covered
1023  }
1024  else
1025  {
1026  if(sl<0) l=1;// Do not list h k 0 twice
1027  else l=0;
1028  }
1029  }
1030  if(stepl==2)
1031  {
1032  if(rpar.mCentering==LATTICE_I) l+=(h+k+l)%2;
1033  if(rpar.mCentering==LATTICE_A) l+=(k+l)%2;// Start at hk1 if k odd
1034  if( (rpar.mCentering==LATTICE_B)
1035  ||(rpar.mCentering==LATTICE_F)) l+=(h+l)%2;// Start at hk1 if h odd
1036  }
1037  for(;;l+=stepl)
1038  {
1039  const float d2=rpar.hkl2d(h,sk*k,sl*l);
1040  if(d2>dmax)
1041  {
1042  //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;
1043  // Only break if d is increasing with l
1044  if((sl*rpar.hkl2d(h,sk*k,sl*l,NULL,3))>=0) break;
1045  else continue;
1046  }
1047  nbCalc++;nbCalcK++;nbCalcH++;
1048  if(storePredictedHKL)
1049  {
1050  dhkl.mvPredictedHKL.push_back(PeakList::hkl(0,0,0,0,h,sk*k,sl*l,d2));
1051  //continue;
1052  }
1053  for(pos=first;pos!=last;++pos)
1054  {
1055  const float tmp=d2-pos->d2obs;
1056  if(tmp<.1)
1057  {
1058  if(tmp<-.1) break;
1059  if(fabs(tmp)<fabs(pos->d2diff))
1060  {
1061  pos->d2diff=tmp;
1062  if(storehkl)
1063  {
1064  pos->h=h;
1065  pos->k=sk*k;
1066  pos->l=sl*l;
1067  pos->isIndexed=true;
1068  pos->d2calc=d2;
1069  }
1070  }
1071  /*
1072  if((verbose)&&(fabs(tmp)<.01))
1073  cout<<__FILE__<<":"<<__LINE__<<" hkl: "<<h<<" "<<k<<" "<<l
1074  <<"#"<<i<<": calc="<<sqrt(d2)<<", obs="<<sqrt(*pd2)<<", min_epsilon="<<*pdq2<<", tmp="<<tmp<<endl;
1075  */
1076  }
1077  }
1078  }
1079  }
1080  if(nbCalcK==0) //d(hk0)>dmax
1081  {
1082  //cout<<__FILE__<<":"<<__LINE__<<" hkl: "<<h<<" "<<sk*k<<" "<<0<<" deriv="<<sk*rpar.hkl2d(h,sk*k,0,NULL,2)<<endl;
1083  if((sk*rpar.hkl2d(h,sk*k,0,NULL,2))>=0) break;
1084  }
1085  }
1086  }
1087  if(nbCalcH==0) break;//h00 beyond limit
1088  }
1089  float epsilon=0.0,zero=0.0;
1090  if(autozero)
1091  {
1092  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos) zero+=pos->d2diff;
1093  zero/=nb;
1094  }
1095  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
1096  {
1097  epsilon +=fabs(pos->d2diff-zero);
1098  }
1099  if(nbSpurious>0)
1100  {// find worst fitting lines and remove them from epsilon calculation
1101  list<pair<float,unsigned int> > vdiff_idx;
1102  unsigned int i=0;
1103  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
1104  vdiff_idx.push_back(make_pair(fabs(pos->d2diff),i++));
1105  vdiff_idx.sort(comparePairFirst<float,unsigned int>);
1106  i=0;
1107  for(list<pair<float,unsigned int> >::reverse_iterator rpos=vdiff_idx.rbegin();rpos!=vdiff_idx.rend();++rpos)
1108  {// :TODO: correct zero after removing spurious lines
1109  epsilon -= fabs(rpos->first-zero);
1110  if(storehkl) dhkl.GetPeakList()[rpos->second].isIndexed=false;
1111  dhkl.GetPeakList()[rpos->second].stats++;
1112  if(++i==nbSpurious) break;
1113  }
1114  }
1115  if(verbose)
1116  {
1117  float epstmp=0;
1118  //unsigned long i=0;
1119  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
1120  {
1121  epstmp+=fabs(pos->d2diff-zero);
1122  //cout<<"Line #"<<i++<<": obs="<<pos->d2obs<<", diff="<<pos->d2diff<<" -> epsilon="<<epstmp<<endl;
1123  }
1124  cout<<"epsilon="<<epstmp<<", dmax="<<dmax<<" ,nb="<<nb<<" ,nbcalc="<<nbCalc<<endl;
1125  }
1126  /*
1127  else
1128  {//Only stat+ the worst
1129  float max=-1;
1130  unsigned int worst=0;
1131  unsigned int i=0;
1132  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
1133  if(abs(pos->d2diff)>max) {worst=i++;max=abs(pos->d2diff);}
1134  else i++;
1135  dhkl.GetPeakList()[worst].stats++;
1136  }
1137  */
1138  if(nbCalc==0) return 0;
1139  const float score=(dmax/predict_coeff)*nb/(2*epsilon*nbCalc);
1140  if(verbose)
1141  {
1142  dhkl.Print(cout);
1143  cout<<"Final score:"<<score<<", nbCalc="<<nbCalc<<" ,<epsilon>="<<epsilon<<" nb="<<nb<<" Qn="<<sqrt(dmax)<<endl;
1144  }
1145  return score;
1146 }
1147 
1149 
1150 CellExplorer::CellExplorer(const PeakList &dhkl, const CrystalSystem lattice, const unsigned int nbSpurious):
1151 mnpar(3),mpPeakList(&dhkl),
1152 mLengthMin(4),mLengthMax(25),
1153 mAngleMin(M_PI),mAngleMax(2*M_PI/3),
1154 mVolumeMin(0),mVolumeMax(1600),
1155 mZeroShiftMin(0),mZeroShiftMax(0),
1156 mlattice(lattice),mCentering(LATTICE_P),mNbSpurious(nbSpurious),
1157 mObs(0),mCalc(0),mWeight(0),mDeriv(0),mBestScore(0.0),
1158 mMinScoreReport(10),mMaxDicVolDepth(6),mDicVolDepthReport(6),
1159 mNbLSQExcept(0)
1160 {
1161  this->Init();
1162 }
1163 
1164 void CellExplorer::Evolution(unsigned int ng,const bool randomize,const float f,const float cr,unsigned int np)
1165 {
1166  this->Init();
1167  const bool autozero=true;
1168  //cout<<__FILE__<<":"<<__LINE__<<"<CellExplorer::Evolution(...): randomizing,ng="<<ng
1169  // <<"random="<<randomize<<"f="<<f<<"cr="<<cr<<"np="<<np<<endl;
1170  vector<pair<RecUnitCell,float> > vRUC(np);
1171  vector<pair<RecUnitCell,float> > vTrial(np);
1172  float bestScore=-1e20;
1173  vector<pair<RecUnitCell,float> >::iterator bestpos=vRUC.begin();
1174 
1175  const clock_t mTime0=clock();
1176 
1177  if(randomize)
1178  {
1179  for(unsigned int i=0;i<vRUC.size();++i)
1180  {
1181  vRUC[i].first.mlattice=mlattice;
1182  vTrial[i].first.mlattice=mlattice;
1183  for(unsigned int k=0;k<mnpar;++k) vRUC[i].first.par[k]=mMin[k]+mAmp[k]*rand()/(float)RAND_MAX;
1184  vRUC[i].second=Score(*mpPeakList,vRUC[i].first,mNbSpurious);
1185  }
1186  }
1187 
1188  for(unsigned long i=ng;i>0;--i)
1189  {
1190  for(unsigned j=0;j<np;j++)
1191  {
1192  if(true)
1193  {// DE/rand/1/exp
1194  unsigned int r1=j,r2=j,r3=j;
1195  while(r1==j)r1=rand()%np;
1196  while((r2==j)||(r1==r2))r2=rand()%np;
1197  while((r3==j)||(r3==r1)||(r3==r2))r3=rand()%np;
1198  unsigned int ncr=1+(int)(cr*mnpar*rand()/(float)RAND_MAX);
1199  unsigned int ncr0=rand()%mnpar;
1200  RecUnitCell *t0=&(vTrial[j].first);
1201  RecUnitCell *c0=&(vRUC[j].first);
1202  RecUnitCell *c1=&(vRUC[r1].first);
1203  RecUnitCell *c2=&(vRUC[r2].first);
1204  RecUnitCell *c3=&(vRUC[r3].first);
1205  for(unsigned int k=0;k<mnpar;++k)t0->par[k] = c0->par[k];
1206  for(unsigned int k=0;k<ncr;++k)
1207  {
1208  const unsigned l=(ncr0+k)%mnpar;
1209  const float v1=c1->par[l]-mMin[l];
1210  const float v2=c2->par[l]-mMin[l];
1211  const float v3=c3->par[l]-mMin[l];
1212  t0->par[l]=mMin[l]+fmod(v1+f*(v2-v3)+3*mAmp[l],mAmp[l]);
1213  }
1214  }
1215  if(false)
1216  {// DE/rand-to-best/1/exp
1217  unsigned int r1=j,r2=j,r3=j;
1218  while(r1==j)r1=rand()%np;
1219  while((r2==j)||(r1==r2))r2=rand()%np;
1220  while((r3==j)||(r3==r1)||(r3==r2))r3=rand()%np;
1221  unsigned int ncr=1+(int)(cr*(mnpar-1)*rand()/(float)RAND_MAX);
1222  unsigned int ncr0=rand()%mnpar;
1223  RecUnitCell *t0=&(vTrial[j].first);
1224  RecUnitCell *c0=&(vRUC[j].first);
1225  //RecUnitCell *c1=&(vRUC[r1].first);
1226  RecUnitCell *c2=&(vRUC[r2].first);
1227  RecUnitCell *c3=&(vRUC[r3].first);
1228  RecUnitCell *best=&(bestpos->first);
1229  for(unsigned int k=0;k<6;++k)t0->par[k] = c0->par[k];//mMin[k]+mAmp[k]*rand()/(float)RAND_MAX;
1230  for(unsigned int k=0;k<ncr;++k)
1231  {
1232  const unsigned l=(ncr0+k)%mnpar;
1233  const float v0=c0->par[l]-mMin[l];
1234  //const float v1=c1->par[l]-mMin[l];
1235  const float v2=c2->par[l]-mMin[l];
1236  const float v3=c3->par[l]-mMin[l];
1237  const float vb=best->par[l]-mMin[l];
1238  t0->par[l]=mMin[l]+fmod(vb+f*(vb-v0)+f*(v2-v3)+5*mAmp[l],mAmp[l]);
1239  }
1240  }
1241  if(false)
1242  {// MC
1243  const float amp=.05/(1+i*.01);
1244  RecUnitCell *t0=&(vTrial[j].first);
1245  for(unsigned int k=0;k<6;++k)
1246  {
1247 
1248  t0->par[k] = mMin[k]+ fmod((float)(amp*mAmp[k]*(rand()/(float)RAND_MAX-0.5)+5*mAmp[k]),(float)mAmp[k]);
1249  }
1250  }
1251  }
1252  // Compute cost for all trials and select best
1253  vector<pair<RecUnitCell,float> >::iterator posTrial,pos;
1254  posTrial=vTrial.begin();
1255  pos=vRUC.begin();
1256  for(;posTrial!=vTrial.end();)
1257  {
1258  // If using auto-zero, fix zero parameter
1259  if(autozero) posTrial->first.par[0]=0;
1260  // Did we go beyond allowed volume ?
1261  switch(mlattice)
1262  {
1263  case TRICLINIC:
1264  break;
1265  case MONOCLINIC:
1266  {
1267  float v0=posTrial->first.par[1]*posTrial->first.par[2]*posTrial->first.par[3];
1268  while(v0<1/mVolumeMax)
1269  {
1270  const unsigned int i=rand()%3+1;
1271  posTrial->first.par[i]*=1/(mVolumeMax*v0)+1e-4;
1272  if(posTrial->first.par[i]>(mMin[i]+mAmp[i])) posTrial->first.par[i]=mMin[i]+mAmp[i];
1273  v0=posTrial->first.par[1]*posTrial->first.par[2]*posTrial->first.par[3];
1274  }
1275  break;
1276  }
1277  case ORTHOROMBIC:
1278  break;
1279  case HEXAGONAL:
1280  break;
1281  case RHOMBOEDRAL:
1282  break;
1283  case TETRAGONAL:
1284  break;
1285  case CUBIC:
1286  break;
1287  }
1288 
1289  const float score=Score(*mpPeakList,posTrial->first,mNbSpurious);
1290  if(score > pos->second)
1291  {
1292  pos->second=score;
1293  const REAL *p0=posTrial->first.par;
1294  REAL *p1=pos->first.par;
1295  for(unsigned int k=0;k<mnpar;++k) *p1++ = *p0++;
1296  if(score>bestScore)
1297  {
1298  bestScore=score;
1299  bestpos=pos;
1300  }
1301  }
1302  /*
1303  else
1304  {
1305  if(log(rand()/(float)RAND_MAX)>(-(score-pos->second)))
1306  {
1307  pos->second=score;
1308  const float *p0=posTrial->first.par;
1309  float *p1=pos->first.par;
1310  for(unsigned int k=0;k<mnpar;++k) *p1++ = *p0++;
1311  }
1312  }
1313  */
1314  ++pos;++posTrial;
1315  }
1316  if((i%100000)==0)
1317  {
1318  vector<float> par=bestpos->first.DirectUnitCell();
1319  cout<<"Generation #"<<ng-i<<", Best score="<<bestScore
1320  <<" Trial: a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]<<", alpha="
1321  <<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG<<", V="<<par[6]
1322  <<" "<<(ng-i)*np/((clock()-mTime0)/(float)CLOCKS_PER_SEC)<<" trials/s"
1323  <<endl;
1324  }
1325  if(false)//((i%10000)==0)
1326  {// Randomize periodically
1327  for(vector<pair<RecUnitCell,float> >::iterator pos=vRUC.begin();pos!=vRUC.end();++pos)
1328  {
1329  if(pos==bestpos) continue;
1330  for(unsigned int k=0;k<mnpar;++k) pos->first.par[k]=mMin[k]+mAmp[k]*rand()/(float)RAND_MAX;
1331  }
1332  }
1333  }
1334  /*
1335  for(vector<pair<RecUnitCell,float> >::iterator pos=vRUC.begin();pos!=vRUC.end();++pos)
1336  {
1337  // Final cost
1338  vector<float> par=pos->first.DirectUnitCell();
1339  cout<<__FILE__<<":"<<__LINE__<<" Trial: a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]<<", alpha="
1340  <<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG<<", V="<<par[6]
1341  <<", score="<<pos->second<<endl;
1342  }
1343  Score(*mpPeakList,bestpos->first,mNbSpurious,true);
1344  */
1345 
1346  //this->ReduceSolutions(true);
1347 
1348  mRecUnitCell=bestpos->first;
1349  float score=Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true);
1350  vector<float> par=mRecUnitCell.DirectUnitCell();
1351  cout<<__FILE__<<":"<<__LINE__<<" Best-DE : a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]<<", alpha="
1352  <<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG<<", V="<<par[6]
1353  <<", score="<<bestpos->second
1354  <<" ("<<ng*np/((clock()-mTime0)/(float)CLOCKS_PER_SEC)<<" trials/s)"<<endl;
1355  if(score>mMinScoreReport*.5)
1356  {
1357  // Now, do a least-squares refinement on best
1358  mRecUnitCell=bestpos->first;
1359  this->LSQRefine(10,true,true);
1360  par=mRecUnitCell.DirectUnitCell();
1361  score=Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true);
1362  cout<<__FILE__<<":"<<__LINE__<<" Best-LSQ: 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="<<score<<endl;
1365  if((score>mMinScoreReport)&&(score>(mBestScore/3)))
1366  {
1367  if(score>mBestScore) mBestScore=score;
1368  mvSolution.push_back(make_pair(mRecUnitCell,score));
1369  mvSolution.back().first.mNbSpurious = mNbSpurious;
1370  this->ReduceSolutions(true);// We may have solutions from previous runs
1371  }
1372  }
1373 }
1374 
1375 void CellExplorer::SetLengthMinMax(const float min,const float max)
1376 {
1377  mLengthMin=min;
1378  mLengthMax=max;
1379 }
1380 void CellExplorer::SetAngleMinMax(const float min,const float max)
1381 {
1382  mAngleMin=min;
1383  mAngleMax=max;
1384 }
1385 void CellExplorer::SetVolumeMinMax(const float min,const float max)
1386 {
1387  mVolumeMin=min;
1388  mVolumeMax=max;
1389 }
1390 void CellExplorer::SetNbSpurious(const unsigned int nb)
1391 {
1392  mNbSpurious=nb;
1393 }
1394 void CellExplorer::SetMinMaxZeroShift(const float min,const float max)
1395 {
1396  mZeroShiftMin=min;
1397  mZeroShiftMax=max;
1398 }
1399 
1400 void CellExplorer::SetCrystalSystem(const CrystalSystem system)
1401 {
1402  mlattice=system;
1403 }
1404 
1405 void CellExplorer::SetCrystalCentering(const CrystalCentering cent)
1406 {
1407  mCentering=cent;
1408 }
1409 
1410 void CellExplorer::SetD2Error(const float err){mD2Error=err;}
1411 
1412 const string& CellExplorer::GetClassName() const
1413 {
1414  const static string className="CellExplorer";
1415  return className;
1416 }
1417 const string& CellExplorer::GetName() const
1418 {
1419  const static string name="Some CellExplorer Object";
1420  return name;
1421 }
1422 void CellExplorer::Print() const
1423 {
1424  this->RefinableObj::Print();
1425 }
1426 unsigned int CellExplorer::GetNbLSQFunction() const
1427 {return 1;}
1428 
1429 const CrystVector_REAL& CellExplorer::GetLSQCalc(const unsigned int) const
1430 {
1431  VFN_DEBUG_ENTRY("CellExplorer::GetLSQCalc()",2)
1432  unsigned int j=0;
1433  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();++pos)
1434  {
1435  if(pos->isIndexed)
1436  mCalc(j++)=mRecUnitCell.hkl2d(pos->h,pos->k,pos->l);
1437  }
1438  //cout<<__FILE__<<":"<<__LINE__<<"LSQCalc : Score:"<<Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true,false)<<endl;
1439  VFN_DEBUG_EXIT("CellExplorer::GetLSQCalc()",2)
1440  return mCalc;
1441 }
1442 const CrystVector_REAL& CellExplorer::GetLSQObs(const unsigned int) const
1443 {
1444  VFN_DEBUG_MESSAGE("CellExplorer::GetLSQObs()",2)
1445  return mObs;
1446 }
1447 const CrystVector_REAL& CellExplorer::GetLSQWeight(const unsigned int) const
1448 {
1449  VFN_DEBUG_MESSAGE("CellExplorer::GetLSQWeight()",2)
1450  //:TODO: exclude the worst points (user-chosen number)
1451  return mWeight;
1452 }
1453 const CrystVector_REAL& CellExplorer::GetLSQDeriv(const unsigned int, RefinablePar &refpar)
1454 {
1455  VFN_DEBUG_ENTRY("CellExplorer::GetLSQDeriv()",2)
1456  REAL *par=NULL;
1457  if(refpar.GetName()=="Reciprocal unit cell par #0") par=mRecUnitCell.par+1;
1458  else
1459  if(refpar.GetName()=="Reciprocal unit cell par #1") par=mRecUnitCell.par+2;
1460  else
1461  if(refpar.GetName()=="Reciprocal unit cell par #2") par=mRecUnitCell.par+3;
1462  else
1463  if(refpar.GetName()=="Reciprocal unit cell par #3") par=mRecUnitCell.par+4;
1464  else
1465  if(refpar.GetName()=="Reciprocal unit cell par #4") par=mRecUnitCell.par+5;
1466  else
1467  if(refpar.GetName()=="Reciprocal unit cell par #5") par=mRecUnitCell.par+6;
1468  else
1469  if(refpar.GetName()=="Zero") par=mRecUnitCell.par+0;
1470  else cout<<__FILE__<<":"<<__LINE__<<":Parameter not found:"<<refpar.GetName()<<endl;
1471  unsigned int j=0;
1472  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();++pos)
1473  {
1474  VFN_DEBUG_MESSAGE("CellExplorer::GetLSQDeriv():"<<j<<"/"<<mpPeakList->GetPeakList().size(),2)
1475  VFN_DEBUG_MESSAGE("CellExplorer::GetLSQDeriv():"<<pos->h<<","<<pos->k<<","<<pos->l,2)
1476  if(pos->isIndexed)
1477  mDeriv(j++)=mRecUnitCell.hkl2d(pos->h,pos->k,pos->l,par);
1478  }
1479  VFN_DEBUG_EXIT("CellExplorer::GetLSQDeriv()",2)
1480  return mDeriv;
1481 }
1482 
1483 void CellExplorer::BeginOptimization(const bool allowApproximations, const bool enableRestraints)
1484 {
1485  VFN_DEBUG_ENTRY("CellExplorer::BeginOptimization()",5)
1486  Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true,false);
1487  const unsigned int nb=mpPeakList->GetPeakList().size();
1488  mCalc.resize(nb-mNbSpurious);
1489  mObs.resize(nb-mNbSpurious);
1490  mWeight.resize(nb-mNbSpurious);
1491  mDeriv.resize(nb-mNbSpurious);
1492  int j=0;
1493  float thres=0.0;
1494  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();++pos)
1495  if(thres<pos->iobs) thres=pos->iobs;
1496  thres/=10;// weight=1 for intensities up to Imax/10
1497 
1498  //cout <<"Beginning optimization with reflexions:"<<endl;
1499  //char buf[100];
1500  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();++pos)
1501  {
1502  if(pos->isIndexed)
1503  {
1504  mObs(j)=pos->d2obs;
1505  if(mObs(j)>thres) mWeight(j)=1;
1506  else mWeight(j)=mObs(j)/thres;
1507  /*
1508  sprintf(buf,"#%2d (%3d %3d %3d) dobs=%6.3f dcalc=%6.3f iobs=%6.3f weight=%6.4f",
1509  i,mpPeakList->mvHKL[i].h,mpPeakList->mvHKL[i].k,mpPeakList->mvHKL[i].l,
1510  1/mpPeakList->mvdobs[i],1/sqrt(mRecUnitCell.hkl2d(mpPeakList->mvHKL[i].h,mpPeakList->mvHKL[i].k,mpPeakList->mvHKL[i].l)),
1511  mpPeakList->mviobs[i],mWeight(j));
1512  */
1513  j++;
1514  }
1515  /*
1516  else
1517  {
1518  sprintf(buf,"#%2d (%3d %3d %3d) dobs=%6.3f dcalc=%6.3f iobs=%6.3f SPURIOUS",
1519  i,mpPeakList->mvHKL[i].h,mpPeakList->mvHKL[i].k,mpPeakList->mvHKL[i].l,
1520  1/mpPeakList->mvdobs[i],1/sqrt(mRecUnitCell.hkl2d(mpPeakList->mvHKL[i].h,mpPeakList->mvHKL[i].k,mpPeakList->mvHKL[i].l)),
1521  mpPeakList->mviobs[i]);
1522  }
1523  cout<<buf<<endl;
1524  */
1525  }
1526  this->RefinableObj::BeginOptimization(allowApproximations,enableRestraints);
1527  VFN_DEBUG_EXIT("CellExplorer::BeginOptimization()",5)
1528 }
1529 
1530 void CellExplorer::LSQRefine(int nbCycle, bool useLevenbergMarquardt, const bool silent)
1531 {
1532  if(mNbLSQExcept>100)
1533  {
1534  if(!silent) cout<<"CellExplorer::LSQRefine(): LSQ was disabled, too many (>100) exceptions caught. Weird unit cell parameters ?";
1535  return;
1536  }
1537  VFN_DEBUG_ENTRY("CellExplorer::LSQRefine()",5)
1538  mLSQObj.SetRefinedObj(*this);
1539  mLSQObj.PrepareRefParList(true);
1540  //this->BeginOptimization();
1541  //cout<<FormatVertVector<REAL>(this->GetLSQObs(0),this->GetLSQCalc(0),this->GetLSQWeight(0),this->GetLSQDeriv(0,this->GetPar((long)0)))<<endl;
1542  try {mLSQObj.Refine(nbCycle,useLevenbergMarquardt,silent);}
1543  catch(const ObjCrystException &except)
1544  {
1545  if(++mNbLSQExcept>100) cout<<"WARNING CellExplorer::LSQRefine(): LSQ was disabled, too many (>100) exceptions caught. Weird unit cell parameters ?"<<endl ;
1546  }
1547  if(!silent) mpPeakList->Print(cout);
1548  VFN_DEBUG_EXIT("CellExplorer::LSQRefine()",5)
1549 }
1550 
1564 bool DichoIndexed(const PeakList &dhkl, const RecUnitCell &par,const RecUnitCell &dpar,
1565  const unsigned int nbUnindexed=0,const bool verbose=false,unsigned int useStoredHKL=0,
1566  const unsigned int maxNbMissingBelow5=0)
1567 {
1568  const unsigned int nb=dhkl.GetPeakList().size();
1569  int nbIndexed=nb-nbUnindexed;// Number of reflections we require to be indexed
1570  float d5=0;
1571  if(maxNbMissingBelow5>0) d5=dhkl.GetPeakList()[4].d2obs;
1572  // number of missing reflections calculated below 5th observed line
1573  unsigned int nbMissingBelow5=0;
1574  // List of indexed reflections
1575  vector<PeakList::hkl>::const_iterator pos,first,last,end;
1576  if(useStoredHKL==1)
1577  {// We already now possible Miller indices for all reflections
1578  unsigned int nbUnIx = 0;
1579  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
1580  {
1581  pos->isIndexed=false;
1582  for(list<PeakList::hkl0>::const_iterator phkl0=pos->vDicVolHKL.begin();phkl0!=pos->vDicVolHKL.end();++phkl0)
1583  {
1584  float d0,d1;
1585  par.hkl2d_delta(phkl0->h,phkl0->k,phkl0->l,dpar,d0,d1);
1586  if((pos->d2obsmax>=d0) && (d1>=pos->d2obsmin))
1587  {
1588  pos->d2calc=(d0+d1)/2;
1589  pos->isIndexed=true;
1590  if(--nbIndexed==0) return true;
1591  break;
1592  }
1593  }
1594  if(!(pos->isIndexed)) if(++nbUnIx>nbUnindexed) return false;
1595  }
1596  return false;
1597  }
1598  const bool storePossibleHKL=(useStoredHKL==2);
1599 
1600  if(storePossibleHKL)
1601  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos)
1602  {
1603  pos->isIndexed=false;
1604  pos->vDicVolHKL.clear();
1605  }
1606  else
1607  for(pos=dhkl.GetPeakList().begin();pos!=dhkl.GetPeakList().end();++pos) pos->isIndexed=false;
1608 
1609  int h,k,l;
1610  float dmax=dhkl.GetPeakList()[nb-1].d2obs;
1611  float dmin=dhkl.GetPeakList()[0 ].d2obs;
1612 
1613 
1614  int sk0,sl0;// do we need >0 *and* <0 indices for k,l ?
1615  switch(par.mlattice)
1616  {
1617  case TRICLINIC:
1618  sk0=-1;sl0=-1;
1619  break;
1620  case MONOCLINIC:
1621  {
1622  sk0=1;sl0=-1;
1623  break;
1624  }
1625  case ORTHOROMBIC:
1626  sk0=1;sl0=1;
1627  break;
1628  case HEXAGONAL:
1629  sk0=-1;sl0=1;
1630  break;
1631  case RHOMBOEDRAL:
1632  sk0=-1;sl0=-1;
1633  break;
1634  case TETRAGONAL:
1635  sk0=1;sl0=1;
1636  break;
1637  case CUBIC:
1638  sk0=1;sl0=1;
1639  break;
1640  // This should never happen. Avoid using unitialized values.
1641  default:
1642  throw 0;
1643  }
1644  int stepk,stepl;// steps in k,l to use for centered lattices
1645  switch(par.mCentering)
1646  {
1647  case LATTICE_P:stepk=1;stepl=1;break;
1648  case LATTICE_I:stepk=1;stepl=2;break;
1649  case LATTICE_A:stepk=1;stepl=2;break;
1650  case LATTICE_B:stepk=1;stepl=2;break;
1651  case LATTICE_C:stepk=2;stepl=1;break;
1652  case LATTICE_F:stepk=2;stepl=2;break;
1653  // This should never happen. Avoid using unitialized values.
1654  default: throw 0;
1655  }
1656  //RecUnitCell par0(par),par1(par);
1657  //for(unsigned int i=0;i<7;++i) {par0.par[i]-=dpar.par[i];par1.par[i]+=dpar.par[i];}
1658 
1659  //currently first & last unindexed dhkl
1660  first=dhkl.GetPeakList().begin(),last=dhkl.GetPeakList().end(),end=dhkl.GetPeakList().end();
1661 
1662  unsigned long nbCalcH,nbCalcK;// Number of calculated lines below dmax for each h,k
1663  for(h=0;;++h)
1664  {
1665  if(verbose) cout<<"H="<<h<<endl;
1666  nbCalcH=0;
1667  for(int sk=sk0;sk<=1;sk+=2)
1668  {
1669  if(h==0) sk=1;
1670  if(stepk==2) k=(h%2);// For LATTICE_C,LATTICE_F: h odd => k odd
1671  else k=0;
1672  for(;;k+=stepk)
1673  {
1674  if(verbose) cout<<"K="<<k*sk<<endl;
1675  nbCalcK=0;
1676  for(int sl=sl0;sl<=1;sl+=2)
1677  {
1678  int l0=0;
1679  if((h+k)==0)
1680  {
1681  sl=1;// No need to list 0 0 l with l<0
1682  l0=1;
1683  }
1684  else
1685  {
1686  if(h==0)
1687  {
1688  if(par.mlattice==MONOCLINIC) sl=1;// 0 k l and 0 k -l are equivalent
1689  if((sk<0)||(sl<0)) l0=1;// Do not list 0 k 0 with k<0
1690  else l0=0;// h==k==0 already covered
1691  }
1692  else
1693  {
1694  if(sl<0) l0=1;// Do not list h k 0 twice
1695  else l0=0;
1696  }
1697  }
1698  if(stepl==2)
1699  {
1700  if(par.mCentering==LATTICE_I) l0+=(h+k+l0)%2;
1701  if(par.mCentering==LATTICE_A) l0+=(k+l0)%2;// Start at k+l even
1702  if( (par.mCentering==LATTICE_B)
1703  ||(par.mCentering==LATTICE_F)) l0+=(h+l0)%2;// Start at h+l even
1704  }
1705  if(verbose) cout<<"SL="<<sl<<", L0="<<l0<<", STEPL="<<stepl<<", Centering="<<par.mCentering<<endl;
1706  for(l=l0;;l+=stepl)
1707  {
1708  if(verbose) cout<<"L="<<l<<","<<sl<<endl;
1709  float d0,d1;
1710  par.hkl2d_delta(h,sk*k,sl*l,dpar,d0,d1);
1711  if(d0<dmax) {nbCalcH++;nbCalcK++;}
1712  if((d1<dmin)&&(maxNbMissingBelow5==0)) continue;
1713  if(d0>dmax)
1714  {
1715  if(par.mlattice==TRICLINIC)
1716  {
1717  // Must check that d is increasing with l, otherwise we still need to increase it
1718  if(verbose) cout<<"L?="<< par.hkl2d(h,sk*k,sl*l,NULL,3)*sl <<", dmax="<<dmax<<endl;
1719  if((par.hkl2d(h,sk*k,sl*l,NULL,3)*sl)>0) break;
1720  }
1721  else break;
1722  }
1723  bool missing=(d0<d5)&&(maxNbMissingBelow5>0);
1724  for(pos=first;pos!=end;++pos)
1725  {
1726  if(pos==last) break;
1727  if((!storePossibleHKL)&&(pos->isIndexed)&&missing) continue;
1728  const float d2obs=pos->d2obs,d2obsmin=pos->d2obsmin, d2obsmax=pos->d2obsmax;
1729  if((d2obsmax>=d0) && (d1>=d2obsmin))
1730  {
1731  missing=false;
1732  if(!(pos->isIndexed))
1733  {
1734  pos->d2calc=(d0+d1)/2;
1735  --nbIndexed;
1736  pos->isIndexed=true;
1737  }
1738  if(verbose) cout<<d1<<" < ? <"<<d0<<"("<<h<<","<<sk*k<<","<<sl*l<<"): "<<d2obs<<" (remaining to index:"<<nbIndexed<<")"<<endl;
1739  if(storePossibleHKL)
1740  pos->vDicVolHKL.push_back(PeakList::hkl0(h,sk*k,sl*l));
1741  else
1742  {
1743  if((!storePossibleHKL)&&(nbIndexed==0)) return true;
1744  if(pos==first){first++;dmin=first->d2obsmin;}
1745  if(pos==last){last--;dmax=last->d2obsmax;}
1746  }
1747  }
1748  }
1749  if(missing) if(++nbMissingBelow5>=maxNbMissingBelow5)return false;
1750  }
1751  }
1752  if(nbCalcK==0) break;// && ((par.hkl2d(h,sk*k,0,NULL,2)*sk)>0)) break; // k too large
1753  }
1754  }
1755  if(nbCalcH==0) break;//h too large
1756  }
1757  if(verbose)
1758  {
1759  dhkl.Print(cout);
1760  }
1761  return nbIndexed<=0;
1762 }
1763 
1764 float CellExplorer::GetBestScore()const{return mBestScore;}
1765 const std::list<std::pair<RecUnitCell,float> >& CellExplorer::GetSolutions()const {return mvSolution;}
1766 std::list<std::pair<RecUnitCell,float> >& CellExplorer::GetSolutions() {return mvSolution;}
1767 
1768 unsigned int CellExplorer::RDicVol(RecUnitCell par0,RecUnitCell dpar, unsigned int depth,unsigned long &nbCalc,const float minV,const float maxV,vector<unsigned int> vdepth)
1769 {
1770  static bool localverbose=false;
1771  if(mlattice==TRICLINIC)
1772  {
1773  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];
1774  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];
1775  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];
1776 
1777  // a*<b*<c*
1778  if((p1m>p2p)||(p2m>p3p)) return 0;
1779 
1780  // max/min absolute values for p4,p5,p6
1781  if((p4m>p1p)||(-p4p>p1p)) return 0;//abs(p4)<p1 <=> b* < b*+/-a*
1782  if((p5m>p2p)||(-p5p>p2p)) return 0;//abs(p5)<p2 <=> c* < c*+/-b*
1783  if((p6m>p1p)||(-p6p>p1p)) return 0;//abs(p6)<p1 <=> c* < c*+/-a*
1784 
1785  const float max6=p1p+p2p-p4m-p5m;
1786  if((p6m>max6)||(-p6p>max6)) return 0;//abs(p6)<p1+p2-p4-p5 <=> c* < c*+/-a*+/-b*
1787 
1788  float p6mm,p6pp,p5mm,p5pp,p4mm,p4pp; // p6pp: smaller V*, larger V, etc..
1789  // derivative of V*^2 with p6
1790  if((p4*p5-2*p2*p6)>0) {p6pp=p6p;p6mm=p6m;}
1791  else{p6pp=p6m;p6mm=p6p;}
1792  // derivative of V*^2 with p5
1793  if((p4*p6-2*p1*p5)>0) {p5pp=p5p;p5mm=p5m;}
1794  else{p5pp=p5m;p5mm=p5p;}
1795  // derivative of V*^2 with p5
1796  if((p5*p6-2*p3*p4)>0) {p4pp=p4p;p4mm=p4m;}
1797  else{p4pp=p4m;p4mm=p4p;}
1798 
1799  //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))));
1800  //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))));
1801  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))));
1802  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))));
1803  if((vmin0>maxV)||(vmax0<minV)) return 0;
1804  }
1805  else
1806  if((depth>0)&&(depth<=2))// test if volume is within range
1807  {
1808  RecUnitCell parm=par0,parp=par0;
1809  for(unsigned int i=0;i<6;++i) {parm.par[i]-=dpar.par[i];parp.par[i]+=dpar.par[i];}
1810  vector<float> parmd=parm.DirectUnitCell();
1811  vector<float> parpd=parp.DirectUnitCell();
1812  if((parpd[6]>maxV)||(parmd[6]<minV))return 0;
1813  }
1814  unsigned int useStoredHKL=1;//Use already stored hkl
1815  if(depth==0) useStoredHKL=2; //Store possible hkl for all observed lines
1816 
1817  unsigned int maxMissingBelow5=0;
1818  // In the triclinic case, accept a maximum of 5 missing reflections below the 5th observed line
1819  if(mlattice==TRICLINIC) maxMissingBelow5=5;
1820 
1821  bool indexed=DichoIndexed(*mpPeakList,par0,dpar,mNbSpurious,localverbose,useStoredHKL,maxMissingBelow5);
1822 
1823  #if 0
1824  // If indexation failed but depth>=4, try adding a zero ?
1825  if( (!indexed) && (depth>=4))
1826  {//:TODO: Check if this is OK ! Vary value with depth
1827  dpar.par[0]=.0001;
1828  indexed=DichoIndexed(*mpPeakList,par0,dpar,mNbSpurious,false,useStoredHKL,maxMissingBelow5);
1829  //if(indexed) cout<<"Added zero - SUCCESS !"<<endl;
1830  }
1831  #endif
1832 
1833  if((indexed)&&(useStoredHKL==2))
1834  {
1835  // Test if two successive lines have been indexed exclusively with the same hkl
1836  unsigned int nbident=0;
1837  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();)
1838  {
1839  if(pos->vDicVolHKL.size()==1)
1840  {
1841  const PeakList::hkl0 h0=pos->vDicVolHKL.front();
1842  if(++pos==mpPeakList->GetPeakList().end()) break;
1843  if(pos->vDicVolHKL.size()==1)
1844  {
1845  const PeakList::hkl0 h1=pos->vDicVolHKL.front();
1846  if((h0.h==h1.h)&&(h0.k==h1.k)&&(h0.l==h1.l)) nbident++;
1847  if(nbident>mNbSpurious) {indexed=false;break;}
1848  }
1849  }
1850  else ++pos;
1851  }
1852  }
1853 
1854  // if we can zoom in for one parameter directly, we need per-parameter depth
1855  if(vdepth.size()==0)
1856  {
1857  vdepth.resize(mnpar-1);
1858  for(vector<unsigned int>::iterator pos=vdepth.begin();pos!=vdepth.end();) *pos++=depth;
1859  }
1860  else
1861  for(vector<unsigned int>::iterator pos=vdepth.begin();pos!=vdepth.end();++pos) if(*pos<depth)*pos=depth;
1862  #if 1
1863  if(false)//((useStoredHKL==2)&&(mNbSpurious==0)&&indexed)
1864  { // If high-d lines have been associated to a single reflection which is either h00, 0k0 or 00l,
1865  // jump the corresponding parameter to higher depth (mDicVolDepthReport, lowest depth report) immediately
1866  vector<pair<unsigned int,float> > vq0(3);
1867  for(unsigned int i=0;i<3;++i) {vq0[i].first=0;vq0[i].second=0.0;}
1868  RecUnitCell par0orig=par0,dparorig=dpar;
1869  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();++pos)
1870  {
1871  if(pos->vDicVolHKL.size()==1)
1872  {
1873  const PeakList::hkl0 h0=pos->vDicVolHKL.front();
1874  if((h0.k==0)&&(h0.l==0)) {vq0[0].first+=1 ; vq0[0].second+=pos->dobs/h0.h;}
1875  else
1876  if((h0.h==0)&&(h0.l==0)) {vq0[1].first+=1 ; vq0[1].second+=pos->dobs/h0.k;}
1877  else
1878  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;}
1879  }
1880  }
1881  switch(mlattice)
1882  {
1883  case TRICLINIC:
1884  {// In the triclinic case we may start with p1 and p2 already at depth>0 (triclinic quick tests)
1885  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;}
1886  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;}
1887  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;}
1888  break;
1889  }
1890  case MONOCLINIC:
1891  {
1892  if(vq0[0].first>0) {par0.par[1]=vq0[0].second/vq0[0].first ; vdepth[0]=mDicVolDepthReport;dpar.par[1]*=.0625;}
1893  if(vq0[1].first>0) {par0.par[2]=vq0[1].second/vq0[1].first ; vdepth[1]=mDicVolDepthReport;dpar.par[2]*=.0625;}
1894  if(vq0[2].first>0) {par0.par[3]=vq0[2].second/vq0[2].first ; vdepth[2]=mDicVolDepthReport;dpar.par[3]*=.0625;}
1895  break;
1896  }
1897  case ORTHOROMBIC:
1898  {
1899  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));}
1900  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));}
1901  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));}
1902  break;
1903  }
1904  case HEXAGONAL:break;
1905  case RHOMBOEDRAL:break;
1906  case TETRAGONAL:break;
1907  case CUBIC:break;
1908  }
1909  // If all parameters are at a higher depth, jump the global depth immediately
1910  unsigned int newdepth=40;
1911  for(vector<unsigned int>::iterator pos=vdepth.begin();pos!=vdepth.end();++pos) if(*pos<newdepth) newdepth=*pos;
1912  if(newdepth>depth) depth=newdepth;
1913  if((vq0[0].first>0)||(vq0[1].first>0)||(vq0[2].first>0))
1914  {
1915  indexed=DichoIndexed(*mpPeakList,par0,dpar,mNbSpurious,false,1,maxMissingBelow5);
1916  if(false)
1917  {
1918  {
1919  RecUnitCell parm=par0orig,parp=par0;
1920  for(unsigned int i=0;i<6;++i) {parm.par[i]-=dparorig.par[i];parp.par[i]+=dparorig.par[i];}
1921  vector<float> parmd=parm.DirectUnitCell();
1922  vector<float> parpd=parp.DirectUnitCell();
1923  char buf[200];
1924  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",
1925  parpd[0],parmd[0],parpd[1],parmd[1],parpd[2],parmd[2],parpd[3]*RAD2DEG,parmd[3]*RAD2DEG,
1926  parpd[4]*RAD2DEG,parmd[4]*RAD2DEG,parpd[5]*RAD2DEG,parmd[5]*RAD2DEG,parpd[6],parmd[6]);
1927  for(unsigned int i = 0; i < depth; ++i) cout << " ";
1928  cout<<buf<<"level="<<depth<<", indexed="<<indexed<<endl;
1929  }
1930  {
1931  RecUnitCell parm=par0,parp=par0;
1932  for(unsigned int i=0;i<6;++i) {parm.par[i]-=dpar.par[i];parp.par[i]+=dpar.par[i];}
1933  vector<float> parmd=parm.DirectUnitCell();
1934  vector<float> parpd=parp.DirectUnitCell();
1935  char buf[200];
1936  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",
1937  parpd[0],parmd[0],parpd[1],parmd[1],parpd[2],parmd[2],parpd[3]*RAD2DEG,parmd[3]*RAD2DEG,
1938  parpd[4]*RAD2DEG,parmd[4]*RAD2DEG,parpd[5]*RAD2DEG,parmd[5]*RAD2DEG,parpd[6],parmd[6]);
1939  for(unsigned int i = 0; i < depth; ++i) cout << " ";
1940  cout<<buf<<"level="<<depth<<", indexed="<<indexed<<endl;
1941  }
1942  }
1943  }
1944  }
1945  #endif
1946  if(false)//(depth==1)&&(rand()%10==0))
1947  {
1948  RecUnitCell parm=par0,parp=par0;
1949  for(unsigned int i=0;i<4;++i) {parm.par[i]-=dpar.par[i];parp.par[i]+=dpar.par[i];}
1950  for(unsigned int i=4;i<7;++i) {parm.par[i]+=dpar.par[i];parp.par[i]-=dpar.par[i];}
1951  vector<float> parmd=parm.DirectUnitCell();
1952  vector<float> parpd=parp.DirectUnitCell();
1953  char buf[200];
1954  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",
1955  parpd[0],parmd[0],parpd[1],parmd[1],parpd[2],parmd[2],parpd[3]*RAD2DEG,parmd[3]*RAD2DEG,
1956  parpd[4]*RAD2DEG,parmd[4]*RAD2DEG,parpd[5]*RAD2DEG,parmd[5]*RAD2DEG,parpd[6],parmd[6]);
1957  for(unsigned int i = 0; i < depth; ++i) cout << " ";
1958  cout<<buf<<"level="<<depth<<", indexed="<<indexed<<"("<<mvSolution.size()<<" sol.)"<<endl;
1959  }
1960  nbCalc++;
1961  // :TODO: if we failed the dichotomy and reached some depth, try guessing a zero shift from the indexed reflections
1962  /*
1963  if((!indexed)&&(depth>=2))
1964  {
1965  vector<float> shifts(mpPeakList->GetPeakList().size());
1966  vector<PeakList::hkl>::const_iterator peakpos=mpPeakList->GetPeakList().begin();
1967  for(vector<float>::iterator spos=shifts.begin();spos!=shifts.end();)
1968  { *spos++ = peakpos->d2diff * (float)(peakpos->isIndexed&&(!peakpos->isSpurious));peakpos++;}
1969  sort(shifts.begin(),shifts.end());
1970  par0.par[0]=shifts[mpPeakList->GetPeakList().size()/2];//use median value
1971  indexed=DichoIndexed(*mpPeakList,par0,dpar,mNbSpurious);
1972  if(indexed) cout<<"Failed Dicho ? Trying auto-zero shifting :Worked !"<<endl;
1973  }
1974  */
1975  if(indexed)
1976  {
1977  unsigned int deeperSolutions=0;
1978  if(depth<mMaxDicVolDepth)
1979  {
1980  if(false)//depth>=5)
1981  {
1982  RecUnitCell parm=par0,parp=par0;
1983  for(unsigned int i=0;i<6;++i) {parm.par[i]-=dpar.par[i];parp.par[i]+=dpar.par[i];}
1984  vector<float> parmd=parm.DirectUnitCell();
1985  vector<float> parpd=parp.DirectUnitCell();
1986  char buf[200];
1987  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",
1988  parpd[0],parmd[0],parpd[1],parmd[1],parpd[2],parmd[2],parpd[3]*RAD2DEG,parmd[3]*RAD2DEG,
1989  parpd[4]*RAD2DEG,parmd[4]*RAD2DEG,parpd[5]*RAD2DEG,parmd[5]*RAD2DEG,parpd[6],parmd[6]);
1990  for(unsigned int i=0;i<depth;++i) cout<<" ";
1991  cout<<"Depth="<<depth<<" "<<buf<<endl;
1992  for(int i=0;i<=6;++i)cout<<par0.par[i]<<",";
1993  for(int i=0;i<=6;++i)cout<<dpar.par[i]<<",";
1994  cout<<endl;
1995  }
1996  RecUnitCell par=par0;
1997  // zero (if used...)
1998  dpar.par[0]=0.5*dpar.par[0];
1999  // Divide interval by 2, except if this parameter is already at a higher depth
2000  // because a main axis has been indexed already.
2001  for(unsigned int i=1;i<mnpar;++i) dpar.par[i]*=(0.5+0.5*(vdepth[i-1]>depth));
2002 
2003  for(int i0=-1;i0<=1;i0+=2)
2004  {
2005  //:TODO: dichotomy on zero shift ?
2006  if(localverbose) cout<<__FILE__<<":"<<__LINE__<<":"<<par.par[3]<<" +/- "<<dpar.par[3]<<" ("<<vdepth[2]<<")"<<endl;
2007  // Don't change parameter if it is already determined at a higher depth
2008  if(vdepth[0]==depth) {par.par[1]=par0.par[1]+i0*dpar.par[1];}
2009  else {i0=2;}// no need to dicho this parameter which is already at higher depth
2010  if(mnpar==2)
2011  deeperSolutions+=RDicVol(par,dpar, depth+1,nbCalc,minV,maxV,vdepth);
2012  else
2013  for(int i1=-1;i1<=1;i1+=2)
2014  {
2015  if(vdepth[1]==depth) {par.par[2]=par0.par[2]+i1*dpar.par[2];}
2016  else {i1=2;}// no need to dicho this parameter which is already at higher depth
2017  if(mnpar==3)
2018  deeperSolutions+=RDicVol(par,dpar, depth+1,nbCalc,minV,maxV,vdepth);
2019  else
2020  for(int i2=-1;i2<=1;i2+=2)
2021  {
2022  if(vdepth[2]==depth) {par.par[3]=par0.par[3]+i2*dpar.par[3];}
2023  else {i2=2;}// no need to dicho this parameter which is already at higher depth
2024  if(mnpar==4)
2025  deeperSolutions+=RDicVol(par,dpar, depth+1,nbCalc,minV,maxV,vdepth);
2026  else
2027  for(int i3=-1;i3<=1;i3+=2)
2028  {
2029  if(vdepth[3]==depth)par.par[4]=par0.par[4]+i3*dpar.par[4];
2030  else i3=2;
2031  if(mnpar==5)
2032  deeperSolutions+=RDicVol(par,dpar, depth+1,nbCalc,minV,maxV,vdepth);
2033  else
2034  for(int i4=-1;i4<=1;i4+=2)
2035  {
2036  par.par[5]=par0.par[5]+i4*dpar.par[5];
2037  //if(mnpar==7)
2038  // deeperSolutions+=RDicVol(par,dpar, depth+1,nbCalc,minV,maxV,vdepth);
2039  //else
2040  for(int i5=-1;i5<=1;i5+=2)
2041  {
2042  par.par[6]=par0.par[6]+i5*dpar.par[6];
2043  //if(localverbose) cout<<__FILE__<<":"<<__LINE__<<":"<<par.par[3]<<" +/- "<<dpar.par[3]<<" ("<<vdepth[2]<<")"<<endl;
2044  deeperSolutions+=RDicVol(par,dpar, depth+1,nbCalc,minV,maxV,vdepth);
2045  }
2046  }
2047  }
2048  }
2049  }
2050  }
2051  }
2052  if((deeperSolutions==0) &&(depth>=mDicVolDepthReport))
2053  {
2054  mRecUnitCell=par0;
2055  vector<float> par=mRecUnitCell.DirectUnitCell();
2056  float score=Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true,false);
2057  // If we already have enough reports at higher depths (depth+2), don't bother record this one
2058  bool report=true;
2059  if(depth<(mMaxDicVolDepth-1))
2060  if(mvNbSolutionDepth[depth+2]>100)report=false;
2061  if(report && (((score>(mMinScoreReport*.5))&&(depth>=mDicVolDepthReport)) || (depth>=mMaxDicVolDepth)))
2062  {
2063  if(false)//score>) mBestScore//((score>mMinScoreReport)||(depth>=mDicVolDepthReport))
2064  cout<<__FILE__<<":"<<__LINE__<<" Depth="<<depth<<" (DIC) ! a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]<<", alpha="
2065  <<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG<<", V="<<par[6]
2066  <<", score="<<score<<endl;
2067  this->LSQRefine(5,true,true);
2068 
2069  // Re-score (may change to a better hkl indexing), and refine again
2070  score=Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true,false);
2071  this->LSQRefine(5,true,true);
2072 
2073  par=mRecUnitCell.DirectUnitCell();
2074  score=Score(*mpPeakList,mRecUnitCell,mNbSpurious,false,true,false);
2075  if( ((score>mMinScoreReport)||(depth>=mDicVolDepthReport))
2076  &&((mvSolution.size()<50)||(score>(mBestScore/3)))
2077  &&((mvSolution.size()<50)||(score>mMinScoreReport)))
2078  {
2079  if((score>(mBestScore))||((score>(mBestScore*0.8))&&(mvSolution.size()<50)))//||(rand()%100==0))
2080  {
2081  char buf[200];
2082  {
2083  RecUnitCell parm=par0,parp=par0;
2084  for(unsigned int i=0;i<4;++i) {parm.par[i]-=dpar.par[i];parp.par[i]+=dpar.par[i];}
2085  for(unsigned int i=4;i<7;++i) {parm.par[i]+=dpar.par[i];parp.par[i]-=dpar.par[i];}
2086  vector<float> parmd=parm.DirectUnitCell();
2087  vector<float> parpd=parp.DirectUnitCell();
2088  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",
2089  parpd[0],parmd[0],parpd[1],parmd[1],parpd[2],parmd[2],parpd[3]*RAD2DEG,parmd[3]*RAD2DEG,
2090  parpd[4]*RAD2DEG,parmd[4]*RAD2DEG,parpd[5]*RAD2DEG,parmd[5]*RAD2DEG,parpd[6],parmd[6]);
2091  for(unsigned int i = 0; i < depth; ++i) cout << " ";
2092 
2093  cout<<buf<<"level="<<depth<<", indexed="<<indexed<<"("<<mvSolution.size()<<" sol.)"<<endl;
2094  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",
2095  parp.par[1],parm.par[1],parp.par[2],parm.par[2],parp.par[3],parm.par[3],parp.par[4],parm.par[4],
2096  parp.par[5],parm.par[5],parp.par[6],parm.par[6]);
2097  for(unsigned int i = 0; i < depth; ++i) cout << " ";
2098  cout<<buf<<"level="<<depth<<", indexed="<<indexed<<"("<<mvSolution.size()<<" sol.)"<<endl;
2099  }
2100  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",
2101  par[0],par[1],par[2],par[3]*RAD2DEG,par[4]*RAD2DEG,par[5]*RAD2DEG,par[6],score,mvSolution.size());
2102  cout<<buf<<endl;
2103  mBestScore=score;
2104  }
2105  mvSolution.push_back(make_pair(mRecUnitCell,score));
2106  mvSolution.back().first.mNbSpurious = mNbSpurious;
2107  mvNbSolutionDepth[depth]+=1;
2108  if((mvSolution.size()>1100)&&(rand()%1000==0))
2109  {
2110  cout<<mvSolution.size()<<" solutions ! Redparing..."<<endl;
2111  this->ReduceSolutions(true);// This will update the min report score
2112  cout<<"-> "<<mvSolution.size()<<" remaining"<<endl;
2113  }
2114  }
2115  }
2116  }
2117  return deeperSolutions+1;
2118  }
2119  return 0;
2120 }
2121 
2122 vector<float> linspace(float min, float max,unsigned int nb)
2123 {
2124  vector<float> v(nb);
2125  for(unsigned int i=0;i<nb;++i) v[i]=min+(max-min)*i/(nb-1);
2126  return v;
2127 }
2128 
2129 void CellExplorer::DicVol(const float minScore,const unsigned int minDepth,const float stopOnScore,const unsigned int stopOnDepth)
2130 {
2131  mNbLSQExcept=0;
2132  mDicVolDepthReport=minDepth;
2133  mMinScoreReport=minScore;
2134  this->Init();
2135  if(minDepth>mMaxDicVolDepth) mMaxDicVolDepth=minDepth;
2136  mvNbSolutionDepth.resize(mMaxDicVolDepth+1);
2137  for(unsigned int i=0;i<=mMaxDicVolDepth;++i) mvNbSolutionDepth[i]=0;
2138 
2139  float latstep=0.5,
2140  vstep=(mVolumeMax-mVolumeMin)/(ceil((mVolumeMax-mVolumeMin)/500)-0.0001);
2141  mCosAngMax=abs(cos(mAngleMax));
2142  const float cosangstep=mCosAngMax/(ceil(mCosAngMax/.08)-.0001);
2143  if(((mVolumeMax-mVolumeMin)/vstep)>10) vstep=(mVolumeMax-mVolumeMin)/9.999;
2144  if(((mLengthMax-mLengthMin)/latstep)>25) latstep=(mLengthMax-mLengthMin)/24.9999;
2145 
2146  cout<<mLengthMin<<"->"<<mLengthMax<<","<<latstep<<","<<(mLengthMax-mLengthMin)/latstep<<endl;
2147  cout<<mAngleMin<<"->"<<mAngleMax<<","<<cosangstep<<","<<mCosAngMax<<","<<(mAngleMax-mAngleMin)/cosangstep<<endl;
2148  cout<<mVolumeMin<<"->"<<mVolumeMax<<","<<vstep<<","<<(mVolumeMax-mVolumeMin)/vstep<<endl;
2149  RecUnitCell par0,dpar;
2150  par0.mlattice=mlattice;
2151  dpar.mlattice=mlattice;
2152  par0.mCentering=mCentering;
2153  dpar.mCentering=mCentering;
2154  //Zero shift parameter - not used for dicvol right now ? :TODO:
2155  par0.par[0]=0.0;
2156  dpar.par[0]=0.0;
2157  unsigned long nbCalc=0;
2158  Chronometer chrono;
2159  float bestscore=0;
2160  list<pair<RecUnitCell,float> >::iterator bestpos;
2161  bool breakDepth=false;
2162  // In the triclinic case, first try assigning a* and b* from the first reflections
2163  if(false) //mlattice==TRICLINIC)
2164  for(float minv=mVolumeMin;minv<mVolumeMax;minv+=vstep)
2165  {
2166  float maxv=minv+vstep;
2167  if(maxv>mVolumeMax)maxv=mVolumeMax;
2168  cout<<"Starting: V="<<minv<<"->"<<maxv<<endl;
2169  const float minr=1/(mLengthMax*mLengthMax);
2170  const float maxr=1/(mLengthMin*mLengthMin);
2171  const float stepr=(maxr-minr)/24.999;
2172  float p1,p2;
2173  for(unsigned int i=0;i<=5;i++)
2174  {
2175  switch(i)
2176  {// Try to find a and b from the first observed reflections
2177  case 0: p1=mpPeakList->GetPeakList()[0].d2obs ;p2=mpPeakList->GetPeakList()[1].d2obs ; break;
2178  case 1: p1=mpPeakList->GetPeakList()[0].d2obs ;p2=mpPeakList->GetPeakList()[2].d2obs ; break;
2179  case 2: p1=mpPeakList->GetPeakList()[1].d2obs/2;p2=mpPeakList->GetPeakList()[0].d2obs ; break;
2180  case 3: p1=mpPeakList->GetPeakList()[1].d2obs/2;p2=mpPeakList->GetPeakList()[2].d2obs ; break;
2181  case 4: p1=mpPeakList->GetPeakList()[2].d2obs/2;p2=mpPeakList->GetPeakList()[0].d2obs ; break;
2182  case 5: p1=mpPeakList->GetPeakList()[2].d2obs/2;p2=mpPeakList->GetPeakList()[1].d2obs ; break;
2183  }
2184  //if(i>0) exit(0);
2185  if(p1>p2) continue;
2186  cout<<"Trying #"<<i<<": a*="<<p1<<", b*="<<p2<<endl;
2187  float min3r=p2,
2188  max3r=maxr;//:TODO: use larger value to take angles into account ?
2189  const float step3r=(max3r-min3r)/(ceil((max3r-min3r)/stepr)-.001);
2190  vector<unsigned int> vdepth(mnpar-1);
2191  for(vector<unsigned int>::iterator pos=vdepth.begin();pos!=vdepth.end();) *pos++=0;
2192  vdepth[0]=3;
2193  vdepth[1]=3;
2194  for(float p3=min3r;p3<max3r;p3+=step3r)
2195  {
2196  //cout<<" p3="<<p3<<endl;
2197  const float max4r=p3+step3r;
2198  const float step4r=max4r/(ceil(max4r/stepr)-.001);
2199  for(float p4=0;p4<max4r;p4+=step4r)
2200  {
2201  //cout<<" p4="<<p4<<endl;
2202  float max5r=(p2+stepr);
2203  const float step5r=max5r/(ceil(max5r/stepr)-.001);
2204  for(float p5=0;p5<max5r;p5+=step5r)
2205  {
2206  float max6r=(p1+stepr);
2207  const float step6r=max6r/(ceil(max6r/stepr)-.001);
2208  for(float p6=-max6r;p6<max6r;p6+=step6r)
2209  {
2210  //cout<<" p6="<<p6<<"/"<<p1<<"/"<<p3<<endl;
2211  dpar.par[1]=stepr*pow(float(0.51),int(vdepth[0]));
2212  dpar.par[2]=stepr*pow(float(0.51),int(vdepth[1]));
2213  dpar.par[3]=step3r*0.51;
2214  dpar.par[4]=step4r*0.51;
2215  dpar.par[5]=step5r*0.51;
2216  dpar.par[6]=step6r*0.51;
2217 
2218  par0.par[0]=0;
2219  par0.par[1]=p1;
2220  par0.par[2]=p2;
2221  par0.par[3]=p3+step3r/2;
2222  par0.par[4]=p4+step4r/2;
2223  par0.par[5]=p5+step5r/2;
2224  par0.par[6]=p6+step6r/2;
2225  //for(int i=0;i<=6;++i)cout<<par0.par[i]<<",";
2226  //cout<<endl;
2227  //for(int i=0;i<=6;++i)cout<<dpar.par[i]<<",";
2228  //cout<<endl;
2229  RDicVol(par0,dpar,0,nbCalc,minv,maxv,vdepth);
2230  }
2231  }
2232  }
2233  }
2234  cout<<"Finished trying: a*="<<p1<<" A, b*="<<p2<<" A, "<<nbCalc
2235  <<" unit cells tested, "<<nbCalc/chrono.seconds()<<" tests/s, Elapsed time="
2236  <<chrono.seconds()<<"s, Best score="<<mBestScore<<", "<<stopOnScore<<", "<<breakDepth<<endl;
2237  breakDepth=false;
2238  if(stopOnDepth>0)
2239  for(unsigned int i=stopOnDepth; i<mvNbSolutionDepth.size();++i)
2240  if(mvNbSolutionDepth[i]>1) {breakDepth=true;break;}
2241  if((mBestScore>stopOnScore)&&(breakDepth)) break;
2242  }//cases
2243  cout<<"Finished triclinic QUICK tests for: V="<<minv<<"->"<<maxv<<" A^3, "<<nbCalc
2244  <<" unit cells tested, "<<nbCalc/chrono.seconds()<<" tests/s, Elapsed time="
2245  <<chrono.seconds()<<"s, Best score="<<mBestScore<<endl;
2246  if((mBestScore>stopOnScore)&&(breakDepth)) break;
2247  }//volume
2248  if((mBestScore<stopOnScore)||(!breakDepth))
2249  for(float minv=mVolumeMin;minv<mVolumeMax;minv+=vstep)
2250  {
2251  float maxv=minv+vstep;
2252  if(maxv>mVolumeMax)maxv=mVolumeMax;
2253  cout<<"Starting: V="<<minv<<"->"<<maxv<<endl;
2254  switch(mlattice)
2255  {
2256  case TRICLINIC:
2257  {
2258  const unsigned int nbstep=25;
2259  vector<float> v1=linspace(mLengthMin,mLengthMax,nbstep);
2260  const float lstep=v1[1]-v1[0];
2261  for(unsigned int i1=0;i1<(nbstep-1);++i1)
2262  {
2263  const float p1 =(1/(v1[i1]*v1[i1])+1/(v1[i1+1]*v1[i1+1]))/2;
2264  const float dp1=(1/(v1[i1]*v1[i1])-1/(v1[i1+1]*v1[i1+1]))/2;
2265  //cout<<"p1="<<p1<<endl;
2266  //cout<<(v1[i1+1]-mLengthMin)/lstep+2.001<<endl;
2267  const unsigned int nb2=int((v1[i1+1]-mLengthMin)/lstep+2.001);
2268  vector<float> v2=linspace(mLengthMin,v1[i1+1],nb2);
2269  //for(unsigned int i2=0;i2<nb2;++i2) cout<<1/v2[i2]/v2[i2]<<" ";
2270  //cout<<endl;
2271  for(unsigned int i2=0;i2<(nb2-1);++i2)
2272  {
2273  const float p2 =(1/(v2[i2]*v2[i2])+1/(v2[i2+1]*v2[i2+1]))/2;
2274  const float dp2=(1/(v2[i2]*v2[i2])-1/(v2[i2+1]*v2[i2+1]))/2;
2275  //cout<<" p2="<<p2<<endl;
2276  const unsigned int nb3=int((v2[i2+1]-mLengthMin)/lstep+2.001);
2277  vector<float> v3=linspace(mLengthMin,v2[i2+1],nb3);
2278  for(unsigned int i3=0;i3<(nb3-1);++i3)
2279  {
2280  const float p3 =(1/(v3[i3]*v3[i3])+1/(v3[i3+1]*v3[i3+1]))/2;
2281  const float dp3=(1/(v3[i3]*v3[i3])-1/(v3[i3+1]*v3[i3+1]))/2;
2282 
2283  //const float vmax3=v1[i1+1]*v2[i2+1]*v3[i3+1];
2284  //const float vmin3=v1[i1]*v2[i2]*v3[i3]/5;
2285  const float vmin3=1/sqrt((p1+dp1)*(p2+dp2)*(p3+dp3));
2286  const float vmax3=1/sqrt((p1-dp1)*(p2-dp2)*(p3-dp3))*2; // *2 - sufficient margin to take into account angles
2287  if(vmax3<minv) continue;
2288  if(vmin3>maxv) continue;
2289 
2290  //cout<<" p3="<<p3<<endl;
2291  //char buf[200];
2292  //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",
2293  // p1,dp1,1/sqrt(p1+dp1),1/sqrt(p1-dp1),p2,dp2,1/sqrt(p2+dp2),1/sqrt(p2-dp2),
2294  // p3,dp3,1/sqrt(p3+dp3),1/sqrt(p3-dp3),vmin3,vmax3);
2295  //cout<<buf<<endl;
2296  const unsigned int nb4=int((p1+dp1)/(4*dp1)+2.001);
2297  vector<float> v4=linspace(0,p1+dp1,nb4);
2298  for(unsigned int i4=0;i4<(nb4-1);++i4)
2299  {
2300  const float p4 =(v4[i4+1]+v4[i4])/2;
2301  const float dp4=(v4[i4+1]-v4[i4])/2;
2302  //cout<<" p4="<<p4<<endl;
2303  const unsigned int nb5=int((p2+dp2)/(4*dp2)+2.001);
2304  vector<float> v5=linspace(0,p2+dp2,nb5);
2305  for(unsigned int i5=0;i5<(nb5-1);++i5)
2306  {
2307  const float p5 =(v5[i5+1]+v5[i5])/2;
2308  const float dp5=(v5[i5+1]-v5[i5])/2;
2309  //cout<<" p5="<<p5<<endl;
2310 
2311  float vmax6=(p1+dp1)+(p2+dp2)-(p4-dp4)-(p5-dp5);
2312  if(vmax6>(p1+dp1)) vmax6=p1+dp1;
2313  if(vmax6<0) continue;
2314  const unsigned int nb6=int((2*vmax6)/(4*dp1)+2.001);
2315  vector<float> v6=linspace(-vmax6,vmax6,nb6);
2316  //cout<<" p6="<<nb6<<","<<vmax6<<endl;
2317  for(unsigned int i6=0;i6<(nb6-1);++i6)
2318  {
2319  const float p6 =(v6[i6+1]+v6[i6])/2;
2320  const float dp6=(v6[i6+1]-v6[i6])/2;
2321  //cout<<" p6="<<p6<<endl;
2322 
2323  //char buf[200];
2324  //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 ",
2325  // p1-dp1,p1+dp1,p2-dp2,p2+dp2,p3-dp3,p3+dp3,p4-dp4,p4+dp4,p5-dp5,p5+dp5,p6-dp6,p6+dp6);
2326  //cout<<"Testing: "<<buf<<endl;
2327  //cout<<i1<<","<<i2<<","<<i3<<","<<i4<<","<<i5<<","<<i6<<endl;
2328  dpar.par[1]=dp1;
2329  dpar.par[2]=dp2;
2330  dpar.par[3]=dp3;
2331  dpar.par[4]=dp4;
2332  dpar.par[5]=dp5;
2333  dpar.par[6]=dp6;
2334 
2335  par0.par[0]=0;
2336  par0.par[1]=p1;
2337  par0.par[2]=p2;
2338  par0.par[3]=p3;
2339  par0.par[4]=p4;
2340  par0.par[5]=p5;
2341  par0.par[6]=p6;
2342 
2343  RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2344  }
2345  }
2346  }
2347  }
2348  }
2349  }
2350  break;
2351  }
2352  case MONOCLINIC:
2353  {
2354  RecUnitCell parlarge,//parm: smallest reciprocal, largest direct cell
2355  parsmall;//parp: largest reciprocal, smallest direct cell
2356  vector<float> parlarged,parsmalld;
2357  latstep=(mLengthMax-mLengthMin)/24.999;
2358  for(float x4=0;x4<mCosAngMax+cosangstep;x4+=cosangstep)
2359  {
2360  const float sinbeta=sqrt(abs(1-x4*x4));
2361  float x1=mLengthMin;
2362  for(;x1<mLengthMax;x1+=latstep)
2363  {
2364  float x2=mLengthMin;
2365  for(;x2<mLengthMax;x2+=latstep)
2366  {
2367  float x3=x1;
2368  const float x3step=(mLengthMax-x1)/(ceil((mLengthMax-x1)/latstep)-0.001);
2369  for(;x3<mLengthMax;x3+=x3step) //x3+=(latstep+x3*sin4)
2370  {
2371  if((x3*x4)>x1) break;// | c * cos(beta) | <a
2372  dpar.par[1]=(1/(x1)-1/(x1+latstep))*0.5/sinbeta;
2373  dpar.par[2]=(1/(x2)-1/(x2+latstep))*0.5/sinbeta;
2374  dpar.par[3]=(1/(x3)-1/(x3+x3step ))*0.5/sinbeta;
2375  dpar.par[4]=cosangstep*0.5;
2376 
2377  par0.par[0]=0;
2378  par0.par[1]=(1/(x1)+1/(x1+latstep))*0.5/sinbeta;
2379  par0.par[2]=(1/(x2)+1/(x2+latstep))*0.5/sinbeta;
2380  par0.par[3]=(1/(x3)+1/(x3+x3step ))*0.5/sinbeta;
2381  par0.par[4]=x4+cosangstep*.5;
2382 
2383  const float smallv=x1*x2*x3*sinbeta;
2384  if(smallv>maxv) break;
2385  const float largev=(x1+latstep)*(x2+latstep)*(x3+latstep)*(sinbeta+cosangstep);
2386  if(largev<minv) continue;
2387  /*
2388  char buf[200];
2389  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",
2390  parsmalld[0],parlarged[0],parsmalld[1],parlarged[1],parsmalld[2],parlarged[2],parsmalld[3]*RAD2DEG,parlarged[3]*RAD2DEG,
2391  parsmalld[4]*RAD2DEG,parlarged[4]*RAD2DEG,parsmalld[5]*RAD2DEG,parlarged[5]*RAD2DEG,parsmalld[6],parlarged[6]);
2392  cout<<buf<<" VM="<<maxv<<", x3="<<x3<<endl;
2393  */
2394  RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2395  }//x3
2396  //if(((parsmalld[6]>maxv)&&(x3==x1))||(parlarged[1]>mLengthMax)) break;
2397  }//x2
2398  }//x1
2399  // Test if we have one solution before going to the next angle range
2400  for(list<pair<RecUnitCell,float> >::iterator pos=mvSolution.begin();pos!=mvSolution.end();++pos)
2401  {
2402  const float score=pos->second;//Score(*mpPeakList,pos->first,mNbSpurious);
2403  if(score>bestscore) {bestscore=score;bestpos=pos;}
2404  }
2405  bool breakDepth=false;
2406  if(stopOnDepth>0)
2407  for(unsigned int i=stopOnDepth; i<mvNbSolutionDepth.size();++i)
2408  if(mvNbSolutionDepth[i]>1) {breakDepth=true;break;}
2409  if((bestscore>stopOnScore)&&(breakDepth)) break;
2410  }//x4
2411  break;
2412  }
2413  case ORTHOROMBIC:
2414  {
2415  if(false)
2416  {
2417  // Test 7.677350 5.803770 10.313160 V=480
2418  //const float a=7.75,b=5.75,c=10.25;
2419  // Test 6.062000 16.779400 16.8881 v=1750
2420  const float a=6.25,b=16.75,c=16.75;
2421  dpar.par[1]=(1/(a-.25)-1/(a+.25))*0.5;
2422  dpar.par[2]=(1/(b-.25)-1/(b+.25))*0.5;
2423  dpar.par[3]=(1/(c-.25)-1/(c+.25))*0.5;
2424  par0.par[0]=0;
2425  par0.par[1]=1/a;
2426  par0.par[2]=1/b;
2427  par0.par[3]=1/c;
2428  RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2429  break;
2430  }
2431  latstep=(mLengthMax-mLengthMin)/24.999;
2432  for(float x1=mLengthMin;x1<mLengthMax;x1+=latstep)
2433  {
2434  for(float x2=x1;x2<mLengthMax;x2+=latstep)
2435  {
2436  for(float x3=x2;x3<mLengthMax;x3+=latstep)
2437  {
2438  dpar.par[1]=(1/(x1)-1/(x1+latstep))*0.5;
2439  dpar.par[2]=(1/(x2)-1/(x2+latstep))*0.5;
2440  dpar.par[3]=(1/(x3)-1/(x3+latstep))*0.5;
2441 
2442  par0.par[0]=0;
2443  par0.par[1]=(1/(x1)+1/(x1+latstep))*0.5;
2444  par0.par[2]=(1/(x2)+1/(x2+latstep))*0.5;
2445  par0.par[3]=(1/(x3)+1/(x3+latstep))*0.5;
2446 
2447  const float vmin=x1*x2*x3,vmax=(x1+latstep)*(x2+latstep)*(x3+latstep);
2448  if(vmin>maxv) break;
2449  if(vmax>=minv) RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2450  }
2451  if((x1*x2*x2)>maxv) break;
2452  }
2453  if((x1*x1*x1)>maxv) break;
2454  }
2455  break;
2456  }
2457  case HEXAGONAL:
2458  {
2459  vector<float> parlarged,parsmalld;// Small & large UC in direct space
2460  latstep=(mLengthMax-mLengthMin)/24.999;
2461  for(float x1=mLengthMin;;x1+=latstep)
2462  {
2463  for(float x2=mLengthMin;x2<(mLengthMax+latstep);x2+=latstep)
2464  {
2465  dpar.par[1]=(1/(x1)-1/(x1+latstep))*0.5;
2466  dpar.par[2]=(1/(x2)-1/(x2+latstep))*0.5;
2467 
2468  par0.par[0]=0;
2469  par0.par[1]=(1/(x1)+1/(x1+latstep))*0.5;
2470  par0.par[2]=(1/(x2)+1/(x2+latstep))*0.5;
2471 
2472  RecUnitCell parlarge=par0,parsmall=par0;
2473  for(unsigned int i=0;i<6;++i) {parlarge.par[i]-=dpar.par[i];parsmall.par[i]+=dpar.par[i];}
2474  parlarged=parlarge.DirectUnitCell();
2475  parsmalld=parsmall.DirectUnitCell();
2476  /*
2477  char buf[200];
2478  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",
2479  parsmalld[0],parlarged[0],parsmalld[1],parlarged[1],parsmalld[2],parlarged[2],parsmalld[3]*RAD2DEG,parlarged[3]*RAD2DEG,
2480  parsmalld[4]*RAD2DEG,parlarged[4]*RAD2DEG,parsmalld[5]*RAD2DEG,parlarged[5]*RAD2DEG,parsmalld[6],parlarged[6]);
2481  */
2482  if((parsmalld[6]<maxv)&&(parlarged[6]>minv))
2483  {
2484  //cout<<buf<<endl;
2485  RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2486  }
2487  //else cout<<buf<<" BREAK"<<endl;
2488  }
2489  if(parlarged[0]>mLengthMax) break;
2490  }
2491  break;
2492  }
2493  case RHOMBOEDRAL: //:TODO:
2494  {
2495  latstep=(mLengthMax-mLengthMin)/24.999;
2496  for(float x1=mLengthMin;x1<(mLengthMax+latstep);x1+=latstep)
2497  {
2498  for(float x2=0;x2<mCosAngMax+cosangstep;x2+=cosangstep)
2499  {
2500  dpar.par[1]=latstep/2*1.1;
2501  dpar.par[2]=cosangstep/2*1.1;
2502 
2503  par0.par[0]=0;
2504  par0.par[1]=x1-latstep/2*1.1;
2505  par0.par[2]=x2-cosangstep/2*1.1;
2506  vector<float> par=par0.DirectUnitCell();
2507  if((par[6]<maxv)&&(par[6]>minv))
2508  {
2509  RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2510  }
2511  }
2512  }
2513  break;
2514  }
2515  case TETRAGONAL:
2516  {
2517  vector<float> parlarged,parsmalld;// Small & large UC in direct space
2518  latstep=(mLengthMax-mLengthMin)/24.999;
2519  for(float x1=mLengthMin;x1<mLengthMax;x1+=latstep)
2520  {
2521  for(float x2=mLengthMin;x2<mLengthMax;x2+=latstep)
2522  {
2523  dpar.par[1]=(1/(x1)-1/(x1+latstep))*0.5;
2524  dpar.par[2]=(1/(x2)-1/(x2+latstep))*0.5;
2525 
2526  par0.par[0]=0;
2527  par0.par[1]=(1/(x1)+1/(x1+latstep))*0.5;
2528  par0.par[2]=(1/(x2)+1/(x2+latstep))*0.5;
2529 
2530  RecUnitCell parlarge=par0,parsmall=par0;
2531  for(unsigned int i=0;i<6;++i) {parlarge.par[i]-=dpar.par[i];parsmall.par[i]+=dpar.par[i];}
2532  parlarged=parlarge.DirectUnitCell();
2533  parsmalld=parsmall.DirectUnitCell();
2534  /*
2535  char buf[200];
2536  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",
2537  parsmalld[0],parlarged[0],parsmalld[1],parlarged[1],parsmalld[2],parlarged[2],parsmalld[3]*RAD2DEG,parlarged[3]*RAD2DEG,
2538  parsmalld[4]*RAD2DEG,parlarged[4]*RAD2DEG,parsmalld[5]*RAD2DEG,parlarged[5]*RAD2DEG,parsmalld[6],parlarged[6]);
2539  */
2540  if((parsmalld[6]<maxv)&&(parlarged[6]>minv))
2541  {
2542  RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2543  }
2544  if(parsmalld[6]>maxv) break;
2545  }
2546  if((x1*mLengthMin*mLengthMin)>maxv) break;
2547  }
2548  break;
2549  }
2550  case CUBIC:
2551  {
2552  latstep=(mLengthMax-mLengthMin)/24.999;
2553  cout<<mLengthMax<<","<<mLengthMin<<","<<latstep<<endl;
2554  for(float x1=mLengthMin;x1<(mLengthMax+latstep);x1+=latstep)
2555  {
2556  dpar.par[1]=(1/(x1)-1/(x1+latstep))*0.5;
2557 
2558  par0.par[0]=0;
2559  par0.par[1]=(1/(x1)+1/(x1+latstep))*0.5;
2560 
2561  const float vmin=x1*x1*x1,vmax=(x1+latstep)*(x1+latstep)*(x1+latstep);
2562  if(vmin>maxv)break;
2563  if(vmax>minv) RDicVol(par0,dpar,0,nbCalc,minv,maxv);
2564  }
2565  break;
2566  }
2567  }
2568  cout<<"Finished: V="<<minv<<"->"<<maxv<<" A^3, "<<nbCalc
2569  <<" unit cells tested, "<<nbCalc/chrono.seconds()<<" tests/s, Elapsed time="
2570  <<chrono.seconds()<<"s"<<endl;
2571  for(list<pair<RecUnitCell,float> >::iterator pos=mvSolution.begin();pos!=mvSolution.end();++pos)
2572  {
2573  const float score=pos->second;//Score(*mpPeakList,pos->first,mNbSpurious);
2574  if(score>bestscore) {bestscore=score;bestpos=pos;}
2575  }
2576  bool breakDepth=false;
2577  if(stopOnDepth>0)
2578  for(unsigned int i=stopOnDepth; i<mvNbSolutionDepth.size();++i)
2579  if(mvNbSolutionDepth[i]>1) {breakDepth=true;break;}
2580  if((bestscore>stopOnScore)&&(breakDepth)) break;
2581  }
2582  /*
2583  {// Tag spurious lines
2584  vector<int> vSpuriousScore;
2585  for(vector<PeakList::hkl>::const_iterator pos=mpPeakList->GetPeakList().begin();pos!=mpPeakList->GetPeakList().end();++pos)
2586  vSpuriousScore.push_back(pos->stats);
2587  sort(vSpuriousScore.begin(),vSpuriousScore.end());
2588  const int threshold=vSpuriousScore[vSpuriousScore.size()/2]*5;
2589  for(vector<PeakList::hkl>::iterator pos=mpPeakList->mvHKL.begin();pos!=mpPeakList->mvHKL.end();++pos)
2590  if(pos->stats > threshold) pos->isSpurious=true;
2591  else pos->isSpurious=false;
2592  mpPeakList->Print(cout);
2593  }
2594  */
2595  this->ReduceSolutions(true);
2596  bestscore=0;bestpos=mvSolution.end();
2597  for(list<pair<RecUnitCell,float> >::iterator pos=mvSolution.begin();pos!=mvSolution.end();++pos)
2598  {
2599  const float score=Score(*mpPeakList,pos->first,mNbSpurious);
2600  if(score>bestscore) {bestpos=pos;bestscore=score;}
2601  vector<float> par=pos->first.DirectUnitCell();
2602  cout<<__FILE__<<":"<<__LINE__<<" Solution ? a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]
2603  <<", alpha="<<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG
2604  <<", V="<<par[6]<<", score="<<score<<endl;
2605  }
2606  if(bestpos!=mvSolution.end())
2607  {
2608  vector<float> par=bestpos->first.DirectUnitCell();
2609  cout<<__FILE__<<":"<<__LINE__<<" BEST ? a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]
2610  <<", alpha="<<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG
2611  <<", V="<<par[6]<<", score="<<bestscore<<endl;
2612  cout<<nbCalc<<"unit cells tested, "<<nbCalc/chrono.seconds()<<" tests/s, Elapsed time="
2613  <<chrono.seconds()<<"s"<<endl;
2614  }
2615 }
2616 
2617 bool SimilarRUC(const RecUnitCell &c0,const RecUnitCell &c1, const float delta=0.005)
2618 {
2619  if(c0.mNbSpurious != c1.mNbSpurious) return false;
2620  vector<float> par0=c0.DirectUnitCell();
2621  vector<float> par1=c1.DirectUnitCell();
2622  float diff=0;
2623  for(unsigned int i=0;i<6;++i) diff += abs(par0[i]-par1[i]);
2624  return (diff/6)<delta;
2625 }
2626 
2627 bool compareRUCScore(std::pair<RecUnitCell,float> &p1, std::pair<RecUnitCell,float> &p2)
2628 {
2629  return p1.second > p2.second;
2630 }
2631 
2632 void CellExplorer::ReduceSolutions(const bool updateReportThreshold)
2633 {
2634  const bool verbose=false;
2635  std::list<std::pair<RecUnitCell,float> > vSolution2;
2636  // TODO: take into account number of spurious lines for cutoff value.
2637  // keep only solutions above mBestScore/5
2638  for(list<pair<RecUnitCell,float> >::iterator pos=mvSolution.begin();pos!=mvSolution.end();)
2639  {
2640  if(pos->second<(mBestScore/5)) pos=mvSolution.erase(pos);
2641  else ++pos;
2642  }
2643  if(updateReportThreshold&& ((mBestScore/5)>mMinScoreReport))
2644  {
2645  cout<<"CellExplorer::ReduceSolutions(): update threshold for report from "
2646  <<mMinScoreReport<<" to "<<mBestScore/5<<endl;
2647  mMinScoreReport=mBestScore/5;
2648  }
2649 
2650  while(mvSolution.size()>0)
2651  {
2652  vSolution2.push_back(mvSolution.front());
2653  mvSolution.pop_front();
2654  vector<float> par=vSolution2.back().first.DirectUnitCell();
2655  if(verbose)
2656  cout<<__FILE__<<":"<<__LINE__<<" SOLUTION: a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]
2657  <<", alpha="<<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG
2658  <<", V="<<par[6]<<", score="<<vSolution2.back().second<<", SIMILAR TO:"<<endl;
2659  for(list<pair<RecUnitCell,float> >::iterator pos=mvSolution.begin();pos!=mvSolution.end();)
2660  {
2661  if(SimilarRUC(pos->first,vSolution2.back().first))
2662  {
2663  par=pos->first.DirectUnitCell();
2664  if(verbose)
2665  cout<<__FILE__<<":"<<__LINE__<<" 1: a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]
2666  <<", alpha="<<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG
2667  <<", V="<<par[6]<<", score="<<pos->second<<" ("<<mvSolution.size()<<")"<<endl;
2668  if(vSolution2.back().first.mlattice==pos->first.mlattice)
2669  {
2670  if(pos->second>vSolution2.back().second) vSolution2.back()=*pos;
2671  }
2672  else if(vSolution2.back().first.mlattice<pos->first.mlattice) vSolution2.back()=*pos;
2673  pos=mvSolution.erase(pos);
2674  }
2675  else
2676  {
2677  par=pos->first.DirectUnitCell();
2678  if(verbose)
2679  cout<<__FILE__<<":"<<__LINE__<<" 0: a="<<par[0]<<", b="<<par[1]<<", c="<<par[2]
2680  <<", alpha="<<par[3]*RAD2DEG<<", beta="<<par[4]*RAD2DEG<<", gamma="<<par[5]*RAD2DEG
2681  <<", V="<<par[6]<<", score="<<pos->second<<" ("<<mvSolution.size()<<")"<<endl;
2682  ++pos;
2683  }
2684  }
2685  }
2686  mvSolution=vSolution2;
2687  mvSolution.sort(compareRUCScore);
2688 
2689  // keep at most 100 solutions, update mDicVolDepthReport and mMinScoreReport if necessary
2690  if(mvSolution.size()>100)
2691  {
2692  mvSolution.resize(100);
2693  if(updateReportThreshold && (mvSolution.back().second>mMinScoreReport))
2694  {
2695  cout<<"CellExplorer::ReduceSolutions(): update threshold for report from "
2696  <<mMinScoreReport<<" to "<<mvSolution.back().second<<endl;
2697  mMinScoreReport=mvSolution.back().second;
2698  }
2699  }
2700 }
2701 
2702 
2703 void CellExplorer::Init()
2704 {
2705  // Prepare global optimisation
2706  //for(unsigned int i=0;i<mpPeakList->nb;++i)
2707  // cout<<__FILE__<<":"<<__LINE__<<":d*="<<mpPeakList->mvdobs[i]<<", d*^2="<<mpPeakList->mvd2obs[i]<<endl;
2708  srand(time(NULL));
2709  vector<pair<RecUnitCell,float> >::iterator pos;
2710  const float min_latt=1./mLengthMax;
2711  const float max_latt=1./mLengthMin;
2712  const float amp_crossp=abs(cos(mAngleMax));
2713  //mMin[0]=-.002;mAmp[0]=.004;
2714  mMin[0]=.00;mAmp[0]=.00;
2715  switch(mlattice)
2716  {
2717  case TRICLINIC:
2718  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2719  mMin[2]=min_latt;mAmp[2]=max_latt-min_latt;
2720  mMin[3]=min_latt;mAmp[3]=max_latt-min_latt;
2721  mMin[4]=0;mAmp[4]=amp_crossp;
2722  mMin[5]=0;mAmp[5]=amp_crossp;
2723  mMin[6]=0;mAmp[6]=amp_crossp;
2724  mnpar=7;
2725  break;
2726  case MONOCLINIC:
2727  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2728  mMin[2]=min_latt;mAmp[2]=max_latt-min_latt;
2729  mMin[3]=min_latt;mAmp[3]=max_latt-min_latt;
2730  mMin[4]=0;mAmp[4]=amp_crossp;
2731  mnpar=5;
2732  break;
2733  case ORTHOROMBIC:
2734  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2735  mMin[2]=min_latt;mAmp[2]=max_latt-min_latt;
2736  mMin[3]=min_latt;mAmp[3]=max_latt-min_latt;
2737  mnpar=4;
2738  break;
2739  case HEXAGONAL:
2740  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2741  mMin[2]=min_latt;mAmp[2]=max_latt-min_latt;
2742  mnpar=3;
2743  break;
2744  case RHOMBOEDRAL:
2745  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2746  mMin[2]=-amp_crossp;mAmp[2]=2*amp_crossp;
2747  mnpar=3;
2748  break;
2749  case TETRAGONAL:
2750  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2751  mMin[2]=min_latt;mAmp[2]=max_latt-min_latt;
2752  mnpar=3;
2753  break;
2754  case CUBIC:
2755  mMin[1]=min_latt;mAmp[1]=max_latt-min_latt;
2756  mnpar=2;
2757  break;
2758  }
2759  //for(unsigned int k=0;k<mnpar;++k) cout<<"par["<<k<<"]: "<<mMin[k]<<"->"<<mMin[k]+mAmp[k]<<endl;
2760 
2761  unsigned int nb1=0,nb2=0;
2762  switch(mlattice)
2763  {
2764  case TRICLINIC:
2765  {
2766  nb1=3;
2767  nb2=3;
2768  break;
2769  }
2770  case MONOCLINIC:
2771  {
2772  nb1=3;
2773  nb2=1;
2774  break;
2775  }
2776  case ORTHOROMBIC:
2777  {
2778  nb1=3;
2779  break;
2780  }
2781  case HEXAGONAL:
2782  {
2783  nb1=2;
2784  break;
2785  }
2786  case RHOMBOEDRAL:
2787  {
2788  nb1=2;
2789  break;
2790  }
2791  case TETRAGONAL:
2792  {
2793  nb1=2;
2794  break;
2795  }
2796  case CUBIC:
2797  {
2798  nb1=1;
2799  break;
2800  }
2801  }
2802  this->ResetParList();
2803  {
2804  RefinablePar tmp("Zero",mRecUnitCell.par+0,-0.01,0.01,gpRefParTypeObjCryst,REFPAR_DERIV_STEP_ABSOLUTE,
2805  true,false,true,false);
2806  tmp.SetDerivStep(1e-4);
2807  this->AddPar(tmp);
2808  }
2809  char buf [50];
2810  string str="Reciprocal unit cell par #";
2811  for(unsigned int i=0;i<nb1;++i)
2812  {
2813  sprintf(buf,"%i",i);
2814  RefinablePar tmp(str+(string)buf,mRecUnitCell.par+i+1,
2815  0.01,1.,gpRefParTypeObjCryst,REFPAR_DERIV_STEP_ABSOLUTE,
2816  false,false,true,false);
2817  tmp.SetDerivStep(1e-4);
2818  this->AddPar(tmp);
2819  }
2820  for(unsigned int i=nb1;i<(nb1+nb2);++i)
2821  {
2822  sprintf(buf,"%i",i);
2823  RefinablePar tmp(str+(string)buf,mRecUnitCell.par+i+1,
2824  0.0,0.5,gpRefParTypeObjCryst,REFPAR_DERIV_STEP_ABSOLUTE,
2825  false,false,true,false);
2826  tmp.SetDerivStep(1e-4);
2827  this->AddPar(tmp);
2828  }
2829 }
2830 
2831 }//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:1564
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:937
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:60
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:109
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:105
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:117
ObjCryst::PeakList::GetPeakList
vector< hkl > & GetPeakList()
Get peak list.
Definition: Indexing.cpp:931
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:203
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:200
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:887
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)
Generate a list of simulated peak positions, from given lattice parameters.
Definition: Indexing.cpp:807
ObjCryst::PeakList::hkl0
One set of Miller indices, a possible indexation for a reflection.
Definition: Indexing.h:152
ObjCryst::PeakList::hkl
One observed diffraction line, to be indexed.
Definition: Indexing.h:161
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