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