ProSHADE  0.6.5 (NOV 2018)
Protein Shape Descriptors and Symmetry Detection
ProSHADE_SO3transform.cpp
Go to the documentation of this file.
1 
21 //============================================ FFTW3 + SOFT
22 #ifdef __cplusplus
23 extern "C" {
24 #endif
25 #include <fftw3.h>
26 #include <wrap_fftw.h>
27 #include <makeweights.h>
28 #include <s2_primitive.h>
29 #include <s2_cospmls.h>
30 #include <s2_legendreTransforms.h>
31 #include <s2_semi_fly.h>
32 #include <rotate_so3_utils.h>
33 #include <utils_so3.h>
34 #include <soft_fftw.h>
35 #include <rotate_so3_fftw.h>
36 #ifdef __cplusplus
37 }
38 #endif
39 
40 //============================================ RVAPI
41 #include <rvapi_interface.h>
42 
43 //============================================ ProSHADE
44 #include "ProSHADE.h"
45 #include "ProSHADE_internal.h"
46 #include "ProSHADE_misc.h"
47 
59 {
60  //======================================== Require E matrices
61  if ( !this->_trSigmaPreComputed )
62  {
63  std::cerr << "!!! ProSHADE ERROR !!! Finding the optimal overlay requires the data from pre-computation of the TrSigma descriptor. This is not really logical for the user, but it saves a lot of time in case both overlay and TrSigma are to be computed. Therefore, you need to call the precomputeTrSigmaDescriptor function. Sorry about the inconvenience." << std::endl;
64 
65  if ( settings->htmlReport )
66  {
67  std::stringstream hlpSS;
68  hlpSS << "<font color=\"red\">" << "Finding the optimal overlay requires the data from pre-computation of the integral over the spheres (E matrices), which has not been done. This looks like an internal bug, please report this case." << "</font>";
69  rvapi_set_text ( hlpSS.str().c_str(),
70  "ProgressSection",
71  settings->htmlReportLineProgress,
72  1,
73  1,
74  1 );
75  settings->htmlReportLineProgress += 1;
76  rvapi_flush ( );
77  }
78 
79  exit ( -1 );
80  }
81 
82  //======================================== Initialise internal values
83  this->_so3Coeffs = reinterpret_cast<fftw_complex*> ( fftw_malloc ( sizeof ( fftw_complex ) * ( ( 4 * pow( static_cast<double> ( this->_bandwidthLimit ),3.0 ) - static_cast<double> ( this->_bandwidthLimit ) ) / 3.0 ) ) );
84 
85  //======================================== Set all coeffs to 0
86  memset ( this->_so3Coeffs,
87  0.0,
88  sizeof ( fftw_complex ) * totalCoeffs_so3 ( this->_bandwidthLimit ) );
89 
90  //======================================== Initialise local variables
91  double wigNorm = 0.0;
92  unsigned int indexO = 0;
93  double signValue = 1.0;
94  bool toBeIgnored = false;
95  unsigned int bandUsageIndex = 0;
96 
97  //======================================== Combine the coefficients into SO(3) Array
98  for ( unsigned int bandIter = 0; bandIter < this->_bandwidthLimit; bandIter++ )
99  {
100  //==================================== Ignore odd l's as they are 0 anyway
101  if ( !this->_keepOrRemove ) { if ( ( bandIter % 2 ) != 0 ) { continue; } }
102 
103  //==================================== Ignore l's designated to be ignored
104  toBeIgnored = false;
105  for ( unsigned int igIt = 0; igIt < static_cast<unsigned int> ( this->_lsToIgnore.size() ); igIt++ ) { if ( bandIter == static_cast<unsigned int> ( this->_lsToIgnore.at(igIt) ) ) { toBeIgnored = true; } }
106  if ( toBeIgnored ) { continue; }
107 
108  //==================================== Get wigner norm factor
109  wigNorm = 2.0 * M_PI * sqrt ( 2.0 / (2.0 * bandIter + 1.0 ) ) ;
110 
111  //==================================== For each order
112  for ( unsigned int ordIter = 0; ordIter < static_cast<unsigned int> ( ( bandIter * 2 ) + 1 ); ordIter++ )
113  {
114  //================================ Set the sign
115  if ( ( ordIter-bandIter + bandIter ) % 2 ) { signValue = -1.0 ; }
116  else { signValue = 1.0 ; }
117 
118  //================================ For each order2
119  for ( unsigned int ordIter2 = 0; ordIter2 < static_cast<unsigned int> ( ( bandIter * 2 ) + 1 ); ordIter2++ )
120  {
121  //============================ Find output index
122  indexO = so3CoefLoc ( ordIter-bandIter,
123  ordIter2-bandIter,
124  bandIter,
125  this->_bandwidthLimit );
126 
127  //============================ Use E Matrix value (it is the integral over radii of c x c*, which is exactly what is needed here)
128  this->_so3Coeffs[indexO][0] = this->_trSigmaEMatrix[bandUsageIndex][ordIter][ordIter2][0] * wigNorm * signValue;
129  this->_so3Coeffs[indexO][1] = this->_trSigmaEMatrix[bandUsageIndex][ordIter][ordIter2][1] * wigNorm * signValue;
130 
131  //============================ Iterate sign value
132  signValue *= -1.0;
133  }
134  }
135  bandUsageIndex += 1;
136  }
137 
138  //======================================== Allocate memory for the inverse SO(3) transform
139  this->_so3InvCoeffs = reinterpret_cast<fftw_complex*> ( fftw_malloc ( sizeof ( fftw_complex ) * 8 * pow( static_cast<double> ( this->_bandwidthLimit ), 3.0 ) ) );
140  this->_so3Workspace1 = reinterpret_cast<fftw_complex*> ( fftw_malloc ( sizeof ( fftw_complex ) * 8 * pow( static_cast<double> ( this->_bandwidthLimit ), 3.0 ) ) );
141  this->_so3Workspace2 = reinterpret_cast<fftw_complex*> ( fftw_malloc ( sizeof ( fftw_complex ) * 14 * pow( static_cast<double> ( this->_bandwidthLimit ), 2.0 ) + 48 * this->_bandwidthLimit ) );
142  this->_so3Workspace3 = new double [static_cast<unsigned int> ( 24 * this->_bandwidthLimit + 2 * pow( static_cast<double> ( this->_bandwidthLimit ), 2.0 ) )];
143 
144  //======================================== Inverse SO(3) FFTW plan
145  fftw_plan inverseSO3;
146 
147  int howmany = 4 * this->_bandwidthLimit * this->_bandwidthLimit;
148  int idist = 2 * this->_bandwidthLimit;
149  int odist = 2 * this->_bandwidthLimit;
150  int rank = 2;
151 
152  int inembed[2], onembed[2];
153  inembed[0] = 2 * this->_bandwidthLimit;
154  inembed[1] = 4 * this->_bandwidthLimit * this->_bandwidthLimit;
155  onembed[0] = 2 * this->_bandwidthLimit;
156  onembed[1] = 4 * this->_bandwidthLimit * this->_bandwidthLimit;
157 
158  int istride = 1;
159  int ostride = 1;
160 
161  int na[2];
162  na[0] = 1;
163  na[1] = 2 * this->_bandwidthLimit;;
164 
165  inverseSO3 = fftw_plan_many_dft ( rank,
166  na,
167  howmany,
168  this->_so3Workspace1,
169  inembed,
170  istride,
171  idist,
172  this->_so3InvCoeffs,
173  onembed,
174  ostride,
175  odist,
176  FFTW_FORWARD,
177  FFTW_ESTIMATE );
178 
179  //======================================== Inverse SO(3) Transform computation
180  Inverse_SO3_Naive_fftw ( this->_bandwidthLimit,
181  this->_so3Coeffs,
182  this->_so3InvCoeffs,
183  this->_so3Workspace1,
184  this->_so3Workspace2,
185  this->_so3Workspace3,
186  &inverseSO3,
187  0 ) ;
188 
189  //======================================== Free memory
190  fftw_free ( this->_so3Coeffs );
191  fftw_free ( this->_so3Workspace1 );
192  fftw_free ( this->_so3Workspace2 );
193  delete this->_so3Workspace3;
194 
195  this->_so3Coeffs = nullptr;
196  this->_so3Workspace1 = nullptr;
197  this->_so3Workspace2 = nullptr;
198  this->_so3Workspace3 = nullptr;
199 
200  fftw_destroy_plan ( inverseSO3 );
201 
202  //============== Done
203  this->_so3InvMapComputed = true;
204 
205 }
206 
218 {
219  //======================================== Sanity check
220  if ( !this->_trSigmaPreComputed )
221  {
222  std::cerr << "!!! ProSHADE ERROR !!! Finding the optimal overlay requires the data from pre-computation of the TrSigma descriptor. This is not really logical for the user, but it saves a lot of time in case both overlay and TrSigma are to be computed. Therefore, you need to call the precomputeTrSigmaDescriptor function. Sorry about the inconvenience." << std::endl;
223 
224  if ( settings->htmlReport )
225  {
226  std::stringstream hlpSS;
227  hlpSS << "<font color=\"red\">" << "Finding the optimal overlay requires the data from pre-computation of the integral over the spheres (E matrices), which has not been done. This looks like an internal bug, please report this case." << "</font>";
228  rvapi_set_text ( hlpSS.str().c_str(),
229  "ProgressSection",
230  settings->htmlReportLineProgress,
231  1,
232  1,
233  1 );
234  settings->htmlReportLineProgress += 1;
235  rvapi_flush ( );
236  }
237 
238  exit ( -1 );
239  }
240 
241  //======================================== Initialise internal values
242  this->_so3InverseCoeffs = std::vector<fftw_complex*> ( this->all->size() );
243  std::vector<fftw_complex*> _so3Coeffs = std::vector<fftw_complex*> ( this->all->size() );
244 
245  //======================================== Initialise local variables
246  double wigNorm = 0.0;
247  unsigned int indexO = 0;
248  double signValue = 1.0;
249  bool toBeIgnored = false;
250  unsigned int bandUsageIndex = 0;
251  unsigned int fileUsageIndex = 0;
252  unsigned int localComparisonBandLim = 0;
253 
254  //======================================== Combine the coefficients into SO(3) Array
255  for ( unsigned int strIt = 0; strIt < static_cast<unsigned int> ( this->all->size() ); strIt++ )
256  {
257  //==================================== If not needed, do not compute
258  if ( this->_enLevelsDoNotFollow.at(strIt) == 1 ) { continue; }
259  if ( this->_trSigmaDoNotFollow.at(strIt) == 1 ) { fileUsageIndex += 1; continue; }
260 
261  //==================================== Find the appropriate bandwidth
262  localComparisonBandLim = std::min( this->one->_bandwidthLimit, this->all->at(strIt)->_bandwidthLimit );
263 
264  //==================================== Allocate memory for results
265  _so3Coeffs.at(strIt) = reinterpret_cast<fftw_complex*> ( fftw_malloc ( sizeof ( fftw_complex ) * ( ( 4 * pow( static_cast<double> ( localComparisonBandLim ), 3.0 ) - static_cast<double> ( localComparisonBandLim ) ) / 3.0 ) ) );
266 
267  //==================================== Set all coeffs to 0
268  memset ( _so3Coeffs.at(strIt),
269  0.0,
270  sizeof ( fftw_complex ) * totalCoeffs_so3 ( localComparisonBandLim ) );
271 
272  bandUsageIndex = 0;
273  for ( unsigned int bandIter = 0; bandIter < localComparisonBandLim; bandIter++ )
274  {
275  //================================ Ignore odd l's as they are 0 anyway
276  if ( !this->_keepOrRemove ) { if ( ( bandIter % 2 ) != 0 ) { continue; } }
277 
278  //================================ Ignore l's designated to be ignored
279  toBeIgnored = false;
280  for ( unsigned int igIt = 0; igIt < static_cast<unsigned int> ( this->_lsToIgnore.size() ); igIt++ ) { if ( bandIter == static_cast<unsigned int> ( this->_lsToIgnore.at(igIt) ) ) { toBeIgnored = true; } }
281  if ( toBeIgnored ) { continue; }
282 
283  //================================ Get wigner norm factor
284  wigNorm = 2.0 * M_PI * sqrt ( 2.0 / (2.0 * bandIter + 1.0 ) ) ;
285 
286  //================================ For each order
287  for ( unsigned int ordIter = 0; ordIter < static_cast<unsigned int> ( ( bandIter * 2 ) + 1 ); ordIter++ )
288  {
289  //============================ Set the sign
290  if ( ( ordIter-bandIter + bandIter ) % 2 ) { signValue = -1.0 ; }
291  else { signValue = 1.0 ; }
292 
293  //============================ For each order2
294  for ( unsigned int ordIter2 = 0; ordIter2 < static_cast<unsigned int> ( ( bandIter * 2 ) + 1 ); ordIter2++ )
295  {
296  //======================== Find output index
297  indexO = so3CoefLoc ( ordIter-bandIter, ordIter2-bandIter, bandIter, localComparisonBandLim );
298 
299  //======================== Use E Matrix value (it is the integral over radii of c x c*, which is exactly what is needed here)
300  _so3Coeffs.at(strIt)[indexO][0] = this->_EMatrices.at(fileUsageIndex).at(bandUsageIndex).at(ordIter).at(ordIter2)[0] * wigNorm * signValue;
301  _so3Coeffs.at(strIt)[indexO][1] = this->_EMatrices.at(fileUsageIndex).at(bandUsageIndex).at(ordIter).at(ordIter2)[1] * wigNorm * signValue;
302 
303  //======================== Iterate sign value
304  signValue *= -1.0;
305  }
306  }
307  bandUsageIndex += 1;
308  }
309  fileUsageIndex += 1;
310  }
311 
312  //======================================== Compute the inverse SO3 transform
313  for ( unsigned int strIt = 0; strIt < static_cast<unsigned int> ( this->all->size() ); strIt++ )
314  {
315  //==================================== If not needed, do not compute
316  if ( this->_enLevelsDoNotFollow.at(strIt) == 1 ) { continue; }
317  if ( this->_trSigmaDoNotFollow.at(strIt) == 1 ) { continue; }
318 
319  //==================================== Find the appropriate bandwidth
320  localComparisonBandLim = std::min( this->one->_bandwidthLimit, this->all->at(strIt)->_bandwidthLimit );
321 
322  //==================================== Allocate memory for the inverse SO(3) transform
323  fftw_complex *_so3Workspace1 = reinterpret_cast<fftw_complex*> ( fftw_malloc ( sizeof ( fftw_complex ) * 8 * pow( static_cast<double> ( localComparisonBandLim ), 3.0 ) ) );
324  fftw_complex *_so3Workspace2 = reinterpret_cast<fftw_complex*> ( fftw_malloc ( sizeof ( fftw_complex ) * 14 * pow( static_cast<double> ( localComparisonBandLim ), 2.0 ) + 48 * localComparisonBandLim ) );
325  double *_so3Workspace3 = new double [static_cast<unsigned int> ( 24 * localComparisonBandLim + 2 * pow( static_cast<double> ( localComparisonBandLim ), 2.0 ) )];
326 
327 
328  this->_so3InverseCoeffs.at(strIt) = reinterpret_cast<fftw_complex*> ( fftw_malloc ( sizeof ( fftw_complex ) * 8 * pow( static_cast<double> ( localComparisonBandLim ), 3.0 ) ) );
329 
330  //==================================== Inverse SO(3) FFTW plan
331  fftw_plan inverseSO3;
332 
333  int howmany = 4 * localComparisonBandLim * localComparisonBandLim;
334  int idist = 2 * localComparisonBandLim;
335  int odist = 2 * localComparisonBandLim;
336  int rank = 2;
337 
338  int inembed[2], onembed[2];
339  inembed[0] = 2 * localComparisonBandLim;
340  inembed[1] = 4 * localComparisonBandLim * localComparisonBandLim;
341  onembed[0] = 2 * localComparisonBandLim;
342  onembed[1] = 4 * localComparisonBandLim * localComparisonBandLim;
343 
344  int istride = 1;
345  int ostride = 1;
346 
347  int na[2];
348  na[0] = 1;
349  na[1] = 2 * localComparisonBandLim;
350 
351  inverseSO3 = fftw_plan_many_dft ( rank,
352  na,
353  howmany,
354  _so3Workspace1,
355  inembed,
356  istride,
357  idist,
358  this->_so3InverseCoeffs.at(strIt),
359  onembed,
360  ostride,
361  odist,
362  FFTW_FORWARD,
363  FFTW_ESTIMATE );
364 
365  //==================================== Inverse SO(3) Transform computation
366  Inverse_SO3_Naive_fftw ( localComparisonBandLim,
367  _so3Coeffs.at(strIt),
368  this->_so3InverseCoeffs.at(strIt),
369  _so3Workspace1,
370  _so3Workspace2,
371  _so3Workspace3,
372  &inverseSO3,
373  0 ) ;
374 
375  //==================================== Free memory
376  fftw_destroy_plan ( inverseSO3 );
377  fftw_free ( _so3Workspace1 );
378  fftw_free ( _so3Workspace2 );
379  fftw_free ( _so3Coeffs.at(strIt) );
380  delete[] _so3Workspace3;
381  }
382 
383  //======================================== Done
384  this->_so3InvMapComputed = true;
385 
386 }
387 
401 void ProSHADE_internal::ProSHADE_comparePairwise::setEulerAngles ( double alpha, double beta, double gamma )
402 {
403  //======================================== Set
404  this->_eulerAngles[0] = alpha;
405  this->_eulerAngles[1] = beta;
406  this->_eulerAngles[2] = gamma;
407 
408  //======================================== Angles set
409  this->_eulerAnglesFound = true;
410 
411  //======================================== Done
412  return;
413 
414 }
415 
430 {
431  //======================================== If done, do not repeat
432  if ( this->_eulerAnglesFound )
433  {
434  return ( this->_eulerAngles );
435  }
436 
437  //======================================== Sanity check
438  if ( !this->_so3InvMapComputed )
439  {
440  std::cerr << "!!! ProSHADE ERROR !!! Error in file compatison !!! Attempted to get Euler angles without pre-computing the required data (SO3 inverse transform). Call the getSO3InverseMap function BEFORE calling the getEulerAngles function." << std::endl;
441 
442  if ( settings->htmlReport )
443  {
444  std::stringstream hlpSS;
445  hlpSS << "<font color=\"red\">" << "Attempted to get Euler angles without pre-computing the required data (SO3 inverse transform). This looks like an internal bug, please report this case." << "</font>";
446  rvapi_set_text ( hlpSS.str().c_str(),
447  "ProgressSection",
448  settings->htmlReportLineProgress,
449  1,
450  1,
451  1 );
452  settings->htmlReportLineProgress += 1;
453  rvapi_flush ( );
454  }
455 
456  exit ( -1 );
457  }
458 
459  //======================================== Initialise global variable
460  this->_eulerAngles = std::array< double, 3 > { { 0.0, 0.0, 0.0 } };
461 
462  //======================================== Find the maximum value and determine the Euler angles
463  unsigned int maxLoc = 0;
464  double maxVal = 0.0;
465  double tmpVal = 0.0;
466 
467  for ( unsigned int iter = 0; iter < ( 8 * pow( static_cast<unsigned int> ( this->_bandwidthLimit ), 3 ) ); iter++ )
468  {
469  tmpVal = pow( this->_so3InvCoeffs[iter][0], 2.0 ) + pow( this->_so3InvCoeffs[iter][1], 2.0 );
470  if ( maxVal < tmpVal )
471  {
472  maxVal = tmpVal;
473  maxLoc = iter;
474  }
475  }
476 
477  //======================================== Save max peak height if required
478  *correlation = maxVal;
479 
480  //======================================== Determine the Euler angles
481  int ii, jj, kk, tmp;
482  ii = std::floor ( maxLoc / ( 4.0 * this->_bandwidthLimit * this->_bandwidthLimit ) );
483  tmp = maxLoc - ( ii * ( 4.0 * this->_bandwidthLimit * this->_bandwidthLimit ) );
484  jj = std::floor ( tmp / ( 2.0 * this->_bandwidthLimit ) );
485  kk = maxLoc - ( ii * ( 4.0 * this->_bandwidthLimit * this->_bandwidthLimit ) ) -
486  jj * ( 2.0 * this->_bandwidthLimit );
487 
488  this->_eulerAngles[2] = ( M_PI * jj / ( static_cast<double> ( this->_bandwidthLimit ) ) );
489  this->_eulerAngles[1] = ( M_PI * ( 2.0 * ii + 1.0 ) / ( 4.0 * this->_bandwidthLimit ) ) ;
490  this->_eulerAngles[0] = ( M_PI * kk / ( static_cast<double> ( this->_bandwidthLimit ) ) );
491 
492  //==================================== Free memory
493  fftw_free ( this->_so3InvCoeffs );
494  this->_so3InvCoeffs = nullptr;
495 
496  //======================================== Done
497  this->_eulerAnglesFound = true;
498 
499  //======================================== Return
500  return ( this->_eulerAngles );
501 
502 }
503 
517 {
518  //======================================== If done, do not repeat
519  if ( this->_eulerAnglesFound )
520  {
521  return ( this->_eulerAngles );
522  }
523 
524  //======================================== Sanity check
525  if ( !this->_so3InvMapComputed )
526  {
527  std::cerr << "!!! ProSHADE ERROR !!! Error in file compatison !!! Attempted to get Euler angles without pre-computing the required data (SO3 inverse transform). Call the getSO3InverseMap function BEFORE calling the getEulerAngles function." << std::endl;
528 
529  if ( settings->htmlReport )
530  {
531  std::stringstream hlpSS;
532  hlpSS << "<font color=\"red\">" << "Attempted to get Euler angles without pre-computing the required data (SO3 inverse transform). This looks like an internal bug, please report this case." << "</font>";
533  rvapi_set_text ( hlpSS.str().c_str(),
534  "ProgressSection",
535  settings->htmlReportLineProgress,
536  1,
537  1,
538  1 );
539  settings->htmlReportLineProgress += 1;
540  rvapi_flush ( );
541  }
542 
543  exit ( -1 );
544  }
545 
546  //======================================== Initialise global variable
547  this->_eulerAngles = std::vector< std::array<double,3> > ( this->all->size() );
548 
549  //======================================== Find the maximum value and determine the Euler angles
550  unsigned int maxLoc = 0;
551  double maxVal = 0.0;
552  double tmpVal = 0.0;
553  unsigned int localComparisonBandLim = 0;
554 
555  for ( unsigned int strIt = 0; strIt < static_cast<unsigned int> ( this->all->size() ); strIt++ )
556  {
557  //==================================== If not needed, do not compute
558  if ( this->_enLevelsDoNotFollow.at(strIt) == 1 ) { continue; }
559  if ( this->_trSigmaDoNotFollow.at(strIt) == 1 ) { continue; }
560 
561  //==================================== Find maximum peak
562  maxLoc = 0;
563  maxVal = 0.0;
564  tmpVal = 0.0;
565 
566  //==================================== Find the appropriate bandwidth
567  localComparisonBandLim = std::min ( this->one->_bandwidthLimit, this->all->at(strIt)->_bandwidthLimit );
568 
569  for ( unsigned int iter = 0; iter < ( 8 * pow( static_cast<unsigned int> ( localComparisonBandLim ), 3 ) ); iter++ )
570  {
571  tmpVal = pow( this->_so3InverseCoeffs.at(strIt)[iter][0], 2.0 ) + pow( this->_so3InverseCoeffs.at(strIt)[iter][1], 2.0 );
572  if ( maxVal < tmpVal )
573  {
574  maxVal = tmpVal;
575  maxLoc = iter;
576  }
577  }
578 
579  //==================================== Determine the Euler angles
580  int ii, jj, kk, tmp;
581  ii = std::floor ( maxLoc / ( 4.0 * localComparisonBandLim * localComparisonBandLim ) );
582  tmp = maxLoc - ( ii * ( 4.0 * localComparisonBandLim * localComparisonBandLim ) );
583  jj = std::floor ( tmp / ( 2.0 * localComparisonBandLim ) );
584  kk = maxLoc - ( ii * ( 4.0 * localComparisonBandLim * localComparisonBandLim ) ) -
585  jj * ( 2.0 * localComparisonBandLim );
586 
587  this->_eulerAngles.at(strIt)[2] = ( M_PI * jj / ( static_cast<double> ( localComparisonBandLim ) ) );
588  this->_eulerAngles.at(strIt)[1] = ( M_PI * ( 2.0 * ii + 1.0 ) / ( 4.0 * localComparisonBandLim ) ) ;
589  this->_eulerAngles.at(strIt)[0] = ( M_PI * kk / ( static_cast<double> ( localComparisonBandLim ) ) );
590 
591  //==================================== Free memory
592  fftw_free ( this->_so3InverseCoeffs.at(strIt) );
593  this->_so3InverseCoeffs.at(strIt) = nullptr;
594  }
595 
596  //======================================== Done
597  this->_eulerAnglesFound = true;
598 
599  //======================================== Return
600  return ( this->_eulerAngles );
601 
602 }
603 
622  double noIQRs,
623  bool freeMem,
624  int peakSize,
625  double realDist,
626  int verbose )
627 {
628  //======================================== Sanity check
629  if ( !this->_so3InvMapComputed )
630  {
631  std::cerr << "!!! ProSHADE ERROR !!! Error in file compatison !!! Attempted to get Euler angles without pre-computing the required data (SO3 inverse transform). Call the getSO3InverseMap function BEFORE calling the getSO3Peaks function." << std::endl;
632 
633  if ( settings->htmlReport )
634  {
635  std::stringstream hlpSS;
636  hlpSS << "<font color=\"red\">" << "Attempted to get Euler angles without pre-computing the required data (SO3 inverse transform). This looks like an internal bug, please report this case." << "</font>";
637  rvapi_set_text ( hlpSS.str().c_str(),
638  "ProgressSection",
639  settings->htmlReportLineProgress,
640  1,
641  1,
642  1 );
643  settings->htmlReportLineProgress += 1;
644  rvapi_flush ( );
645  }
646 
647  exit ( -1 );
648  }
649 
650  if ( this->_so3InvCoeffs == nullptr )
651  {
652  std::cerr << "!!! ProSHADE ERROR !!! Error in function getSO3Peaks. The pointer to the SO3 data is empty. Common cause is when you call this function more than once, or when you call the getEulerAngles function before this one. If either is the case, please add \"false\" as an additional parameter to all previous calls of both functions." << std::endl;
653 
654  if ( settings->htmlReport )
655  {
656  std::stringstream hlpSS;
657  hlpSS << "<font color=\"red\">" << "The pointer to the SO3 data is empty. This looks like an internal bug, please report this case." << "</font>";
658  rvapi_set_text ( hlpSS.str().c_str(),
659  "ProgressSection",
660  settings->htmlReportLineProgress,
661  1,
662  1,
663  1 );
664  settings->htmlReportLineProgress += 1;
665  rvapi_flush ( );
666  }
667 
668  exit ( -1 );
669  }
670 
671  //======================================== Initialise local variables
672  int x = 0;
673  int y = 0;
674  int z = 0;
675  int newIter = 0;
676  int xDim = 4 * pow ( this->_bandwidthLimit, 2 );
677  int yDim = 2 * this->_bandwidthLimit;
678 
679  std::array<double,8> hlpArr;
680  std::vector< std::array<double,8> > peakList;
681  std::array<double,4> hlpInterp;
682  std::vector< std::array<double,4> > peakInterpVals;
683  std::vector< std::vector< std::array<double,4> > > peakInterpValsList;
684  std::vector< std::vector<double> > tMat;
685  std::vector< std::vector<double> > avgMat;
686  std::vector< std::vector< std::vector<double> > > uAndVMat;
687  std::vector< std::array<double,4> > rotAxes;
688  std::array<double,4> curAxis;
689  std::vector<double> hlpVec;
690  std::vector< std::array<double,8> > ret;
691  std::vector<double> peakHeights;
692 
693  double tmpVal = 0.0;
694  double angAlpha = 0.0;
695  double angBeta = 0.0;
696  double angGamma = 0.0;
697  double singularityErrTolerance = 0.05;
698  double avgMatWeight = 0.0;
699  double rotVal = 0.0;
700  double median = 0.0;
701  double interQuartileRange = 0.0;
702 
703  int peakX = 0;
704  int peakY = 0;
705  int peakZ = 0;
706 
707  unsigned int medianPos = 0;
708 
709  bool breakPeak = false;
710 
711  //======================================== Find peak heights
712  for ( unsigned int iter = 0; iter < ( 8 * pow( static_cast<unsigned int> ( this->_bandwidthLimit ), 3 ) ); iter++ )
713  {
714  //==================================== Find point value
715  tmpVal = pow( this->_so3InvCoeffs[iter][0], 2.0 ) + pow( this->_so3InvCoeffs[iter][1], 2.0 );
716 
717  //==================================== Find the x, y and z
718  x = std::floor ( iter / xDim );
719  y = std::floor ( ( iter - x * xDim ) / yDim );
720  z = iter - x * xDim - y * yDim;
721 
722  //==================================== Get the X surrounding values as well as the peak
723  breakPeak = false;
724  for ( int xCh = -peakSize; xCh <= +peakSize; xCh++ )
725  {
726  if ( breakPeak ) { break; }
727  for ( int yCh = -peakSize; yCh <= +peakSize; yCh++ )
728  {
729  if ( breakPeak ) { break; }
730  for ( int zCh = -peakSize; zCh <= +peakSize; zCh++ )
731  {
732  if ( breakPeak ) { break; }
733  if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) { continue; }
734 
735  peakX = x + xCh; if ( peakX >= yDim ) { peakX -= yDim; }; if ( peakX < 0 ) { peakX += yDim; }
736  peakY = y + yCh; if ( peakY >= yDim ) { peakY -= yDim; }; if ( peakY < 0 ) { peakY += yDim; }
737  peakZ = z + zCh; if ( peakZ >= yDim ) { peakZ -= yDim; }; if ( peakZ < 0 ) { peakZ += yDim; }
738  newIter = peakX * xDim + peakY * yDim + peakZ;
739 
740  if ( (pow( this->_so3InvCoeffs[newIter][0], 2.0 ) + pow( this->_so3InvCoeffs[newIter][1], 2.0 )) > tmpVal ) { breakPeak = true; break; }
741  }
742  }
743  }
744  if ( breakPeak ) { continue; }
745 
746  //==================================== Find Euler angles
747  angGamma = ( M_PI * y / ( static_cast<double> ( this->_bandwidthLimit ) ) );
748  angBeta = ( M_PI * ( 2.0 * x + 1.0 ) / ( 4.0 * this->_bandwidthLimit ) ) ;
749  angAlpha = ( M_PI * z / ( static_cast<double> ( this->_bandwidthLimit ) ) );
750 
751  //==================================== Convert ZXZ Euler angle ranges
752  if ( angAlpha > M_PI ) { angAlpha = angAlpha - M_PI * 2.0; }
753  if ( angBeta > M_PI ) { angBeta = angBeta - M_PI * 2.0; }
754  if ( angGamma > M_PI ) { angGamma = angGamma - M_PI * 2.0; }
755 
756  //==================================== Remove 0 angle peaks
757  ProSHADE_internal_misc::getAxisAngleFromEuler ( angAlpha, angBeta, angGamma, &median, &median, &median, &rotVal, true );
758  if ( ( rotVal < ( 0.0 + singularityErrTolerance ) ) && ( rotVal > ( 0.0 - singularityErrTolerance ) ) ) { continue; }
759 
760  //==================================== Save the peak heights to determine threshold
761  peakHeights.emplace_back ( tmpVal );
762  }
763 
764  //======================================== Find the threshold for peaks
765  medianPos = static_cast<unsigned int> ( peakHeights.size() / 2 );
766  std::sort ( peakHeights.begin(), peakHeights.end() );
767  median = peakHeights.at(medianPos);
768 
769  interQuartileRange = peakHeights.at(medianPos+(medianPos/2)) - peakHeights.at(medianPos-(medianPos/2));
770  this->_peakHeightThr = peakHeights.at(medianPos+(medianPos/2)) + ( interQuartileRange * noIQRs );
771 
772  //======================================== Find peaks
773  for ( unsigned int iter = 0; iter < ( 8 * pow( static_cast<unsigned int> ( this->_bandwidthLimit ), 3 ) ); iter++ )
774  {
775  //==================================== Find point value
776  tmpVal = pow( this->_so3InvCoeffs[iter][0], 2.0 ) + pow( this->_so3InvCoeffs[iter][1], 2.0 );
777 
778  //==================================== If below threshold, go to next
779  if ( tmpVal < this->_peakHeightThr ) { continue; }
780 
781  //==================================== Find the x, y and z
782  x = std::floor ( iter / xDim );
783  y = std::floor ( ( iter - x * xDim ) / yDim );
784  z = iter - x * xDim - y * yDim;
785 
786  //==================================== Initialise interpolation values
787  peakInterpVals.clear ( );
788 
789  //==================================== Get the X surrounding values as well as the peak
790  breakPeak = false;
791  for ( int xCh = -peakSize; xCh <= +peakSize; xCh++ )
792  {
793  if ( breakPeak ) { break; }
794  for ( int yCh = -peakSize; yCh <= +peakSize; yCh++ )
795  {
796  if ( breakPeak ) { break; }
797  for ( int zCh = -peakSize; zCh <= +peakSize; zCh++ )
798  {
799  if ( breakPeak ) { break; }
800  if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) { continue; }
801 
802  hlpArr[0] = x + xCh; if ( hlpArr[0] >= yDim ) { hlpArr[0] -= yDim; }; if ( hlpArr[0] < 0 ) { hlpArr[0] += yDim; }
803  hlpArr[1] = y + yCh; if ( hlpArr[1] >= yDim ) { hlpArr[1] -= yDim; }; if ( hlpArr[1] < 0 ) { hlpArr[1] += yDim; }
804  hlpArr[2] = z + zCh; if ( hlpArr[2] >= yDim ) { hlpArr[2] -= yDim; }; if ( hlpArr[2] < 0 ) { hlpArr[2] += yDim; }
805  newIter = hlpArr[0] * xDim + hlpArr[1] * yDim + hlpArr[2];
806  hlpArr[3] = pow( this->_so3InvCoeffs[newIter][0], 2.0 ) + pow( this->_so3InvCoeffs[newIter][1], 2.0 );
807  if ( hlpArr[3] > tmpVal ) { breakPeak = true; break; }
808  hlpInterp[0] = hlpArr[0]; hlpInterp[1] = hlpArr[1]; hlpInterp[2] = hlpArr[2]; hlpInterp[3] = hlpArr[3];
809  peakInterpVals.emplace_back ( hlpInterp );
810  }
811  }
812  }
813  if ( breakPeak ) { continue; }
814 
815  //==================================== This value is surrounded by smaller values! Check if it is not a 0 angle
816  angGamma = ( M_PI * y / ( static_cast<double> ( this->_bandwidthLimit ) ) );
817  angBeta = ( M_PI * ( 2.0 * x + 1.0 ) / ( 4.0 * this->_bandwidthLimit ) ) ;
818  angAlpha = ( M_PI * z / ( static_cast<double> ( this->_bandwidthLimit ) ) );
819 
820  ProSHADE_internal_misc::getAxisAngleFromEuler ( angAlpha, angBeta, angGamma, &rotVal, &rotVal, &rotVal, &rotVal, true );
821  if ( ( rotVal < ( 0.0 + singularityErrTolerance ) ) && ( rotVal > ( 0.0 - singularityErrTolerance ) ) ) { continue; }
822 
823  //==================================== OK, now check if this peak actually leads to something similar
824  if ( this->checkPeakExistence ( angAlpha, angBeta, angGamma, settings, realDist ) )
825  {
826  hlpArr[0] = x + 0; if ( hlpArr[0] >= yDim ) { hlpArr[0] -= yDim; }; if ( hlpArr[0] < 0 ) { hlpArr[0] += yDim; }
827  hlpArr[1] = y + 0; if ( hlpArr[1] >= yDim ) { hlpArr[1] -= yDim; }; if ( hlpArr[1] < 0 ) { hlpArr[1] += yDim; }
828  hlpArr[2] = z + 0; if ( hlpArr[2] >= yDim ) { hlpArr[2] -= yDim; }; if ( hlpArr[2] < 0 ) { hlpArr[2] += yDim; }
829  newIter = hlpArr[0] * xDim + hlpArr[1] * yDim + hlpArr[2];
830  hlpArr[3] = pow( this->_so3InvCoeffs[newIter][0], 2.0 ) + pow( this->_so3InvCoeffs[newIter][1], 2.0 );
831  peakList.emplace_back ( hlpArr );
832  hlpInterp[0] = hlpArr[0];
833  hlpInterp[1] = hlpArr[1];
834  hlpInterp[2] = hlpArr[2];
835  hlpInterp[3] = hlpArr[3];
836  peakInterpVals.emplace_back ( hlpInterp );
837  peakInterpValsList.emplace_back ( peakInterpVals );
838  }
839  }
840  if ( verbose > 3 )
841  {
842  std::cout << ">>>>> " << peakList.size() << " peaks detected." << std::endl;
843  }
844 
845  //======================================== Free memory
846  if ( freeMem )
847  {
848  fftw_free ( this->_so3InvCoeffs );
849  this->_so3InvCoeffs = nullptr;
850  }
851 
852  //======================================== Interpolate angles
853  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( peakList.size() ); iter++ )
854  {
855  //==================================== Init
856  avgMat.clear ();
857  hlpVec.clear ();
858  hlpVec.emplace_back ( 0.0 );
859  hlpVec.emplace_back ( 0.0 );
860  hlpVec.emplace_back ( 0.0 );
861  avgMat.emplace_back ( hlpVec );
862  avgMat.emplace_back ( hlpVec );
863  avgMat.emplace_back ( hlpVec );
864 
865  avgMatWeight = 0.0;
866 
867  //==================================== For each surrounding value (including the peak value)
868  for ( unsigned int valIter = 0; valIter < static_cast<unsigned int> ( peakInterpValsList.at(iter).size() ); valIter++ )
869  {
870  //================================ Get angles
871  angGamma = ( M_PI * peakInterpValsList.at(iter).at(valIter)[1] / ( static_cast<double> ( this->_bandwidthLimit ) ) );
872  angBeta = ( M_PI * ( 2.0 * peakInterpValsList.at(iter).at(valIter)[0] + 1.0 ) / ( 4.0 * this->_bandwidthLimit ) ) ;
873  angAlpha = ( M_PI * peakInterpValsList.at(iter).at(valIter)[2] / ( static_cast<double> ( this->_bandwidthLimit ) ) );
874 
875  //================================ Convert ZXZ Euler angle ranges
876  if ( angAlpha > M_PI ) { angAlpha = angAlpha - M_PI * 2.0; }
877  if ( angBeta > M_PI ) { angBeta = angBeta - M_PI * 2.0; }
878  if ( angGamma > M_PI ) { angGamma = angGamma - M_PI * 2.0; }
879 
880  //================================ Init
881  tMat.clear ( );
882 
883  //================================ Get rotation matrices
884  hlpVec.clear ( );
885  hlpVec.emplace_back ( cos ( angAlpha ) * cos ( angBeta ) * cos ( angGamma ) - sin ( angAlpha ) * sin ( angGamma ) );
886  hlpVec.emplace_back ( sin ( angAlpha ) * cos ( angBeta ) * cos ( angGamma ) + cos ( angAlpha ) * sin ( angGamma ) );
887  hlpVec.emplace_back ( -sin ( angBeta ) * cos ( angGamma ) );
888  tMat.emplace_back ( hlpVec );
889 
890  hlpVec.clear ( );
891  hlpVec.emplace_back ( -cos ( angAlpha ) * cos ( angBeta ) * sin ( angGamma ) - sin ( angAlpha ) * cos ( angGamma ) );
892  hlpVec.emplace_back ( -sin ( angAlpha ) * cos ( angBeta ) * sin ( angGamma ) + cos ( angAlpha ) * cos ( angGamma ) );
893  hlpVec.emplace_back ( sin ( angBeta ) * sin ( angGamma ) );
894  tMat.emplace_back ( hlpVec );
895 
896  hlpVec.clear ( );
897  hlpVec.emplace_back ( cos ( angAlpha ) * sin ( angBeta ) );
898  hlpVec.emplace_back ( sin ( angAlpha ) * sin ( angBeta ) );
899  hlpVec.emplace_back ( cos ( angBeta ) );
900  tMat.emplace_back ( hlpVec );
901 
902  //================================ Add to average matrix
903  avgMat.at(0).at(0) += tMat.at(0).at(0) * peakInterpValsList.at(iter).at(valIter)[3];
904  avgMat.at(0).at(1) += tMat.at(0).at(1) * peakInterpValsList.at(iter).at(valIter)[3];
905  avgMat.at(0).at(2) += tMat.at(0).at(2) * peakInterpValsList.at(iter).at(valIter)[3];
906 
907  avgMat.at(1).at(0) += tMat.at(1).at(0) * peakInterpValsList.at(iter).at(valIter)[3];
908  avgMat.at(1).at(1) += tMat.at(1).at(1) * peakInterpValsList.at(iter).at(valIter)[3];
909  avgMat.at(1).at(2) += tMat.at(1).at(2) * peakInterpValsList.at(iter).at(valIter)[3];
910 
911  avgMat.at(2).at(0) += tMat.at(2).at(0) * peakInterpValsList.at(iter).at(valIter)[3];
912  avgMat.at(2).at(1) += tMat.at(2).at(1) * peakInterpValsList.at(iter).at(valIter)[3];
913  avgMat.at(2).at(2) += tMat.at(2).at(2) * peakInterpValsList.at(iter).at(valIter)[3];
914 
915  //================================ and to average weight
916  avgMatWeight += peakInterpValsList.at(iter).at(valIter)[3];
917  }
918 
919  //==================================== Normalise average mats by weights
920  avgMat.at(0).at(0) /= avgMatWeight;
921  avgMat.at(0).at(1) /= avgMatWeight;
922  avgMat.at(0).at(2) /= avgMatWeight;
923 
924  avgMat.at(1).at(0) /= avgMatWeight;
925  avgMat.at(1).at(1) /= avgMatWeight;
926  avgMat.at(1).at(2) /= avgMatWeight;
927 
928  avgMat.at(2).at(0) /= avgMatWeight;
929  avgMat.at(2).at(1) /= avgMatWeight;
930  avgMat.at(2).at(2) /= avgMatWeight;
931 
932  //==================================== SVD average mat
933  uAndVMat = this->getSingularValuesUandVMatrix ( avgMat, 3, settings );
934 
935  //==================================== Init matrix for multiplication
936  tMat.clear ( );
937  hlpVec.clear ( );
938  hlpVec.emplace_back ( 0.0 );
939  hlpVec.emplace_back ( 0.0 );
940  hlpVec.emplace_back ( 0.0 );
941  tMat.emplace_back ( hlpVec ); tMat.emplace_back ( hlpVec ); tMat.emplace_back ( hlpVec );
942 
943  //==================================== Compute U * V^T
944  for ( unsigned int row = 0; row < 3; row++ )
945  {
946  for ( unsigned int col = 0; col < 3; col++ )
947  {
948  for ( unsigned int inner = 0; inner < 3; inner++ )
949  {
950  tMat.at(row).at(col) += uAndVMat.at(0).at(inner).at(row) * uAndVMat.at(1).at(col).at(inner);
951  }
952  }
953  }
954 
955  //==================================== Find alpha, beta and gamma from matrix
956  //=! Deal with beta angle first
957  if ( ( angBeta >= 0.0 ) ) { angBeta = acos ( tMat.at(2).at(2) ); }
958  else { angBeta = -acos ( tMat.at(2).at(2) ); }
959 
960 
961  //=! Now alpha
962  if ( angAlpha > ( M_PI / 2.0 ) ) { angAlpha = M_PI - asin ( tMat.at(2).at(1) / sin ( angBeta ) ); }
963  else if ( angAlpha < ( -M_PI / 2.0 ) ) { angAlpha = -M_PI - asin ( tMat.at(2).at(1) / sin ( angBeta ) ); }
964  else { angAlpha = asin ( tMat.at(2).at(1) / sin ( angBeta ) ); }
965 
966  //=! Now gamma
967  if ( angGamma > ( M_PI / 2.0 ) ) { angGamma = M_PI - asin ( tMat.at(1).at(2) / sin ( angBeta ) ); }
968  else if ( angGamma < ( -M_PI / 2.0 ) ) { angGamma = -M_PI - asin ( tMat.at(1).at(2) / sin ( angBeta ) ); }
969  else { angGamma = asin ( tMat.at(1).at(2) / sin ( angBeta ) ); }
970 
971  //==================================== Now, change the peak angles to interpolated angles
972  peakList.at(iter)[0] = angAlpha;
973  peakList.at(iter)[1] = angBeta;
974  peakList.at(iter)[2] = angGamma;
975  }
976  if ( verbose > 3 )
977  {
978  std::cout << ">>>>> Peaks angles interpolated." << std::endl;
979  }
980 
981  //======================================== Get axis-angle representation values
982  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( peakList.size() ); iter++ )
983  {
984  //==================================== Get angles
985  angGamma = peakList.at(iter)[2];
986  angBeta = peakList.at(iter)[1];
987  angAlpha = peakList.at(iter)[0];
988 
989  //==================================== Find the axis-angle representation axis
990  ProSHADE_internal_misc::getAxisAngleFromEuler ( angAlpha, angBeta, angGamma, &curAxis[0], &curAxis[1], &curAxis[2], &curAxis[3] );
991 
992  //==================================== Save axis-angle representation results
993  rotAxes.emplace_back ( curAxis );
994  }
995 
996  //======================================== Save
997  for ( unsigned int rot = 0; rot < static_cast<unsigned int> ( peakList.size() ); rot++ )
998  {
999  hlpArr[0] = peakList.at(rot)[0];
1000  hlpArr[1] = peakList.at(rot)[1];
1001  hlpArr[2] = peakList.at(rot)[2];
1002  hlpArr[3] = peakList.at(rot)[3];
1003  hlpArr[4] = rotAxes.at(rot)[0];
1004  hlpArr[5] = rotAxes.at(rot)[1];
1005  hlpArr[6] = rotAxes.at(rot)[2];
1006  hlpArr[7] = rotAxes.at(rot)[3];
1007  ret.emplace_back ( hlpArr );
1008  }
1009  if ( verbose > 3 )
1010  {
1011  std::cout << ">>>>> Peaks converted to angle-axis." << std::endl;
1012  }
1013 
1014  //======================================== Sort (descending order)
1015  std::sort ( ret.begin(), ret.end(), [](const std::array<double,8>& a, const std::array<double,8>& b) { return a[3] > b[3]; });
1016 
1017  //======================================== Return
1018  return ( ret );
1019 
1020 }
1021 
1037 double ProSHADE_internal::ProSHADE_comparePairwise::maxPeakHeightFromAngleAxis ( double X,
1038  double Y,
1039  double Z,
1040  double Angle,
1041  ProSHADE::ProSHADE_settings* settings )
1042 {
1043  //======================================== Sanity check
1044  if ( this->_so3InvCoeffs == nullptr )
1045  {
1046  std::cerr << "!!! ProSHADE ERROR !!! Error in function checkPeakFromAngleAxis. The pointer to the SO3 data is empty. Common cause is when you call this function more than once, or when you call the getEulerAngles function before this one. If either is the case, please add \"false\" as an additional parameter to all previous calls of both functions." << std::endl;
1047 
1048  if ( settings->htmlReport )
1049  {
1050  std::stringstream hlpSS;
1051  hlpSS << "<font color=\"red\">" << "The pointer to the SO3 data is empty. This looks like an internal bug, please report this case." << "</font>";
1052  rvapi_set_text ( hlpSS.str().c_str(),
1053  "ProgressSection",
1054  settings->htmlReportLineProgress,
1055  1,
1056  1,
1057  1 );
1058  settings->htmlReportLineProgress += 1;
1059  rvapi_flush ( );
1060  }
1061 
1062  exit ( -1 );
1063  }
1064 
1065  //======================================== Initialise
1066  double peakHeight = 0.0;
1067  double x = 0.0;
1068  double y = 0.0;
1069  double z = 0.0;
1070  double angAlpha = 0.0;
1071  double angBeta = 0.0;
1072  double angGamma = 0.0;
1073  double axX = 0.0;
1074  double axY = 0.0;
1075  double axZ = 0.0;
1076  double angle = 0.0;
1077 
1078  int xDim = 4 * pow ( this->_bandwidthLimit, 2 );
1079  int yDim = 2 * this->_bandwidthLimit;
1080 
1081  double ret = 0.0;
1082 
1083  //======================================== Iterate through SO3 Inv data and check all coordinates which give the correct Angle-Axis representation
1084  for ( unsigned int iter = 0; iter < ( 8 * pow( static_cast<unsigned int> ( this->_bandwidthLimit ), 3 ) ); iter++ )
1085  {
1086  //==================================== Find point value
1087  peakHeight = pow( this->_so3InvCoeffs[iter][0], 2.0 ) + pow( this->_so3InvCoeffs[iter][1], 2.0 );
1088 
1089  //==================================== Quick check
1090  if ( peakHeight <= ret ) { continue; }
1091 
1092  //==================================== Find the x, y and z
1093  x = std::floor ( iter / xDim );
1094  y = std::floor ( ( iter - x * xDim ) / yDim );
1095  z = iter - x * xDim - y * yDim;
1096 
1097  //==================================== Find Euler angles
1098  angGamma = ( M_PI * y / ( static_cast<double> ( this->_bandwidthLimit ) ) );
1099  angBeta = ( M_PI * ( 2.0 * x + 1.0 ) / ( 4.0 * this->_bandwidthLimit ) ) ;
1100  angAlpha = ( M_PI * z / ( static_cast<double> ( this->_bandwidthLimit ) ) );
1101 
1102  //==================================== Convert ZXZ Euler angle ranges
1103  if ( angAlpha > M_PI ) { angAlpha = angAlpha - M_PI * 2.0; }
1104  if ( angBeta > M_PI ) { angBeta = angBeta - M_PI * 2.0; }
1105  if ( angGamma > M_PI ) { angGamma = angGamma - M_PI * 2.0; }
1106 
1107  //==================================== Convert to Angle-Axis representation and check compatibility to tested values
1108  ProSHADE_internal_misc::getAxisAngleFromEuler ( angAlpha, angBeta, angGamma, &axX, &axY, &axZ, &angle );
1109  if ( !( ( angle < ( Angle + 0.1 ) ) && ( ( angle > ( Angle - 0.1 ) ) ) ) )
1110  {
1111  continue;
1112  }
1113 
1114  if ( !( ( axX < ( X + 0.1 ) ) && ( ( axX > ( X - 0.1 ) ) ) ) ||
1115  !( ( axY < ( Y + 0.1 ) ) && ( ( axY > ( Y - 0.1 ) ) ) ) ||
1116  !( ( axZ < ( Z + 0.1 ) ) && ( ( axZ > ( Z - 0.1 ) ) ) ) )
1117  {
1118  continue;
1119  }
1120 
1121  //==================================== Save the peak height
1122  if ( ret < peakHeight ) { ret = peakHeight; }
1123  }
1124 
1125  //======================================== Done
1126  return ( ret );
1127 
1128 }
1129 
1144  double Y,
1145  double Z,
1146  double Angle,
1147  ProSHADE::ProSHADE_settings* settings )
1148 {
1149  //======================================== Sanity check
1150  if ( this->_so3InvCoeffs == nullptr )
1151  {
1152  std::cerr << "!!! ProSHADE ERROR !!! Error in function maxAvgPeakForSymmetry. The pointer to the SO3 data is empty. Common cause is when you call this function more than once, or when you call the getEulerAngles function before this one. If either is the case, please add \"false\" as an additional parameter to all previous calls of both functions." << std::endl;
1153 
1154  if ( settings->htmlReport )
1155  {
1156  std::stringstream hlpSS;
1157  hlpSS << "<font color=\"red\">" << "The pointer to the SO3 data is empty. This looks like an internal bug, please report this case." << "</font>";
1158  rvapi_set_text ( hlpSS.str().c_str(),
1159  "ProgressSection",
1160  settings->htmlReportLineProgress,
1161  1,
1162  1,
1163  1 );
1164  settings->htmlReportLineProgress += 1;
1165  rvapi_flush ( );
1166  }
1167 
1168  exit ( -1 );
1169  }
1170 
1171  //======================================== Initialise
1172  double x = 0.0;
1173  double y = 0.0;
1174  double z = 0.0;
1175  double angAlpha = 0.0;
1176  double angBeta = 0.0;
1177  double angGamma = 0.0;
1178  double axX = 0.0;
1179  double axY = 0.0;
1180  double axZ = 0.0;
1181  double angle = 0.0;
1182  int xDim = 4 * pow ( this->_bandwidthLimit, 2 );
1183  int yDim = 2 * this->_bandwidthLimit;
1184  double ret = 0.0;
1185  std::vector< std::array<double,2> > peaks;
1186  std::array<double,2> hlpArr;
1187  bool xPos = true; if ( X < 0 ) { xPos = false; }
1188  bool yPos = true; if ( Y < 0 ) { yPos = false; }
1189  bool zPos = true; if ( Z < 0 ) { zPos = false; }
1190 
1191  //======================================== Iterate through SO3 Inv data and check all coordinates which give the correct Angle-Axis representation
1192  for ( unsigned int iter = 0; iter < ( 8 * pow( static_cast<unsigned int> ( this->_bandwidthLimit ), 3 ) ); iter++ )
1193  {
1194  //==================================== Find the x, y and z
1195  x = std::floor ( iter / xDim );
1196  y = std::floor ( ( iter - x * xDim ) / yDim );
1197  z = iter - x * xDim - y * yDim;
1198 
1199  //==================================== Find Euler angles
1200  angGamma = ( M_PI * y / ( static_cast<double> ( this->_bandwidthLimit ) ) );
1201  angBeta = ( M_PI * ( 2.0 * x + 1.0 ) / ( 4.0 * this->_bandwidthLimit ) ) ;
1202  angAlpha = ( M_PI * z / ( static_cast<double> ( this->_bandwidthLimit ) ) );
1203 
1204  //==================================== Convert to Angle-Axis representation and check compatibility to tested values
1205  ProSHADE_internal_misc::getAxisAngleFromEuler ( angAlpha, angBeta, angGamma, &axX, &axY, &axZ, &angle );
1206  if ( !( ( std::abs( axX ) < ( std::abs( X ) + 0.1 ) ) && ( ( std::abs( axX ) > ( std::abs( X ) - 0.1 ) ) ) ) ||
1207  !( ( std::abs( axY ) < ( std::abs( Y ) + 0.1 ) ) && ( ( std::abs( axY ) > ( std::abs( Y ) - 0.1 ) ) ) ) ||
1208  !( ( std::abs( axZ ) < ( std::abs( Z ) + 0.1 ) ) && ( ( std::abs( axZ ) > ( std::abs( Z ) - 0.1 ) ) ) ) )
1209  {
1210  continue;
1211  }
1212  else
1213  {
1214  if ( ( ( xPos && ( axX > 0.0 ) ) || ( !xPos && ( axX < 0.0 ) ) ) &&
1215  ( ( yPos && ( axY > 0.0 ) ) || ( !yPos && ( axY < 0.0 ) ) ) &&
1216  ( ( zPos && ( axZ > 0.0 ) ) || ( !zPos && ( axZ < 0.0 ) ) ) )
1217  {
1218  //============================ Same signs
1219  ;
1220  }
1221  else
1222  {
1223  if ( ( ( xPos && ( axX < 0.0 ) ) || ( !xPos && ( axX > 0.0 ) ) ) &&
1224  ( ( yPos && ( axY < 0.0 ) ) || ( !yPos && ( axY > 0.0 ) ) ) &&
1225  ( ( zPos && ( axZ < 0.0 ) ) || ( !zPos && ( axZ > 0.0 ) ) ) )
1226  {
1227  //======================== All opposite signs
1228  angle *= -1.0;
1229  }
1230  else
1231  {
1232  continue;
1233  }
1234  }
1235  }
1236 
1237  //==================================== Save vals
1238  hlpArr[0] = angle + M_PI;
1239  hlpArr[1] = pow( this->_so3InvCoeffs[iter][0], 2.0 ) + pow( this->_so3InvCoeffs[iter][1], 2.0 );
1240  peaks.emplace_back ( hlpArr );
1241  }
1242 
1243  //======================================== Sort the vector
1244  std::sort ( peaks.begin(), peaks.end(), [](const std::array<double,2>& a, const std::array<double,2>& b) { return a[0] < b[0]; } );
1245 
1246  //======================================== Find the best X peaks with correct distances
1247  double angStep = 0.1;
1248  std::vector<double> sums ( std::floor ( (2.0*M_PI/angStep) / Angle ) );
1249  double curSum = 0.0;
1250  double maxVal = 0.0;
1251  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( std::floor ( (2.0*M_PI/angStep) / Angle ) ); iter++ )
1252  {
1253  curSum = 0.0;
1254  for ( unsigned int angCmb = 0; angCmb < static_cast<unsigned int> ( Angle ); angCmb++ )
1255  {
1256  maxVal = 0.0;
1257  for ( unsigned int angIt = 0; angIt < static_cast<unsigned int> ( peaks.size() ); angIt++ )
1258  {
1259  if ( peaks.at(angIt)[0] < ( ( iter*angStep ) + ( ( 2.0 * M_PI / Angle ) * angCmb ) ) ) { continue; }
1260  if ( peaks.at(angIt)[0] > ( ( (iter+1)*angStep ) + ( ( 2.0 * M_PI / Angle ) * angCmb ) ) ) { break; }
1261 
1262  if ( peaks.at(angIt)[1] > maxVal ) { maxVal = peaks.at(angIt)[1]; }
1263  }
1264  curSum += maxVal;
1265  }
1266  curSum /= Angle;
1267  if ( ret < curSum ) { ret = curSum; }
1268  }
1269 
1270  //======================================== Done
1271  return ( ret );
1272 
1273 }
1274 
1292 std::vector< std::array<double,5> > ProSHADE_internal::ProSHADE_comparePairwise::findCnSymmetry ( std::vector< std::array<double,8> > peaks,
1293  ProSHADE::ProSHADE_settings* settings,
1294  double axisErrorTolerance,
1295  bool freeMem,
1296  double percAllowedToMiss,
1297  int verbose )
1298 {
1299  //======================================== Sanity check
1300  if ( peaks.size() < 1 )
1301  {
1302  std::cerr << "!!! ProSHADE ERROR !!! No peaks in the data - cannot search for symmetry. Either loosen up the criteria for the getSO3Peaks function, or conclude no symmetry." << std::endl;
1303 
1304  if ( settings->htmlReport )
1305  {
1306  std::stringstream hlpSS;
1307  hlpSS << "<font color=\"red\">" << "Found no peaks in the rotation function map. This could signify no symmetry being present, or very low sampling of the map. Please try increasing resolution and bandwidth if you believe symmetry should be present." << "</font>";
1308  rvapi_set_text ( hlpSS.str().c_str(),
1309  "ProgressSection",
1310  settings->htmlReportLineProgress,
1311  1,
1312  1,
1313  1 );
1314  settings->htmlReportLineProgress += 1;
1315  rvapi_flush ( );
1316  }
1317 
1318  exit ( -1 );
1319  }
1320 
1321  //======================================== Initialise local variables
1322  std::vector< std::array<double,5> > ret;
1323  std::vector< std::array<double,7> > grpPeaks;
1324  std::vector< std::array<double,7> > additionalPeaks;
1325  std::vector<double> angVals;
1326  std::array<double,7> hlpArr;
1327  std::array<double,5> retArr;
1328  bool sameAxisFound = false;
1329  bool groupContinuous = false;
1330  bool oneDifferent = false;
1331  bool matchedThisPeak = false;
1332  bool additionalPeakFound = false;
1333  double errTolerance = 0.0;
1334  double angTolerance = 0.0;
1335  double angDivisionRemainder = 0.0;
1336  double angDivisionBasis = 0.0;
1337  double groupAngle = 0.0;
1338  double meanX = 0.0;
1339  double meanY = 0.0;
1340  double meanZ = 0.0;
1341  double totVals = 0.0;
1342  double meanPeak = 0.0;
1343  double nextPeakError = ( M_PI * 2.0 ) / ( static_cast<double> ( this->_bandwidthLimit ) * 2.0 );
1344  double nextSymmetryError = 0.0;
1345  double peakVal = 0.0;
1346  double minAngDist = 0.0;
1347  double currentAverage = 0.0;
1348  int peaksNeededMinimum = 0;
1349  int peaksNeededRemaining = 0;
1350  int peaksNeededFound = 0;
1351  std::vector<double> groupAngleHits;
1352  std::vector< std::array<double,2> > groupMembers;
1353  std::array<double,2> hlpArr2;
1354  std::vector<double> missingPeaks;
1355  std::vector<double> angsToTry;
1356  std::vector<unsigned int> thisAxisSyms;
1357  std::vector<unsigned int> delVecElems;
1358 
1359  std::vector< std::vector<unsigned int> > uniqueAxisGroups;
1360 
1361  if ( axisErrorTolerance == 0.0 ) { errTolerance = 0.1; }
1362  else { errTolerance = axisErrorTolerance; }
1363 
1364  //======================================== Find unique axes
1365  for ( unsigned int peakIter = 0; peakIter < static_cast<unsigned int> ( peaks.size() ); peakIter++ )
1366  {
1367  //==================================== If axis is 0 ; 0 ; 0 then ignore, this will be dealt with later
1368  if ( ( peaks.at(peakIter)[4] == 0.0 ) && ( peaks.at(peakIter)[5] == 0.0 ) && ( peaks.at(peakIter)[6] == 0.0 ) ) { continue; }
1369 
1370  //==================================== Initialise variables for next peak
1371  sameAxisFound = false;
1372 
1373  //==================================== Give same direction to maximum vector elements
1374  if ( ( std::abs(peaks.at(peakIter)[4]) > std::abs(peaks.at(peakIter)[5]) ) && ( std::abs(peaks.at(peakIter)[4]) > std::abs(peaks.at(peakIter)[6]) ) )
1375  {
1376  if ( peaks.at(peakIter)[4] < 0 )
1377  {
1378  peaks.at(peakIter)[4] *= -1.0;
1379  peaks.at(peakIter)[5] *= -1.0;
1380  peaks.at(peakIter)[6] *= -1.0;
1381  peaks.at(peakIter)[7] *= -1.0;
1382  }
1383  }
1384  else if ( ( std::abs(peaks.at(peakIter)[5]) > std::abs(peaks.at(peakIter)[4]) ) && ( std::abs(peaks.at(peakIter)[5]) > std::abs(peaks.at(peakIter)[6]) ) )
1385  {
1386  if ( peaks.at(peakIter)[5] < 0 )
1387  {
1388  peaks.at(peakIter)[4] *= -1.0;
1389  peaks.at(peakIter)[5] *= -1.0;
1390  peaks.at(peakIter)[6] *= -1.0;
1391  peaks.at(peakIter)[7] *= -1.0;
1392  }
1393  }
1394  else if ( ( std::abs(peaks.at(peakIter)[6]) > std::abs(peaks.at(peakIter)[4]) ) && ( std::abs(peaks.at(peakIter)[6]) > std::abs(peaks.at(peakIter)[5]) ) )
1395  {
1396  if ( peaks.at(peakIter)[6] < 0 )
1397  {
1398  peaks.at(peakIter)[4] *= -1.0;
1399  peaks.at(peakIter)[5] *= -1.0;
1400  peaks.at(peakIter)[6] *= -1.0;
1401  peaks.at(peakIter)[7] *= -1.0;
1402  }
1403  }
1404 
1405  //==================================== If angle is 0, ignore
1406  if ( ( peaks.at(peakIter)[7] < ( 0.0 + errTolerance ) ) && ( peaks.at(peakIter)[7] > ( 0.0 - errTolerance ) ) )
1407  {
1408  continue;
1409  }
1410 
1411  //==================================== Compare to all already present axes
1412  for ( unsigned int sameAxisGrp = 0; sameAxisGrp < static_cast<unsigned int> ( uniqueAxisGroups.size() ); sameAxisGrp++ )
1413  {
1414  //====== Try matching
1415  for ( unsigned int sameAxis = 0; sameAxis < static_cast<unsigned int> ( uniqueAxisGroups.at(sameAxisGrp).size() ); sameAxis++ )
1416  {
1417  //== Is there identical axis?
1418  if ( ( ( (peaks.at(uniqueAxisGroups.at(sameAxisGrp).at(sameAxis))[4]-errTolerance) <= peaks.at(peakIter)[4] ) &&
1419  ( (peaks.at(uniqueAxisGroups.at(sameAxisGrp).at(sameAxis))[4]+errTolerance) > peaks.at(peakIter)[4] ) ) &&
1420  ( ( (peaks.at(uniqueAxisGroups.at(sameAxisGrp).at(sameAxis))[5]-errTolerance) <= peaks.at(peakIter)[5] ) &&
1421  ( (peaks.at(uniqueAxisGroups.at(sameAxisGrp).at(sameAxis))[5]+errTolerance) > peaks.at(peakIter)[5] ) ) &&
1422  ( ( (peaks.at(uniqueAxisGroups.at(sameAxisGrp).at(sameAxis))[6]-errTolerance) <= peaks.at(peakIter)[6] ) &&
1423  ( (peaks.at(uniqueAxisGroups.at(sameAxisGrp).at(sameAxis))[6]+errTolerance) > peaks.at(peakIter)[6] ) ) )
1424  {
1425  sameAxisFound = true;
1426  uniqueAxisGroups.at(sameAxisGrp).emplace_back ( peakIter );
1427  break;
1428  }
1429  }
1430  }
1431 
1432  //==================================== If same axis was found, do nothing
1433  if ( sameAxisFound ) { continue; }
1434 
1435  //==================================== No similar axis was found
1436  std::vector<unsigned int> hlpVec;
1437  hlpVec.emplace_back ( peakIter );
1438  uniqueAxisGroups.emplace_back ( hlpVec );
1439  }
1440 
1441  //======================================== Now deal with the axis = 0 ; 0 ; 0
1442  for ( unsigned int peakIter = 0; peakIter < static_cast<unsigned int> ( peaks.size() ); peakIter++ )
1443  {
1444  //==================================== If angle is 0, ignore
1445  if ( ( peaks.at(peakIter)[7] < ( 0.0 + errTolerance ) ) && ( peaks.at(peakIter)[7] > ( 0.0 - errTolerance ) ) )
1446  {
1447  continue;
1448  }
1449 
1450  //==================================== Now, deal with the axis = 0 ; 0 ; 0
1451  if ( ( peaks.at(peakIter)[4] == 0.0 ) && ( peaks.at(peakIter)[5] == 0.0 ) && ( peaks.at(peakIter)[6] == 0.0 ) )
1452  {
1453  for ( unsigned int sameAxisGrp = 0; sameAxisGrp < static_cast<unsigned int> ( uniqueAxisGroups.size() ); sameAxisGrp++ )
1454  {
1455  uniqueAxisGroups.at(sameAxisGrp).emplace_back ( peakIter );
1456  }
1457  }
1458  }
1459 
1460  //======================================== Check each unique axis group for pairs of identical angle-axis representation with different Euler angles and similar peak height.
1461  for ( unsigned int axisGroup = 0; axisGroup < static_cast<unsigned int> ( uniqueAxisGroups.size() ); axisGroup++ )
1462  {
1463  for ( unsigned int ex1 = 0; ex1 < static_cast<unsigned int> ( uniqueAxisGroups.at(axisGroup).size() ); ex1++ )
1464  {
1465  for ( unsigned int ex2 = 1; ex2 < static_cast<unsigned int> ( uniqueAxisGroups.at(axisGroup).size() ); ex2++ )
1466  {
1467  //============================ Keep only unique pairs
1468  if ( ex1 >= ex2 ) { continue; }
1469 
1470  //============================ Check if the pair has the same angle representation (it has to have the same axis)
1471  if ( ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[7] < ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[7] + errTolerance ) ) &&
1472  ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[7] > ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[7] - errTolerance ) ) )
1473  {
1474  //======================== Same angles. Check if different Euler angles
1475  oneDifferent = false;
1476  if ( ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[0] < ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[0] + errTolerance ) ) &&
1477  ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[0] > ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[0] - errTolerance ) ) )
1478  {
1479  ;
1480  }
1481  else
1482  {
1483  oneDifferent = true;
1484  }
1485 
1486  if ( ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[1] < ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[1] + errTolerance ) ) &&
1487  ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[1] > ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[1] - errTolerance ) ) )
1488  {
1489  ;
1490  }
1491  else
1492  {
1493  oneDifferent = true;
1494  }
1495 
1496  if ( ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[2] < ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[2] + errTolerance ) ) &&
1497  ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[2] > ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[2] - errTolerance ) ) )
1498  {
1499  ;
1500  }
1501  else
1502  {
1503  oneDifferent = true;
1504  }
1505 
1506  if ( !oneDifferent )
1507  {
1508  continue;
1509  }
1510  }
1511  }
1512  }
1513  }
1514  if ( verbose > 3 )
1515  {
1516  std::cout << ">>>>> Peaks clustered by symmetry axis." << std::endl;
1517  }
1518 
1519  //======================================== For each group of peaks with the same axis
1520  for ( unsigned int axisGroup = 0; axisGroup < static_cast<unsigned int> ( uniqueAxisGroups.size() ); axisGroup++ )
1521  {
1522  //==================================== Re-organise the peak data
1523  thisAxisSyms.clear ( );
1524  additionalPeaks.clear ( );
1525 
1526  ROUGH_RESTART:
1527 
1528  additionalPeakFound = false;
1529  grpPeaks.clear ( );
1530  for ( unsigned int grpIter = 0; grpIter < static_cast<unsigned int> ( uniqueAxisGroups.at(axisGroup).size() ); grpIter++ )
1531  {
1532  hlpArr[0] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpIter))[4];
1533  hlpArr[1] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpIter))[5];
1534  hlpArr[2] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpIter))[6];
1535  hlpArr[3] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpIter))[7];
1536  hlpArr[4] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpIter))[3];
1537  hlpArr[5] = static_cast<double> ( grpIter );
1538  hlpArr[6] = 0.0;
1539 
1540  grpPeaks.emplace_back ( hlpArr );
1541  }
1542 
1543  if ( additionalPeaks.size () > 0 )
1544  {
1545  for ( unsigned int adIt = 0; adIt < static_cast<unsigned int> ( additionalPeaks.size() ); adIt++ )
1546  {
1547  grpPeaks.emplace_back ( additionalPeaks.at(adIt) );
1548  }
1549  }
1550 
1551  if ( verbose > 3 )
1552  {
1553  printf ( ">>>>> Symmetry axis group:\n" );
1554  printf ( " Group\tPeak\t x\t y\t z\t Angle\t Peak\n" );
1555  printf ( " Index\tIndex\t elem.\t elem.\t elem.\t (rad)\theight\n" );
1556  for ( int i = 0; i < static_cast<int> ( grpPeaks.size() ); i++ )
1557  {
1558  std::sort ( grpPeaks.begin(), grpPeaks.end(), [](const std::array<double,7>& a, const std::array<double,7>& b) { return a[4] > b[4]; } );
1559 
1560  printf ( " %d\t %d\t%+.3f\t%+.3f\t%+.3f\t%+.3f\t%+.3f\n", axisGroup, i, grpPeaks.at(i)[0], grpPeaks.at(i)[1], grpPeaks.at(i)[2], grpPeaks.at(i)[3], grpPeaks.at(i)[4] );
1561  }
1562  }
1563 
1564  while ( true )
1565  {
1566  //================================ Init
1567  angsToTry.clear ( );
1568 
1569  //================================ Find the peak with highest value
1570  std::sort ( grpPeaks.begin(), grpPeaks.end(), [](const std::array<double,7>& a, const std::array<double,7>& b) { return a[4] > b[4]; });
1571 
1572  //================================ If no peak is over err, exit loop
1573  if ( grpPeaks.at(0)[4] < ( 0.0 + this->_peakHeightThr ) )
1574  {
1575  break;
1576  }
1577 
1578  //================================ Find smallest distance to all other peaks
1579  minAngDist = 999.9;
1580  for ( unsigned int pAD = 1; pAD < static_cast<unsigned int> ( grpPeaks.size() ); pAD++ )
1581  {
1582  if ( ( std::abs( std::abs ( grpPeaks.at(0)[3] ) - std::abs ( grpPeaks.at(pAD)[3] ) ) < minAngDist ) && std::abs( std::abs ( grpPeaks.at(0)[3] ) - std::abs ( grpPeaks.at(pAD)[3] ) ) > 0.1 )
1583  {
1584  minAngDist = std::abs( std::abs ( grpPeaks.at(0)[3] ) - std::abs ( grpPeaks.at(pAD)[3] ) );
1585  }
1586  }
1587 
1588  //================================ Should this symmetry be tested?
1589  angDivisionRemainder = std::modf ( static_cast<double> ( (2.0*M_PI) / std::abs ( minAngDist ) ), &angDivisionBasis );
1590  if ( angDivisionRemainder > 0.8 )
1591  {
1592  angDivisionRemainder -= 1.0;
1593  angDivisionBasis += 1.0;
1594  }
1595  angDivisionRemainder /= angDivisionBasis;
1596  nextSymmetryError = ( M_PI * 2.0 / angDivisionBasis ) - ( M_PI * 2.0 / ( angDivisionBasis + 1.0 ) );
1597  angTolerance = ( nextPeakError / nextSymmetryError );
1598 
1599  if ( ( angDivisionRemainder < ( 0.0 + angTolerance ) ) && ( angDivisionRemainder > ( 0.0 - angTolerance ) ) )
1600  {
1601  if ( angTolerance < 0.5 )
1602  {
1603  angsToTry.emplace_back ( angDivisionBasis );
1604  }
1605  else
1606  {
1607  angsToTry.emplace_back ( angDivisionBasis - 1.0 );
1608  angsToTry.emplace_back ( angDivisionBasis );
1609  angsToTry.emplace_back ( angDivisionBasis + 1.0 );
1610  }
1611  }
1612 
1613  //================================ Is there symmetry with the highest peak angle?
1614  angDivisionRemainder = std::modf ( static_cast<double> ( (2.0*M_PI) / std::abs ( grpPeaks.at(0)[3] ) ), &angDivisionBasis );
1615 
1616  if ( angDivisionRemainder > 0.8 )
1617  {
1618  angDivisionRemainder -= 1.0;
1619  angDivisionBasis += 1.0;
1620  }
1621  angDivisionRemainder /= angDivisionBasis;
1622 
1623  //================================ Find appropriate error tolerance
1624  nextSymmetryError = ( M_PI * 2.0 / angDivisionBasis ) - ( M_PI * 2.0 / ( angDivisionBasis + 1.0 ) );
1625  angTolerance = ( nextPeakError / nextSymmetryError );
1626 
1627  if ( ( angDivisionRemainder < ( 0.0 + angTolerance ) ) && ( angDivisionRemainder > ( 0.0 - angTolerance ) ) )
1628  {
1629  if ( angTolerance < 0.5 )
1630  {
1631  angsToTry.emplace_back ( angDivisionBasis );
1632  }
1633  else
1634  {
1635  angsToTry.emplace_back ( angDivisionBasis - 1.0 );
1636  angsToTry.emplace_back ( angDivisionBasis );
1637  angsToTry.emplace_back ( angDivisionBasis + 1.0 );
1638  }
1639  }
1640 
1641  //================================ Use only unique values
1642  std::sort ( angsToTry.begin(), angsToTry.end(), std::greater<double>() );
1643  angsToTry.erase ( std::unique( angsToTry.begin(), angsToTry.end() ), angsToTry.end() );
1644 
1645  //================================ Check if the symmetry has not already been found for this axis
1646  delVecElems.clear();
1647  for ( unsigned int tAS = 0; tAS < static_cast<unsigned int> ( thisAxisSyms.size() ); tAS++ )
1648  {
1649  for ( unsigned int nTS = 0; nTS < static_cast<unsigned int> ( angsToTry.size() ); nTS++ )
1650  {
1651  if ( thisAxisSyms.at(tAS) == static_cast<unsigned int> ( angsToTry.at(nTS) ) )
1652  {
1653  delVecElems.emplace_back ( nTS );
1654  }
1655  }
1656  }
1657 
1658  std::sort ( delVecElems.begin(), delVecElems.end(), std::greater<double>() );
1659  for ( unsigned int remEl = 0; remEl < static_cast<unsigned int> ( delVecElems.size() ); remEl++ )
1660  {
1661  angsToTry.erase ( angsToTry.begin() + delVecElems.at(remEl) );
1662  }
1663 
1664  if ( verbose > 3 )
1665  {
1666  printf ( ">>>>>>>>> Testing the following symmetries: " );
1667  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( angsToTry.size() ); iter++ )
1668  {
1669  printf ( "C%.0f\t", angsToTry.at(iter) );
1670  }
1671  printf ( "\n" );
1672  }
1673 
1674  if ( angsToTry.size () == 0 )
1675  {
1676  grpPeaks.at(0)[4] = 0.0;
1677  continue;
1678  }
1679 
1680  for ( unsigned int angToTry = 0; angToTry < static_cast<unsigned int> ( angsToTry.size() ); angToTry++ )
1681  {
1682  //============================ Yes! Determine other peaks expected positions
1683  groupAngle = ( 2.0 * M_PI ) / angsToTry.at(angToTry);
1684  angVals.clear ( );
1685 
1686  for ( int iter = static_cast<int> ( -(angsToTry.at(angToTry)/2.0 + 1.0) ); iter <= static_cast<int> ( angsToTry.at(angToTry)/2.0 + 1.0 ); iter++ )
1687  {
1688  angVals.emplace_back ( static_cast<double> ( iter * groupAngle ) );
1689  }
1690 
1691  //============================ Now check for existence of these peaks
1692  groupAngleHits.clear ( );
1693  groupMembers.clear ( );
1694  missingPeaks.clear ( );
1695  for ( unsigned int grpAngIt = 0; grpAngIt < static_cast<unsigned int> ( angVals.size() ); grpAngIt++ )
1696  {
1697  matchedThisPeak = false;
1698  for ( unsigned int peakIt = 0; peakIt < static_cast<unsigned int> ( grpPeaks.size() ); peakIt++ )
1699  {
1700  if ( angVals.at(grpAngIt) == 0.0 )
1701  {
1702  groupAngleHits.emplace_back ( angVals.at(grpAngIt) );
1703  hlpArr2[0] = static_cast<double> ( peakIt );
1704  hlpArr2[1] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpPeaks.at(peakIt)[5]))[3];
1705  groupMembers.emplace_back ( hlpArr2 );
1706  matchedThisPeak = true;
1707  break;
1708  }
1709 
1710  if ( ( angVals.at(grpAngIt) < ( grpPeaks.at(peakIt)[3] + errTolerance ) ) &&
1711  ( angVals.at(grpAngIt) > ( grpPeaks.at(peakIt)[3] - errTolerance ) ) )
1712  {
1713  groupAngleHits.emplace_back ( angVals.at(grpAngIt) );
1714  hlpArr2[0] = static_cast<double> ( peakIt );
1715  hlpArr2[1] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpPeaks.at(peakIt)[5]))[3];
1716  groupMembers.emplace_back ( hlpArr2 );
1717  matchedThisPeak = true;
1718  break;
1719  }
1720  }
1721 
1722  if ( !matchedThisPeak )
1723  {
1724  missingPeaks.emplace_back ( angVals.at(grpAngIt) );
1725  }
1726  }
1727 
1728  //============================ And finally check for continuity
1729  groupContinuous = false;
1730  if ( groupAngleHits.size () > 1 )
1731  {
1732  groupContinuous = true;
1733  for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( groupAngleHits.size () ); iter++ )
1734  {
1735  if ( ( ( groupAngleHits.at(iter-1) + groupAngle ) < ( groupAngleHits.at(iter) + errTolerance ) ) &&
1736  ( ( groupAngleHits.at(iter-1) + groupAngle ) > ( groupAngleHits.at(iter) - errTolerance ) ) )
1737  {
1738  ;
1739  }
1740  else
1741  {
1742  groupContinuous = false;
1743  break;
1744  }
1745  }
1746  }
1747 
1748  //============================ Could it be that there are a few peaks missing (not highest around, but high enough) and so the symmetry fails?
1749  if ( !groupContinuous || ( static_cast<int> ( groupAngleHits.size() ) < static_cast<int> ( angsToTry.at(angToTry) ) ) )
1750  {
1751  if ( verbose > 3 )
1752  {
1753  printf ( ">>>>>>>>>>>>>> Checking for missing peaks...\n" );
1754  }
1755 
1756  //======================== Is it few or many?
1757  if ( ( ( static_cast<double> ( angsToTry.at(angToTry) ) - static_cast<double> ( groupAngleHits.size() ) ) / static_cast<double> ( angsToTry.at(angToTry) ) ) < percAllowedToMiss )
1758  {
1759  //==================== Few. Initialise the values
1760  peaksNeededMinimum = static_cast<int> ( angsToTry.at(angToTry) ) - static_cast<int> ( groupAngleHits.size() );
1761  peaksNeededFound = 0;
1762  peaksNeededRemaining = static_cast<int> ( missingPeaks.size() );
1763  currentAverage = 0.0;
1764  for ( unsigned int grpHlpIt = 0; grpHlpIt < static_cast<unsigned int> ( grpPeaks.size() ); grpHlpIt++ )
1765  {
1766  currentAverage += grpPeaks.at(grpHlpIt)[4];
1767  }
1768  currentAverage /= static_cast<double> ( grpPeaks.size() );
1769  currentAverage = std::max ( currentAverage - this->_peakHeightThr, this->_peakHeightThr );
1770  if ( angsToTry.at(0) < 9 )
1771  {
1772  currentAverage = std::max ( this->_peakHeightThr, currentAverage / 4.0 );
1773  }
1774 
1775  //==================== Now, for each missing, check if it exists
1776  for ( unsigned int misPeak = 0; misPeak < static_cast<unsigned int> ( missingPeaks.size() ); misPeak++ )
1777  {
1778  if ( peaksNeededMinimum > ( peaksNeededRemaining + peaksNeededFound ) ) { break; }
1779 
1780  if ( missingPeaks.at(misPeak) > M_PI ) { peaksNeededRemaining -= 1; continue; }
1781  if ( missingPeaks.at(misPeak) < -M_PI ) { peaksNeededRemaining -= 1; continue; }
1782  if ( missingPeaks.at(misPeak) == 0.0 ) { peaksNeededRemaining -= 1; continue; }
1783 
1784  double newX = 0.0;
1785  double newY = 0.0;
1786  double newZ = 0.0;
1787  bool pExists = false;
1788 
1789  for ( unsigned int grpHlpIt = 0; grpHlpIt < static_cast<unsigned int> ( grpPeaks.size() ); grpHlpIt++ )
1790  {
1791  newX += grpPeaks.at(grpHlpIt)[0];
1792  newY += grpPeaks.at(grpHlpIt)[1];
1793  newZ += grpPeaks.at(grpHlpIt)[2];
1794  }
1795  newX /= static_cast<double> ( grpPeaks.size() );
1796  newY /= static_cast<double> ( grpPeaks.size() );
1797  newZ /= static_cast<double> ( grpPeaks.size() );
1798 
1799  if ( missingPeaks.at(misPeak) > 0.0 )
1800  {
1801  peakVal = this->maxPeakHeightFromAngleAxis ( newX, newY, newZ, missingPeaks.at(misPeak), settings );
1802  if ( peakVal >= currentAverage )
1803  {
1804  peaksNeededFound += 1;
1805 
1806  hlpArr[0] = newX;
1807  hlpArr[1] = newY;
1808  hlpArr[2] = newZ;
1809  hlpArr[3] = missingPeaks.at(misPeak);
1810  hlpArr[4] = 0.0;
1811  hlpArr[5] = 0.0;
1812  hlpArr[6] = 1.0;
1813 
1814  additionalPeaks.emplace_back ( hlpArr );
1815  additionalPeakFound = true;
1816  pExists = true;
1817  }
1818 
1819  if ( !pExists )
1820  {
1821  peakVal = this->maxPeakHeightFromAngleAxis ( -newX, -newY, -newZ, -missingPeaks.at(misPeak), settings );
1822  if ( peakVal >= currentAverage )
1823  {
1824  peaksNeededFound += 1;
1825 
1826  hlpArr[0] = newX;
1827  hlpArr[1] = newY;
1828  hlpArr[2] = newZ;
1829  hlpArr[3] = missingPeaks.at(misPeak);
1830  hlpArr[4] = 0.0;
1831  hlpArr[5] = 0.0;
1832  hlpArr[6] = 1.0;
1833 
1834  additionalPeaks.emplace_back ( hlpArr );
1835  additionalPeakFound = true;
1836  pExists = true;
1837  }
1838  }
1839  }
1840  else
1841  {
1842  peakVal = this->maxPeakHeightFromAngleAxis ( -newX, -newY, -newZ, -missingPeaks.at(misPeak), settings );
1843  if ( peakVal >= currentAverage )
1844  {
1845  peaksNeededFound += 1;
1846 
1847  hlpArr[0] = newX;
1848  hlpArr[1] = newY;
1849  hlpArr[2] = newZ;
1850  hlpArr[3] = missingPeaks.at(misPeak);
1851  hlpArr[4] = 0.0;
1852  hlpArr[5] = 0.0;
1853  hlpArr[6] = 1.0;
1854 
1855  additionalPeaks.emplace_back ( hlpArr );
1856  additionalPeakFound = true;
1857  pExists = true;
1858  }
1859 
1860  if ( !pExists )
1861  {
1862  peakVal = this->maxPeakHeightFromAngleAxis ( newX, newY, newZ, missingPeaks.at(misPeak), settings );
1863  if ( peakVal >= currentAverage )
1864  {
1865  peaksNeededFound += 1;
1866 
1867  hlpArr[0] = newX;
1868  hlpArr[1] = newY;
1869  hlpArr[2] = newZ;
1870  hlpArr[3] = missingPeaks.at(misPeak);
1871  hlpArr[4] = 0.0;
1872  hlpArr[5] = 0.0;
1873  hlpArr[6] = 1.0;
1874 
1875  additionalPeaks.emplace_back ( hlpArr );
1876  additionalPeakFound = true;
1877  pExists = true;
1878  }
1879  }
1880  }
1881 
1882  peaksNeededRemaining -= 1;
1883  }
1884  }
1885  }
1886 
1887  //============================ I really hate to do this... SORRY!
1888  if ( additionalPeakFound )
1889  {
1890  if ( verbose > 3 )
1891  {
1892  std::cout << ">>>>>>>>>>>>>>>>> Missing peaks discovered. Redoing the analysis." << std::endl;
1893  }
1894  goto ROUGH_RESTART;
1895  }
1896 
1897  //============================ if not continuous, check if 0 angle helps
1898  if ( !groupContinuous )
1899  {
1900  groupAngleHits.emplace_back ( 0.0 );
1901  std::sort ( groupAngleHits.begin(), groupAngleHits.end() );
1902 
1903  if ( groupAngleHits.size () > 1 )
1904  {
1905  groupContinuous = true;
1906  for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( groupAngleHits.size () ); iter++ )
1907  {
1908  if ( ( ( groupAngleHits.at(iter-1) + groupAngle ) < ( groupAngleHits.at(iter) + errTolerance ) ) &&
1909  ( ( groupAngleHits.at(iter-1) + groupAngle ) > ( groupAngleHits.at(iter) - errTolerance ) ) )
1910  {
1911  ;
1912  }
1913  else
1914  {
1915  groupContinuous = false;
1916  break;
1917  }
1918  }
1919  }
1920  }
1921 
1922  //============================ Remove this symmetry for this axis as it has been tried
1923  thisAxisSyms.emplace_back ( angsToTry.at(angToTry) );
1924 
1925  //============================ If symmetry found, save it and set group peak heights to 0
1926  if ( groupContinuous )
1927  {
1928  // If ALL peaks are found
1929  if ( static_cast<int> ( groupAngleHits.size() ) >= static_cast<int> ( angsToTry.at(angToTry) ) )
1930  {
1931  // Init
1932  std::sort ( groupMembers.begin(), groupMembers.end(), [](const std::array<double,2>& a, const std::array<double,2>& b) { return a[1] > b[1]; } );
1933  groupMembers.erase ( std::unique( groupMembers.begin(), groupMembers.end() ), groupMembers.end() );
1934 
1935  if ( verbose > 3 )
1936  {
1937  printf ( ">>>>>>>>>>> Found symmetry C%+.0f on lines ", angsToTry.at(angToTry) );
1938  for ( unsigned int i = 0; i < static_cast<unsigned int> ( groupMembers.size() ); i++ )
1939  {
1940  printf ( "%+.0f ", groupMembers.at(i)[0] );
1941  }
1942  printf ( "\n" );
1943  }
1944 
1945  // Save
1946  retArr[0] = angsToTry.at(angToTry);
1947 
1948  meanX = 0.0;
1949  meanY = 0.0;
1950  meanZ = 0.0;
1951  meanPeak = 0.0;
1952  totVals = 0.0;
1953 
1954  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( groupMembers.size() ); iter++ )
1955  {
1956  if ( ( grpPeaks.at(groupMembers.at(iter)[0])[0] == 0.0 ) && ( grpPeaks.at(groupMembers.at(iter)[0])[1] == 0.0 ) &&
1957  ( grpPeaks.at(groupMembers.at(iter)[0])[2] == 0.0 ) ) { continue; }
1958  if ( groupMembers.at(iter)[1] == 0.0 ) { continue; }
1959 
1960  meanX += grpPeaks.at(groupMembers.at(iter)[0])[0];
1961  meanY += grpPeaks.at(groupMembers.at(iter)[0])[1];
1962  meanZ += grpPeaks.at(groupMembers.at(iter)[0])[2];
1963  if ( grpPeaks.at(groupMembers.at(iter)[0])[6] == 0.0 )
1964  {
1965  meanPeak += groupMembers.at(iter)[1];
1966  }
1967  totVals += 1.0;
1968  }
1969 
1970  if ( ( meanX == 0.0 ) && ( meanY == 0.0 ) && ( meanZ == 0.0 ) )
1971  {
1972  retArr[1] = 0.0;
1973  retArr[2] = 0.0;
1974  retArr[3] = 1.0;
1975  retArr[4] = groupMembers.at(0)[1];
1976  }
1977  else
1978  {
1979  meanX /= totVals;
1980  meanY /= totVals;
1981  meanZ /= totVals;
1982  meanPeak /= totVals;
1983 
1984  retArr[1] = meanX;
1985  retArr[2] = meanY;
1986  retArr[3] = meanZ;
1987  retArr[4] = meanPeak;
1988  }
1989 
1990  ret.emplace_back ( retArr );
1991 
1992  grpPeaks.at(0)[4] = 0.0;
1993  }
1994  else
1995  {
1996  //==================== The group is continuous, but does not include ALL angles and so symmetry is NOT found
1997  grpPeaks.at(0)[4] = 0.0;
1998  }
1999  }
2000 
2001  //============================ Set the main peak height to 0 and try again
2002  if ( !groupContinuous )
2003  {
2004  grpPeaks.at(0)[4] = 0.0;
2005  }
2006  }
2007  }
2008  }
2009  if ( verbose > 3 )
2010  {
2011  std::cout << ">>>>> C symmetries detection complete." << std::endl;
2012  }
2013 
2014  //======================================== Check for identical symmetries (coming from finding two peaks with almost identical angle)
2015  if ( ret.size() > 1 )
2016  {
2017  std::vector<unsigned int> removeThese;
2018  for ( unsigned int symIt = 0; symIt < static_cast<unsigned int> ( ret.size() ); symIt++ )
2019  {
2020  for ( unsigned int existSym = 0; existSym < static_cast<unsigned int> ( ret.size() ); existSym++ )
2021  {
2022  if ( symIt >= existSym ) { continue; }
2023 
2024  //============================ Is axis the same?
2025  if ( ( ( ret.at(existSym)[1] < ( ret.at(symIt)[1] + errTolerance ) ) && ( ret.at(existSym)[1] > ( ret.at(symIt)[1] - errTolerance ) ) ) &&
2026  ( ( ret.at(existSym)[2] < ( ret.at(symIt)[2] + errTolerance ) ) && ( ret.at(existSym)[2] > ( ret.at(symIt)[2] - errTolerance ) ) ) &&
2027  ( ( ret.at(existSym)[3] < ( ret.at(symIt)[3] + errTolerance ) ) && ( ret.at(existSym)[3] > ( ret.at(symIt)[3] - errTolerance ) ) ) )
2028  {
2029  //======================== Is symmetry the same?
2030  if ( ret.at(existSym)[0] == ret.at(symIt)[0] )
2031  {
2032  //==================== Mark for removal
2033  if ( ret.at(existSym)[4] > ret.at(symIt)[4] )
2034  {
2035  removeThese.emplace_back ( symIt );
2036  }
2037  else
2038  {
2039  removeThese.emplace_back ( existSym );
2040  }
2041  }
2042 
2043  }
2044  }
2045  }
2046 
2047  std::sort ( removeThese.begin(), removeThese.end(), std::greater<unsigned int>() );
2048  removeThese.erase ( std::unique( removeThese.begin(), removeThese.end() ), removeThese.end() );
2049 
2050  for ( unsigned int remIt = 0; remIt < static_cast<unsigned int> ( removeThese.size() ); remIt++ )
2051  {
2052  ret.erase ( ret.begin() + removeThese.at(remIt) );
2053  }
2054  }
2055  if ( verbose > 3 )
2056  {
2057  std::cout << ">>>>> Duplicate C symmetries removal complete." << std::endl;
2058  }
2059 
2060  //======================================== Free memory
2061  if ( freeMem )
2062  {
2063  fftw_free ( this->_so3InvCoeffs );
2064  this->_so3InvCoeffs = nullptr;
2065  }
2066 
2067  //======================================== Sort by mean peak height
2068  std::sort ( ret.begin(), ret.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
2069 
2070  //======================================== Done
2071  this->_CSymmsFound = true;
2072  return ( ret );
2073 
2074 }
2075 
2085 {
2086  //======================================== Free SO3 Inv Map
2087  if ( this->_so3InvCoeffs != nullptr )
2088  {
2089  fftw_free ( this->_so3InvCoeffs );
2090  this->_so3InvCoeffs = nullptr;
2091  }
2092 
2093  //============== Done
2094  return ;
2095 
2096 }
2097 
2109 {
2110  return ( this->_peakHeightThr );
2111 
2112 }
2113 
2128 std::vector< std::array<double,5> > ProSHADE_internal::ProSHADE_comparePairwise::findCnSymmetryClear ( std::vector< std::array<double,5> > CnSymm,
2129  ProSHADE::ProSHADE_settings* settings,
2130  double maxGap,
2131  bool *pf )
2132 {
2133  //======================================== Sanity check
2134  if ( !this->_CSymmsFound )
2135  {
2136  std::cerr << "!!! ProSHADE ERROR !!! Tried to simplify the C symmetry peaks before detecting them. Terminating..." << std::endl;
2137 
2138  if ( settings->htmlReport )
2139  {
2140  std::stringstream hlpSS;
2141  hlpSS << "<font color=\"red\">" << "Tried to simplify the C symmetry peaks before detecting them. This looks like internal bug, please report this case." << "</font>";
2142  rvapi_set_text ( hlpSS.str().c_str(),
2143  "ProgressSection",
2144  settings->htmlReportLineProgress,
2145  1,
2146  1,
2147  1 );
2148  settings->htmlReportLineProgress += 1;
2149  rvapi_flush ( );
2150  }
2151 
2152  exit ( -1 );
2153  }
2154  if ( CnSymm.size() < 1 )
2155  {
2156  std::cout << "!!! ProSHADE WARNING !!! Tried to simplify the C-symmetries list, but the list has less than 1 entry." << std::endl;
2157  if ( pf != nullptr )
2158  {
2159  *pf = true;
2160  }
2161 
2162  return ( CnSymm );
2163  }
2164 
2165  //======================================== Initialise
2166  std::vector< std::array<double,5> > ret;
2167  std::vector< std::array<double,2> > pHeights;
2168  std::array<double,2> hlpArr;
2169 
2170  //======================================== Get all the peak heights
2171  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( CnSymm.size() ); iter++ )
2172  {
2173  hlpArr[0] = static_cast<double> ( iter );
2174  hlpArr[1] = CnSymm.at(iter)[4];
2175  pHeights.emplace_back ( hlpArr );
2176  }
2177 
2178  //============== Find gaps
2179  std::sort ( pHeights.begin(), pHeights.end(), [](const std::array<double,2>& a, const std::array<double,2>& b) { return a[1] > b[1]; });
2180  double maxPeakThres = pHeights.at(0)[1] - ( pHeights.at(0)[1] * maxGap );
2181 
2182  ret.emplace_back ( CnSymm.at(pHeights.at(0)[0]) );
2183  for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( pHeights.size() ); iter++ )
2184  {
2185  if ( pHeights.at(iter)[1] >= maxPeakThres )
2186  {
2187  ret.emplace_back ( CnSymm.at(pHeights.at(iter)[0]) );
2188  }
2189  }
2190 
2191  //======================================== Give the highest symmetry first, if small, else sort by peak
2192  std::sort ( ret.begin(), ret.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[0] > b[0]; });
2193 
2194  //======================================== Done
2195  return ( ret );
2196 
2197 }
2198 
2211 std::vector< std::vector< std::array<double,6> > > ProSHADE_internal::ProSHADE_comparePairwise::findDnSymmetry ( std::vector< std::array<double,5> > CnSymm,
2212  double axisErrorTolerance )
2213 {
2214  //======================================== Initialise local values
2215  double dotProduct = 0.0;
2216  bool changeFound = true;
2217  bool matchedAll = true;
2218  bool identical = true;
2219  std::vector< std::vector<std::array<double,6> > > ret;
2220  std::vector< std::vector<std::array<double,6> > > retHlp;
2221  std::vector< std::array<double,6> > hlpVec;
2222  std::array<double,6> hlpArr;
2223 
2224  //======================================== Create all pairs
2225  for ( unsigned int cSym1 = 0; cSym1 < static_cast<unsigned int> ( CnSymm.size() ); cSym1++ )
2226  {
2227  for ( unsigned int cSym2 = 1; cSym2 < static_cast<unsigned int> ( CnSymm.size() ); cSym2++ )
2228  {
2229  //================================ Use unique only
2230  if ( cSym1 >= cSym2 ) { continue; }
2231 
2232  //================================ Initialise
2233  hlpVec.clear ( );
2234 
2235  //================================ Compare
2236  dotProduct = ( CnSymm.at(cSym1)[1] * CnSymm.at(cSym2)[1] +
2237  CnSymm.at(cSym1)[2] * CnSymm.at(cSym2)[2] +
2238  CnSymm.at(cSym1)[3] * CnSymm.at(cSym2)[3] );
2239 
2240  //================================ Check if it is a pair
2241  if ( std::abs ( dotProduct ) < axisErrorTolerance )
2242  {
2243  hlpArr[0] = CnSymm.at(cSym1)[0];
2244  hlpArr[1] = CnSymm.at(cSym1)[1];
2245  hlpArr[2] = CnSymm.at(cSym1)[2];
2246  hlpArr[3] = CnSymm.at(cSym1)[3];
2247  hlpArr[4] = CnSymm.at(cSym1)[4];
2248  hlpVec.emplace_back ( hlpArr );
2249 
2250  hlpArr[0] = CnSymm.at(cSym2)[0];
2251  hlpArr[1] = CnSymm.at(cSym2)[1];
2252  hlpArr[2] = CnSymm.at(cSym2)[2];
2253  hlpArr[3] = CnSymm.at(cSym2)[3];
2254  hlpArr[4] = CnSymm.at(cSym2)[4];
2255  hlpVec.emplace_back ( hlpArr );
2256  retHlp.emplace_back ( hlpVec );
2257  }
2258  }
2259  }
2260 
2261  //======================================== If nothing found, done
2262  if ( retHlp.size() == 0 ) { this->_DSymmsFound = true; return ( ret ); }
2263 
2264  //======================================== Keep adding matches until no more are found
2265  changeFound = true;
2266  while ( changeFound )
2267  {
2268  //==================================== Initialise
2269  changeFound = false;
2270 
2271  //==================================== Compare each CSym to all existing DSyms
2272  for ( unsigned int cSym1 = 0; cSym1 < static_cast<unsigned int> ( CnSymm.size() ); cSym1++ )
2273  {
2274  //================================ Compare to all existing DSyms
2275  for ( unsigned int eSym = 0; eSym < static_cast<unsigned int> ( retHlp.size() ); eSym++ )
2276  {
2277  //============================ Initialise
2278  matchedAll = true;
2279  identical = false;
2280 
2281  //============================ Check all members if ALL are orthogonal
2282  for ( unsigned int mem = 0; mem < static_cast<unsigned int> ( retHlp.at(eSym).size() ); mem++ )
2283  {
2284  //======================== Are identical
2285  if ( ( CnSymm.at(cSym1)[0] == retHlp.at(eSym).at(mem)[0] ) &&
2286  ( CnSymm.at(cSym1)[1] == retHlp.at(eSym).at(mem)[1] ) &&
2287  ( CnSymm.at(cSym1)[2] == retHlp.at(eSym).at(mem)[2] ) &&
2288  ( CnSymm.at(cSym1)[3] == retHlp.at(eSym).at(mem)[3] ) )
2289  {
2290  identical = true;
2291  }
2292 
2293  //======================== Compare
2294  dotProduct = ( CnSymm.at(cSym1)[1] * retHlp.at(eSym).at(mem)[1] +
2295  CnSymm.at(cSym1)[2] * retHlp.at(eSym).at(mem)[2] +
2296  CnSymm.at(cSym1)[3] * retHlp.at(eSym).at(mem)[3] );
2297 
2298  if ( std::abs ( dotProduct ) < axisErrorTolerance ) { ; }
2299  else { matchedAll = false; }
2300  }
2301 
2302  //============================ If identical, go to next
2303  if ( identical ) { continue; }
2304 
2305  //============================ If ALL are, then add
2306  if ( matchedAll )
2307  {
2308  hlpArr[0] = CnSymm.at(cSym1)[0];
2309  hlpArr[1] = CnSymm.at(cSym1)[1];
2310  hlpArr[2] = CnSymm.at(cSym1)[2];
2311  hlpArr[3] = CnSymm.at(cSym1)[3];
2312  hlpArr[4] = CnSymm.at(cSym1)[4];
2313  retHlp.at(eSym).emplace_back ( hlpArr );
2314 
2315  changeFound = true;
2316  }
2317  }
2318  }
2319  }
2320 
2321  //======================================== Remove identical groups
2322  for ( unsigned int eSym1 = 0; eSym1 < static_cast<unsigned int> ( retHlp.size() ); eSym1++ )
2323  {
2324  //==================================== Already removed?
2325  if ( retHlp.at(eSym1).at(0)[0] == -1.0 ) { continue; }
2326 
2327  for ( unsigned int eSym2 = 1; eSym2 < static_cast<unsigned int> ( retHlp.size() ); eSym2++ )
2328  {
2329  //================================ Already removed?
2330  if ( retHlp.at(eSym2).at(0)[0] == -1.0 ) { continue; }
2331 
2332  //================================ Use unique only
2333  if ( eSym1 >= eSym2 ) { continue; }
2334 
2335  //================================ Are identical?
2336  identical = true;
2337  for ( unsigned int mem1 = 0; mem1 < static_cast<unsigned int> ( retHlp.at(eSym1).size() ); mem1++ )
2338  {
2339  matchedAll = false;
2340  for ( unsigned int mem2 = 0; mem2 < static_cast<unsigned int> ( retHlp.at(eSym2).size() ); mem2++ )
2341  {
2342  if ( ( retHlp.at(eSym1).at(mem1)[0] == retHlp.at(eSym2).at(mem2)[0] ) &&
2343  ( retHlp.at(eSym1).at(mem1)[1] == retHlp.at(eSym2).at(mem2)[1] ) &&
2344  ( retHlp.at(eSym1).at(mem1)[2] == retHlp.at(eSym2).at(mem2)[2] ) &&
2345  ( retHlp.at(eSym1).at(mem1)[3] == retHlp.at(eSym2).at(mem2)[3] ) )
2346  {
2347  matchedAll = true;
2348  }
2349  }
2350 
2351  if ( !matchedAll ) { identical = false; break; }
2352  }
2353 
2354  if ( identical )
2355  {
2356  retHlp.at(eSym2).at(0)[0] = -1.0;
2357  }
2358  }
2359  }
2360 
2361  //======================================== Re-save only non-identical
2362  double avgPeak = 0.0;
2363  double sumSym = 0.0;
2364 
2365  for ( unsigned int eSym1 = 0; eSym1 < static_cast<unsigned int> ( retHlp.size() ); eSym1++ )
2366  {
2367  if ( retHlp.at(eSym1).at(0)[0] == -1.0 ) { continue; }
2368  else
2369  {
2370  //================================ This is a true result. Now find the average peak
2371  avgPeak = retHlp.at(eSym1).at(0)[4] * retHlp.at(eSym1).at(0)[0];
2372  sumSym = retHlp.at(eSym1).at(0)[0];
2373  for ( unsigned int mem = 1; mem < static_cast<unsigned int> ( retHlp.at(eSym1).size() ); mem++ )
2374  {
2375  avgPeak += retHlp.at(eSym1).at(mem)[4] * retHlp.at(eSym1).at(mem)[0];
2376  sumSym += retHlp.at(eSym1).at(mem)[0];
2377  }
2378  avgPeak /= sumSym;
2379 
2380  hlpVec.clear ( );
2381  for ( unsigned int mem = 0; mem < static_cast<unsigned int> ( retHlp.at(eSym1).size() ); mem++ )
2382  {
2383  hlpArr[0] = retHlp.at(eSym1).at(mem)[0];
2384  hlpArr[1] = retHlp.at(eSym1).at(mem)[1];
2385  hlpArr[2] = retHlp.at(eSym1).at(mem)[2];
2386  hlpArr[3] = retHlp.at(eSym1).at(mem)[3];
2387  hlpArr[4] = retHlp.at(eSym1).at(mem)[4];
2388  hlpArr[5] = avgPeak;
2389  hlpVec.emplace_back ( hlpArr );
2390  }
2391  ret.emplace_back ( hlpVec );
2392  }
2393  }
2394 
2395  //======================================== Sort by average peak height
2396  std::sort ( ret.begin(), ret.end(), [](const std::vector< std::array<double,6> >& a, const std::vector< std::array<double,6> >& b) { return a.at(0)[5] > b.at(0)[5]; });
2397 
2398  //======================================== Done
2399  this->_DSymmsFound = true;
2400  return ( ret );
2401 
2402 }
2403 
2418 std::vector< std::vector< std::array<double,6> > > ProSHADE_internal::ProSHADE_comparePairwise::findDnSymmetryClear ( std::vector< std::vector< std::array<double,6> > > DnSymm,
2419  ProSHADE::ProSHADE_settings* settings,
2420  double maxGap,
2421  bool *pf )
2422 {
2423  //======================================== Sanity check
2424  if ( !this->_DSymmsFound )
2425  {
2426  std::cerr << "!!! ProSHADE ERROR !!! Tried to simplify the D symmetry peaks before detecting them. Terminating..." << std::endl;
2427 
2428  if ( settings->htmlReport )
2429  {
2430  std::stringstream hlpSS;
2431  hlpSS << "<font color=\"red\">" << "Tried to simplify the D symmetry peaks before detecting them. This looks like internal bug, please report this case." << "</font>";
2432  rvapi_set_text ( hlpSS.str().c_str(),
2433  "ProgressSection",
2434  settings->htmlReportLineProgress,
2435  1,
2436  1,
2437  1 );
2438  settings->htmlReportLineProgress += 1;
2439  rvapi_flush ( );
2440  }
2441 
2442  exit ( -1 );
2443  }
2444  if ( DnSymm.size() < 1 )
2445  {
2446  std::cout << "!!! ProSHADE WARNING !!! Tried to simplify the D-symmetries list, but the list has less than 1 entry." << std::endl;
2447  if ( pf != nullptr )
2448  {
2449  *pf = true;
2450  }
2451  std::vector< std::vector< std::array<double,6> > > ret;
2452  return ( ret );
2453  }
2454 
2455  //======================================== Initialise
2456  std::vector< std::vector<std::array<double,6> > > ret;
2457  std::vector< std::array<double,6> > hVec;
2458  std::vector< std::vector<std::array<double,4> > > pHeights;
2459  std::vector< std::array<double,4> > hlpVec;
2460  std::array<double,4> hlpArr;
2461 
2462  //======================================== Get all the average peak heights
2463  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( DnSymm.size() ); iter++ )
2464  {
2465  hlpVec.clear ( );
2466  hlpArr[0] = static_cast<double> ( iter );
2467  hlpArr[3] = 0;
2468 
2469  for ( unsigned int it = 0; it < static_cast<unsigned int> ( DnSymm.at(iter).size() ); it++ )
2470  {
2471  hlpArr[1] = static_cast<double> ( it );
2472  hlpArr[2] = DnSymm.at(iter).at(it)[0];
2473  hlpArr[3] = DnSymm.at(iter).at(it)[4];
2474  hlpVec.emplace_back ( hlpArr );
2475  }
2476 
2477  pHeights.emplace_back ( hlpVec );
2478  }
2479 
2480  //======================================== Sort peak groups by C symmetry and then by peak height
2481  double lSym = 0.0;
2482  double mPeak = 0.0;
2483  double firstIt = 0.0;
2484  double secondIt = 0.0;
2485  double avgPeak = 0.0;
2486  std::vector< std::array<double,2> > combsEx;
2487  std::array<double,2> hA;
2488  bool uniqCombo = true;
2489  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( pHeights.size() ); iter++ )
2490  {
2491  //==================================== Sort by symmetry
2492  std::sort ( pHeights.at(iter).begin(), pHeights.at(iter).end(), [](const std::array<double,4>& a, const std::array<double,4>& b) { return a[2] > b[2]; });
2493 
2494  //==================================== For the highest symmetry, find highest peak
2495  lSym = pHeights.at(iter).at(0)[2];
2496  mPeak = 0.0;
2497  for ( unsigned int it = 0; it < static_cast<unsigned int> ( pHeights.at(iter).size() ); it++ )
2498  {
2499  if ( pHeights.at(iter).at(it)[2] == lSym )
2500  {
2501  if ( mPeak < pHeights.at(iter).at(it)[3] )
2502  {
2503  mPeak = pHeights.at(iter).at(it)[3];
2504  firstIt = static_cast<double> ( it );
2505  }
2506  }
2507  else
2508  {
2509  secondIt = std::max ( secondIt, pHeights.at(iter).at(it)[2] );
2510  }
2511  }
2512  if ( secondIt == 0.0 ) { ; }
2513  else { lSym = secondIt; }
2514  mPeak = 0.0;
2515  for ( unsigned int it = 0; it < static_cast<unsigned int> ( pHeights.at(iter).size() ); it++ )
2516  {
2517  if ( pHeights.at(iter).at(it)[2] == lSym )
2518  {
2519  if ( mPeak < pHeights.at(iter).at(it)[3] )
2520  {
2521  if ( static_cast<double> ( it ) != firstIt )
2522  {
2523  mPeak = pHeights.at(iter).at(it)[3];
2524  secondIt = static_cast<double> ( it );
2525  }
2526  }
2527  }
2528  }
2529 
2530  if ( iter == 0 )
2531  {
2532  hVec.clear ( );
2533  hVec.emplace_back ( DnSymm.at(iter).at(pHeights.at(iter).at(firstIt)[1]) );
2534  hVec.emplace_back ( DnSymm.at(iter).at(pHeights.at(iter).at(secondIt)[1]) );
2535  avgPeak = ( DnSymm.at(iter).at(pHeights.at(iter).at(firstIt)[1])[4] + DnSymm.at(iter).at(pHeights.at(iter).at(secondIt)[1])[4] ) / 2.0;
2536  hVec.at(0)[5] = avgPeak;
2537  hVec.at(1)[5] = avgPeak;
2538  ret.emplace_back ( hVec );
2539  hA[0] = DnSymm.at(iter).at(pHeights.at(iter).at(firstIt)[1])[0];
2540  hA[1] = DnSymm.at(iter).at(pHeights.at(iter).at(secondIt)[1])[0];
2541  combsEx.emplace_back ( hA );
2542  }
2543  else
2544  {
2545  uniqCombo = true;
2546  for ( unsigned int it = 0; it < static_cast<unsigned int> ( combsEx.size() ); it++ )
2547  {
2548  if ( ( combsEx.at(it)[0] == DnSymm.at(iter).at(pHeights.at(iter).at(firstIt)[1])[0] ) &&
2549  ( combsEx.at(it)[1] == DnSymm.at(iter).at(pHeights.at(iter).at(secondIt)[1])[0] ) )
2550  {
2551  uniqCombo = false;
2552  }
2553  }
2554 
2555  if ( uniqCombo )
2556  {
2557  hVec.clear ( );
2558  hVec.emplace_back ( DnSymm.at(iter).at(pHeights.at(iter).at(firstIt)[1]) );
2559  hVec.emplace_back ( DnSymm.at(iter).at(pHeights.at(iter).at(secondIt)[1]) );
2560  ret.emplace_back ( hVec );
2561  hA[0] = DnSymm.at(iter).at(pHeights.at(iter).at(firstIt)[1])[0];
2562  hA[1] = DnSymm.at(iter).at(pHeights.at(iter).at(secondIt)[1])[0];
2563  combsEx.emplace_back ( hA );
2564  }
2565  }
2566  }
2567 
2568  //======================================== Sort by symmetry size if possible
2569  if ( ret.size () < 6 )
2570  {
2571  std::sort ( ret.begin(), ret.end(), [](const std::vector< std::array<double,6> >& a, const std::vector< std::array<double,6> >& b) { return a.at(0)[0] > b.at(0)[0]; });
2572  }
2573 
2574  //======================================== Done
2575  return ( ret );
2576 
2577 }
2578 
2595  double beta,
2596  double gamma,
2597  ProSHADE::ProSHADE_settings* settings,
2598  double passThreshold )
2599 {
2600  //======================================== Initialise
2601  bool ret = false;
2602  double dist = 0.0;
2603 
2604  //======================================== Get Wigner D matrices for the given angles
2605  this->setEulerAngles ( alpha, beta, gamma );
2606  this->generateWignerMatrices ( settings );
2607 
2608  //======================================== Compute full rotation function distance
2609  dist = this->getRotCoeffDistance ( settings );
2610 
2611  //======================================== Decide if rotation is real or ghost
2612  if ( dist >= passThreshold ) { ret = true; }
2613  else { ret = false; }
2614 
2615  //======================================== Done
2616  return ( ret );
2617 
2618 }
2619 
2633 std::vector< std::vector< std::array<double,5> > > ProSHADE_internal::ProSHADE_comparePairwise::findTetrSymmetry ( std::vector< std::array<double,5> > CnSymm,
2634  double *tetrSymmPeakAvg,
2635  double axisErrorTolerance )
2636 {
2637  //======================================== Intialise local variables
2638  std::vector< std::vector<std::array<double,5> > > ret;
2639  std::vector< std::array<double,5> > hlpVec;
2640  std::array<double,5> hlpArr;
2641  std::vector<unsigned int> C3s;
2642  double dotProduct = 0.0;
2643  *tetrSymmPeakAvg = 0.0;
2644  double totPairs = 0.0;
2645 
2646  //======================================== Find all C3 symmetries
2647  for ( unsigned int cSym = 0; cSym < static_cast<unsigned int> ( CnSymm.size() ); cSym++ )
2648  {
2649  if ( CnSymm.at(cSym)[0] == 3 )
2650  {
2651  C3s.emplace_back ( cSym );
2652  }
2653  }
2654 
2655  //======================================== Check for any C3 symmetry with axis angle acos ( 1.0 / 3.0 )
2656  for ( unsigned int c3Sym = 0; c3Sym < static_cast<unsigned int> ( C3s.size() ); c3Sym++ )
2657  {
2658  for ( unsigned int cSym = 0; cSym < static_cast<unsigned int> ( CnSymm.size() ); cSym++ )
2659  {
2660  if ( CnSymm.at(cSym)[0] == 3 )
2661  {
2662  if ( cSym == C3s.at(c3Sym) ) { continue; }
2663 
2664  dotProduct = ( CnSymm.at(C3s.at(c3Sym))[1] * CnSymm.at(cSym)[1] +
2665  CnSymm.at(C3s.at(c3Sym))[2] * CnSymm.at(cSym)[2] +
2666  CnSymm.at(C3s.at(c3Sym))[3] * CnSymm.at(cSym)[3] );
2667 
2668  if ( ( ( 1.0 / 3.0 ) > ( dotProduct - axisErrorTolerance ) ) &&
2669  ( ( 1.0 / 3.0 ) < ( dotProduct + axisErrorTolerance ) ) )
2670  {
2671  //======================== Found tetrahedral symmetry
2672  hlpVec.clear ( );
2673 
2674  hlpArr[0] = CnSymm.at(C3s.at(c3Sym))[0];
2675  hlpArr[1] = CnSymm.at(C3s.at(c3Sym))[1];
2676  hlpArr[2] = CnSymm.at(C3s.at(c3Sym))[2];
2677  hlpArr[3] = CnSymm.at(C3s.at(c3Sym))[3];
2678  hlpArr[4] = ( CnSymm.at(C3s.at(c3Sym))[4] + CnSymm.at(cSym)[4] ) / 2.0;
2679  hlpVec.emplace_back ( hlpArr );
2680 
2681  hlpArr[0] = CnSymm.at(cSym)[0];
2682  hlpArr[1] = CnSymm.at(cSym)[1];
2683  hlpArr[2] = CnSymm.at(cSym)[2];
2684  hlpArr[3] = CnSymm.at(cSym)[3];
2685  hlpVec.emplace_back ( hlpArr );
2686  ret.emplace_back ( hlpVec );
2687  *tetrSymmPeakAvg += hlpArr[4];
2688  totPairs += 1.0;
2689  }
2690  }
2691  }
2692  }
2693 
2694  //======================================== Find overall average peak
2695  *tetrSymmPeakAvg /= totPairs;
2696 
2697  //======================================== Done
2698  return ( ret );
2699 
2700 }
2701 
2715 std::vector< std::vector< std::array<double,5> > > ProSHADE_internal::ProSHADE_comparePairwise::findOctaSymmetry ( std::vector< std::array<double,5> > CnSymm,
2716  double *octaSymmPeakAvg,
2717  double axisErrorTolerance )
2718 {
2719  //======================================== Intialise local variables
2720  std::vector< std::vector< std::array<double,5 > > > ret;
2721  std::vector< std::array<double,5> > hlpVec;
2722  std::array<double,5> hlpArr;
2723  std::vector<unsigned int> C4s;
2724  double dotProduct = 0.0;
2725  *octaSymmPeakAvg = 0.0;
2726  double totPairs = 0.0;
2727 
2728  //======================================== Find all C4 symmetries
2729  for ( unsigned int cSym = 0; cSym < static_cast<unsigned int> ( CnSymm.size() ); cSym++ )
2730  {
2731  if ( CnSymm.at(cSym)[0] == 4 )
2732  {
2733  C4s.emplace_back ( cSym );
2734  }
2735  }
2736 
2737  //======================================== Check for any C3 symmetry with axis angle PI - acos ( 1.0 / 3.0 )
2738  for ( unsigned int c4Sym = 0; c4Sym < static_cast<unsigned int> ( C4s.size() ); c4Sym++ )
2739  {
2740  for ( unsigned int cSym = 0; cSym < static_cast<unsigned int> ( CnSymm.size() ); cSym++ )
2741  {
2742  if ( CnSymm.at(cSym)[0] == 3 )
2743  {
2744  dotProduct = ( CnSymm.at(C4s.at(c4Sym))[1] * CnSymm.at(cSym)[1] +
2745  CnSymm.at(C4s.at(c4Sym))[2] * CnSymm.at(cSym)[2] +
2746  CnSymm.at(C4s.at(c4Sym))[3] * CnSymm.at(cSym)[3] );
2747 
2748  if ( ( ( 1.0 / sqrt ( 3.0 ) ) > ( dotProduct - axisErrorTolerance ) ) &&
2749  ( ( 1.0 / sqrt ( 3.0 ) ) < ( dotProduct + axisErrorTolerance ) ) )
2750  {
2751  //======================== Found octahedral symmetry
2752  hlpVec.clear ( );
2753 
2754  hlpArr[0] = CnSymm.at(C4s.at(c4Sym))[0];
2755  hlpArr[1] = CnSymm.at(C4s.at(c4Sym))[1];
2756  hlpArr[2] = CnSymm.at(C4s.at(c4Sym))[2];
2757  hlpArr[3] = CnSymm.at(C4s.at(c4Sym))[3];
2758  hlpArr[4] = ( CnSymm.at(C4s.at(c4Sym))[4] + CnSymm.at(cSym)[4] ) / 2.0;
2759  hlpVec.emplace_back ( hlpArr );
2760 
2761  hlpArr[0] = CnSymm.at(cSym)[0];
2762  hlpArr[1] = CnSymm.at(cSym)[1];
2763  hlpArr[2] = CnSymm.at(cSym)[2];
2764  hlpArr[3] = CnSymm.at(cSym)[3];
2765  hlpVec.emplace_back ( hlpArr );
2766  ret.emplace_back ( hlpVec );
2767  *octaSymmPeakAvg += hlpArr[4];
2768  totPairs += 1.0;
2769  }
2770  }
2771  }
2772  }
2773 
2774  //======================================== Find overall average peak
2775  *octaSymmPeakAvg /= totPairs;
2776 
2777  //======================================== Done
2778  return ( ret );
2779 
2780 }
2781 
2795 std::vector< std::vector< std::array<double,5> > > ProSHADE_internal::ProSHADE_comparePairwise::findIcosSymmetry ( std::vector< std::array<double,5> > CnSymm,
2796  double *icosSymmPeakAvg,
2797  double axisErrorTolerance )
2798 {
2799  //======================================== Intialise local variables
2800  std::vector< std::vector< std::array<double,5> > > ret;
2801  std::vector< std::array<double,5> > hlpVec;
2802  std::array<double,5> hlpArr;
2803  std::vector<unsigned int> C5s;
2804  double dotProduct = 0.0;
2805  *icosSymmPeakAvg = 0.0;
2806  double totPairs = 0.0;
2807 
2808  //======================================== Find all C5 symmetries
2809  for ( unsigned int cSym = 0; cSym < static_cast<unsigned int> ( CnSymm.size() ); cSym++ )
2810  {
2811  if ( CnSymm.at(cSym)[0] == 5 )
2812  {
2813  C5s.emplace_back ( cSym );
2814  }
2815  }
2816 
2817  //======================================== Check for any C3 symmetry with axis angle acos ( sqrt ( 5.0 ) / 3.0 )
2818  for ( unsigned int c5Sym = 0; c5Sym < static_cast<unsigned int> ( C5s.size() ); c5Sym++ )
2819  {
2820  for ( unsigned int cSym = 0; cSym < static_cast<unsigned int> ( CnSymm.size() ); cSym++ )
2821  {
2822  if ( CnSymm.at(cSym)[0] == 3 )
2823  {
2824  dotProduct = ( CnSymm.at(C5s.at(c5Sym))[1] * CnSymm.at(cSym)[1] +
2825  CnSymm.at(C5s.at(c5Sym))[2] * CnSymm.at(cSym)[2] +
2826  CnSymm.at(C5s.at(c5Sym))[3] * CnSymm.at(cSym)[3] );
2827 
2828  if ( ( ( sqrt ( 5.0 ) / 3.0 ) > ( dotProduct - axisErrorTolerance ) ) &&
2829  ( ( sqrt ( 5.0 ) / 3.0 ) < ( dotProduct + axisErrorTolerance ) ) )
2830  {
2831  //======================== Found icosahedral symmetry
2832  hlpVec.clear ( );
2833 
2834  hlpArr[0] = CnSymm.at(C5s.at(c5Sym))[0];
2835  hlpArr[1] = CnSymm.at(C5s.at(c5Sym))[1];
2836  hlpArr[2] = CnSymm.at(C5s.at(c5Sym))[2];
2837  hlpArr[3] = CnSymm.at(C5s.at(c5Sym))[3];
2838  hlpArr[4] = ( CnSymm.at(C5s.at(c5Sym))[4] + CnSymm.at(cSym)[4] ) / 2.0;
2839  hlpVec.emplace_back ( hlpArr );
2840 
2841  hlpArr[0] = CnSymm.at(cSym)[0];
2842  hlpArr[1] = CnSymm.at(cSym)[1];
2843  hlpArr[2] = CnSymm.at(cSym)[2];
2844  hlpArr[3] = CnSymm.at(cSym)[3];
2845  hlpVec.emplace_back ( hlpArr );
2846  ret.emplace_back ( hlpVec );
2847  *icosSymmPeakAvg += hlpArr[4];
2848  totPairs += 1.0;
2849  }
2850  }
2851  }
2852  }
2853 
2854  //======================================== Find overall average peak
2855  *icosSymmPeakAvg /= totPairs;
2856 
2857  //======================================== Done
2858  return ( ret );
2859 
2860 }
2861 
2880  ProSHADE::ProSHADE_settings* settings,
2881  std::vector< std::array<double,5> > tetrSymms,
2882  std::vector< std::array<double,5> > allCs,
2883  double axisErrorTolerance,
2884  int verbose )
2885 {
2886  //======================================== Initialise
2887  std::vector< std::array<double,5> > ret;
2888  std::array<double,5> hlpArr;
2889  bool allAnglesMet = false;
2890  double dotProduct = 0.0;
2891 
2892  //======================================== Push the two already known C3 to ret. This is to allow different starting points, should the first one fail
2893  ret.emplace_back ( tetrSymms.at(0) );
2894  ret.emplace_back ( tetrSymms.at(1) );
2895 
2896  if ( verbose > 1 )
2897  {
2898  std::cout << ">> Searching for tetrahederal symmetry axes." << std::endl;
2899  }
2900  if ( verbose > 2 )
2901  {
2902  std::cout << ">>>>> Searching for the four C3 symmetry axes." << std::endl;
2903  }
2904 
2905  //======================================== Find four C3 symmetries with acos(1/3) angle between them
2906  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
2907  {
2908  //==================================== If already found four, do not bother
2909  if ( ret.size() > 2 )
2910  {
2911  break;
2912  }
2913 
2914  //==================================== Search only using C3s
2915  if ( allCs.at(iter)[0] != 3.0 )
2916  {
2917  continue;
2918  }
2919 
2920  //==================================== If this is the first C3, then just save it
2921  if ( ret.size() == 0 )
2922  {
2923  ret.emplace_back ( allCs.at(iter) );
2924  continue;
2925  }
2926 
2927  //==================================== If second or more C3, check if it has the correct angle to all other already found C3s
2928  allAnglesMet = true;
2929  for ( unsigned int i = 0; i < static_cast<unsigned int> ( ret.size() ); i++ )
2930  {
2931  dotProduct = ( ret.at(i)[1] * allCs.at(iter)[1] +
2932  ret.at(i)[2] * allCs.at(iter)[2] +
2933  ret.at(i)[3] * allCs.at(iter)[3] );
2934 
2935  if ( ( ( 1.0 / 3.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
2936  ( ( 1.0 / 3.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
2937  {
2938  ;
2939  }
2940  else
2941  {
2942  allAnglesMet = false;
2943  }
2944  }
2945 
2946  if ( allAnglesMet )
2947  {
2948  ret.emplace_back ( allCs.at(iter) );
2949  }
2950  }
2951 
2952  //======================================== Not enough C3's. Report empty - no tetrahedral elements found
2953  if ( ret.size() != 3 )
2954  {
2955  if ( verbose > 3 )
2956  {
2957  std::cout << ">>>>>>>> Found only " << ret.size() << " C3 symmetries. Searching for the missing ones." << std::endl;
2958  }
2959 
2960  //==================================== If just one is missing, try to find it
2961  if ( ret.size() == 2 )
2962  {
2963  //================================ Find the missing axes
2964  std::vector< std::array<double,3> > axesVec;
2965  std::array<double,3> hArr;
2966  double axNorm, axX, axY, axZ;
2967  for ( double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
2968  {
2969  for ( double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
2970  {
2971  for ( double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
2972  {
2973  if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) { continue; }
2974 
2975  axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
2976  axX = xIt / axNorm;
2977  axY = yIt / axNorm;
2978  axZ = zIt / axNorm;
2979 
2980  allAnglesMet = true;
2981  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
2982  {
2983  dotProduct = ( ret.at(it)[1] * axX +
2984  ret.at(it)[2] * axY +
2985  ret.at(it)[3] * axZ );
2986 
2987  if ( ( ( 1.0 / 3.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
2988  ( ( 1.0 / 3.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
2989  {
2990  ;
2991  }
2992  else
2993  {
2994  allAnglesMet = false;
2995  continue;
2996  }
2997  }
2998 
2999  if ( allAnglesMet )
3000  {
3001  hArr[0] = axX;
3002  hArr[1] = axY;
3003  hArr[2] = axZ;
3004 
3005  allAnglesMet = true;
3006  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3007  {
3008  if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
3009  ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
3010  ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
3011  {
3012  allAnglesMet = false;
3013  }
3014  }
3015 
3016  if ( allAnglesMet )
3017  {
3018  axesVec.emplace_back ( hArr );
3019  }
3020  }
3021  }
3022  }
3023  }
3024 
3025  std::vector< std::array<double,5> > hRes;
3026  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3027  {
3028  //============================ Check if this axis solves the missing problem
3029  double res = cmpObj->maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 3.0, settings );
3030 
3031  hlpArr[0] = 3.0;
3032  hlpArr[1] = axesVec.at(axIt)[0];
3033  hlpArr[2] = axesVec.at(axIt)[1];
3034  hlpArr[3] = axesVec.at(axIt)[2];
3035  hlpArr[4] = res;
3036  hRes.emplace_back ( hlpArr );
3037  }
3038 
3039  std::sort ( hRes.begin(), hRes.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
3040 
3041  if ( hRes.at(0)[4] > ( cmpObj->getPeakThreshold ( ) ) )
3042  {
3043  ret.emplace_back ( hRes.at(0) );
3044  }
3045 
3046  if ( ret.size() != 4 )
3047  {
3048  ret.clear ( );
3049  return ( ret );
3050  }
3051  }
3052  else
3053  {
3054  ret.clear ( );
3055  return ( ret );
3056  }
3057  }
3058  if ( verbose > 3 )
3059  {
3060  std::cout << ">>>>>>>> All C3 symmetries located." << std::endl;
3061  }
3062 
3063  if ( verbose > 2 )
3064  {
3065  std::cout << ">>>>> Searching for the three perpendicular C2 symmetry axes." << std::endl;
3066  }
3067 
3068  //======================================== Search for the D222 required; find three C2s which are perpendicular to each other and have angle acos(1/6) to the first C3
3069  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
3070  {
3071  //==================================== Search only using C2s
3072  if ( allCs.at(iter)[0] != 2.0 )
3073  {
3074  continue;
3075  }
3076 
3077  //==================================== C2 having angle acos(0.5) to the first C3?
3078  dotProduct = ( ret.at(0)[1] * allCs.at(iter)[1] +
3079  ret.at(0)[2] * allCs.at(iter)[2] +
3080  ret.at(0)[3] * allCs.at(iter)[3] );
3081 
3082  if ( ( ( 1.0 / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3083  ( ( 1.0 / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3084  {
3085  //================================ If first C2, just save it
3086  if ( ret.size() == 4 )
3087  {
3088  ret.emplace_back ( allCs.at(iter) );
3089  continue;
3090  }
3091 
3092  //================================ Else, check if perpendicular to other C2s
3093  allAnglesMet = true;
3094  for ( unsigned int i = 4; i < static_cast<unsigned int> ( ret.size() ); i++ )
3095  {
3096  dotProduct = ( ret.at(i)[1] * allCs.at(iter)[1] +
3097  ret.at(i)[2] * allCs.at(iter)[2] +
3098  ret.at(i)[3] * allCs.at(iter)[3] );
3099 
3100  if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3101  ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3102  {
3103  ;
3104  }
3105  else
3106  {
3107  allAnglesMet = false;
3108  }
3109  }
3110 
3111  if ( allAnglesMet )
3112  {
3113  ret.emplace_back ( allCs.at(iter) );
3114  }
3115  }
3116  }
3117 
3118  //======================================== If not found all, report so
3119  if ( ret.size() != 6 )
3120  {
3121  if ( verbose > 3 )
3122  {
3123  std::cout << ">>>>>>>> Found only " << ret.size() - 4 << " C3 symmetries. Searching for the missing ones." << std::endl;
3124  }
3125 
3126  //==================================== If just one is missing, try to find it
3127  if ( ret.size() == 6 )
3128  {
3129  //================================ Find the missing axes
3130  std::vector< std::array<double,3> > axesVec;
3131  std::array<double,3> hArr;
3132  double axNorm, axX, axY, axZ;
3133  for ( double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
3134  {
3135  for ( double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
3136  {
3137  for ( double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
3138  {
3139  if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) { continue; }
3140 
3141  axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
3142  axX = xIt / axNorm;
3143  axY = yIt / axNorm;
3144  axZ = zIt / axNorm;
3145 
3146  allAnglesMet = true;
3147  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3148  {
3149  dotProduct = ( ret.at(it)[1] * axX +
3150  ret.at(it)[2] * axY +
3151  ret.at(it)[3] * axZ );
3152 
3153  if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3154  ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3155  {
3156  ;
3157  }
3158  else
3159  {
3160  allAnglesMet = false;
3161  continue;
3162  }
3163  }
3164 
3165  if ( allAnglesMet )
3166  {
3167  hArr[0] = axX;
3168  hArr[1] = axY;
3169  hArr[2] = axZ;
3170 
3171  allAnglesMet = true;
3172  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3173  {
3174  if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
3175  ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
3176  ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
3177  {
3178  allAnglesMet = false;
3179  }
3180  }
3181 
3182  if ( allAnglesMet )
3183  {
3184  axesVec.emplace_back ( hArr );
3185  }
3186  }
3187  }
3188  }
3189  }
3190 
3191  std::vector< std::array<double,5> > hRes;
3192  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3193  {
3194  //============================ Check if this axis solves the missing problem
3195  double res = cmpObj->maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 2.0, settings );
3196 
3197  hlpArr[0] = 2.0;
3198  hlpArr[1] = axesVec.at(axIt)[0];
3199  hlpArr[2] = axesVec.at(axIt)[1];
3200  hlpArr[3] = axesVec.at(axIt)[2];
3201  hlpArr[4] = res;
3202  hRes.emplace_back ( hlpArr );
3203  }
3204 
3205  std::sort ( hRes.begin(), hRes.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
3206 
3207  if ( hRes.at(0)[4] > ( cmpObj->getPeakThreshold ( ) ) )
3208  {
3209  ret.emplace_back ( hRes.at(0) );
3210  }
3211 
3212  if ( ret.size() != 7 )
3213  {
3214  ret.clear ( );
3215  return ( ret );
3216  }
3217  }
3218  else
3219  {
3220  ret.clear ( );
3221  return ( ret );
3222  }
3223  }
3224  if ( verbose > 3 )
3225  {
3226  std::cout << ">>>>>>>> All C2 symmetries located." << std::endl;
3227  }
3228 
3229  //======================================== Done
3230  return ( ret );
3231 
3232 }
3233 
3243 std::vector< std::array<double,5> > ProSHADE_internal::ProSHADE_symmetry::generateTetrElements ( std::vector< std::array<double,5> > symmAxes,
3244  ProSHADE::ProSHADE_settings* settings,
3245  int verbose )
3246 {
3247  //======================================== Initialise
3248  std::vector< std::array<double,5> > ret;
3249  std::array<double,5> hlpArr;
3250  if ( verbose > 1 )
3251  {
3252  std::cout << ">> Searching for tetrahederal symmetry elements." << std::endl;
3253  }
3254 
3255  //============== Identity element
3256  hlpArr[0] = 1.0;
3257  hlpArr[1] = 1.0;
3258  hlpArr[2] = 0.0;
3259  hlpArr[3] = 0.0;
3260  hlpArr[4] = 0.0;
3261  ret.emplace_back ( hlpArr );
3262 
3263  if ( verbose > 3 )
3264  {
3265  std::cout << ">>>>>>>> Identity element." << std::endl;
3266  }
3267 
3268  //======================================== 3 times 120 degrees by each vertex axes clockwise
3269  for ( unsigned int iter = 0; iter < 3; iter++ )
3270  {
3271  hlpArr[0] = 3.0;
3272  hlpArr[1] = symmAxes.at(iter)[1];
3273  hlpArr[2] = symmAxes.at(iter)[2];
3274  hlpArr[3] = symmAxes.at(iter)[3];
3275  hlpArr[4] = 120.0;
3276  ret.emplace_back ( hlpArr );
3277  }
3278  if ( verbose > 3 )
3279  {
3280  std::cout << ">>>>>>>> 4 clockwise C3 rotations." << std::endl;
3281  }
3282 
3283  //======================================== 4 times 120 degrees by each vertex axes anti-clockwise
3284  for ( unsigned int iter = 0; iter < 3; iter++ )
3285  {
3286  hlpArr[0] = 3.0;
3287  hlpArr[1] = symmAxes.at(iter)[1];
3288  hlpArr[2] = symmAxes.at(iter)[2];
3289  hlpArr[3] = symmAxes.at(iter)[3];
3290  hlpArr[4] = -120.0;
3291  ret.emplace_back ( hlpArr );
3292  }
3293  if ( verbose > 3 )
3294  {
3295  std::cout << ">>>>>>>> 4 anti-clockwise C3 rotations." << std::endl;
3296  }
3297 
3298  //======================================== 3 times 180 degrees by each D222 axis
3299  for ( unsigned int iter = 3; iter < 6; iter++ )
3300  {
3301  hlpArr[0] = 2.0;
3302  hlpArr[1] = symmAxes.at(iter)[1];
3303  hlpArr[2] = symmAxes.at(iter)[2];
3304  hlpArr[3] = symmAxes.at(iter)[3];
3305  hlpArr[4] = 180.0;
3306  ret.emplace_back ( hlpArr );
3307  }
3308  if ( verbose > 3 )
3309  {
3310  std::cout << ">>>>>>>> 3 clockwise D222 rotations." << std::endl;
3311  }
3312  if ( verbose > 2 )
3313  {
3314  std::cout << ">>>>> Tetrahedral symmetry elements generated." << std::endl;
3315  }
3316 
3317  //======================================== Done
3318  return ( ret );
3319 
3320 }
3321 
3340  ProSHADE::ProSHADE_settings* settings,
3341  std::vector< std::array<double,5> > octaSymms,
3342  std::vector< std::array<double,5> > allCs,
3343  double axisErrorTolerance,
3344  int verbose )
3345 {
3346  //======================================== Initialise
3347  std::vector< std::array<double,5> > ret;
3348  std::array<double,5> hlpArr;
3349  bool allAnglesMet = false;
3350  double dotProduct = 0.0;
3351 
3352  //======================================== Push the already known C4 to ret. This is to allow different starting points, should the first one fail
3353  if ( octaSymms.at(0)[0] == 4.0 )
3354  {
3355  ret.emplace_back ( octaSymms.at(0) );
3356  }
3357  else
3358  {
3359  ret.emplace_back ( octaSymms.at(1) );
3360  }
3361 
3362  if ( verbose > 1 )
3363  {
3364  std::cout << ">> Searching for (cub)octahedral symmetry elements." << std::endl;
3365  }
3366  if ( verbose > 2 )
3367  {
3368  std::cout << ">>>>> Searching for the three C4 symmetry elements." << std::endl;
3369  }
3370 
3371  //======================================== Search for the other two C4's perpendicular to the starting C4
3372  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
3373  {
3374  //==================================== Search only using C4s
3375  if ( allCs.at(iter)[0] != 4.0 )
3376  {
3377  continue;
3378  }
3379 
3380  //==================================== If first C4, check for it being perpendicular to the one already present
3381  if ( ret.size() == 1 )
3382  {
3383  dotProduct = ( ret.at(0)[1] * allCs.at(iter)[1] +
3384  ret.at(0)[2] * allCs.at(iter)[2] +
3385  ret.at(0)[3] * allCs.at(iter)[3] );
3386 
3387  if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3388  ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3389  {
3390  ret.emplace_back ( allCs.at(iter) );
3391  continue;
3392  }
3393  }
3394 
3395  //======================================== If two C4s already exist, check both against this one
3396  if ( ret.size() == 2 )
3397  {
3398  allAnglesMet = true;
3399  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3400  {
3401  dotProduct = ( ret.at(it)[1] * allCs.at(iter)[1] +
3402  ret.at(it)[2] * allCs.at(iter)[2] +
3403  ret.at(it)[3] * allCs.at(iter)[3] );
3404 
3405  if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3406  ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3407  {
3408  ;
3409  }
3410  else
3411  {
3412  allAnglesMet = false;
3413  }
3414  }
3415 
3416  if ( allAnglesMet )
3417  {
3418  ret.emplace_back ( allCs.at(iter) );
3419  }
3420  }
3421  }
3422 
3423  //======================================== Not enough C4's. Try to find the missing C4
3424  if ( ret.size() != 3 )
3425  {
3426  if ( verbose > 3 )
3427  {
3428  std::cout << ">>>>>>>> Found only " << ret.size() << " C4 symmetries. Searching for the missing ones." << std::endl;
3429  }
3430 
3431  //==================================== If just one is missing, try to find it
3432  if ( ret.size() == 2 )
3433  {
3434  //================================ Find the missing axes
3435  std::vector< std::array<double,3> > axesVec;
3436  std::array<double,3> hArr;
3437  double axNorm, axX, axY, axZ;
3438  for ( double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
3439  {
3440  for ( double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
3441  {
3442  for ( double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
3443  {
3444  if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) { continue; }
3445 
3446  axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
3447  axX = xIt / axNorm;
3448  axY = yIt / axNorm;
3449  axZ = zIt / axNorm;
3450 
3451  allAnglesMet = true;
3452  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3453  {
3454  dotProduct = ( ret.at(it)[1] * axX +
3455  ret.at(it)[2] * axY +
3456  ret.at(it)[3] * axZ );
3457 
3458  if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3459  ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3460  {
3461  ;
3462  }
3463  else
3464  {
3465  allAnglesMet = false;
3466  continue;
3467  }
3468  }
3469 
3470  if ( allAnglesMet )
3471  {
3472  hArr[0] = axX;
3473  hArr[1] = axY;
3474  hArr[2] = axZ;
3475 
3476  allAnglesMet = true;
3477  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3478  {
3479  if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
3480  ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
3481  ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
3482  {
3483  allAnglesMet = false;
3484  }
3485  }
3486 
3487  if ( allAnglesMet )
3488  {
3489  axesVec.emplace_back ( hArr );
3490  }
3491  }
3492  }
3493  }
3494  }
3495 
3496  std::vector< std::array<double,5> > hRes;
3497  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3498  {
3499  //============================ Check if this axis solves the missing problem
3500  double res = cmpObj->maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 4.0, settings );
3501 
3502  hlpArr[0] = 4.0;
3503  hlpArr[1] = axesVec.at(axIt)[0];
3504  hlpArr[2] = axesVec.at(axIt)[1];
3505  hlpArr[3] = axesVec.at(axIt)[2];
3506  hlpArr[4] = res;
3507  hRes.emplace_back ( hlpArr );
3508  }
3509 
3510  std::sort ( hRes.begin(), hRes.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
3511 
3512  if ( hRes.at(0)[4] > ( cmpObj->getPeakThreshold ( ) ) )
3513  {
3514  ret.emplace_back ( hRes.at(0) );
3515  }
3516 
3517  if ( ret.size() != 3 )
3518  {
3519  ret.clear ( );
3520  return ( ret );
3521  }
3522  }
3523  else
3524  {
3525  ret.clear ( );
3526  return ( ret );
3527  }
3528  }
3529  if ( verbose > 3 )
3530  {
3531  std::cout << ">>>>>>>> All C4 symmetries located." << std::endl;
3532  }
3533 
3534  if ( verbose > 2 )
3535  {
3536  std::cout << ">>>>> Searching for the four C3 symmetry elements." << std::endl;
3537  }
3538 
3539  //======================================== Search for the 4 unique C3s with angles 1/sqrt(3) to the C4s
3540  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
3541  {
3542  //==================================== Search only using C3s
3543  if ( allCs.at(iter)[0] != 3.0 )
3544  {
3545  continue;
3546  }
3547 
3548  //==================================== Check that the prospective C3 has angle 1/sqrt(3) (dihedral angle) to the C4s
3549  allAnglesMet = true;
3550  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3551  {
3552  //================================ Do not check the C3s against themselves
3553  if ( ret.at(it)[0] == 3.0 )
3554  {
3555  continue;
3556  }
3557 
3558  dotProduct = ( ret.at(it)[1] * allCs.at(iter)[1] +
3559  ret.at(it)[2] * allCs.at(iter)[2] +
3560  ret.at(it)[3] * allCs.at(iter)[3] );
3561 
3562  if ( ( ( 1.0 / sqrt ( 3.0 ) ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3563  ( ( 1.0 / sqrt ( 3.0 ) ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3564  {
3565  ;
3566  }
3567  else
3568  {
3569  allAnglesMet = false;
3570  }
3571  }
3572 
3573  if ( allAnglesMet )
3574  {
3575  ret.emplace_back ( allCs.at(iter) );
3576  }
3577  }
3578 
3579  //======================================== Not enough C3's. Try to find it
3580  if ( ret.size() != 7 )
3581  {
3582  if ( verbose > 3 )
3583  {
3584  std::cout << ">>>>>>>> Found only " << ret.size() - 3 << " C3 symmetries. Searching for the missing ones." << std::endl;
3585  }
3586 
3587  //==================================== If just one is missing, try to find it
3588  if ( ret.size() == 6 )
3589  {
3590  //================================ Find the missing axes
3591  std::vector< std::array<double,3> > axesVec;
3592  std::array<double,3> hArr;
3593  double axNorm, axX, axY, axZ;
3594  for ( double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
3595  {
3596  for ( double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
3597  {
3598  for ( double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
3599  {
3600  if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) { continue; }
3601 
3602  axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
3603  axX = xIt / axNorm;
3604  axY = yIt / axNorm;
3605  axZ = zIt / axNorm;
3606 
3607  allAnglesMet = true;
3608  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3609  {
3610  //================ Do not check the C3s against themselves
3611  if ( ret.at(it)[0] == 3.0 )
3612  {
3613  continue;
3614  }
3615 
3616  dotProduct = ( ret.at(it)[1] * axX +
3617  ret.at(it)[2] * axY +
3618  ret.at(it)[3] * axZ );
3619 
3620  if ( ( ( 1.0 / sqrt ( 3.0 ) ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3621  ( ( 1.0 / sqrt ( 3.0 ) ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3622  {
3623  ;
3624  }
3625  else
3626  {
3627  allAnglesMet = false;
3628  continue;
3629  }
3630  }
3631 
3632  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3633  {
3634  //================ Do check the C3s against themselves
3635  if ( ret.at(it)[0] == 4.0 )
3636  {
3637  continue;
3638  }
3639 
3640  dotProduct = ( ret.at(it)[1] * axX +
3641  ret.at(it)[2] * axY +
3642  ret.at(it)[3] * axZ );
3643 
3644  if ( ( 1.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3645  ( 1.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3646  {
3647  allAnglesMet = false;
3648  continue;
3649  }
3650  else
3651  {
3652  ;
3653  }
3654  }
3655 
3656  if ( allAnglesMet )
3657  {
3658  hArr[0] = axX;
3659  hArr[1] = axY;
3660  hArr[2] = axZ;
3661 
3662  allAnglesMet = true;
3663  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3664  {
3665  if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
3666  ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
3667  ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
3668  {
3669  allAnglesMet = false;
3670  }
3671  }
3672 
3673  if ( allAnglesMet )
3674  {
3675  axesVec.emplace_back ( hArr );
3676  }
3677  }
3678  }
3679  }
3680  }
3681 
3682  std::vector< std::array<double,5> > hRes;
3683  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3684  {
3685  //============================ Check if this axis solves the missing problem
3686  double res = cmpObj->maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 3.0, settings );
3687 
3688  hlpArr[0] = 3.0;
3689  hlpArr[1] = axesVec.at(axIt)[0];
3690  hlpArr[2] = axesVec.at(axIt)[1];
3691  hlpArr[3] = axesVec.at(axIt)[2];
3692  hlpArr[4] = res;
3693  hRes.emplace_back ( hlpArr );
3694  }
3695 
3696  std::sort ( hRes.begin(), hRes.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
3697 
3698  if ( hRes.at(0)[4] > ( cmpObj->getPeakThreshold ( ) ) )
3699  {
3700  ret.emplace_back ( hRes.at(0) );
3701  }
3702 
3703  if ( ret.size() != 7 )
3704  {
3705  ret.clear ( );
3706  return ( ret );
3707  }
3708  }
3709  else
3710  {
3711  ret.clear ( );
3712  return ( ret );
3713  }
3714  }
3715  if ( verbose > 3 )
3716  {
3717  std::cout << ">>>>>>>> All C3 symmetries located." << std::endl;
3718  }
3719 
3720  if ( verbose > 2 )
3721  {
3722  std::cout << ">>>>> Searching for the six C2 symmetry elements." << std::endl;
3723  }
3724 
3725  //======================================== Search for the 6 unique C2s with angles 1/sqrt(2) or 0.0 to the C4s
3726  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
3727  {
3728  //==================================== Search only using C2s
3729  if ( allCs.at(iter)[0] != 2.0 )
3730  {
3731  continue;
3732  }
3733 
3734  //==================================== Check that the prospective C2 has angle 1/sqrt(2) or 0.0 to the C4s
3735  allAnglesMet = true;
3736  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3737  {
3738  //====== Take only the C4s
3739  if ( ret.at(it)[0] != 4.0 )
3740  {
3741  continue;
3742  }
3743 
3744  dotProduct = ( ret.at(it)[1] * allCs.at(iter)[1] +
3745  ret.at(it)[2] * allCs.at(iter)[2] +
3746  ret.at(it)[3] * allCs.at(iter)[3] );
3747 
3748  if ( ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3749  ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) ) ||
3750  ( ( ( 1.0 / sqrt ( 2.0 ) ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3751  ( ( 1.0 / sqrt ( 2.0 ) ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) ) )
3752  {
3753  ;
3754  }
3755  else
3756  {
3757  allAnglesMet = false;
3758  }
3759  }
3760 
3761  if ( allAnglesMet )
3762  {
3763  ret.emplace_back ( allCs.at(iter) );
3764  }
3765  }
3766 
3767  //======================================== Not enough C2's. Try to find them
3768  if ( ret.size() != 13 )
3769  {
3770  if ( verbose > 3 )
3771  {
3772  std::cout << ">>>>>>>> Found only " << ret.size() - 13 << " C2 symmetries. Searching for the missing ones." << std::endl;
3773  }
3774 
3775  //==================================== If just one is missing, try to find it
3776  if ( ret.size() > 10 )
3777  {
3778  //================================ Find the missing axes
3779  std::vector< std::array<double,3> > axesVec;
3780  std::array<double,3> hArr;
3781  double axNorm, axX, axY, axZ;
3782  for ( double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
3783  {
3784  for ( double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
3785  {
3786  for ( double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
3787  {
3788  if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) { continue; }
3789 
3790  axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
3791  axX = xIt / axNorm;
3792  axY = yIt / axNorm;
3793  axZ = zIt / axNorm;
3794 
3795  allAnglesMet = true;
3796  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3797  {
3798  //================ Do not check the C3s and C4s
3799  if ( ( ret.at(it)[0] == 3.0 ) || ( ret.at(it)[0] == 4.0 ) )
3800  {
3801  continue;
3802  }
3803 
3804  dotProduct = ( ret.at(it)[1] * axX +
3805  ret.at(it)[2] * axY +
3806  ret.at(it)[3] * axZ );
3807 
3808  if ( ( ( ( 1.0 / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3809  ( ( 1.0 / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) ) ||
3810  ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3811  ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) ) )
3812  {
3813  ;
3814  }
3815  else
3816  {
3817  allAnglesMet = false;
3818  continue;
3819  }
3820  }
3821 
3822  if ( allAnglesMet )
3823  {
3824  hArr[0] = axX;
3825  hArr[1] = axY;
3826  hArr[2] = axZ;
3827 
3828  allAnglesMet = true;
3829  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3830  {
3831  if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
3832  ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
3833  ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
3834  {
3835  allAnglesMet = false;
3836  }
3837  }
3838 
3839  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( ret.size() ); axIt++ )
3840  {
3841  if ( ( ret.at(axIt)[0] == 3.0 ) || ( ret.at(axIt)[0] == 4.0 ) )
3842  {
3843  continue;
3844  }
3845 
3846  if ( ( ( ret.at(axIt)[1] > ( axX - axisErrorTolerance ) ) && ( ret.at(axIt)[1] < ( axX + axisErrorTolerance ) ) ) &&
3847  ( ( ret.at(axIt)[2] > ( axY - axisErrorTolerance ) ) && ( ret.at(axIt)[2] < ( axY + axisErrorTolerance ) ) ) &&
3848  ( ( ret.at(axIt)[3] > ( axZ - axisErrorTolerance ) ) && ( ret.at(axIt)[3] < ( axZ + axisErrorTolerance ) ) ) )
3849  {
3850  allAnglesMet = false;
3851  }
3852  }
3853 
3854  if ( allAnglesMet )
3855  {
3856  axesVec.emplace_back ( hArr );
3857  }
3858  }
3859  }
3860  }
3861  }
3862 
3863  std::vector< std::array<double,5> > hRes;
3864  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3865  {
3866  //============================ Check if this axis solves the missing problem
3867  double res = cmpObj->maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 2.0, settings );
3868 
3869  hlpArr[0] = 2.0;
3870  hlpArr[1] = axesVec.at(axIt)[0];
3871  hlpArr[2] = axesVec.at(axIt)[1];
3872  hlpArr[3] = axesVec.at(axIt)[2];
3873  hlpArr[4] = res;
3874  hRes.emplace_back ( hlpArr );
3875  }
3876 
3877  std::sort ( hRes.begin(), hRes.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
3878 
3879  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( hRes.size() ); axIt++ )
3880  {
3881  if ( ret.size () != 13 )
3882  {
3883  if ( hRes.at(axIt)[4] > ( cmpObj->getPeakThreshold ( ) ) )
3884  {
3885  ret.emplace_back ( hRes.at(axIt) );
3886  }
3887  }
3888  }
3889 
3890  if ( ret.size() != 13 )
3891  {
3892  ret.clear ( );
3893  return ( ret );
3894  }
3895  }
3896  else
3897  {
3898  ret.clear ( );
3899  return ( ret );
3900  }
3901  }
3902  if ( verbose > 3 )
3903  {
3904  std::cout << ">>>>>>>> All C2 symmetries located." << std::endl;
3905  }
3906 
3907  //======================================== Done
3908  return ( ret );
3909 
3910 }
3911 
3921 std::vector< std::array<double,5> > ProSHADE_internal::ProSHADE_symmetry::generateOctaElements ( std::vector< std::array<double,5> > symmAxes,
3922  ProSHADE::ProSHADE_settings* settings,
3923  int verbose )
3924 {
3925  //======================================== Initialise
3926  std::vector< std::array<double,5> > ret;
3927  std::array<double,5> hlpArr;
3928  if ( verbose > 1 )
3929  {
3930  std::cout << ">> Searching for octahedral symmetry elements." << std::endl;
3931  }
3932 
3933  //======================================== Identity element
3934  hlpArr[0] = 1.0;
3935  hlpArr[1] = 1.0;
3936  hlpArr[2] = 0.0;
3937  hlpArr[3] = 0.0;
3938  hlpArr[4] = 0.0;
3939  ret.emplace_back ( hlpArr );
3940  if ( verbose > 3 )
3941  {
3942  std::cout << ">>>>>>>> Identity element." << std::endl;
3943  }
3944 
3945  //======================================== 3 x 180 around each C4 axis clockwise
3946  for ( unsigned int iter = 0; iter < 3; iter++ )
3947  {
3948  hlpArr[0] = 4.0;
3949  hlpArr[1] = symmAxes.at(iter)[1];
3950  hlpArr[2] = symmAxes.at(iter)[2];
3951  hlpArr[3] = symmAxes.at(iter)[3];
3952  hlpArr[4] = 180.0;
3953  ret.emplace_back ( hlpArr );
3954  }
3955  if ( verbose > 3 )
3956  {
3957  std::cout << ">>>>>>>> 3 clockwise 180 C4 rotations." << std::endl;
3958  }
3959 
3960  //======================================== 6 x +-90 around each C4 axis clockwise
3961  for ( unsigned int iter = 0; iter < 3; iter++ )
3962  {
3963  hlpArr[0] = 4.0;
3964  hlpArr[1] = symmAxes.at(iter)[1];
3965  hlpArr[2] = symmAxes.at(iter)[2];
3966  hlpArr[3] = symmAxes.at(iter)[3];
3967  hlpArr[4] = 90.0;
3968  ret.emplace_back ( hlpArr );
3969 
3970  hlpArr[0] = 4.0;
3971  hlpArr[1] = symmAxes.at(iter)[1];
3972  hlpArr[2] = symmAxes.at(iter)[2];
3973  hlpArr[3] = symmAxes.at(iter)[3];
3974  hlpArr[4] = -90.0;
3975  ret.emplace_back ( hlpArr );
3976  }
3977  if ( verbose > 3 )
3978  {
3979  std::cout << ">>>>>>>> 6 (anti)clockwise +-90 C4 rotations." << std::endl;
3980  }
3981 
3982  //======================================== 8 x +-120 around each C3 axis clockwise
3983  for ( unsigned int iter = 3; iter < 7; iter++ )
3984  {
3985  hlpArr[0] = 3.0;
3986  hlpArr[1] = symmAxes.at(iter)[1];
3987  hlpArr[2] = symmAxes.at(iter)[2];
3988  hlpArr[3] = symmAxes.at(iter)[3];
3989  hlpArr[4] = 120.0;
3990  ret.emplace_back ( hlpArr );
3991 
3992  hlpArr[0] = 3.0;
3993  hlpArr[1] = symmAxes.at(iter)[1];
3994  hlpArr[2] = symmAxes.at(iter)[2];
3995  hlpArr[3] = symmAxes.at(iter)[3];
3996  hlpArr[4] = -120.0;
3997  ret.emplace_back ( hlpArr );
3998  }
3999  if ( verbose > 3 )
4000  {
4001  std::cout << ">>>>>>>> 8 (anti)clockwise +-120 C3 rotations." << std::endl;
4002  }
4003 
4004  //======================================== 6 x 180 around each C2 axis clockwise
4005  for ( unsigned int iter = 7; iter < 13; iter++ )
4006  {
4007  hlpArr[0] = 2.0;
4008  hlpArr[1] = symmAxes.at(iter)[1];
4009  hlpArr[2] = symmAxes.at(iter)[2];
4010  hlpArr[3] = symmAxes.at(iter)[3];
4011  hlpArr[4] = 180.0;
4012  ret.emplace_back ( hlpArr );
4013  }
4014  if ( verbose > 3 )
4015  {
4016  std::cout << ">>>>>>>> 6 clockwise 120 C2 rotations." << std::endl;
4017  }
4018 
4019  if ( verbose > 2 )
4020  {
4021  std::cout << ">>>>> Octahederal symmetry elements generated." << std::endl;
4022  }
4023 
4024  //======================================== Done
4025  return ( ret );
4026 
4027 }
4028 
4048  ProSHADE::ProSHADE_settings* settings,
4049  std::vector< std::array<double,5> > icosSymms,
4050  std::vector< std::array<double,5> > allCs,
4051  double axisErrorTolerance,
4052  int verbose )
4053 {
4054  //======================================== Initialise
4055  std::vector< std::array<double,5> > ret;
4056  std::array<double,5> hlpArr;
4057  bool allAnglesMet = false;
4058  double dotProduct = 0.0;
4059 
4060  //======================================== Push the already known C5 to ret. This is to allow different starting points, should the first one fail
4061  if ( icosSymms.at(0)[0] == 5.0 )
4062  {
4063  ret.emplace_back ( icosSymms.at(0) );
4064  }
4065  else
4066  {
4067  ret.emplace_back ( icosSymms.at(1) );
4068  }
4069 
4070  if ( verbose > 1 )
4071  {
4072  std::cout << ">> Searching for icosahedral symmetry elements." << std::endl;
4073  }
4074  if ( verbose > 2 )
4075  {
4076  std::cout << ">>>>> Searching for the six C5 symmetry elements." << std::endl;
4077  }
4078  //======================================== Search for the other five C5's with the angle cos(0.5) to each other
4079  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
4080  {
4081  //==================================== Search only using C5s
4082  if ( allCs.at(iter)[0] != 5.0 )
4083  {
4084  continue;
4085  }
4086 
4087  //==================================== Check the known C5's against the prospective one
4088  allAnglesMet = true;
4089  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
4090  {
4091  dotProduct = ( ret.at(it)[1] * allCs.at(iter)[1] +
4092  ret.at(it)[2] * allCs.at(iter)[2] +
4093  ret.at(it)[3] * allCs.at(iter)[3] );
4094 
4095  if ( ( ( 1.0 / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4096  ( ( 1.0 / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4097  {
4098  ;
4099  }
4100  else
4101  {
4102  allAnglesMet = false;
4103  }
4104  }
4105 
4106  if ( allAnglesMet )
4107  {
4108 
4109  bool alreadyExists = false;
4110  for ( unsigned int i = 0; i < static_cast<unsigned int> ( ret.size() ); i++ )
4111  {
4112  if ( ( ( ( ret.at(i)[1] + axisErrorTolerance ) < ( allCs.at(iter)[1] ) ) && ( ( ret.at(i)[1] - axisErrorTolerance ) > ( allCs.at(iter)[1] ) ) ) &&
4113  ( ( ( ret.at(i)[2] + axisErrorTolerance ) < ( allCs.at(iter)[2] ) ) && ( ( ret.at(i)[2] - axisErrorTolerance ) > ( allCs.at(iter)[2] ) ) ) &&
4114  ( ( ( ret.at(i)[3] + axisErrorTolerance ) < ( allCs.at(iter)[3] ) ) && ( ( ret.at(i)[3] - axisErrorTolerance ) > ( allCs.at(iter)[3] ) ) ) )
4115  {
4116  alreadyExists = true;
4117  }
4118  }
4119 
4120  if ( !alreadyExists )
4121  {
4122  ret.emplace_back ( allCs.at(iter) );
4123  }
4124 
4125  }
4126  }
4127 
4128  //======================================== Not enough C5's. Try to find the missing C5
4129  if ( ret.size() != 6 )
4130  {
4131  if ( verbose > 3 )
4132  {
4133  std::cout << ">>>>>>>> Found only " << ret.size() << " C5 symmetries. Searching for the missing ones." << std::endl;
4134  }
4135 
4136  //==================================== If few are missing, try to find them
4137  if ( ret.size() > 2 )
4138  {
4139  //================================ Find the missing axes
4140  std::vector< std::array<double,3> > axesVec;
4141  std::array<double,3> hArr;
4142  double axNorm, axX, axY, axZ;
4143  for ( double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
4144  {
4145  for ( double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
4146  {
4147  for ( double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
4148  {
4149  if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) { continue; }
4150 
4151  axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
4152  axX = xIt / axNorm;
4153  axY = yIt / axNorm;
4154  axZ = zIt / axNorm;
4155 
4156  allAnglesMet = true;
4157  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
4158  {
4159  dotProduct = ( ret.at(it)[1] * axX +
4160  ret.at(it)[2] * axY +
4161  ret.at(it)[3] * axZ );
4162 
4163  if ( ( ( 1.0 / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4164  ( ( 1.0 / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4165  {
4166  ;
4167  }
4168  else
4169  {
4170  allAnglesMet = false;
4171  continue;
4172  }
4173  }
4174 
4175  if ( allAnglesMet )
4176  {
4177  hArr[0] = axX;
4178  hArr[1] = axY;
4179  hArr[2] = axZ;
4180 
4181  allAnglesMet = true;
4182  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
4183  {
4184  if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
4185  ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
4186  ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
4187  {
4188  allAnglesMet = false;
4189  }
4190  }
4191 
4192  if ( allAnglesMet )
4193  {
4194  axesVec.emplace_back ( hArr );
4195  }
4196  }
4197  }
4198  }
4199  }
4200 
4201  std::vector< std::array<double,5> > hRes;
4202  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
4203  {
4204  //============================ Check if this axis solves the missing problem
4205  double res = cmpObj->maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 5.0, settings );
4206 
4207  hlpArr[0] = 5.0;
4208  hlpArr[1] = axesVec.at(axIt)[0];
4209  hlpArr[2] = axesVec.at(axIt)[1];
4210  hlpArr[3] = axesVec.at(axIt)[2];
4211  hlpArr[4] = res;
4212  hRes.emplace_back ( hlpArr );
4213  }
4214 
4215  std::sort ( hRes.begin(), hRes.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
4216 
4217  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( hRes.size() ); axIt++ )
4218  {
4219  if ( ret.size() != 6 )
4220  {
4221  if ( hRes.at(axIt)[4] > ( cmpObj->getPeakThreshold ( ) ) )
4222  {
4223  ret.emplace_back ( hRes.at(axIt) );
4224  }
4225  }
4226  }
4227 
4228  if ( ret.size() != 6 )
4229  {
4230  ret.clear ( );
4231  return ( ret );
4232  }
4233  }
4234  else
4235  {
4236  ret.clear ( );
4237  return ( ret );
4238  }
4239  }
4240  if ( verbose > 3 )
4241  {
4242  std::cout << ">>>>>>>> All C5 symmetries located." << std::endl;
4243  }
4244 
4245  if ( verbose > 2 )
4246  {
4247  std::cout << ">>>>> Searching for the ten C3 symmetry elements." << std::endl;
4248  }
4249 
4250  //======================================== Search for the ten C3's with the angle cos(-sqrt(5)/3) and 1.0 - cos(-sqrt(5)/3) to all the C5's
4251  double closeC5s = 0.0;
4252  double awayC5s = 0.0;
4253  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
4254  {
4255  //==================================== Search only using C3s
4256  if ( allCs.at(iter)[0] != 3.0 )
4257  {
4258  continue;
4259  }
4260 
4261  //==================================== Check the known C5's against the new C3
4262  allAnglesMet = true;
4263  closeC5s = 0.0;
4264  awayC5s = 0.0;
4265  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
4266  {
4267  //================================ Use only the C5's
4268  if ( ret.at(it)[0] != 5.0 )
4269  {
4270  continue;
4271  }
4272 
4273  dotProduct = ( ret.at(it)[1] * allCs.at(iter)[1] +
4274  ret.at(it)[2] * allCs.at(iter)[2] +
4275  ret.at(it)[3] * allCs.at(iter)[3] );
4276 
4277  if ( ( ( sqrt ( 5.0 ) / 3.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4278  ( ( sqrt ( 5.0 ) / 3.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4279  {
4280  closeC5s += 1.0;
4281  continue;
4282  }
4283  if ( ( ( 1.0 - ( sqrt ( 5.0 ) / 3.0 ) ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4284  ( ( 1.0 - ( sqrt ( 5.0 ) / 3.0 ) ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4285  {
4286  awayC5s += 1.0;
4287  continue;
4288  }
4289  }
4290 
4291  if ( ( closeC5s == 3.0 ) && ( awayC5s == 3.0 ) )
4292  {
4293  ret.emplace_back ( allCs.at(iter) );
4294  }
4295  }
4296 
4297  //======================================== Not enough C3's. Try to find the missing C3
4298  if ( ret.size() != 16 )
4299  {
4300  if ( verbose > 3 )
4301  {
4302  std::cout << ">>>>>>>> Found only " << ret.size() - 6 << " C3 symmetries. Searching for the missing ones." << std::endl;
4303  }
4304 
4305  //==================================== If few are missing, try to find them
4306  if ( ret.size() > 5 )
4307  {
4308  //====== Find the missing axes
4309  std::vector< std::array<double,3 >> axesVec;
4310  std::array<double,3> hArr;
4311  double axNorm, axX, axY, axZ;
4312  for ( double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
4313  {
4314  for ( double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
4315  {
4316  for ( double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
4317  {
4318  if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) { continue; }
4319 
4320  axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
4321  axX = xIt / axNorm;
4322  axY = yIt / axNorm;
4323  axZ = zIt / axNorm;
4324 
4325  closeC5s = 0.0;
4326  awayC5s = 0.0;
4327  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
4328  {
4329  //====== Use only the C5's
4330  if ( ret.at(it)[0] != 5.0 )
4331  {
4332  continue;
4333  }
4334 
4335  dotProduct = ( ret.at(it)[1] * axX +
4336  ret.at(it)[2] * axY +
4337  ret.at(it)[3] * axZ );
4338 
4339  if ( ( ( sqrt ( 5.0 ) / 3.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4340  ( ( sqrt ( 5.0 ) / 3.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4341  {
4342  closeC5s += 1.0;
4343  continue;
4344  }
4345  if ( ( ( 1.0 - ( sqrt ( 5.0 ) / 3.0 ) ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4346  ( ( 1.0 - ( sqrt ( 5.0 ) / 3.0 ) ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4347  {
4348  awayC5s += 1.0;
4349  continue;
4350  }
4351  }
4352 
4353  if ( ( closeC5s == 3.0 ) && ( awayC5s == 3.0 ) )
4354  {
4355  hArr[0] = axX;
4356  hArr[1] = axY;
4357  hArr[2] = axZ;
4358 
4359  allAnglesMet = true;
4360  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
4361  {
4362  if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
4363  ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
4364  ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
4365  {
4366  allAnglesMet = false;
4367  }
4368  }
4369 
4370  if ( allAnglesMet )
4371  {
4372  axesVec.emplace_back ( hArr );
4373  }
4374  }
4375  }
4376  }
4377  }
4378 
4379  std::vector< std::array<double,5> > hRes;
4380  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
4381  {
4382  //============================ Check if this axis solves the missing problem
4383  double res = cmpObj->maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 3.0, settings );
4384 
4385  hlpArr[0] = 3.0;
4386  hlpArr[1] = axesVec.at(axIt)[0];
4387  hlpArr[2] = axesVec.at(axIt)[1];
4388  hlpArr[3] = axesVec.at(axIt)[2];
4389  hlpArr[4] = res;
4390  hRes.emplace_back ( hlpArr );
4391  }
4392 
4393  std::sort ( hRes.begin(), hRes.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
4394 
4395  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( hRes.size() ); axIt++ )
4396  {
4397  if ( ret.size() != 16 )
4398  {
4399  if ( hRes.at(axIt)[4] > ( cmpObj->getPeakThreshold ( ) ) )
4400  {
4401  ret.emplace_back ( hRes.at(axIt) );
4402  }
4403  }
4404  }
4405 
4406  if ( ret.size() != 16 )
4407  {
4408  ret.clear ( );
4409  return ( ret );
4410  }
4411  }
4412  else
4413  {
4414  ret.clear ( );
4415  return ( ret );
4416  }
4417  }
4418  if ( verbose > 3 )
4419  {
4420  std::cout << ">>>>>>>> All C3 symmetries located." << std::endl;
4421  }
4422 
4423  if ( verbose > 2 )
4424  {
4425  std::cout << ">>>>> Searching for the fifteen C2 symmetry elements." << std::endl;
4426  }
4427 
4428  //======================================== Search for the fifteen C2's with the angle cos(0), cos(0.5) and cos(sqrt(3)/2) to all the C5's
4429  double perpC2s = 0.0;
4430  double sixthC2s = 0.0;
4431  double halfSixtbC2s = 0.0;
4432  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
4433  {
4434  //==================================== Search only using C2s
4435  if ( allCs.at(iter)[0] != 2.0 )
4436  {
4437  continue;
4438  }
4439 
4440  //==================================== Check the known C5's against the new C2
4441  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
4442  {
4443  //================================ Use only the C5's
4444  if ( ret.at(it)[0] != 5.0 )
4445  {
4446  continue;
4447  }
4448 
4449  dotProduct = ( ret.at(it)[1] * allCs.at(iter)[1] +
4450  ret.at(it)[2] * allCs.at(iter)[2] +
4451  ret.at(it)[3] * allCs.at(iter)[3] );
4452 
4453  if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4454  ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4455  {
4456  perpC2s += 1.0;
4457  continue;
4458  }
4459  if ( ( ( 1.0 / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4460  ( ( 1.0 / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4461  {
4462  sixthC2s += 1.0;
4463  continue;
4464  }
4465  if ( ( ( sqrt ( 3.0 ) / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4466  ( ( sqrt ( 3.0 ) / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4467  {
4468  halfSixtbC2s += 1.0;
4469  continue;
4470  }
4471  }
4472 
4473  if ( ( perpC2s == 2.0 ) && ( sixthC2s == 2.0 ) && ( halfSixtbC2s == 2.0 ) )
4474  {
4475  ret.emplace_back ( allCs.at(iter) );
4476  }
4477  }
4478 
4479  //======================================== Not enough C3's. Try to find the missing C3
4480  if ( ret.size() != 31 )
4481  {
4482  if ( verbose > 3 )
4483  {
4484  std::cout << ">>>>>>>> Found only " << ret.size() - 16 << " C2 symmetries. Searching for the missing ones." << std::endl;
4485  }
4486 
4487  //==================================== If few are missing, try to find them
4488  if ( ret.size() > 15 )
4489  {
4490  //====== Find the missing axes
4491  std::vector< std::array<double,3> > axesVec;
4492  std::array<double,3> hArr;
4493  double axNorm, axX, axY, axZ;
4494  for ( double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
4495  {
4496  for ( double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
4497  {
4498  for ( double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
4499  {
4500  if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) { continue; }
4501 
4502  axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
4503  axX = xIt / axNorm;
4504  axY = yIt / axNorm;
4505  axZ = zIt / axNorm;
4506 
4507  perpC2s = 0.0;
4508  sixthC2s = 0.0;
4509  halfSixtbC2s = 0.0;
4510  for ( unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
4511  {
4512  //================ Use only the C5's
4513  if ( ret.at(it)[0] != 5.0 )
4514  {
4515  continue;
4516  }
4517 
4518  dotProduct = ( ret.at(it)[1] * axX +
4519  ret.at(it)[2] * axY +
4520  ret.at(it)[3] * axZ );
4521 
4522  if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4523  ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4524  {
4525  perpC2s += 1.0;
4526  continue;
4527  }
4528  if ( ( ( 1.0 / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4529  ( ( 1.0 / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4530  {
4531  sixthC2s += 1.0;
4532  continue;
4533  }
4534  if ( ( ( sqrt ( 3.0 ) / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4535  ( ( sqrt ( 3.0 ) / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4536  {
4537  halfSixtbC2s += 1.0;
4538  continue;
4539  }
4540  }
4541 
4542  if ( ( perpC2s == 2.0 ) && ( sixthC2s == 2.0 ) && ( halfSixtbC2s == 2.0 ) )
4543  {
4544  hArr[0] = axX;
4545  hArr[1] = axY;
4546  hArr[2] = axZ;
4547 
4548  allAnglesMet = true;
4549  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
4550  {
4551  if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
4552  ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
4553  ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
4554  {
4555  allAnglesMet = false;
4556  }
4557  }
4558 
4559  if ( allAnglesMet )
4560  {
4561  axesVec.emplace_back ( hArr );
4562  }
4563  }
4564  }
4565  }
4566  }
4567 
4568  std::vector< std::array<double,5> > hRes;
4569  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
4570  {
4571  //============================ Check if this axis solves the missing problem
4572  double res = cmpObj->maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 2.0, settings );
4573 
4574  hlpArr[0] = 2.0;
4575  hlpArr[1] = axesVec.at(axIt)[0];
4576  hlpArr[2] = axesVec.at(axIt)[1];
4577  hlpArr[3] = axesVec.at(axIt)[2];
4578  hlpArr[4] = res;
4579  hRes.emplace_back ( hlpArr );
4580  }
4581 
4582  std::sort ( hRes.begin(), hRes.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
4583 
4584  for ( unsigned int axIt = 0; axIt < static_cast<unsigned int> ( hRes.size() ); axIt++ )
4585  {
4586  if ( ret.size() != 31 )
4587  {
4588  if ( hRes.at(axIt)[4] > ( cmpObj->getPeakThreshold ( ) ) )
4589  {
4590  ret.emplace_back ( hRes.at(axIt) );
4591  }
4592  }
4593  }
4594 
4595  if ( ret.size() != 31 )
4596  {
4597  ret.clear ( );
4598  return ( ret );
4599  }
4600  }
4601  else
4602  {
4603  ret.clear ( );
4604  return ( ret );
4605  }
4606  }
4607  if ( verbose > 3 )
4608  {
4609  std::cout << ">>>>>>>> All C2 symmetries located." << std::endl;
4610  }
4611 
4612  //======================================== Done
4613  return ( ret );
4614 
4615 }
4616 
4626 std::vector< std::array<double,5> > ProSHADE_internal::ProSHADE_symmetry::generateIcosElements ( std::vector< std::array<double,5> > symmAxes,
4627  ProSHADE::ProSHADE_settings* settings,
4628  int verbose )
4629 {
4630  //======================================== Initialise
4631  std::vector< std::array<double,5> > ret;
4632  std::array<double,5> hlpArr;
4633  if ( verbose > 1 )
4634  {
4635  std::cout << ">> Searching for icosahedral symmetry elements." << std::endl;
4636  }
4637 
4638  //======================================== Identity element
4639  hlpArr[0] = 1.0;
4640  hlpArr[1] = 1.0;
4641  hlpArr[2] = 0.0;
4642  hlpArr[3] = 0.0;
4643  hlpArr[4] = 0.0;
4644  ret.emplace_back ( hlpArr );
4645  if ( verbose > 3 )
4646  {
4647  std::cout << ">>>>>>>> Identity element." << std::endl;
4648  }
4649 
4650  //======================================== 12 x +-72 around each C5 axis (anti)clockwise
4651  for ( unsigned int iter = 0; iter < 6; iter++ )
4652  {
4653  hlpArr[0] = 5.0;
4654  hlpArr[1] = symmAxes.at(iter)[1];
4655  hlpArr[2] = symmAxes.at(iter)[2];
4656  hlpArr[3] = symmAxes.at(iter)[3];
4657  hlpArr[4] = 72.0;
4658  ret.emplace_back ( hlpArr );
4659 
4660  hlpArr[0] = 5.0;
4661  hlpArr[1] = symmAxes.at(iter)[1];
4662  hlpArr[2] = symmAxes.at(iter)[2];
4663  hlpArr[3] = symmAxes.at(iter)[3];
4664  hlpArr[4] = -72.0;
4665  ret.emplace_back ( hlpArr );
4666  }
4667  if ( verbose > 3 )
4668  {
4669  std::cout << ">>>>>>>> 12 (anti)clockwise 144 C5 rotations." << std::endl;
4670  }
4671 
4672  //======================================== 12 x +-72 around each C5 axis (anti)clockwise
4673  for ( unsigned int iter = 0; iter < 6; iter++ )
4674  {
4675  hlpArr[0] = 5.0;
4676  hlpArr[1] = symmAxes.at(iter)[1];
4677  hlpArr[2] = symmAxes.at(iter)[2];
4678  hlpArr[3] = symmAxes.at(iter)[3];
4679  hlpArr[4] = 144.0;
4680  ret.emplace_back ( hlpArr );
4681 
4682  hlpArr[0] = 5.0;
4683  hlpArr[1] = symmAxes.at(iter)[1];
4684  hlpArr[2] = symmAxes.at(iter)[2];
4685  hlpArr[3] = symmAxes.at(iter)[3];
4686  hlpArr[4] = -144.0;
4687  ret.emplace_back ( hlpArr );
4688  }
4689  if ( verbose > 3 )
4690  {
4691  std::cout << ">>>>>>>> 12 (anti)clockwise 144 C5 rotations." << std::endl;
4692  }
4693 
4694  //======================================== 20 x +-120 around each C3 axis (anti)clockwise
4695  for ( unsigned int iter = 6; iter < 16; iter++ )
4696  {
4697  hlpArr[0] = 3.0;
4698  hlpArr[1] = symmAxes.at(iter)[1];
4699  hlpArr[2] = symmAxes.at(iter)[2];
4700  hlpArr[3] = symmAxes.at(iter)[3];
4701  hlpArr[4] = 120.0;
4702  ret.emplace_back ( hlpArr );
4703 
4704  hlpArr[0] = 3.0;
4705  hlpArr[1] = symmAxes.at(iter)[1];
4706  hlpArr[2] = symmAxes.at(iter)[2];
4707  hlpArr[3] = symmAxes.at(iter)[3];
4708  hlpArr[4] = -120.0;
4709  ret.emplace_back ( hlpArr );
4710  }
4711  if ( verbose > 3 )
4712  {
4713  std::cout << ">>>>>>>> 20 (anti)clockwise 120 C3 rotations." << std::endl;
4714  }
4715 
4716  //======================================== 15 x 180 around each C2 axis clockwise
4717  for ( unsigned int iter = 16; iter < 31; iter++ )
4718  {
4719  hlpArr[0] = 2.0;
4720  hlpArr[1] = symmAxes.at(iter)[1];
4721  hlpArr[2] = symmAxes.at(iter)[2];
4722  hlpArr[3] = symmAxes.at(iter)[3];
4723  hlpArr[4] = 180.0;
4724  ret.emplace_back ( hlpArr );
4725  }
4726  if ( verbose > 3 )
4727  {
4728  std::cout << ">>>>>>>> 15 clockwise 180 C2 rotations." << std::endl;
4729  }
4730 
4731  if ( verbose > 2 )
4732  {
4733  std::cout << ">>>>> Icosahedral symmetry elements generated." << std::endl;
4734  }
4735 
4736  //======================================== Done
4737  return ( ret );
4738 
4739 }
4740 
4755  ProSHADE_data *obj2,
4756  double *ob2XMov,
4757  double *ob2YMov,
4758  double *ob2ZMov )
4759 {
4760  //======================================== Initalise
4761  std::array<double,4> ret = std::array<double,4> { { 0.0, 0.0, 0.0, 0.0 } };
4762  int arrPos = 0;
4763  int hlpPos = 0;
4764 
4765  double newXSampling = std::min ( obj1->_xSamplingRate, obj2->_xSamplingRate );
4766  double newYSampling = std::min ( obj1->_ySamplingRate, obj2->_ySamplingRate );
4767  double newZSampling = std::min ( obj1->_zSamplingRate, obj2->_zSamplingRate );
4768 
4769  int newU = 0;
4770  int newV = 0;
4771  int newW = 0;
4772 
4773  //======================================== For structure 1, get the sampling right
4774  {
4775  //================================ Find cell parameters
4776  newU = std::floor ( obj1->_xRange / newXSampling );
4777  newV = std::floor ( obj1->_yRange / newYSampling );
4778  newW = std::floor ( obj1->_zRange / newZSampling );
4779 
4780  double newXRange = static_cast<double> ( ( newU ) * newXSampling );
4781  double newYRange = static_cast<double> ( ( newV ) * newYSampling );
4782  double newZRange = static_cast<double> ( ( newW ) * newZSampling );
4783 
4784  int uDiff = newU - obj1->_maxMapU;
4785  int vDiff = newV - obj1->_maxMapV;
4786  int wDiff = newW - obj1->_maxMapW;
4787 
4788  int newXTo = 0;
4789  int newXFrom = 0;
4790  int newYTo = 0;
4791  int newYFrom = 0;
4792  int newZTo = 0;
4793  int newZFrom = 0;
4794 
4795  if ( uDiff % 2 == 0 )
4796  {
4797  newXTo = obj1->_xTo + (uDiff/2);
4798  newXFrom = obj1->_xFrom - (uDiff/2);
4799  }
4800  else
4801  {
4802  newXTo = obj1->_xTo + ((uDiff+1)/2);
4803  newXFrom = obj1->_xFrom - ((uDiff+1)/2);
4804  newU += 1;
4805  }
4806  if ( vDiff % 2 == 0 )
4807  {
4808  newYTo = obj1->_yTo + (vDiff/2);
4809  newYFrom = obj1->_yFrom - (vDiff/2);
4810  }
4811  else
4812  {
4813  newYTo = obj1->_yTo + ((vDiff+1)/2);
4814  newYFrom = obj1->_yFrom - ((vDiff+1)/2);
4815  newV += 1;
4816  }
4817  if ( wDiff % 2 == 0 )
4818  {
4819  newZTo = obj1->_zTo + (wDiff/2);
4820  newZFrom = obj1->_zFrom - (wDiff/2);
4821  }
4822  else
4823  {
4824  newZTo = obj1->_zTo + ((wDiff+1)/2);
4825  newZFrom = obj1->_zFrom - ((wDiff+1)/2);
4826  newW += 1;
4827  }
4828 
4829  if ( ( ( obj1->_xFrom >= 0 ) && ( newXFrom < 0 ) ) &&
4830  ( ( obj1->_xTo >= 0 ) && ( newXTo > 0 ) ) )
4831  {
4832  newU += 1;
4833  }
4834 
4835  if ( ( ( obj1->_yFrom >= 0 ) && ( newYFrom < 0 ) ) &&
4836  ( ( obj1->_yTo >= 0 ) && ( newYTo > 0 ) ) )
4837  {
4838  newV += 1;
4839  }
4840 
4841  if ( ( ( obj1->_zFrom >= 0 ) && ( newZFrom < 0 ) ) &&
4842  ( ( obj1->_zTo >= 0 ) && ( newZTo > 0 ) ) )
4843  {
4844  newW += 1;
4845  }
4846 
4847  //================================ Make sure map dimmensions are even
4848  if ( newU % 2 != 0 )
4849  {
4850  newU += 1;
4851  newXTo += 1;
4852  }
4853 
4854  if ( newV % 2 != 0 )
4855  {
4856  newV += 1;
4857  newYTo += 1;
4858  }
4859 
4860  if ( newW % 2 != 0 )
4861  {
4862  newW += 1;
4863  newZTo += 1;
4864  }
4865 
4866  //================================ Interpolate map
4867  double *hlpMap = new double [(newU+1) * (newV+1) * (newW+1)];
4868  double newX = 0.0;
4869  double newY = 0.0;
4870  double newZ = 0.0;
4871  int oldX = 0;
4872  int oldY = 0;
4873  int oldZ = 0;
4874  bool xOver = false;
4875  bool yOver = false;
4876  bool zOver = false;
4877  bool xFound = false;
4878  bool yFound = false;
4879  bool zFound = false;
4880  double xd, yd, zd;
4881  std::array<double,4> p000, p001, p010, p011, p100, p101, p110, p111;
4882  std::array<double,3> p00, p01, p10, p11;
4883  std::array<double,2> p0, p1;
4884  for ( int uIt = 0; uIt < (newU+1); uIt++ )
4885  {
4886  for ( int vIt = 0; vIt < (newV+1); vIt++ )
4887  {
4888  for ( int wIt = 0; wIt < (newW+1); wIt++ )
4889  {
4890  //==================== Init
4891  arrPos = wIt + (newW+1) * ( vIt + (newV+1) * uIt );
4892 
4893  xFound = false;
4894  yFound = false;
4895  zFound = false;
4896 
4897  //==================== Get new X, Y and Z in angstroms
4898  newX = static_cast<double> ( uIt ) * newXSampling;
4899  newY = static_cast<double> ( vIt ) * newYSampling;
4900  newZ = static_cast<double> ( wIt ) * newZSampling;
4901 
4902  //==================== Find corresponding old X, Y and Z indices
4903  for ( int iter = 0; iter < static_cast<int> (obj1->_maxMapU+1); iter++ )
4904  {
4905  if ( ( newX >= ( static_cast<double> ( iter ) * obj1->_xSamplingRate ) ) && ( newX < ( static_cast<double> ( iter + 1 ) * obj1->_xSamplingRate ) ) )
4906  {
4907  oldX = iter;
4908  xFound = true;
4909  break;
4910  }
4911  }
4912  xOver = false;
4913  if ( ( oldX + 1 ) >= static_cast<int> (obj1->_maxMapU+1) ) { xOver = true; }
4914  if ( !xFound ) { hlpMap[arrPos] = 0.0; continue; }
4915 
4916  for ( int iter = 0; iter < static_cast<int> (obj1->_maxMapV+1); iter++ )
4917  {
4918  if ( ( newY >= ( static_cast<double> ( iter ) * obj1->_ySamplingRate ) ) && ( newY < ( static_cast<double> ( iter + 1 ) * obj1->_ySamplingRate ) ) )
4919  {
4920  oldY = iter;
4921  yFound = true;
4922  break;
4923  }
4924  }
4925  yOver = false;
4926  if ( ( oldY + 1 ) >= static_cast<int> (obj1->_maxMapV+1) ) { yOver = true; }
4927  if ( !yFound ) { hlpMap[arrPos] = 0.0; continue; }
4928 
4929  for ( int iter = 0; iter < static_cast<int> (obj1->_maxMapW+1); iter++ )
4930  {
4931  if ( ( newZ >= ( static_cast<double> ( iter ) * obj1->_zSamplingRate ) ) && ( newZ < ( static_cast<double> ( iter + 1 ) * obj1->_zSamplingRate ) ) )
4932  {
4933  oldZ = iter;
4934  zFound = true;
4935  break;
4936  }
4937  }
4938  zOver = false;
4939  if ( ( oldZ + 1 ) >= static_cast<int> (obj1->_maxMapW+1) ) { zOver = true; }
4940  if ( !zFound ) { hlpMap[arrPos] = 0.0; continue; }
4941 
4942  //==================== Get the surrounding values and distances
4943  p000[0] = oldX;
4944  p000[1] = oldY;
4945  p000[2] = oldZ;
4946  hlpPos = p000[2] + (obj1->_maxMapW+1) * ( p000[1] + (obj1->_maxMapV+1) * p000[0] );
4947  p000[3] = obj1->_densityMapCor[hlpPos];
4948 
4949  p001[0] = oldX;
4950  p001[1] = oldY;
4951  p001[2] = oldZ + 1;
4952  if ( zOver ) { p001[2] = 0.0; }
4953  hlpPos = p001[2] + (obj1->_maxMapW+1) * ( p001[1] + (obj1->_maxMapV+1) * p001[0] );
4954  p001[3] = obj1->_densityMapCor[hlpPos];
4955 
4956  p010[0] = oldX;
4957  p010[1] = oldY + 1;
4958  p010[2] = oldZ;
4959  if ( yOver ) { p010[1] = 0.0; }
4960  hlpPos = p010[2] + (obj1->_maxMapW+1) * ( p010[1] + (obj1->_maxMapV+1) * p010[0] );
4961  p010[3] = obj1->_densityMapCor[hlpPos];
4962 
4963  p011[0] = oldX;
4964  p011[1] = oldY + 1;
4965  p011[2] = oldZ + 1;
4966  if ( yOver ) { p011[1] = 0.0; }
4967  if ( zOver ) { p011[2] = 0.0; }
4968  hlpPos = p011[2] + (obj1->_maxMapW+1) * ( p011[1] + (obj1->_maxMapV+1) * p011[0] );
4969  p011[3] = obj1->_densityMapCor[hlpPos];
4970 
4971  p100[0] = oldX + 1;
4972  p100[1] = oldY;
4973  p100[2] = oldZ;
4974  if ( xOver ) { p100[0] = 0.0; }
4975  hlpPos = p100[2] + (obj1->_maxMapW+1) * ( p100[1] + (obj1->_maxMapV+1) * p100[0] );
4976  p100[3] = obj1->_densityMapCor[hlpPos];
4977 
4978  p101[0] = oldX + 1;
4979  p101[1] = oldY;
4980  p101[2] = oldZ + 1;
4981  if ( xOver ) { p101[0] = 0.0; }
4982  if ( zOver ) { p101[2] = 0.0; }
4983  hlpPos = p101[2] + (obj1->_maxMapW+1) * ( p101[1] + (obj1->_maxMapV+1) * p101[0] );
4984  p101[3] = obj1->_densityMapCor[hlpPos];
4985 
4986  p110[0] = oldX + 1;
4987  p110[1] = oldY + 1;
4988  p110[2] = oldZ;
4989  if ( xOver ) { p110[0] = 0.0; }
4990  if ( yOver ) { p110[1] = 0.0; }
4991  hlpPos = p110[2] + (obj1->_maxMapW+1) * ( p110[1] + (obj1->_maxMapV+1) * p110[0] );
4992  p110[3] = obj1->_densityMapCor[hlpPos];
4993 
4994  p111[0] = oldX + 1;
4995  p111[1] = oldY + 1;
4996  p111[2] = oldZ + 1;
4997  if ( xOver ) { p111[0] = 0.0; }
4998  if ( yOver ) { p111[1] = 0.0; }
4999  if ( zOver ) { p111[2] = 0.0; }
5000  hlpPos = p111[2] + (obj1->_maxMapW+1) * ( p111[1] + (obj1->_maxMapV+1) * p111[0] );
5001  p111[3] = obj1->_densityMapCor[hlpPos];
5002 
5003  //==================== Interpolate
5004  xd = 1.0 - ( ( newX - ( static_cast<double> ( oldX ) * obj1->_xSamplingRate ) ) / obj1->_xSamplingRate );
5005  p00[0] = p000[1]; p00[1] = p000[2]; p00[2] = ( xd * p000[3] ) + ( (1.0 - xd) * p100[3] );
5006  p01[0] = p001[1]; p01[1] = p001[2]; p01[2] = ( xd * p001[3] ) + ( (1.0 - xd) * p101[3] );
5007  p10[0] = p010[1]; p10[1] = p010[2]; p10[2] = ( xd * p010[3] ) + ( (1.0 - xd) * p110[3] );
5008  p11[0] = p011[1]; p11[1] = p011[2]; p11[2] = ( xd * p011[3] ) + ( (1.0 - xd) * p111[3] );
5009 
5010  yd = 1.0 - ( ( newY - ( static_cast<double> ( oldY ) * obj1->_ySamplingRate ) ) / obj1->_ySamplingRate );
5011  p0[0] = p00[1]; p0[1] = ( yd * p00[2] ) + ( (1.0 - yd) * p10[2] );
5012  p1[0] = p01[1]; p1[1] = ( yd * p01[2] ) + ( (1.0 - yd) * p11[2] );
5013 
5014  zd = 1.0 - ( ( newZ - ( static_cast<double> ( oldZ ) * obj1->_zSamplingRate ) ) / obj1->_zSamplingRate );
5015  hlpMap[arrPos] = ( zd * p0[1] ) + ( (1.0 - zd) * p1[1] );
5016  }
5017  }
5018  }
5019 
5020  //================================ Copy map
5021  delete[] obj1->_densityMapCor;
5022  obj1->_densityMapCor = new double [(newU+1) * (newV+1) * (newW+1)];
5023  for ( int iter = 0; iter < static_cast<int> ( (newU+1) * (newV+1) * (newW+1) ); iter++ )
5024  {
5025  obj1->_densityMapCor[iter] = hlpMap[iter];
5026  }
5027 
5028  //================================ Free memory
5029  delete[] hlpMap;
5030 
5031  //================================ Set new cell parameters
5032  int xMov1 = std::round ( ( ( obj1->_xFrom * obj1->_xSamplingRate ) - ( newXFrom * newXSampling ) ) / newXSampling );
5033  int yMov1 = std::round ( ( ( obj1->_yFrom * obj1->_ySamplingRate ) - ( newYFrom * newYSampling ) ) / newYSampling );
5034  int zMov1 = std::round ( ( ( obj1->_zFrom * obj1->_zSamplingRate ) - ( newZFrom * newZSampling ) ) / newZSampling );
5035 
5036  obj1->_maxMapU = newU;
5037  obj1->_maxMapV = newV;
5038  obj1->_maxMapW = newW;
5039 
5040  obj1->_xFrom = newXFrom + xMov1;
5041  obj1->_yFrom = newYFrom + yMov1;
5042  obj1->_zFrom = newZFrom + zMov1;
5043 
5044  obj1->_xTo = newXTo + xMov1;
5045  obj1->_yTo = newYTo + yMov1;
5046  obj1->_zTo = newZTo + zMov1;
5047 
5048  obj1->_xRange = newXRange;
5049  obj1->_yRange = newYRange;
5050  obj1->_zRange = newZRange;
5051 
5052  obj1->_xSamplingRate = newXSampling;
5053  obj1->_ySamplingRate = newYSampling;
5054  obj1->_zSamplingRate = newZSampling;
5055  }
5056 
5057  //======================================== For structure 2, get the sampling right
5058  {
5059  //================================ Find cell parameters
5060  newU = std::floor ( obj2->_xRange / newXSampling );
5061  newV = std::floor ( obj2->_yRange / newYSampling );
5062  newW = std::floor ( obj2->_zRange / newZSampling );
5063 
5064  double newXRange = static_cast<double> ( ( newU + 1 ) * newXSampling );
5065  double newYRange = static_cast<double> ( ( newV + 1 ) * newYSampling );
5066  double newZRange = static_cast<double> ( ( newW + 1 ) * newZSampling );
5067 
5068  int uDiff = newU - obj2->_maxMapU;
5069  int vDiff = newV - obj2->_maxMapV;
5070  int wDiff = newW - obj2->_maxMapW;
5071 
5072  int newXTo = 0;
5073  int newXFrom = 0;
5074  int newYTo = 0;
5075  int newYFrom = 0;
5076  int newZTo = 0;
5077  int newZFrom = 0;
5078 
5079  if ( uDiff % 2 == 0 )
5080  {
5081  newXTo = obj2->_xTo + (uDiff/2);
5082  newXFrom = obj2->_xFrom - (uDiff/2);
5083  }
5084  else
5085  {
5086  newXTo = obj2->_xTo + ((uDiff+1)/2);
5087  newXFrom = obj2->_xFrom - ((uDiff+1)/2);
5088  newU += 1;
5089  }
5090  if ( vDiff % 2 == 0 )
5091  {
5092  newYTo = obj2->_yTo + (vDiff/2);
5093  newYFrom = obj2->_yFrom - (vDiff/2);
5094  }
5095  else
5096  {
5097  newYTo = obj2->_yTo + ((vDiff+1)/2);
5098  newYFrom = obj2->_yFrom - ((vDiff+1)/2);
5099  newV += 1;
5100  }
5101  if ( wDiff % 2 == 0 )
5102  {
5103  newZTo = obj2->_zTo + (wDiff/2);
5104  newZFrom = obj2->_zFrom - (wDiff/2);
5105  }
5106  else
5107  {
5108  newZTo = obj2->_zTo + ((wDiff+1)/2);
5109  newZFrom = obj2->_zFrom - ((wDiff+1)/2);
5110  newW += 1;
5111  }
5112 
5113  if ( ( ( obj2->_xFrom >= 0 ) && ( newXFrom < 0 ) ) &&
5114  ( ( obj2->_xTo >= 0 ) && ( newXTo > 0 ) ) )
5115  {
5116  newU += 1;
5117  }
5118 
5119  if ( ( ( obj2->_yFrom >= 0 ) && ( newYFrom < 0 ) ) &&
5120  ( ( obj2->_yTo >= 0 ) && ( newYTo > 0 ) ) )
5121  {
5122  newV += 1;
5123  }
5124 
5125  if ( ( ( obj2->_zFrom >= 0 ) && ( newZFrom < 0 ) ) &&
5126  ( ( obj2->_zTo >= 0 ) && ( newZTo > 0 ) ) )
5127  {
5128  newW += 1;
5129  }
5130 
5131  //================================ Make sure map dimmensions are even
5132  if ( newU % 2 == 0 )
5133  {
5134  newU += 1;
5135  newXTo += 1;
5136  }
5137 
5138  if ( newV % 2 == 0 )
5139  {
5140  newV += 1;
5141  newYTo += 1;
5142  }
5143 
5144  if ( newW % 2 == 0 )
5145  {
5146  newW += 1;
5147  newZTo += 1;
5148  }
5149 
5150  //================================ Interpolate map
5151  double *hlpMap = new double [(newU+1) * (newV+1) * (newW+1)];
5152  double newX = 0.0;
5153  double newY = 0.0;
5154  double newZ = 0.0;
5155  int oldX = 0;
5156  int oldY = 0;
5157  int oldZ = 0;
5158  bool xOver = false;
5159  bool yOver = false;
5160  bool zOver = false;
5161  bool xFound = false;
5162  bool yFound = false;
5163  bool zFound = false;
5164  double xd, yd, zd;
5165  std::array<double,4> p000, p001, p010, p011, p100, p101, p110, p111;
5166  std::array<double,3> p00, p01, p10, p11;
5167  std::array<double,2> p0, p1;
5168 
5169  for ( int uIt = 0; uIt < (newU+1); uIt++ )
5170  {
5171  for ( int vIt = 0; vIt < (newV+1); vIt++ )
5172  {
5173  for ( int wIt = 0; wIt < (newW+1); wIt++ )
5174  {
5175  //==================== Init
5176  arrPos = wIt + (newW+1) * ( vIt + (newV+1) * uIt );
5177 
5178  xFound = false;
5179  yFound = false;
5180  zFound = false;
5181 
5182  //==================== Get new X, Y and Z in angstroms
5183  newX = static_cast<double> ( uIt ) * newXSampling;
5184  newY = static_cast<double> ( vIt ) * newYSampling;
5185  newZ = static_cast<double> ( wIt ) * newZSampling;
5186 
5187  //==================== Find corresponding old X, Y and Z indices
5188  for ( int iter = 0; iter < static_cast<int> (obj2->_maxMapU+1); iter++ )
5189  {
5190  if ( ( newX >= ( static_cast<double> ( iter ) * obj2->_xSamplingRate ) ) && ( newX < ( static_cast<double> ( iter + 1 ) * obj2->_xSamplingRate ) ) )
5191  {
5192  oldX = iter;
5193  xFound = true;
5194  break;
5195  }
5196  }
5197  xOver = false;
5198  if ( ( oldX + 1 ) >= static_cast<int> (obj2->_maxMapU+1) ) { xOver = true; }
5199  if ( !xFound ) { hlpMap[arrPos] = 0.0; continue; }
5200 
5201  for ( int iter = 0; iter < static_cast<int> (obj2->_maxMapV+1); iter++ )
5202  {
5203  if ( ( newY >= ( static_cast<double> ( iter ) * obj2->_ySamplingRate ) ) && ( newY < ( static_cast<double> ( iter + 1 ) * obj2->_ySamplingRate ) ) )
5204  {
5205  oldY = iter;
5206  yFound = true;
5207  break;
5208  }
5209  }
5210  yOver = false;
5211  if ( ( oldY + 1 ) >= static_cast<int> (obj2->_maxMapV+1) ) { yOver = true; }
5212  if ( !yFound ) { hlpMap[arrPos] = 0.0; continue; }
5213 
5214  for ( int iter = 0; iter < static_cast<int> (obj2->_maxMapW+1); iter++ )
5215  {
5216  if ( ( newZ >= ( static_cast<double> ( iter ) * obj2->_zSamplingRate ) ) && ( newZ < ( static_cast<double> ( iter + 1 ) * obj2->_zSamplingRate ) ) )
5217  {
5218  oldZ = iter;
5219  zFound = true;
5220  break;
5221  }
5222  }
5223  zOver = false;
5224  if ( ( oldZ + 1 ) >= static_cast<int> (obj2->_maxMapW+1) ) { zOver = true; }
5225  if ( !zFound ) { hlpMap[arrPos] = 0.0; continue; }
5226 
5227  //==================== Get the surrounding values and distances
5228  p000[0] = oldX;
5229  p000[1] = oldY;
5230  p000[2] = oldZ;
5231  hlpPos = p000[2] + (obj2->_maxMapW+1) * ( p000[1] + (obj2->_maxMapV+1) * p000[0] );
5232  p000[3] = obj2->_densityMapCor[hlpPos];
5233 
5234  p001[0] = oldX;
5235  p001[1] = oldY;
5236  p001[2] = oldZ + 1;
5237  if ( zOver ) { p001[2] = 0.0; }
5238  hlpPos = p001[2] + (obj2->_maxMapW+1) * ( p001[1] + (obj2->_maxMapV+1) * p001[0] );
5239  p001[3] = obj2->_densityMapCor[hlpPos];
5240 
5241  p010[0] = oldX;
5242  p010[1] = oldY + 1;
5243  p010[2] = oldZ;
5244  if ( yOver ) { p010[1] = 0.0; }
5245  hlpPos = p010[2] + (obj2->_maxMapW+1) * ( p010[1] + (obj2->_maxMapV+1) * p010[0] );
5246  p010[3] = obj2->_densityMapCor[hlpPos];
5247 
5248  p011[0] = oldX;
5249  p011[1] = oldY + 1;
5250  p011[2] = oldZ + 1;
5251  if ( yOver ) { p011[1] = 0.0; }
5252  if ( zOver ) { p011[2] = 0.0; }
5253  hlpPos = p011[2] + (obj2->_maxMapW+1) * ( p011[1] + (obj2->_maxMapV+1) * p011[0] );
5254  p011[3] = obj2->_densityMapCor[hlpPos];
5255 
5256  p100[0] = oldX + 1;
5257  p100[1] = oldY;
5258  p100[2] = oldZ;
5259  if ( xOver ) { p100[0] = 0.0; }
5260  hlpPos = p100[2] + (obj2->_maxMapW+1) * ( p100[1] + (obj2->_maxMapV+1) * p100[0] );
5261  p100[3] = obj2->_densityMapCor[hlpPos];
5262 
5263  p101[0] = oldX + 1;
5264  p101[1] = oldY;
5265  p101[2] = oldZ + 1;
5266  if ( xOver ) { p101[0] = 0.0; }
5267  if ( zOver ) { p101[2] = 0.0; }
5268  hlpPos = p101[2] + (obj2->_maxMapW+1) * ( p101[1] + (obj2->_maxMapV+1) * p101[0] );
5269  p101[3] = obj2->_densityMapCor[hlpPos];
5270 
5271  p110[0] = oldX + 1;
5272  p110[1] = oldY + 1;
5273  p110[2] = oldZ;
5274  if ( xOver ) { p110[0] = 0.0; }
5275  if ( yOver ) { p110[1] = 0.0; }
5276  hlpPos = p110[2] + (obj2->_maxMapW+1) * ( p110[1] + (obj2->_maxMapV+1) * p110[0] );
5277  p110[3] = obj2->_densityMapCor[hlpPos];
5278 
5279  p111[0] = oldX + 1;
5280  p111[1] = oldY + 1;
5281  p111[2] = oldZ + 1;
5282  if ( xOver ) { p111[0] = 0.0; }
5283  if ( yOver ) { p111[1] = 0.0; }
5284  if ( zOver ) { p111[2] = 0.0; }
5285  hlpPos = p111[2] + (obj2->_maxMapW+1) * ( p111[1] + (obj2->_maxMapV+1) * p111[0] );
5286  p111[3] = obj2->_densityMapCor[hlpPos];
5287 
5288  //==================== Interpolate
5289  xd = 1.0 - ( ( newX - ( static_cast<double> ( oldX ) * obj2->_xSamplingRate ) ) / obj2->_xSamplingRate );
5290  p00[0] = p000[1]; p00[1] = p000[2]; p00[2] = ( xd * p000[3] ) + ( (1.0 - xd) * p100[3] );
5291  p01[0] = p001[1]; p01[1] = p001[2]; p01[2] = ( xd * p001[3] ) + ( (1.0 - xd) * p101[3] );
5292  p10[0] = p010[1]; p10[1] = p010[2]; p10[2] = ( xd * p010[3] ) + ( (1.0 - xd) * p110[3] );
5293  p11[0] = p011[1]; p11[1] = p011[2]; p11[2] = ( xd * p011[3] ) + ( (1.0 - xd) * p111[3] );
5294 
5295  yd = 1.0 - ( ( newY - ( static_cast<double> ( oldY ) * obj2->_ySamplingRate ) ) / obj2->_ySamplingRate );
5296  p0[0] = p00[1]; p0[1] = ( yd * p00[2] ) + ( (1.0 - yd) * p10[2] );
5297  p1[0] = p01[1]; p1[1] = ( yd * p01[2] ) + ( (1.0 - yd) * p11[2] );
5298 
5299  zd = 1.0 - ( ( newZ - ( static_cast<double> ( oldZ ) * obj2->_zSamplingRate ) ) / obj2->_zSamplingRate );
5300  hlpMap[arrPos] = ( zd * p0[1] ) + ( (1.0 - zd) * p1[1] );
5301  }
5302  }
5303  }
5304 
5305  //================================ Copy map
5306  delete[] obj2->_densityMapCor;
5307  obj2->_densityMapCor = new double [(newU+1) * (newV+1) * (newW+1)];
5308  for ( int iter = 0; iter < static_cast<int> ( (newU+1) * (newV+1) * (newW+1) ); iter++ )
5309  {
5310  obj2->_densityMapCor[iter] = hlpMap[iter];
5311  }
5312 
5313  //================================ Free memory
5314  delete[] hlpMap;
5315 
5316  //================================ Make sure start is the same as structure 1
5317  *ob2XMov = std::round ( ( ( obj1->_xFrom * obj1->_xSamplingRate ) - ( newXFrom * newXSampling ) ) / newXSampling );
5318  *ob2YMov = std::round ( ( ( obj1->_yFrom * obj1->_ySamplingRate ) - ( newYFrom * newYSampling ) ) / newYSampling );
5319  *ob2ZMov = std::round ( ( ( obj1->_zFrom * obj1->_zSamplingRate ) - ( newZFrom * newZSampling ) ) / newZSampling );
5320 
5321  //================================ Set new cell parameters
5322  obj2->_maxMapU = newU;
5323  obj2->_maxMapV = newV;
5324  obj2->_maxMapW = newW;
5325 
5326  obj2->_xFrom = newXFrom + (*ob2XMov);
5327  obj2->_yFrom = newYFrom + (*ob2YMov);
5328  obj2->_zFrom = newZFrom + (*ob2ZMov);
5329 
5330  obj2->_xTo = newXTo + (*ob2XMov);
5331  obj2->_yTo = newYTo + (*ob2YMov);
5332  obj2->_zTo = newZTo + (*ob2ZMov);
5333 
5334  obj2->_xRange = newXRange;
5335  obj2->_yRange = newYRange;
5336  obj2->_zRange = newZRange;
5337 
5338  obj2->_xSamplingRate = newXSampling;
5339  obj2->_ySamplingRate = newYSampling;
5340  obj2->_zSamplingRate = newZSampling;
5341 
5342  *ob2XMov *= obj2->_xSamplingRate;
5343  *ob2YMov *= obj2->_ySamplingRate;
5344  *ob2ZMov *= obj2->_zSamplingRate;
5345  }
5346 
5347  //======================================== Add zero padding to the smaller structure to make sure number of points is the same
5348  newU = std::max ( obj1->_maxMapU, obj2->_maxMapU );
5349  newV = std::max ( obj1->_maxMapV, obj2->_maxMapV );
5350  newW = std::max ( obj1->_maxMapW, obj2->_maxMapW );
5351 
5352  if ( ( static_cast<int> ( obj1->_maxMapU ) != newU ) || ( static_cast<int> ( obj1->_maxMapV ) != newV ) || ( static_cast<int> ( obj1->_maxMapW ) != newW ) )
5353  {
5354  //==================================== Initialise
5355  double *hlpMap = new double [(newU+1) * (newV+1) * (newW+1)];
5356 
5357  //==================================== Zero padd
5358  for ( int uIt = 0; uIt < ( newU + 1 ); uIt++ )
5359  {
5360  for ( int vIt = 0; vIt < ( newV + 1 ); vIt++ )
5361  {
5362  for ( int wIt = 0; wIt < ( newW + 1 ); wIt++ )
5363  {
5364  arrPos = wIt + (newW+1) * ( vIt + (newV+1) * uIt );
5365  hlpPos = wIt + (obj1->_maxMapW+1) * ( vIt + (obj1->_maxMapV+1) * uIt );
5366 
5367  if ( ( uIt < ( static_cast<int> ( obj1->_maxMapU ) + 1 ) ) &&
5368  ( vIt < ( static_cast<int> ( obj1->_maxMapV ) + 1 ) ) &&
5369  ( wIt < ( static_cast<int> ( obj1->_maxMapW ) + 1 ) ) )
5370  {
5371  hlpMap[arrPos] = obj1->_densityMapCor[hlpPos];
5372  }
5373  else
5374  {
5375  hlpMap[arrPos] = 0.0;
5376  }
5377  }
5378  }
5379  }
5380 
5381  //==================================== Copy map
5382  delete[] obj1->_densityMapCor;
5383  obj1->_densityMapCor = new double [(newU+1) * (newV+1) * (newW+1)];
5384  for ( int iter = 0; iter < static_cast<int> ( (newU+1) * (newV+1) * (newW+1) ); iter++ )
5385  {
5386  obj1->_densityMapCor[iter] = hlpMap[iter];
5387  }
5388 
5389  //=================================== Free memory
5390  delete[] hlpMap;
5391 
5392  //==================================== Set new cell parameters
5393  obj1->_xTo += newU - obj1->_maxMapU;
5394  obj1->_yTo += newV - obj1->_maxMapV;
5395  obj1->_zTo += newW - obj1->_maxMapW;
5396 
5397  obj1->_maxMapU = newU;
5398  obj1->_maxMapV = newV;
5399  obj1->_maxMapW = newW;
5400 
5401  obj1->_xRange = ( newU + 1 ) * obj1->_xSamplingRate;
5402  obj1->_yRange = ( newV + 1 ) * obj1->_ySamplingRate;
5403  obj1->_zRange = ( newW + 1 ) * obj1->_zSamplingRate;
5404  }
5405 
5406  if ( ( static_cast<int> ( obj2->_maxMapU ) != newU ) || ( static_cast<int> ( obj2->_maxMapV ) != newV ) || ( static_cast<int> ( obj2->_maxMapW ) != newW ) )
5407  {
5408  //==================================== Initialise
5409  double *hlpMap = new double [(newU+1) * (newV+1) * (newW+1)];
5410 
5411  //==================================== Zero padd
5412  for ( int uIt = 0; uIt < ( newU + 1 ); uIt++ )
5413  {
5414  for ( int vIt = 0; vIt < ( newV + 1 ); vIt++ )
5415  {
5416  for ( int wIt = 0; wIt < ( newW + 1 ); wIt++ )
5417  {
5418  arrPos = wIt + (newW+1) * ( vIt + (newV+1) * uIt );
5419  hlpPos = wIt + (obj2->_maxMapW+1) * ( vIt + (obj2->_maxMapV+1) * uIt );
5420 
5421  if ( ( uIt < ( static_cast<int> ( obj2->_maxMapU ) + 1 ) ) &&
5422  ( vIt < ( static_cast<int> ( obj2->_maxMapV ) + 1 ) ) &&
5423  ( wIt < ( static_cast<int> ( obj2->_maxMapW ) + 1 ) ) )
5424  {
5425  hlpMap[arrPos] = obj2->_densityMapCor[hlpPos];
5426  }
5427  else
5428  {
5429  hlpMap[arrPos] = 0.0;
5430  }
5431  }
5432  }
5433  }
5434 
5435  //==================================== Copy map
5436  delete[] obj2->_densityMapCor;
5437  obj2->_densityMapCor = new double [(newU+1) * (newV+1) * (newW+1)];
5438  for ( int iter = 0; iter < static_cast<int> ( (newU+1) * (newV+1) * (newW+1) ); iter++ )
5439  {
5440  obj2->_densityMapCor[iter] = hlpMap[iter];
5441  }
5442 
5443  //=================================== Free memory
5444  delete[] hlpMap;
5445 
5446  //==================================== Set new cell parameters
5447  obj2->_xTo += newU - obj2->_maxMapU;
5448  obj2->_yTo += newV - obj2->_maxMapV;
5449  obj2->_zTo += newW - obj2->_maxMapW;
5450 
5451  obj2->_maxMapU = newU;
5452  obj2->_maxMapV = newV;
5453  obj2->_maxMapW = newW;
5454 
5455  obj2->_xRange = ( newU + 1 ) * obj2->_xSamplingRate;
5456  obj2->_yRange = ( newV + 1 ) * obj2->_ySamplingRate;
5457  obj2->_zRange = ( newW + 1 ) * obj2->_zSamplingRate;
5458  }
5459 
5460  //======================================== Initialise FFTWs for Translation function
5461  int minU = std::min ( (obj1->_maxMapU+1), (obj2->_maxMapU+1) );
5462  int minV = std::min ( (obj1->_maxMapV+1), (obj2->_maxMapV+1) );
5463  int minW = std::min ( (obj1->_maxMapW+1), (obj2->_maxMapW+1) );
5464 
5465  fftw_complex *tmpIn1 = new fftw_complex[(obj1->_maxMapU+1) * (obj1->_maxMapV+1) * (obj1->_maxMapW+1)];
5466  fftw_complex *tmpOut1 = new fftw_complex[(obj1->_maxMapU+1) * (obj1->_maxMapV+1) * (obj1->_maxMapW+1)];
5467  fftw_complex *tmpIn2 = new fftw_complex[(obj2->_maxMapU+1) * (obj2->_maxMapV+1) * (obj2->_maxMapW+1)];
5468  fftw_complex *tmpOut2 = new fftw_complex[(obj2->_maxMapU+1) * (obj2->_maxMapV+1) * (obj2->_maxMapW+1)];
5469  fftw_complex *resIn = new fftw_complex[minU * minV * minW];
5470  fftw_complex *resOut = new fftw_complex[minU * minV * minW];
5471 
5472  //======================================== Get Fourier transforms of the maps
5473  fftw_plan forwardFourierObj1;
5474  fftw_plan forwardFourierObj2;
5475  fftw_plan inverseFourierCombo;
5476 
5477  forwardFourierObj1 = fftw_plan_dft_3d ( (obj1->_maxMapU+1), (obj1->_maxMapV+1), (obj1->_maxMapW+1), tmpIn1, tmpOut1, FFTW_FORWARD, FFTW_ESTIMATE );
5478  forwardFourierObj2 = fftw_plan_dft_3d ( (obj2->_maxMapU+1), (obj2->_maxMapV+1), (obj2->_maxMapW+1), tmpIn2, tmpOut2, FFTW_FORWARD, FFTW_ESTIMATE );
5479  inverseFourierCombo = fftw_plan_dft_3d ( minU, minV, minW, resOut, resIn, FFTW_BACKWARD, FFTW_ESTIMATE );
5480 
5481  //======================================== Fill in input data
5482  for ( int iter = 0; iter < static_cast<int> ( (obj1->_maxMapU+1) * (obj1->_maxMapV+1) * (obj1->_maxMapW+1) ); iter++ )
5483  {
5484  tmpIn1[iter][0] = obj1->_densityMapCor[iter];
5485  tmpIn1[iter][1] = 0.0;
5486  }
5487 
5488  for ( int iter = 0; iter < static_cast<int> ( (obj2->_maxMapU+1) * (obj2->_maxMapV+1) * (obj2->_maxMapW+1) ); iter++ )
5489  {
5490  tmpIn2[iter][0] = obj2->_densityMapCor[iter];
5491  tmpIn2[iter][1] = 0.0;
5492  }
5493 
5494  //======================================== Calculate Fourier
5495  fftw_execute ( forwardFourierObj1 );
5496  fftw_execute ( forwardFourierObj2 );
5497 
5498  //======================================== Combine Fourier coeffs
5499  double normFactor = static_cast<double> ( minU * minV * minW );
5500  std::array<double,2> hlpArr = std::array<double,2> { { 2.0, 2.0 } };
5501  int h1, k1, l1;
5502  int uo2 = (obj1->_maxMapU+1) / 2;
5503  int vo2 = (obj1->_maxMapV+1) / 2;
5504  int wo2 = (obj1->_maxMapW+1) / 2;
5505  int muo2 = (minU+1) / 2;
5506  int mvo2 = (minV+1) / 2;
5507  int mwo2 = (minW+1) / 2;
5508 
5509  for ( int uIt = 0; uIt < static_cast<int> ( obj1->_maxMapU+1 ); uIt++ )
5510  {
5511  for ( int vIt = 0; vIt < static_cast<int> ( obj1->_maxMapV+1 ); vIt++ )
5512  {
5513  for ( int wIt = 0; wIt < static_cast<int> ( obj1->_maxMapW+1 ); wIt++ )
5514  {
5515  if ( uIt > uo2 ) { h1 = uIt - (obj1->_maxMapU+1); } else { h1 = uIt; }
5516  if ( vIt > vo2 ) { k1 = vIt - (obj1->_maxMapV+1); } else { k1 = vIt; }
5517  if ( wIt > wo2 ) { l1 = wIt - (obj1->_maxMapW+1); } else { l1 = wIt; }
5518 
5519  if ( ( std::abs ( h1 ) <= muo2 ) && ( std::abs ( k1 ) <= mvo2 ) && ( std::abs ( l1 ) <= mwo2 ) )
5520  {
5521  if ( h1 < 0 ) { h1 += minU; }
5522  if ( k1 < 0 ) { k1 += minV; }
5523  if ( l1 < 0 ) { l1 += minW; }
5524 
5525  hlpPos = l1 + minW * ( k1 + minV * h1 );
5526  arrPos = wIt + (obj1->_maxMapW+1) * ( vIt + (obj1->_maxMapV+1) * uIt );
5527  resIn[hlpPos][0] = tmpOut1[arrPos][0];
5528  resIn[hlpPos][1] = tmpOut1[arrPos][1];
5529  }
5530  }
5531  }
5532  }
5533 
5534  uo2 = (obj2->_maxMapU+1) / 2;
5535  vo2 = (obj2->_maxMapV+1) / 2;
5536  wo2 = (obj2->_maxMapW+1) / 2;
5537  for ( int uIt = 0; uIt < static_cast<int> ( obj2->_maxMapU+1 ); uIt++ )
5538  {
5539  for ( int vIt = 0; vIt < static_cast<int> ( obj2->_maxMapV+1 ); vIt++ )
5540  {
5541  for ( int wIt = 0; wIt < static_cast<int> ( obj2->_maxMapW+1 ); wIt++ )
5542  {
5543  if ( uIt > uo2 ) { h1 = uIt - (obj2->_maxMapU+1); } else { h1 = uIt; }
5544  if ( vIt > vo2 ) { k1 = vIt - (obj2->_maxMapV+1); } else { k1 = vIt; }
5545  if ( wIt > wo2 ) { l1 = wIt - (obj2->_maxMapW+1); } else { l1 = wIt; }
5546 
5547  if ( ( std::abs ( h1 ) <= muo2 ) && ( std::abs ( k1 ) <= mvo2 ) && ( std::abs ( l1 ) <= mwo2 ) )
5548  {
5549  if ( h1 < 0 ) { h1 += minU; }
5550  if ( k1 < 0 ) { k1 += minV; }
5551  if ( l1 < 0 ) { l1 += minW; }
5552 
5553  hlpPos = l1 + minW * ( k1 + minV * h1 );
5554  arrPos = wIt + (obj2->_maxMapW+1) * ( vIt + (obj2->_maxMapV+1) * uIt );
5555  hlpArr = ProSHADE_internal_misc::complexMultiplicationConjug ( &resIn[hlpPos][0],
5556  &resIn[hlpPos][1],
5557  &tmpOut2[arrPos][0],
5558  &tmpOut2[arrPos][1] );
5559  resOut[hlpPos][0] = hlpArr[0] / normFactor;
5560  resOut[hlpPos][1] = hlpArr[1] / normFactor;
5561 
5562  }
5563  }
5564  }
5565  }
5566 
5567  fftw_execute ( inverseFourierCombo );
5568 
5569  //======================================== Find highest peak in translation function
5570  double mapPeak = 0.0;
5571  for ( int uIt = 0; uIt < minU; uIt++ )
5572  {
5573  for ( int vIt = 0; vIt < minV; vIt++ )
5574  {
5575  for ( int wIt = 0; wIt < minW; wIt++ )
5576  {
5577  arrPos = wIt + minW * ( vIt + minV * uIt );
5578  if ( resIn[arrPos][0] > mapPeak )
5579  {
5580  mapPeak = resIn[arrPos][0];
5581  ret[0] = uIt;
5582  ret[1] = vIt;
5583  ret[2] = wIt;
5584  }
5585  }
5586  }
5587  }
5588  ret[3] = mapPeak / ( ( obj2->_maxMapU+1 ) * ( obj2->_maxMapV+1 ) * ( obj2->_maxMapW+1 ) );
5589 
5590  //======================================== Dont translate over half
5591  if ( ret[0] > muo2 ) { ret[0] = ret[0] - minU; }
5592  if ( ret[1] > mvo2 ) { ret[1] = ret[1] - minV; }
5593  if ( ret[2] > mwo2 ) { ret[2] = ret[2] - minW; }
5594 
5595  //======================================== Free memory
5596  fftw_destroy_plan ( forwardFourierObj1 );
5597  fftw_destroy_plan ( forwardFourierObj2 );
5598  fftw_destroy_plan ( inverseFourierCombo );
5599  delete[] tmpIn1;
5600  delete[] tmpIn2;
5601  delete[] tmpOut1;
5602  delete[] tmpOut2;
5603  delete[] resIn;
5604  delete[] resOut;
5605 
5606  //======================================== Convert to angstroms
5607  ret[0] *= obj2->_xSamplingRate;
5608  ret[1] *= obj2->_ySamplingRate;
5609  ret[2] *= obj2->_zSamplingRate;
5610 
5611  //======================================== Do not translate over half the cell size, do minus instead
5612  if ( ret[0] > ( obj2->_xRange / 2.0 ) ) { ret[0] -= obj2->_xRange; }
5613  if ( ret[1] > ( obj2->_yRange / 2.0 ) ) { ret[1] -= obj2->_yRange; }
5614  if ( ret[2] > ( obj2->_zRange / 2.0 ) ) { ret[2] -= obj2->_zRange; }
5615 
5616  if ( ret[0] < ( -obj2->_xRange / 2.0 ) ) { ret[0] += obj2->_xRange; }
5617  if ( ret[1] < ( -obj2->_yRange / 2.0 ) ) { ret[1] += obj2->_yRange; }
5618  if ( ret[2] < ( -obj2->_zRange / 2.0 ) ) { ret[2] += obj2->_zRange; }
5619 
5620  //======================================== Done
5621  return ( ret );
5622 }
5623 
5642  ProSHADE::ProSHADE_settings* settings,
5643  std::string saveName,
5644  int verbose,
5645  std::string axOrd,
5646  bool internalUse )
5647 {
5648  //======================================== Sanity checks
5649  if ( ( this->_eulerAngles[0] == 0.0 ) && ( this->_eulerAngles[1] == 0.0 ) && ( this->_eulerAngles[2] == 0.0 ) && ( !internalUse ) )
5650  {
5651  cmpObj1->writeMap ( saveName, cmpObj1->_densityMapCor, axOrd );
5652  return ;
5653  }
5654 
5655  if ( ( this->_eulerAngles[0] == 0.0 ) && ( this->_eulerAngles[1] == 0.0 ) && ( this->_eulerAngles[2] == 0.0 ) && ( internalUse ) )
5656  {
5657  return ;
5658  }
5659 
5660  //======================================== Initalise
5661  std::array<double,2> hlpArr;
5662  int arrPos;
5663  int arrPos2;
5664  double **realSHCoeffsRes = new double* [cmpObj1->_noShellsWithData];
5665  double **imagSHCoeffsRes = new double* [cmpObj1->_noShellsWithData];
5666  for ( unsigned int i = 0; i < cmpObj1->_noShellsWithData; i++ )
5667  {
5668  realSHCoeffsRes[i] = new double [cmpObj1->_oneDimmension * cmpObj1->_oneDimmension];
5669  imagSHCoeffsRes[i] = new double [cmpObj1->_oneDimmension * cmpObj1->_oneDimmension];
5670 
5671  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( cmpObj1->_oneDimmension * cmpObj1->_oneDimmension ); iter++ )
5672  {
5673  realSHCoeffsRes[i][iter] = 0.0;
5674  imagSHCoeffsRes[i][iter] = 0.0;
5675  }
5676  }
5677 
5678  //======================================== Get Wigner matrices
5679  this->generateWignerMatrices ( settings );
5680 
5681  if ( verbose > 3 )
5682  {
5683  std::cout << ">>>>>>>> Wigner matrices for given Euler angles obtained." << std::endl;
5684  }
5685 
5686  //======================================== Multiply coeffs by Wigner
5687  // ... for each shell
5688  for ( int shell = 0; shell < static_cast<int> ( cmpObj1->_noShellsWithData ); shell++ )
5689  {
5690  //==================================== ... for each band
5691  for ( int bandIter = 0; bandIter < static_cast<int> ( this->_bandwidthLimit ); bandIter++ )
5692  {
5693  //================================ ... for each order
5694  for ( int ord1 = 0; ord1 < (bandIter * 2 + 1); ord1++ )
5695  {
5696  arrPos2 = seanindex ( ord1-bandIter, bandIter, this->_bandwidthLimit );
5697 
5698  //============================ ... again, for each order prime
5699  for ( int ord2 = 0; ord2 < (bandIter * 2 + 1); ord2++ )
5700  {
5701  arrPos = seanindex ( ord2-bandIter, bandIter, this->_bandwidthLimit );
5702  hlpArr = ProSHADE_internal_misc::complexMultiplication ( &cmpObj1->_realSHCoeffs[shell][arrPos],
5703  &cmpObj1->_imagSHCoeffs[shell][arrPos],
5704  &this->_wignerMatrices.at(bandIter).at(ord1).at(ord2)[0],
5705  &this->_wignerMatrices.at(bandIter).at(ord1).at(ord2)[1] );
5706 
5707  realSHCoeffsRes[shell][arrPos2] += hlpArr[0];
5708  imagSHCoeffsRes[shell][arrPos2] += hlpArr[1];
5709  }
5710  }
5711  }
5712  }
5713  if ( verbose > 0 )
5714  {
5715  std::cout << "Rotated all coefficients." << std::endl;
5716  }
5717 
5718  //======================================== Inverse the SH coeffs to shells
5719  // ... Initalise memory
5720  fftw_plan idctPlan, ifftPlan;
5721  int rank, howmany_rank ;
5722  fftw_iodim dims[1], howmany_dims[1];
5723 
5724  double *sigR = new double [cmpObj1->_oneDimmension * cmpObj1->_oneDimmension];
5725  double *sigI = new double [cmpObj1->_oneDimmension * cmpObj1->_oneDimmension];
5726  double *rcoeffs = new double [cmpObj1->_oneDimmension * cmpObj1->_oneDimmension];
5727  double *icoeffs = new double [cmpObj1->_oneDimmension * cmpObj1->_oneDimmension];
5728  double *weights = new double [4 * this->_bandwidthLimit];
5729  double *workspace = new double [( 10 * this->_bandwidthLimit * this->_bandwidthLimit ) + ( 24 * this->_bandwidthLimit )];
5730 
5731  //======================================== Initialise internal values ( and allocate memory )
5732  double **newShellMappedData = new double*[cmpObj1->_noShellsWithData];
5733  for ( unsigned int sh = 0; sh < cmpObj1->_noShellsWithData; sh++ )
5734  {
5735  newShellMappedData[sh] = new double[static_cast<unsigned int> ( 4 * this->_bandwidthLimit * this->_bandwidthLimit )];
5736  }
5737  if ( verbose > 2 )
5738  {
5739  std::cout << ">>>>> Memory allocation complete." << std::endl;
5740  }
5741 
5742  //======================================== Initalise FFTW plans
5743  idctPlan = fftw_plan_r2r_1d ( 2 * this->_bandwidthLimit,
5744  weights,
5745  workspace,
5746  FFTW_REDFT01,
5747  FFTW_ESTIMATE );
5748 
5749  rank = 1;
5750  dims[0].n = 2 * this->_bandwidthLimit;
5751  dims[0].is = 2 * this->_bandwidthLimit;
5752  dims[0].os = 1;
5753  howmany_rank = 1;
5754  howmany_dims[0].n = 2 * this->_bandwidthLimit;
5755  howmany_dims[0].is = 1;
5756  howmany_dims[0].os = 2 * this->_bandwidthLimit;
5757 
5758  //======================================== Inverse FFT
5759  ifftPlan = fftw_plan_guru_split_dft ( rank,
5760  dims,
5761  howmany_rank,
5762  howmany_dims,
5763  sigR,
5764  sigI,
5765  rcoeffs,
5766  icoeffs,
5767  FFTW_ESTIMATE );
5768 
5769  //======================================== ... for each shell
5770  for ( int shell = 0; shell < static_cast<int> ( cmpObj1->_noShellsWithData ); shell++ )
5771  {
5772  //==================================== Load SH coeffs to arrays
5773  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( cmpObj1->_oneDimmension * cmpObj1->_oneDimmension ); iter++ )
5774  {
5775  rcoeffs[iter] = realSHCoeffsRes[shell][iter];
5776  icoeffs[iter] = imagSHCoeffsRes[shell][iter];
5777  }
5778 
5779  //==================================== Get inverse spherical harmonics transform for the shell
5780  InvFST_semi_fly ( rcoeffs,
5781  icoeffs,
5782  sigR,
5783  sigI,
5784  this->_bandwidthLimit,
5785  workspace,
5786  0,
5787  this->_bandwidthLimit,
5788  &idctPlan,
5789  &ifftPlan );
5790 
5791  //==================================== Copy the results to the rotated shells array
5792  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( cmpObj1->_oneDimmension * cmpObj1->_oneDimmension ); iter++ )
5793  {
5794  newShellMappedData[shell][iter] = sigR[iter];
5795  }
5796  }
5797  if ( verbose > 1 )
5798  {
5799  std::cout << ">> Inverted all rotated spherical harmonics coefficients." << std::endl;
5800  }
5801 
5802  //======================================== Clear memory
5803  delete[] sigR;
5804  delete[] sigI;
5805  delete[] rcoeffs;
5806  delete[] icoeffs;
5807  delete[] weights;
5808  delete[] workspace;
5809 
5810  //======================================== Clear memory
5811  for ( unsigned int i = 0; i < cmpObj1->_noShellsWithData; i++ )
5812  {
5813  delete[] realSHCoeffsRes[i];
5814  delete[] imagSHCoeffsRes[i];
5815  }
5816  delete[] realSHCoeffsRes;
5817  delete[] imagSHCoeffsRes;
5818 
5819  //======================================== Convert the shell data back to cartesian grid
5820  // ... Create the Cartesian grid
5821  double *densityMapRotated = new double [(cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1)];
5822 
5823  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ((cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1)); iter++ )
5824  {
5825  densityMapRotated[iter] = 0.0;
5826  }
5827 
5828  //======================================== Find spherical cut-offs
5829  std::vector<double> lonCO ( cmpObj1->_thetaAngle + 1 );
5830  for ( unsigned int iter = 0; iter <= cmpObj1->_thetaAngle; iter++ ) { lonCO.at(iter) = static_cast<double> ( iter ) * ( ( static_cast<double> ( M_PI ) * 2.0 ) / static_cast<double> ( this->_thetaAngle ) ) - ( static_cast<double> ( M_PI ) ); }
5831  std::vector<double> latCO ( cmpObj1->_phiAngle + 1 );
5832  for ( unsigned int iter = 0; iter <= cmpObj1->_phiAngle; iter++ ) { latCO.at(iter) = ( static_cast<double> ( iter ) * ( static_cast<double> ( M_PI ) / static_cast<double> ( this->_phiAngle ) ) - ( static_cast<double> ( M_PI ) / 2.0 ) ); }
5833 
5834  //======================================== Interpolate onto cartesian grid
5835  double rad = 0.0;
5836  double lon = 0.0;
5837  double lat = 0.0;
5838  double newU = 0.0;
5839  double newV = 0.0;
5840  double newW = 0.0;
5841  unsigned int lowerLon = 0;
5842  unsigned int upperLon = 0;
5843  unsigned int lowerLat = 0;
5844  unsigned int upperLat = 0;
5845  unsigned int lowerShell = 0;
5846  unsigned int upperShell = 0;
5847  double x00 = 0.0;
5848  double x01 = 0.0;
5849  double x10 = 0.0;
5850  double x11 = 0.0;
5851  double dist00 = 0.0;
5852  double dist01 = 0.0;
5853  double dist10 = 0.0;
5854  double dist11 = 0.0;
5855  double wCoeff = 0.0;
5856  double lsVal = 0.0;
5857  double upVal = 0.0;
5858 
5859  for ( int uIt = 0; uIt < static_cast<int> (cmpObj1->_maxMapU+1); uIt++ )
5860  {
5861  for ( int vIt = 0; vIt < static_cast<int> (cmpObj1->_maxMapV+1); vIt++ )
5862  {
5863  for ( int wIt = 0; wIt < static_cast<int> (cmpObj1->_maxMapW+1); wIt++ )
5864  {
5865 
5866  //============================ Convert to centered coords
5867  newU = static_cast<double> ( uIt - ( static_cast<int> (cmpObj1->_maxMapU+1) / 2 ) );
5868  newV = static_cast<double> ( vIt - ( static_cast<int> (cmpObj1->_maxMapV+1) / 2 ) );
5869  newW = static_cast<double> ( wIt - ( static_cast<int> (cmpObj1->_maxMapW+1) / 2 ) );
5870 
5871  //============================ Deal with 0 ; 0 ; 0
5872  if ( ( newU == 0.0 ) && ( newV == 0.0 ) && ( newW == 0.0 ) )
5873  {
5874  arrPos = wIt + (cmpObj1->_maxMapW + 1) * ( vIt + (cmpObj1->_maxMapV + 1) * uIt );
5875  densityMapRotated[arrPos] = cmpObj1->_densityMapCor[arrPos];
5876  continue;
5877  }
5878 
5879  //============================ Convert to spherical coords
5880  rad = sqrt ( pow( (newU*cmpObj1->_xSamplingRate), 2.0 ) +
5881  pow( (newV*cmpObj1->_ySamplingRate), 2.0 ) +
5882  pow( (newW*cmpObj1->_zSamplingRate), 2.0 ) );
5883  lon = atan2 ( (newV*cmpObj1->_ySamplingRate) , (newU*cmpObj1->_xSamplingRate) );
5884  lat = asin ( (newW*cmpObj1->_zSamplingRate) / rad );
5885 
5886  //============================ Deal with nan's
5887  if ( rad != rad ) { rad = 0.0; }
5888  if ( lon != lon ) { lon = 0.0; }
5889  if ( lat != lat ) { lat = 0.0; }
5890 
5891  //============================ Find the angle cutoffs around the point
5892  for ( unsigned int iter = 0; iter < cmpObj1->_thetaAngle; iter++ )
5893  {
5894  if ( ( std::trunc(10000. * lonCO.at(iter)) <= std::trunc(10000. * lon) ) && ( std::trunc(10000. * lonCO.at(iter+1)) > std::trunc(10000. * lon) ) )
5895  {
5896  lowerLon = iter;
5897  upperLon = iter+1;
5898  break;
5899  }
5900  }
5901  if ( upperLon == cmpObj1->_thetaAngle ) { upperLon = 0; }
5902 
5903  for ( unsigned int iter = 0; iter < cmpObj1->_phiAngle; iter++ )
5904  {
5905  if ( ( std::trunc(10000. * latCO.at(iter)) <= std::trunc(10000. * lat) ) && ( std::trunc(10000. * latCO.at(iter+1)) > std::trunc(10000. * lat) ) )
5906  {
5907  lowerLat = iter;
5908  upperLat = iter+1;
5909  break;
5910  }
5911  }
5912  if ( upperLat == cmpObj1->_phiAngle ) { upperLat = 0; }
5913 
5914  //============================ Find shells above and below
5915  lowerShell = 0;
5916  upperShell = 0;
5917  for ( unsigned int iter = 0; iter < (cmpObj1->_noShellsWithData-1); iter++ )
5918  {
5919  if ( ( cmpObj1->_shellPlacement.at(iter) <= rad ) && ( cmpObj1->_shellPlacement.at(iter+1) > rad ) )
5920  {
5921  lowerShell = iter;
5922  upperShell = iter+1;
5923  break;
5924  }
5925  }
5926  if ( upperShell == 0 )
5927  {
5928  arrPos = wIt + (cmpObj1->_maxMapW + 1) * ( vIt + (cmpObj1->_maxMapV + 1) * uIt );
5929  densityMapRotated[arrPos] = 0.0;
5930  continue;
5931  }
5932 
5933  //============================ Interpolate lower shell
5934  x00 = newShellMappedData[lowerShell][static_cast<int> ( lowerLat * cmpObj1->_thetaAngle + lowerLon )];
5935  dist00 = sqrt( pow( (lon - lonCO.at(lowerLon)), 2.0 ) +
5936  pow( (lat - latCO.at(lowerLat)), 2.0 ) );
5937  x01 = newShellMappedData[lowerShell][static_cast<int> ( lowerLat * cmpObj1->_thetaAngle + upperLon )];
5938  dist01 = sqrt( pow( (lon - lonCO.at(upperLon)), 2.0 ) +
5939  pow( (lat - latCO.at(lowerLat)), 2.0 ) );
5940  x10 = newShellMappedData[lowerShell][static_cast<int> ( upperLat * cmpObj1->_thetaAngle + lowerLon )];
5941  dist10 = sqrt( pow( (lon - lonCO.at(lowerLon)), 2.0 ) +
5942  pow( (lat - latCO.at(upperLat)), 2.0 ) );
5943  x11 = newShellMappedData[lowerShell][static_cast<int> ( upperLat * cmpObj1->_thetaAngle + upperLon )];
5944  dist11 = sqrt( pow( (lon - lonCO.at(upperLon)), 2.0 ) +
5945  pow( (lat - latCO.at(upperLat)), 2.0 ) );
5946 
5947  if ( !( dist00 == 0.0 ) || !( dist01 == 0.0 ) || !( dist10 == 0.0 ) || !( dist11 == 0.0 ) )
5948  {
5949  if ( ( dist00 == 0.0 ) && ( dist01 == 0.0 ) && ( dist10 == 0.0 ) && ( dist11 == 0.0 ) )
5950  {
5951  arrPos = wIt + (cmpObj1->_maxMapW + 1) * ( vIt + (cmpObj1->_maxMapV + 1) * uIt );
5952  densityMapRotated[arrPos] = 0.0;
5953  continue;
5954  }
5955  else
5956  {
5957  if ( dist00 == 0.0 ) { dist00 = std::max ( dist01, std::max ( dist10, dist11 ) ) * 0.0000001; }
5958  if ( dist01 == 0.0 ) { dist01 = std::max ( dist00, std::max ( dist10, dist11 ) ) * 0.0000001; }
5959  if ( dist10 == 0.0 ) { dist10 = std::max ( dist01, std::max ( dist00, dist11 ) ) * 0.0000001; }
5960  if ( dist11 == 0.0 ) { dist11 = std::max ( dist01, std::max ( dist10, dist00 ) ) * 0.0000001; }
5961  }
5962  }
5963 
5964  wCoeff = 1.0/( 1.0/dist00 + 1.0/dist01 + 1.0/dist10 + 1.0/dist11 );
5965  lsVal = (x00 * (( 1.0/dist00 )*wCoeff)) + (x01 * ((1.0/dist01)*wCoeff)) + (x10 * ((1.0/dist10)*wCoeff)) + (x11 * ((1.0/dist11)*wCoeff));
5966 
5967  //============================ Interpolate upper shell
5968  x00 = newShellMappedData[upperShell][static_cast<int> ( lowerLat * cmpObj1->_thetaAngle + lowerLon )];
5969  x01 = newShellMappedData[upperShell][static_cast<int> ( lowerLat * cmpObj1->_thetaAngle + upperLon )];
5970  x10 = newShellMappedData[upperShell][static_cast<int> ( upperLat * cmpObj1->_thetaAngle + lowerLon )];
5971  x11 = newShellMappedData[upperShell][static_cast<int> ( upperLat * cmpObj1->_thetaAngle + upperLon )];
5972 
5973  upVal = (x00 * (( 1.0/dist00 )*wCoeff)) + (x01 * ((1.0/dist01)*wCoeff)) + (x10 * ((1.0/dist10)*wCoeff)) + (x11 * ((1.0/dist11)*wCoeff));
5974 
5975  //============================ Interpolate radius
5976  dist00 = cmpObj1->_shellPlacement.at(upperShell) - rad;
5977  dist11 = rad - cmpObj1->_shellPlacement.at(lowerShell);
5978 
5979  arrPos = wIt + (cmpObj1->_maxMapW + 1) * ( vIt + (cmpObj1->_maxMapV + 1) * uIt );
5980  densityMapRotated[arrPos] = upVal*dist11 + lsVal*dist00;
5981 
5982  continue;
5983  }
5984  }
5985  }
5986 
5987  //======================================== Normalise
5988  std::vector<double> rMapVals ( (cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1) );
5989  std::vector<double> oMapVals ( (cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1) );
5990 
5991 // cmpObj1->writeMap ( "test.map", cmpObj1->_densityMapCor ); exit (0);
5992 //
5993 // if ( cmpObj1->_densityMapCor != nullptr ) { std::cout << "! NICE !" << std::endl; } else { std::cout << "! FUCK !" << std::endl; }
5994 // std::cout << cmpObj1->_maxMapU << " x " << cmpObj1->_maxMapV << " x " << cmpObj1->_maxMapW << " = " << (cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1) << std::endl;
5995  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( (cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1) ); iter++ )
5996  {
5997 // std::cout << iter << std::endl;
5998 // std::cout << " ... " << densityMapRotated[iter] << std::endl;
5999 // std::cout << " xxx " << rMapVals.at(iter) << std::endl;
6000  if ( densityMapRotated[iter] == densityMapRotated[iter] ) { rMapVals.at(iter) = densityMapRotated[iter]; }
6001  else { rMapVals.at(iter) = 0.0; }
6002 // std::cout << " xxx " << oMapVals.at(iter) << std::endl;
6003 // std::cout << " ... " << cmpObj1->_densityMapCor[iter] << std::endl;
6004  oMapVals.at(iter) = cmpObj1->_densityMapCor[iter];
6005  }
6006 
6007  double rMapMean = std::accumulate ( rMapVals.begin(), rMapVals.end(), 0.0 ) / static_cast<double> ( rMapVals.size() );
6008  double rMapSdev = std::sqrt ( static_cast<double> ( std::inner_product ( rMapVals.begin(), rMapVals.end(), rMapVals.begin(), 0.0 ) ) / static_cast<double> ( rMapVals.size() ) - rMapMean * rMapMean );
6009 
6010  double oMapMean = std::accumulate ( oMapVals.begin(), oMapVals.end(), 0.0 ) / static_cast<double> ( oMapVals.size() );
6011  double oMapSdev = std::sqrt ( static_cast<double> ( std::inner_product ( oMapVals.begin(), oMapVals.end(), oMapVals.begin(), 0.0 ) ) / static_cast<double> ( oMapVals.size() ) - oMapMean * oMapMean );
6012  oMapVals.clear ( );
6013 
6014  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( (cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1) ); iter++ )
6015  {
6016  densityMapRotated[iter] = ( rMapVals.at(iter) - rMapMean ) / rMapSdev;
6017  cmpObj1->_densityMapCor[iter] = ( cmpObj1->_densityMapCor[iter] - oMapMean ) / oMapSdev;
6018  }
6019  rMapVals.clear ( );
6020 
6021  //======================================== Remove all outside the max shell
6022  double dist = std::max ( (cmpObj1->_maxMapU+1)/2.0, std::max ( (cmpObj1->_maxMapV+1)/2.0, (cmpObj1->_maxMapW+1)/2.0 ) );
6023 
6024  for ( int uIt = 0; uIt < static_cast<int> (cmpObj1->_maxMapU+1); uIt++ )
6025  {
6026  for ( int vIt = 0; vIt < static_cast<int> (cmpObj1->_maxMapV+1); vIt++ )
6027  {
6028  for ( int wIt = 0; wIt < static_cast<int> (cmpObj1->_maxMapW+1); wIt++ )
6029  {
6030  arrPos = wIt + (cmpObj1->_maxMapW + 1) * ( vIt + (cmpObj1->_maxMapV + 1) * uIt );
6031 
6032  if ( sqrt ( pow( uIt-((cmpObj1->_maxMapU+1)/2.0), 2.0 ) +
6033  pow( vIt-((cmpObj1->_maxMapV+1)/2.0), 2.0 ) +
6034  pow( wIt-((cmpObj1->_maxMapW+1)/2.0), 2.0 ) ) > dist )
6035  {
6036  densityMapRotated[arrPos] = 0.0;
6037  cmpObj1->_densityMapCor[arrPos] = 0.0;
6038  }
6039  }
6040  }
6041  }
6042 
6043  //======================================== Clear memory
6044  for ( unsigned int sh = 0; sh < cmpObj1->_noShellsWithData; sh++ )
6045  {
6046  delete[] newShellMappedData[sh];
6047  }
6048  delete[] newShellMappedData;
6049 
6050  //======================================== Write out the rotated map and clean up
6051  if ( !internalUse )
6052  {
6053  cmpObj1->writeMap ( saveName, densityMapRotated, axOrd );
6054 
6055  delete[] cmpObj1->_densityMapCor;
6056  cmpObj1->_densityMapCor = nullptr;
6057  }
6058  else
6059  {
6060  //==================================== ... or save the rotated map internally
6061  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( (cmpObj1->_maxMapU+1) * (cmpObj1->_maxMapV+1) * (cmpObj1->_maxMapW+1) ); iter++ )
6062  {
6063  cmpObj1->_densityMapCor[iter] = densityMapRotated[iter];
6064  }
6065  }
6066 
6067  //======================================== Clear memory
6068  delete[] densityMapRotated;
6069 
6070  //======================================== Done
6071  return ;
6072 
6073 }
6074 
6086 {
6087  //======================================== Sanity checks
6088  if ( ( this->all->size() < 1 ) || ( ( this->all->size() % 2 ) != 0 ) )
6089  {
6090  std::cerr << "!!! ProSHADE ERROR !!! Error in file comparison !!! The compare against part has incorrect number of entries." << std::endl;
6091  exit ( -1 );
6092  }
6093 
6094  int checkValue = static_cast<double> ( this->one->_keepOrRemove );
6095  for ( unsigned int strIt = 0; strIt < static_cast<unsigned int> ( this->all->size() ); strIt += 2 )
6096  {
6097  if ( checkValue != static_cast<double> ( this->all->at(strIt)->_keepOrRemove ) )
6098  {
6099  std::cerr << "!!! ProSHADE ERROR !!! Error in file comparison between " << this->one->_inputFileName << " AND " << this->all->at(strIt)->_inputFileName << " !!! The phase treatment is different, use the same phase treatment (i.e. remove all, or keep all) to make the strucutres comparable." << std::endl;
6100  exit ( -1 );
6101  }
6102  }
6103 
6104  checkValue = static_cast<double> ( this->two->_keepOrRemove );
6105  for ( unsigned int strIt = 1; strIt < static_cast<unsigned int> ( this->all->size() ); strIt += 2 )
6106  {
6107  if ( checkValue != static_cast<double> ( this->all->at(strIt)->_keepOrRemove ) )
6108  {
6109  std::cerr << "!!! ProSHADE ERROR !!! Error in file comparison between " << this->one->_inputFileName << " AND " << this->all->at(strIt)->_inputFileName << " !!! The phase treatment is different, use the same phase treatment (i.e. remove all, or keep all) to make the strucutres comparable." << std::endl;
6110  exit ( -1 );
6111  }
6112  }
6113 
6114  //======================================== Initialise
6115  std::vector<ProSHADE_internal::ProSHADE_data*> newAll;
6116  std::vector<double> pattCorrelationVec;
6117  std::vector<std::array<double,4>> translationVec;
6118  int totStrIt = 0;
6119  bool reverseOrder = false;
6120 
6121  //======================================== Find phaseless translation
6122  for ( unsigned int strIt = 0; strIt < static_cast<unsigned int> ( this->all->size() ); strIt += 2 )
6123  {
6124  if ( this->one->_keepOrRemove == false )
6125  {
6126  //================================ Pick the right range of objects (phase-wise)
6127  if ( ( strIt == 0 ) && ( this->all->at(strIt)->_keepOrRemove != false ) )
6128  {
6129  strIt += 1;
6130  reverseOrder = true;
6131  }
6132 
6133  //================================ Create Patterson comparison object
6134  ProSHADE_comparePairwise* cmpObj = new ProSHADE_comparePairwise ( this->one,
6135  this->all->at(strIt),
6136  settings->mPower,
6137  settings->ignoreLs,
6138  settings->glIntegOrder,
6139  settings );
6140 
6141  //================================ Get the angles and correlation
6142  pattCorrelationVec.emplace_back ( 0.0 );
6143  cmpObj->precomputeTrSigmaDescriptor ( );
6144  cmpObj->getSO3InverseMap ( settings );
6145  double xMov1 = 0.0;
6146  double yMov1 = 0.0;
6147  double zMov1 = 0.0;
6148  std::array<double,3> euAngs = cmpObj->getEulerAngles ( settings, &pattCorrelationVec.at(totStrIt) );
6149 
6150  //================================ Clear the memory
6151  delete cmpObj;
6152 
6153  //================================ Create the object for rotation
6154  if ( reverseOrder )
6155  {
6156  cmpObj = new ProSHADE_comparePairwise ( this->all->at(strIt-1),
6157  this->all->at(strIt-1),
6158  settings->mPower,
6159  settings->ignoreLs,
6160  settings->glIntegOrder,
6161  settings );
6162  }
6163  else
6164  {
6165  cmpObj = new ProSHADE_comparePairwise ( this->all->at(strIt+1),
6166  this->all->at(strIt+1),
6167  settings->mPower,
6168  settings->ignoreLs,
6169  settings->glIntegOrder,
6170  settings );
6171  }
6172 
6173  //================================ Find 'opposite' Euler angles and set them
6174  // Get mat from Euler
6175  double mat22 = cos ( euAngs[1] );
6176  double mat02 = -sin ( euAngs[1] ) * cos ( euAngs[2] );
6177  double mat12 = sin ( euAngs[1] ) * sin ( euAngs[2] );
6178  double mat20 = cos ( euAngs[0] ) * sin ( euAngs[1] );
6179  double mat21 = sin ( euAngs[0] ) * sin ( euAngs[1] );
6180 
6181  // Transpose
6182  double rMat02 = mat20;
6183  double rMat20 = mat02;
6184  double rMat21 = mat12;
6185  double rMat12 = mat21;
6186 
6187  // Get Euler from map
6188  euAngs[0] = atan2 ( rMat21, rMat20 );
6189  euAngs[1] = acos ( mat22 );
6190  euAngs[2] = atan2 ( rMat12, -rMat02 );
6191 
6192  if ( euAngs[0] < 0.0 ) { euAngs[0]= 2.0 * M_PI + euAngs[0]; }
6193  if ( euAngs[1] < 0.0 ) { euAngs[1]= M_PI + euAngs[1]; }
6194  if ( euAngs[2] < 0.0 ) { euAngs[2]= 2.0 * M_PI + euAngs[2]; }
6195 
6196  cmpObj->setEulerAngles ( euAngs[0], euAngs[1], euAngs[2] );
6197 
6198  //================================ Rotate by the Euler
6199  if ( reverseOrder )
6200  {
6201  cmpObj->rotateStructure ( this->all->at(strIt-1),
6202  settings,
6203  settings->clearMapFile,
6204  settings->verbose,
6205  settings->axisOrder,
6206  true );
6207 
6208  ProSHADE_data str1Copy = ProSHADE_data ( this->one );
6209  translationVec.emplace_back ( cmpObj->getTranslationFunctionMap ( &str1Copy, this->all->at(strIt-1), &xMov1, &yMov1, &zMov1 ) );
6210 
6211  this->all->at(strIt-1)->translateMap ( -translationVec.at(totStrIt)[0],
6212  -translationVec.at(totStrIt)[1],
6213  -translationVec.at(totStrIt)[2] );
6214  }
6215  else
6216  {
6217  cmpObj->rotateStructure ( this->all->at(strIt+1),
6218  settings,
6219  settings->clearMapFile,
6220  settings->verbose,
6221  settings->axisOrder,
6222  true );
6223 
6224  ProSHADE_data str1Copy = ProSHADE_data ( this->one );
6225  translationVec.emplace_back ( cmpObj->getTranslationFunctionMap ( &str1Copy, this->all->at(strIt+1), &xMov1, &yMov1, &zMov1 ) );
6226 
6227  this->all->at(strIt+1)->translateMap ( -translationVec.at(totStrIt)[0],
6228  -translationVec.at(totStrIt)[1],
6229  -translationVec.at(totStrIt)[2] );
6230  }
6231  }
6232  else
6233  {
6234  //================================ Pick the right range of objects (phase-wise)
6235  if ( ( strIt == 0 ) && ( this->all->at(strIt)->_keepOrRemove != false ) )
6236  {
6237  strIt += 1;
6238  reverseOrder = true;
6239  }
6240 
6241  //================================ Create Patterson comparison object
6242  ProSHADE_comparePairwise* cmpObj = new ProSHADE_comparePairwise ( this->two,
6243  this->all->at(strIt),
6244  settings->mPower,
6245  settings->ignoreLs,
6246  settings->glIntegOrder,
6247  settings );
6248 
6249  //================================ Get the angles and correlation
6250  pattCorrelationVec.emplace_back ( 0.0 );
6251  cmpObj->precomputeTrSigmaDescriptor ( );
6252  cmpObj->getSO3InverseMap ( settings );
6253  double xMov1 = 0.0;
6254  double yMov1 = 0.0;
6255  double zMov1 = 0.0;
6256  std::array<double,3> euAngs = cmpObj->getEulerAngles ( settings, &pattCorrelationVec.at(totStrIt) );
6257 
6258  //================================ Clear the memory
6259  delete cmpObj;
6260 
6261  //================================ Create the object for rotation
6262  if ( reverseOrder )
6263  {
6264  cmpObj = new ProSHADE_comparePairwise ( this->all->at(strIt-1),
6265  this->all->at(strIt-1),
6266  settings->mPower,
6267  settings->ignoreLs,
6268  settings->glIntegOrder,
6269  settings );
6270  }
6271  else
6272  {
6273  cmpObj = new ProSHADE_comparePairwise ( this->all->at(strIt+1),
6274  this->all->at(strIt+1),
6275  settings->mPower,
6276  settings->ignoreLs,
6277  settings->glIntegOrder,
6278  settings );
6279  }
6280 
6281  //================================ Find 'opposite' Euler angles and set them
6282  // Get mat from Euler
6283  double mat22 = cos ( euAngs[1] );
6284  double mat02 = -sin ( euAngs[1] ) * cos ( euAngs[2] );
6285  double mat12 = sin ( euAngs[1] ) * sin ( euAngs[2] );
6286  double mat20 = cos ( euAngs[0] ) * sin ( euAngs[1] );
6287  double mat21 = sin ( euAngs[0] ) * sin ( euAngs[1] );
6288 
6289  // Transpose
6290  double rMat02 = mat20;
6291  double rMat20 = mat02;
6292  double rMat21 = mat12;
6293  double rMat12 = mat21;
6294 
6295  // Get Euler from map
6296  euAngs[0] = atan2 ( rMat21, rMat20 );
6297  euAngs[1] = acos ( mat22 );
6298  euAngs[2] = atan2 ( rMat12, -rMat02 );
6299 
6300  if ( euAngs[0] < 0.0 ) { euAngs[0]= 2.0 * M_PI + euAngs[0]; }
6301  if ( euAngs[1] < 0.0 ) { euAngs[1]= M_PI + euAngs[1]; }
6302  if ( euAngs[2] < 0.0 ) { euAngs[2]= 2.0 * M_PI + euAngs[2]; }
6303 
6304  cmpObj->setEulerAngles ( euAngs[0], euAngs[1], euAngs[2] );
6305 
6306  //================================ Rotate by the Euler
6307  if ( reverseOrder )
6308  {
6309  cmpObj->rotateStructure ( this->all->at(strIt-1),
6310  settings,
6311  settings->clearMapFile,
6312  settings->verbose,
6313  settings->axisOrder,
6314  true );
6315 
6316  ProSHADE_data str1Copy = ProSHADE_data ( this->two );
6317  translationVec.emplace_back ( cmpObj->getTranslationFunctionMap ( &str1Copy, this->all->at(strIt-1), &xMov1, &yMov1, &zMov1 ) );
6318 
6319  this->all->at(strIt-1)->translateMap ( -translationVec.at(totStrIt)[0],
6320  -translationVec.at(totStrIt)[1],
6321  -translationVec.at(totStrIt)[2] );
6322  }
6323  else
6324  {
6325  cmpObj->rotateStructure ( this->all->at(strIt+1),
6326  settings,
6327  settings->clearMapFile,
6328  settings->verbose,
6329  settings->axisOrder,
6330  true );
6331 
6332  ProSHADE_data str1Copy = ProSHADE_data ( this->two );
6333  translationVec.emplace_back ( cmpObj->getTranslationFunctionMap ( &str1Copy, this->all->at(strIt+1), &xMov1, &yMov1, &zMov1 ) );
6334 
6335  this->all->at(strIt+1)->translateMap ( -translationVec.at(totStrIt)[0],
6336  -translationVec.at(totStrIt)[1],
6337  -translationVec.at(totStrIt)[2] );
6338  }
6339  }
6340 
6341  totStrIt += 1;
6342  }
6343 
6344  //======================================== Delete the Patterson maps, they are no longer required
6345  for ( signed int strIt = static_cast<signed int> ( this->all->size() - 1 ); strIt >= 0 ; strIt-- )
6346  {
6347  if ( this->all->at(strIt)->_keepOrRemove == false )
6348  {
6349  delete this->all->at(strIt);
6350  this->all->at(strIt) = nullptr;
6351  this->all->erase ( this->all->begin ( ) + strIt );
6352  }
6353  }
6354 
6355  //======================================== Done
6356  return ;
6357 }
6358 
6359 
6360 
6361 
6362 
6363 
6364 
6365 
This class deals with reading in the data and computing structure specific information including the ...
void generateWignerMatrices(ProSHADE::ProSHADE_settings *settings)
This function is responsible for computing the Wigner D matrices for symmetry detection purposes...
std::vector< std::array< double, 8 > > getSO3Peaks(ProSHADE::ProSHADE_settings *settings, double noIQRs=1.5, bool freeMem=true, int peakSize=1, double realDist=0.4, int verbose=1)
This function detects &#39;true&#39; peaks from the background of the SO3 inverse transform map...
std::vector< std::vector< std::array< double, 5 > > > findTetrSymmetry(std::vector< std::array< double, 5 > > CnSymm, double *tetrSymmPeakAvg, double axisErrorTolerance=0.1)
This function detects the tetrahedral symmetries in the structure.
void rotateStructure(ProSHADE_data *cmpObj1, ProSHADE::ProSHADE_settings *settings, std::string saveName, int verbose=0, std::string axOrd="xyz", bool internalUse=false)
Rotates the density map be given Angle-Axis rotation using the Wigner matrices and inverse spharical ...
std::string clearMapFile
If map features are to be extracted, should the clear map be saved (then give file name here)...
Definition: ProSHADE.h:144
std::vector< std::vector< std::array< double, 5 > > > findIcosSymmetry(std::vector< std::array< double, 5 > > CnSymm, double *icosSymmPeakAvg, double axisErrorTolerance=0.1)
This function detects the icosahedral symmetries in the structure.
void alignDensities(ProSHADE::ProSHADE_settings *settings)
Takes the internal objects with and without phases and aligns them to all the other objects...
bool htmlReport
Should HTML report for the run be created?
Definition: ProSHADE.h:186
std::vector< std::array< double, 5 > > findCnSymmetry(std::vector< std::array< double, 8 > > peaks, ProSHADE::ProSHADE_settings *settings, double axisErrorTolerance=0.0, bool freeMem=true, double percAllowedToMiss=0.33, int verbose=1)
This function attempts to detect the C-symmetries present in the list of overlay peaks.
void setEulerAngles(double alpha, double beta, double gamma)
This function is used to set the Euler angles for further processing.
int verbose
Should the software report on the progress, or just be quiet? Value between 0 (quiet) and 4 (loud) ...
Definition: ProSHADE.h:189
std::vector< std::array< double, 5 > > generateTetrElements(std::vector< std::array< double, 5 > > symmAxes, ProSHADE::ProSHADE_settings *settings, int verbose=0)
This function generates the 12 unique tetrahedral symmetry group elements from the already detected a...
The main header file containing all declarations the user of the library needs.
double getRotCoeffDistance(ProSHADE::ProSHADE_settings *settings)
This function computes the full rotation function descriptor distances.
void precomputeTrSigmaDescriptor()
This function computes the E matrices required for the trace sigma descriptor, the rotation function ...
std::vector< std::array< double, 5 > > generateTetrAxes(ProSHADE_comparePairwise *cmpObj, ProSHADE::ProSHADE_settings *settings, std::vector< std::array< double, 5 > > tetrSymm, std::vector< std::array< double, 5 > > allCs, double axisErrorTolerance=0.1, int verbose=0)
This function generates all the tetrahedral symmetry axes from a pair detected by findTetrSymmetry fu...
std::vector< std::vector< std::array< double, 5 > > > findOctaSymmetry(std::vector< std::array< double, 5 > > CnSymm, double *octaSymmPeakAvg, double axisErrorTolerance=0.1)
This function detects the octahedral symmetries in the structure.
void getSO3InverseMap(ProSHADE::ProSHADE_settings *settings)
This function is responsible for computing the SO3 inverse transform.
double maxAvgPeakForSymmetry(double X, double Y, double Z, double Angle, ProSHADE::ProSHADE_settings *settings)
This function takes angle and axis and checks the SO(3) map for this specific symmetry only...
double getPeakThreshold()
This is an accessor function to the peak height threshold.
std::vector< int > ignoreLs
This vector lists all the bandwidth values which should be ignored and not part of the computations...
Definition: ProSHADE.h:117
std::vector< std::array< double, 5 > > generateOctaElements(std::vector< std::array< double, 5 > > symmAxes, ProSHADE::ProSHADE_settings *settings, int verbose=0)
This function generates the 24 octahedral symmetry group elements from the axes generated by generate...
std::vector< std::vector< std::array< double, 6 > > > findDnSymmetry(std::vector< std::array< double, 5 > > CnSymm, double axisErrorTolerance=0.1)
This function detects dihedral (D) symmetries from the list of already detected C symmetries...
std::vector< std::array< double, 5 > > generateOctaAxes(ProSHADE_comparePairwise *cmpObj, ProSHADE::ProSHADE_settings *settings, std::vector< std::array< double, 5 > > octaSymm, std::vector< std::array< double, 5 > > allCs, double axisErrorTolerance=0.1, int verbose=0)
This function generates all the ctahedral symmetry elements from a pair detected by findOctaSymmetry ...
This is the executive class responsible for comparing strictly two structures.
std::vector< std::array< double, 5 > > findCnSymmetryClear(std::vector< std::array< double, 5 > > CnSymm, ProSHADE::ProSHADE_settings *settings, double maxGap=0.2, bool *pf=nullptr)
This function takes the detected C-symmetries list, removes low probability symmetries and sorts as t...
std::vector< std::array< double, 5 > > generateIcosAxes(ProSHADE_comparePairwise *cmpObj, ProSHADE::ProSHADE_settings *settings, std::vector< std::array< double, 5 > > icosSymm, std::vector< std::array< double, 5 > > allCs, double axisErrorTolerance=0.1, int verbose=0)
This function generates all the icosahedral symmetry elements from a pair detected by findIcosSymmetr...
void getSO3InverseMap(ProSHADE::ProSHADE_settings *settings)
This function is responsible for computing the SO3 inverse transform.
std::vector< std::array< double, 5 > > generateIcosElements(std::vector< std::array< double, 5 > > symmAxes, ProSHADE::ProSHADE_settings *settings, int verbose=0)
This function generates the 60 icosahedral symmetry group elements from the axes generated by generat...
std::string axisOrder
A string specifying the order of the axis. Must have three characters and any permutation of &#39;x&#39;...
Definition: ProSHADE.h:183
The main header file containing all declarations for the innter workings of the library.
bool checkPeakExistence(double alpha, double beta, double gamma, ProSHADE::ProSHADE_settings *settings, double passThreshold=0.3)
This function checks if a peak Euler angles actually produce something the descriptors would agree is...
int htmlReportLineProgress
Iterator for current HTML line in the progress bar.
Definition: ProSHADE.h:188
ProSHADE_comparePairwise(ProSHADE_data *cmpObj1, ProSHADE_data *cmpObj2, double mPower, std::vector< int > ignoreL, unsigned int order, ProSHADE::ProSHADE_settings *settings)
Contructor for the ProSHADE_comparePairwise class.
std::vector< std::vector< std::array< double, 6 > > > findDnSymmetryClear(std::vector< std::vector< std::array< double, 6 > > > DnSymm, ProSHADE::ProSHADE_settings *settings, double maxGap=0.2, bool *pf=nullptr)
This function sorts the D symmetries and removed low probability ones.
std::vector< std::array< double, 3 > > getEulerAngles(ProSHADE::ProSHADE_settings *settings)
This function finds the highest peak in the SO3 inverse transform map and sets it as the optimal over...
This class stores all the settings and is passed to the executive classes instead of multitude of par...
Definition: ProSHADE.h:74
void freeInvMap()
This function frees the SOFT inverse map.
double mPower
This parameter determines the scaling for trace sigma descriptor.
Definition: ProSHADE.h:114
std::array< double, 4 > getTranslationFunctionMap(ProSHADE_data *obj1, ProSHADE_data *obj2, double *ob2XMov=nullptr, double *ob2YMov=nullptr, double *ob2ZMov=nullptr)
Computes the translation function for the two objects and returns the position of the highest peak...
This file contains the ProSHADE_internal_misc namespace and its miscellaneous functions.
std::array< double, 3 > getEulerAngles(ProSHADE::ProSHADE_settings *settings, double *correlation=nullptr)
This function finds the highest peak in the SO3 inverse transform map and sets it as the optimal over...
unsigned int glIntegOrder
This parameter controls the Gauss-Legendre integration order and so the radial resolution.
Definition: ProSHADE.h:82