ProSHADE  0.7.6.1 (AUG 2021)
Protein Shape Detection
ProSHADE_symmetry.cpp
Go to the documentation of this file.
1 
22 //==================================================== ProSHADE
23 #include "ProSHADE_symmetry.hpp"
24 
25 //==================================================== Local functions prototypes
26 proshade_double determinePeakThreshold ( std::vector < proshade_double > inArr, proshade_double noIQRsFromMedian );
27 bool sortProSHADESymmetryByPeak ( proshade_double* a, proshade_double* b );
28 std::vector < std::pair< proshade_unsign, proshade_unsign > > findBestIcosDihedralPair ( std::vector< proshade_double* >* CSymList, proshade_double minPeakHeight, proshade_double axErr );
29 std::pair< proshade_unsign, proshade_unsign > findBestOctaDihedralPair ( std::vector< proshade_double* >* CSymList, proshade_double minPeakHeight, proshade_double axErr );
30 std::pair< proshade_unsign, proshade_unsign > findBestTetraDihedralPair ( std::vector< proshade_double* >* CSymList, proshade_double minPeakHeight, proshade_double axErr );
31 
42 {
43  //================================================ Report progress
44  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting self-rotation function computation." );
45 
46  //================================================ Compute un-weighted E matrices and their weights
47  ProSHADE_internal_distances::computeEMatrices ( this, this, settings );
48 
49  //================================================ Normalise E matrices by the magnitudes
50  ProSHADE_internal_distances::normaliseEMatrices ( this, this, settings );
51 
52  //================================================ Generate SO(3) coefficients
54 
55  //================================================ Compute the inverse SO(3) Fourier Transform (SOFT) on the newly computed coefficients
57 
58  //================================================ Report completion
59  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Self-rotation function obtained." );
60 
61  //================================================ Done
62  return ;
63 
64 }
65 
71 proshade_double determinePeakThreshold ( std::vector < proshade_double > inArr, proshade_double noIQRsFromMedian )
72 {
73  //================================================ Initialise variables
74  proshade_double ret = 0.0;
75  proshade_unsign vecSize = static_cast< proshade_unsign > ( inArr.size() );
76  proshade_double* meadianAndIQR = new proshade_double[2];
77 
78  //================================================ Deal with low number of input cases
79  if ( vecSize == 0 ) { delete[] meadianAndIQR; return ( ret ); } // Return 0
80  if ( vecSize <= 4 ) { ret = std::accumulate ( inArr.begin(), inArr.end(), 0.0 ) / static_cast< proshade_double > ( vecSize ); } // Return mean
81 
82  //================================================ Deal with reasonable number in input cases
83  else
84  {
85  //============================================ Allocate memory for median and IQR computation
86  ProSHADE_internal_misc::checkMemoryAllocation ( meadianAndIQR, __FILE__, __LINE__, __func__ );
87 
88  //============================================ Find median and IQR
89  ProSHADE_internal_maths::vectorMedianAndIQR ( &inArr, meadianAndIQR );
90 
91  //============================================ Get the threshold
92  ret = meadianAndIQR[0] + ( meadianAndIQR[1] * noIQRsFromMedian );
93  }
94 
95  //================================================ Sanity checks
96  if ( ret > *( std::max_element ( inArr.begin(), inArr.end() ) ) )
97  {
98  ret = *( std::max_element ( inArr.begin(), inArr.end() ) );
99  }
100 
101  //================================================ Release memory
102  delete[] meadianAndIQR;
103 
104  //================================================ Done
105  return ( ret );
106 
107 }
108 
121 {
122  //================================================ Report progress
123  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting self-rotation function conversion to angle-axis representation." );
124 
125  //================================================ Initialise variables
126  proshade_double shellSpacing = ( 2.0 * M_PI ) / static_cast<proshade_double> ( this->maxShellBand ) * 2.0;
127  std::vector< proshade_double > allPeakHeights;
128 
129  //================================================ Initialise the spheres
130  for ( proshade_unsign spIt = 1; spIt < ( this->maxShellBand * 2 ); spIt++ )
131  {
132  this->sphereMappedRotFun.emplace_back ( new ProSHADE_internal_spheres::ProSHADE_rotFun_sphere( static_cast<proshade_double> ( spIt ) * shellSpacing,
133  shellSpacing,
134  this->maxShellBand * 2,
135  static_cast<proshade_double> ( spIt ) * shellSpacing,
136  spIt - 1 ) );
137  }
138 
139  //================================================ Interpolate the rotation function onto the spheres
140  for ( proshade_unsign shIt = 0; shIt < static_cast<proshade_unsign> ( sphereMappedRotFun.size() ); shIt++ )
141  {
142  //============================================ Report progress
143  std::stringstream hlpSS;
144  hlpSS << "Interpolating sphere " << shIt << " ( radius: " << this->sphereMappedRotFun.at(shIt)->getRadius() << " ).";
145  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 3, hlpSS.str() );
146 
147  //============================================ Interpolate onto spheres
148  this->sphereMappedRotFun.at(shIt)->interpolateSphereValues ( this->getInvSO3Coeffs ( ) );
149  }
150 
151  //================================================ Report completion
152  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Self-rotation function converted to spherical angle-axis space." );
153  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Started peak detection on the angle-axis spheres." );
154 
155  //================================================ Find all peaks in the sphere grids
156  for ( proshade_unsign shIt = 0; shIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.size() ); shIt++ )
157  {
158  this->sphereMappedRotFun.at(shIt)->findAllPeaks ( static_cast< proshade_signed > ( settings->peakNeighbours ), &allPeakHeights );
159  }
160 
161  //================================================ Report progress
162  std::stringstream hlpSS;
163  hlpSS << "Detected " << allPeakHeights.size() << " peaks with any height.";
164  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 3, hlpSS.str() );
165 
166  //================================================ Compute threshold for small peaks
167  proshade_double peakThres = std::max ( settings->minSymPeak, determinePeakThreshold ( allPeakHeights, settings->noIQRsFromMedianNaivePeak ) );
168 
169  //================================================ Report progress
170  std::stringstream hlpSS2;
171  hlpSS2 << "From these peaks, decided the threshold will be " << peakThres << " peak height.";
172  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 4, hlpSS2.str() );
173 
174  //================================================ Remove too small peaks
175  for ( proshade_unsign shIt = 0; shIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.size() ); shIt++ )
176  {
177  this->sphereMappedRotFun.at(shIt)->removeSmallPeaks ( peakThres );
178  }
179 
180  //================================================ Report progress
181  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 3, "Peaks detected for all spheres." );
182 
183  //================================================ Done
184  return ;
185 
186 }
187 
199 {
200  //================================================ Initialise variables
201  std::vector< proshade_double* > ret;
202 
203  //================================================ Report progress
204  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting C symmetry detection." );
205 
206  //================================================ Get list of peaks in the self-rotation map
207  std::vector< proshade_double* > allPeaks = ProSHADE_internal_peakSearch::getAllPeaksNaive ( this->getInvSO3Coeffs (), this->getMaxBand() * 2,
208  static_cast< proshade_signed > ( settings->peakNeighbours ),
209  settings->noIQRsFromMedianNaivePeak );
210 
211  //================================================ Convert peaks to angle-axis
212  std::vector< proshade_double* > peaksAA = ProSHADE_internal_symmetry::getPeaksAngleAxisPositions ( allPeaks, settings->verbose );
213 
214  //================================================ Sort peaks by height groups
215  std::vector< proshade_double > peakGroupsBoundaries = ProSHADE_internal_symmetry::findPeaksByHeightBoundaries ( peaksAA, settings->smoothingFactor );
216 
217  //================================================ Get symmetry per group
218  std::vector< std::vector< proshade_unsign > > detectedCSymmetries;
219  for ( proshade_signed iter = static_cast< proshade_signed > ( peakGroupsBoundaries.size() - 1 ); iter >= 0; iter-- )
220  {
221  //============================================ Get peaks group peaks only
222  std::vector< proshade_double* > symPeaks;
223  for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( peaksAA.size() ); it++ )
224  {
225  if ( peaksAA.at(it)[4] > peakGroupsBoundaries.at( static_cast< size_t > ( iter ) ) ) { ProSHADE_internal_misc::addToDblPtrVector ( &symPeaks, peaksAA.at(it) ); }
226  }
227 
228  //============================================ Search for symmetry in these peaks
229  detectedCSymmetries = ProSHADE_internal_symmetry::findPeaksCSymmetry ( &symPeaks, settings->verbose,
230  this->getMaxBand(),
231  settings->symMissPeakThres,
232  settings->axisErrTolerance,
233  settings->axisErrToleranceDefault,
234  this );
235 
236  //============================================ Print detected symmetries
237  for ( proshade_unsign detIt = 0; detIt < static_cast<proshade_unsign> ( detectedCSymmetries.size() ); detIt++ ) { ProSHADE_internal_symmetry::printSymmetryGroup ( detectedCSymmetries.at(detIt), symPeaks, settings->verbose ); }
238 
239  //============================================ Save detected
240  ProSHADE_internal_symmetry::saveAllCSymmetries ( detectedCSymmetries, symPeaks, &ret, settings->axisErrTolerance );
241  }
242 
243  //================================================ Release memory
244  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( peaksAA.size() ); iter++ ) { delete[] allPeaks.at(iter); delete[] peaksAA.at(iter); }
245 
246  //================================================ Report completion
247  ProSHADE_internal_symmetry::printSymmetryCompletion ( static_cast<proshade_unsign>( ret.size() ), settings->verbose );
248 
249  //================================================ Done
250  return ( ret );
251 
252 }
253 
263 std::vector< proshade_double* > ProSHADE_internal_symmetry::getPeaksAngleAxisPositions ( std::vector< proshade_double* > allPeaks, proshade_signed verbose )
264 {
265  //================================================ Initialise variables
266  std::vector< proshade_double* > ret;
267  proshade_double* hlpP = nullptr;
268  proshade_double* rotMat = new proshade_double [9];
269  ProSHADE_internal_misc::checkMemoryAllocation ( rotMat, __FILE__, __LINE__, __func__ );
270 
271  //================================================ For each peak
272  for ( proshade_unsign peakIter = 0; peakIter < static_cast<proshade_unsign> ( allPeaks.size() ); peakIter++ )
273  {
274  //============================================ Convert Euler ZXZ angles to rotation matrix
275  ProSHADE_internal_maths::getRotationMatrixFromEulerZXZAngles ( allPeaks.at(peakIter)[0], allPeaks.at(peakIter)[1], allPeaks.at(peakIter)[2], rotMat );
276 
277  //============================================ Allocate pointer to results
278  hlpP = new proshade_double [5];
279  ProSHADE_internal_misc::checkMemoryAllocation ( hlpP, __FILE__, __LINE__, __func__ );
280 
281  //============================================ Convert rotation matrix to Angle-axis reporesentation
282  ProSHADE_internal_maths::getAxisAngleFromRotationMatrix ( rotMat, &hlpP[0], &hlpP[1], &hlpP[2], &hlpP[3] );
283  hlpP[4] = allPeaks.at(peakIter)[3];
284 
285  //============================================ Save results
287  }
288 
289  //================================================ Release memory
290  delete[] rotMat;
291 
292  //================================================ Report progress
293  std::stringstream hlpSSP;
294  hlpSSP << "Found " << ret.size() << " possible peaks.";
295  ProSHADE_internal_messages::printProgressMessage ( verbose, 3, hlpSSP.str() );
296 
297  //================================================ Warning if no peaks!
298  if ( ret.size() < 1 )
299  {
300  ProSHADE_internal_messages::printWarningMessage ( verbose, "!!! ProSHADE WARNING !!! Failed to detect any symmetries. There are no reasonable peaks in the self-rotation map. If you believe there should be some symmetry, you can try decreasing the resolution or changing the peak IQR threshold.", "WS00029" );
301  }
302 
303  //================================================ Done
304  return ( ret );
305 
306 }
307 
317 std::vector< proshade_double > ProSHADE_internal_symmetry::findPeaksByHeightBoundaries ( std::vector< proshade_double* > allPeaks, proshade_double smoothing )
318 {
319  //================================================ Initialise variables
320  std::vector< proshade_double > boundaries;
321  ProSHADE_internal_misc::addToDoubleVector ( &boundaries, 0.0 );
322  proshade_double peakContribution = 0.0;
323 
324  //================================================ Generate Probability Density function (PDF)
325  std::vector< proshade_double > pdf;
326  for ( proshade_double iter = 0.0; iter <= 1.0; iter += 0.01 )
327  {
328  //============================================ Initialise point
329  peakContribution = 0.0;
330 
331  //============================================ Sum peak contributions
332  for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( allPeaks.size() ); it++ )
333  {
334  peakContribution += ProSHADE_internal_maths::normalDistributionValue ( allPeaks.at(it)[4], smoothing, iter );
335  }
336 
337  //============================================ Save result
338  ProSHADE_internal_misc::addToDoubleVector ( &pdf, peakContribution );
339  }
340 
341  //================================================ Find boundaries
342  proshade_double prev = pdf.at(0);
343  for ( proshade_unsign iter = 1; iter < static_cast<proshade_unsign> ( pdf.size() - 1 ); iter ++ )
344  {
345  //============================================ Check for local minima
346  if ( ( prev > pdf.at(iter) ) && ( pdf.at(iter+1) > pdf.at(iter) ) )
347  {
348  ProSHADE_internal_misc::addToDoubleVector ( &boundaries, static_cast< proshade_double > ( iter ) * 0.01 );
349  }
350 
351  //============================================ Prepare next iteration
352  prev = pdf.at(iter);
353  }
354 
355  //================================================ Done
356  return ( boundaries );
357 
358 }
359 
375 std::vector< std::vector< proshade_unsign > > ProSHADE_internal_symmetry::findPeaksCSymmetry ( std::vector< proshade_double* >* peaks, proshade_signed verbose, proshade_unsign band, proshade_double missPeakThres, proshade_double axisErrTolerance, bool axisErrToleranceDef, ProSHADE_internal_data::ProSHADE_data* dataObj )
376 {
377  //======================================== Initialise variables
378  std::vector< std::vector< proshade_unsign > > ret;
379  std::vector< proshade_double > triedAlready;
380  std::vector< proshade_unsign > angsToTry, testedAlready;
381  proshade_double angDist, angDivisionRemainder, angDivisionBasis, nextSymmetryError, nextPeakError = ( M_PI * 2.0 ) / ( static_cast<proshade_double> ( band ) * 2.0 );
382 
383  //================================================ Sanity check
384  if ( peaks->size() < 1 ) { return ( ret ); }
385 
386  //================================================ Group peaks by axes
387  std::vector< std::vector< proshade_unsign > > sameAxesGroups = ProSHADE_internal_symmetry::groupSameAxes ( *peaks, axisErrTolerance );
388 
389  //================================================ For each axis group
390  for ( proshade_unsign grpIt = 0; grpIt < static_cast<proshade_unsign> ( sameAxesGroups.size() ); grpIt++ )
391  {
392  //============================================ Print axis group if need be
393  ProSHADE_internal_symmetry::printSymmetryPeaks ( sameAxesGroups.at(grpIt), *peaks, verbose, grpIt );
394 
395  //============================================ While there are distances between rotation angles in the group
396  triedAlready.clear ( );
397  testedAlready.clear ( );
398  while ( ProSHADE_internal_symmetry::smallestDistanceBetweenAngles ( sameAxesGroups.at(grpIt), *peaks, &triedAlready, &angDist ) )
399  {
400  //======================================== Check if testable fold value exists, otherwise test other distances
401  angsToTry.clear ( );
402  if ( !ProSHADE_internal_symmetry::determineFoldToTry ( angDist, &angDivisionBasis, &angDivisionRemainder, nextPeakError, &nextSymmetryError, &angsToTry ) ) { continue; }
403 
404  //======================================== If reasonable folds are found, test these for being complete symmetries
405  ProSHADE_internal_symmetry::findSymmetryUsingFold ( dataObj, &angsToTry, &sameAxesGroups.at(grpIt), peaks, &ret, &testedAlready, axisErrTolerance, axisErrToleranceDef, missPeakThres, verbose );
406  }
407 
408  }
409 
410  //================================================ Done
411  return ( ret );
412 
413 }
414 
426 std::vector< std::vector< proshade_unsign > > ProSHADE_internal_symmetry::groupSameAxes ( std::vector< proshade_double* >& peaks, proshade_double errTolerance )
427 {
428  //================================================ Initialise variables
429  std::vector< std::vector< proshade_unsign > > ret;
430  bool sameAxisFound = false;
431  proshade_double angTolerance = std::acos ( 1.0 - errTolerance );
432 
433  //================================================ Set all largest axis value to positive (this will make the 0,0,1 and 0,0,-1 axes the same)
435 
436  //================================================ For each axis
437  for ( proshade_unsign peakIter = 0; peakIter < static_cast<proshade_unsign> ( peaks.size() ); peakIter++ )
438  {
439  //============================================ Initialise variables for next peak
440  sameAxisFound = false;
441 
442  //============================================ Ignore zero angle peaks
443  if ( ( peaks.at(peakIter)[3] - angTolerance <= 0.0 ) && ( peaks.at(peakIter)[3] + angTolerance > 0.0 ) ) { continue; }
444 
445  //============================================ Ignore very small axis peaks - the axis may be wrong here.
446  // !! The value of 0.1 is hardcoded, but arbitrary
447  if ( ( ( peaks.at(peakIter)[0] - 0.1 <= 0.0 ) && ( peaks.at(peakIter)[0] + 0.1 > 0.0 ) ) &&
448  ( ( peaks.at(peakIter)[1] - 0.1 <= 0.0 ) && ( peaks.at(peakIter)[1] + 0.1 > 0.0 ) ) &&
449  ( ( peaks.at(peakIter)[2] - 0.1 <= 0.0 ) && ( peaks.at(peakIter)[2] + 0.1 > 0.0 ) ) ) { continue; }
450 
451  //============================================ Compare to all already detected axes groups
452  for ( proshade_unsign sameAxisGrp = 0; sameAxisGrp < static_cast<proshade_unsign> ( ret.size() ); sameAxisGrp++ )
453  {
454  //======================================== and all their members
455  for ( proshade_unsign sameAxis = 0; sameAxis < static_cast<proshade_unsign> ( ret.at(sameAxisGrp).size() ); sameAxis++ )
456  {
457  //==================================== Is this identical axis to the tested one?
458  if ( ProSHADE_internal_maths::vectorOrientationSimilarity ( peaks.at(ret.at(sameAxisGrp).at(sameAxis))[0],
459  peaks.at(ret.at(sameAxisGrp).at(sameAxis))[1],
460  peaks.at(ret.at(sameAxisGrp).at(sameAxis))[2],
461  peaks.at(peakIter)[0],
462  peaks.at(peakIter)[1],
463  peaks.at(peakIter)[2],
464  errTolerance ) )
465  {
466  sameAxisFound = true;
467  ProSHADE_internal_misc::addToUnsignVector ( &ret.at(sameAxisGrp), peakIter );
468  break;
469  }
470  }
471  }
472 
473  //============================================ If same axis was found, do nothing
474  if ( sameAxisFound ) { continue; }
475 
476  //============================================ No similar axis was found
477  std::vector<proshade_unsign> hlpVec;
478  ProSHADE_internal_misc::addToUnsignVector ( &hlpVec, peakIter );
480  }
481 
482  //================================================ Add zero peak to all axes
484 
485  //================================================ Done
486  return ( ret );
487 
488 }
489 
498 void ProSHADE_internal_symmetry::giveOppositeAxesSameDirection ( std::vector< proshade_double* > peaks )
499 {
500  //================================================ Apply to all peaks
501  for ( proshade_unsign i = 0; i < static_cast<proshade_unsign> ( peaks.size() ); i++ )
502  {
503  const FloatingPoint< proshade_double > lhs1 ( std::max ( std::abs ( peaks.at(i)[0] ), std::max( std::abs ( peaks.at(i)[1] ), std::abs ( peaks.at(i)[2] ) ) ) );
504  const FloatingPoint< proshade_double > rhs1 ( std::abs ( peaks.at(i)[0] ) );
505  const FloatingPoint< proshade_double > rhs2 ( std::abs ( peaks.at(i)[1] ) );
506  const FloatingPoint< proshade_double > rhs3 ( std::abs ( peaks.at(i)[2] ) );
507  if ( ( ( lhs1.AlmostEquals ( rhs1 ) ) && ( peaks.at(i)[0] < 0.0 ) ) ||
508  ( ( lhs1.AlmostEquals ( rhs2 ) ) && ( peaks.at(i)[1] < 0.0 ) ) ||
509  ( ( lhs1.AlmostEquals ( rhs3 ) ) && ( peaks.at(i)[2] < 0.0 ) ) )
510  {
511  peaks.at(i)[0] *= -1.0;
512  peaks.at(i)[1] *= -1.0;
513  peaks.at(i)[2] *= -1.0;
514  peaks.at(i)[3] *= -1.0;
515  }
516  }
517 
518  //================================================ Done
519  return ;
520 
521 }
522 
529 void ProSHADE_internal_symmetry::printSymmetryPeaks ( std::vector< proshade_unsign > grp, std::vector< proshade_double* > peaks, proshade_signed verbose, proshade_unsign groupNo )
530 {
531  //================================================ Symmetry group output header
532  std::stringstream hlpSS;
533  hlpSS << "Symmetry axis group " << groupNo;
534  ProSHADE_internal_messages::printProgressMessage ( verbose, 6, hlpSS.str() );
535  ProSHADE_internal_messages::printProgressMessage ( verbose, 6, "Peak index\t\tx\t y\t z\tAngle\tPeak heiht" );
536 
537  //================================================ Print the symmetry group
538  for ( proshade_unsign axIt = 0; axIt < static_cast<proshade_unsign> ( grp.size() ); axIt++ )
539  {
540  std::stringstream SS;
541  SS << " " << axIt << "\t " << static_cast<int>( peaks.at(grp.at(axIt))[0] * 100.0 ) / 100.0 << "\t" << static_cast<int>( peaks.at(grp.at(axIt))[1] * 100.0 ) / 100.0 << "\t" << static_cast<int>( peaks.at(grp.at(axIt))[2] * 100.0 ) / 100.0 << "\t" << static_cast<int>( peaks.at(grp.at(axIt))[3] * 100.0 ) / 100.0 << "\t" << static_cast<int>( peaks.at(grp.at(axIt))[4] * 100.0 ) / 100.0;
542  ProSHADE_internal_messages::printProgressMessage ( verbose, 6, SS.str() );
543  }
544 
545  //================================================ Done
546  return ;
547 
548 }
549 
562 bool ProSHADE_internal_symmetry::smallestDistanceBetweenAngles ( std::vector< proshade_unsign > grp, std::vector< proshade_double* > peaks, std::vector< proshade_double >* tried, proshade_double* dist )
563 {
564  //================================================ Initialise variables
565  bool ret = false;
566  bool skip = false;
567  proshade_unsign g1 = 0, g2 = 0;
568  *dist = 999.9;
569 
570  //================================================ For each group pair
571  for ( proshade_unsign gr1It = 0; gr1It < static_cast<proshade_unsign> ( grp.size() ); gr1It++ )
572  {
573  for ( proshade_unsign gr2It = 1; gr2It < static_cast<proshade_unsign> ( grp.size() ); gr2It++ )
574  {
575  //======================================== Unique pairs only
576  if ( gr1It >= gr2It ) { continue; }
577 
578  //======================================== Have we tried this already?
579  skip = false;
580  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( tried->size() ); iter += 3 )
581  {
582  //==================================== Avoid already tested combinations
583  const FloatingPoint< proshade_double > lhs1 ( static_cast< proshade_double > ( gr2It ) ), rhs1 ( tried->at( iter + 1 ) );
584  const FloatingPoint< proshade_double > lhs2 ( static_cast< proshade_double > ( gr1It ) ), rhs2 ( tried->at( iter ) );
585  if ( ( lhs1.AlmostEquals ( rhs1 ) ) && ( lhs2.AlmostEquals ( rhs2 ) ) ) { skip = true; }
586 
587  //==================================== Also avoid distances very close to already tested distances (no problem until approx C700)
588  if ( !skip &&
589  ( ( std::abs( std::abs ( peaks.at(grp.at(gr1It))[3] ) - std::abs ( peaks.at(grp.at(gr2It))[3] ) ) - 0.01 ) < tried->at( iter + 2 ) ) &&
590  ( ( std::abs( std::abs ( peaks.at(grp.at(gr1It))[3] ) - std::abs ( peaks.at(grp.at(gr2It))[3] ) ) + 0.01 ) > tried->at( iter + 2 ) ) )
591  {
592  skip = true;
593  }
594  }
595  if ( skip ) { continue; }
596 
597  //======================================== Is this the smallest distance?
598  if ( std::abs( std::abs ( peaks.at(grp.at(gr1It))[3] ) - std::abs ( peaks.at(grp.at(gr2It))[3] ) ) < (*dist) )
599  {
600  //==================================== Avoid very small angle distances as they would just take time (the hardcoded value would only be a problem for C700 and larger symmetries...
601  if ( std::abs( std::abs ( peaks.at(grp.at(gr1It))[3] ) - std::abs ( peaks.at(grp.at(gr2It))[3] ) ) > 0.01 )
602  {
603  g1 = gr1It;
604  g2 = gr2It;
605  *dist = std::abs( std::abs ( peaks.at(grp.at(gr1It))[3] ) - std::abs ( peaks.at(grp.at(gr2It))[3] ) );
606  }
607  }
608  }
609  }
610 
611  //================================================ If new dist found, save to tried and return success
612  const FloatingPoint< proshade_double > lhs1 ( *dist ), rhs1 ( 999.9 );
613  if ( !lhs1.AlmostEquals ( rhs1 ) )
614  {
615  ret = true;
616  ProSHADE_internal_misc::addToDoubleVector ( tried, static_cast< proshade_double > ( g1 ) );
617  ProSHADE_internal_misc::addToDoubleVector ( tried, static_cast< proshade_double > ( g2 ) );
619  }
620 
621  //================================================ Done
622  return ( ret );
623 
624 }
625 
635 void ProSHADE_internal_symmetry::addZeroPeakToGroups ( std::vector< std::vector< proshade_unsign > >& grpsVec, std::vector< proshade_double* >& peaks )
636 {
637  //================================================ Do your job
638  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( grpsVec.size() ); iter++ )
639  {
640  proshade_double* hlpP = new proshade_double [5];
641  ProSHADE_internal_misc::checkMemoryAllocation ( hlpP, __FILE__, __LINE__, __func__ );
642  hlpP[0] = peaks.at(grpsVec.at(iter).at(0))[0];
643  hlpP[1] = peaks.at(grpsVec.at(iter).at(0))[1];
644  hlpP[2] = peaks.at(grpsVec.at(iter).at(0))[2];
645  hlpP[3] = 0.0;
646  hlpP[4] = peaks.at(grpsVec.at(iter).at(0))[4];
647  ProSHADE_internal_misc::addToUnsignVector ( &grpsVec.at(iter), static_cast<proshade_unsign> ( peaks.size() ) );
649  }
650 
651  //================================================ Done
652  return ;
653 
654 }
655 
672 bool ProSHADE_internal_symmetry::determineFoldToTry ( proshade_double dist, proshade_double* divBasis, proshade_double* divRem, proshade_double peakErr, proshade_double* symmErr, std::vector< proshade_unsign >* angsToTry )
673 {
674  //================================================ Initialise variables
675  bool ret = false;
676 
677  //================================================ Find the basis and remainder of the 2pi/dist equation
678  *divRem = std::modf ( static_cast<proshade_double> ( ( 2.0 * M_PI ) / std::abs ( dist ) ), divBasis );
679 
680  //================================================ If the remainder would be smaller for larger basis, so this basis
681  if ( *divRem > 0.5 )
682  {
683  *divRem -= 1.0;
684  *divBasis += 1.0;
685  }
686 
687  //================================================ Determine errors on peaks and on folds
688  *symmErr = ( M_PI * 2.0 / *divBasis ) - ( M_PI * 2.0 / ( *divBasis + 1.0 ) );
689  proshade_double angTolerance = ( peakErr / *symmErr );
690 
691  //================================================ Is remainder small enough?
692  if ( ( *divRem < ( 0.0 + angTolerance ) ) && ( *divRem > ( 0.0 - angTolerance ) ) )
693  {
694  //============================================ Are we sure about the fold determination accuracy
695  proshade_signed angTolRound = std::min ( ProSHADE_internal_mapManip::myRound ( angTolerance ), static_cast<proshade_signed> ( 10 ) );
696  for ( proshade_signed iter = -angTolRound; iter <= angTolRound; iter++ )
697  {
698  ProSHADE_internal_misc::addToUnsignVector ( angsToTry, static_cast<proshade_unsign> ( std::max ( *divBasis + static_cast< proshade_double > ( iter ), 2.0 ) ) );
699  }
700  }
701 
702  //================================================ Return indication of whether testable fold value(s) was found.
703  if ( angsToTry->size() == 0 ) { ret = false; }
704  else { ret = true; }
705 
706  //================================================ Done
707  return ( ret );
708 
709 }
710 
719 void ProSHADE_internal_symmetry::findExpectedPeakRotations ( proshade_unsign fold, std::vector< proshade_double >* expAngs )
720 {
721  //================================================ Initialise variables
722  proshade_double groupAngle = ( 2.0 * M_PI ) / static_cast<proshade_double> ( fold );
723 
724  //================================================ Generate expected angles
725  for ( proshade_signed iter = static_cast<proshade_signed> ( -( static_cast<proshade_double> ( fold ) / 2.0 + 1.0) ); iter <= static_cast<proshade_signed> ( static_cast<proshade_double> ( fold )/2.0 + 1.0 ); iter++ )
726  {
727  ProSHADE_internal_misc::addToDoubleVector ( expAngs, static_cast< proshade_double > ( iter ) * groupAngle );
728  }
729 
730  //================================================ Done
731  return ;
732 
733 }
734 
748 proshade_unsign ProSHADE_internal_symmetry::checkExpectedAgainstFound ( std::vector< proshade_unsign > grp, std::vector< proshade_double* > peaks, std::vector< proshade_double >* expAngs, std::vector< proshade_unsign >* matchedAngs, std::vector< proshade_unsign >* missingAngs, proshade_double angTol )
749 {
750  //================================================ Initialise variables
751  proshade_unsign ret = 0;
752  proshade_unsign retHlp = 0;
753  proshade_double groupAngle = expAngs->at(1) - expAngs->at(0);
754  bool matchedThisPeak = false;
755  bool noDoubleMatches = false;
756  std::vector < proshade_unsign > matchedAlready;
757 
758  //================================================ For each expected peak rotation angle value
759  for ( proshade_unsign expAngIt = 0; expAngIt < static_cast<proshade_unsign> ( expAngs->size() ); expAngIt++ )
760  {
761  //============================================ For each peak in the group
762  matchedThisPeak = false;
763  for ( proshade_unsign peakIt = 0; peakIt < static_cast<proshade_unsign> ( grp.size() ); peakIt++ )
764  {
765  if ( ( expAngs->at(expAngIt) < ( peaks.at(grp.at(peakIt))[3] + angTol ) ) &&
766  ( expAngs->at(expAngIt) > ( peaks.at(grp.at(peakIt))[3] - angTol ) ) )
767  {
768  noDoubleMatches = false;
769  for ( proshade_unsign ndm = 0; ndm < static_cast<proshade_unsign> ( matchedAlready.size() ); ndm++ )
770  {
771  if ( matchedAlready.at(ndm) == grp.at(peakIt) ) { noDoubleMatches = true; break; }
772  }
773 
774  if ( !noDoubleMatches )
775  {
776  ProSHADE_internal_misc::addToUnsignVector ( matchedAngs, grp.at(peakIt) );
777  ProSHADE_internal_misc::addToUnsignVector ( &matchedAlready, grp.at(peakIt) );
778  matchedThisPeak = true;
779  break;
780  }
781  }
782  }
783 
784  //============================================ If not matched, add to missing
785  if ( !matchedThisPeak )
786  {
787  ProSHADE_internal_misc::addToUnsignVector ( missingAngs, expAngIt );
788  }
789  }
790 
791  //================================================ Find the number of consecutive matches
792  if ( matchedAngs->size () > 1 )
793  {
794  for ( proshade_unsign iter = 1; iter < static_cast<unsigned int> ( matchedAngs->size () ); iter++ )
795  {
796  if ( ( ( peaks.at(matchedAngs->at(iter-1))[3] + groupAngle ) < ( peaks.at(matchedAngs->at(iter))[3] + angTol ) ) &&
797  ( ( peaks.at(matchedAngs->at(iter-1))[3] + groupAngle ) > ( peaks.at(matchedAngs->at(iter))[3] - angTol ) ) )
798  {
799  retHlp += 1;
800  }
801  else
802  {
803  retHlp = 0;
804  }
805  if ( retHlp > ret ) { ret = retHlp; }
806  }
807  }
808 
809  //================================================ Done
810  return ( ret + 1 ); // This is because the count of matches is the count of intervals between numbers, so +1 to get the count of matched numbers.
811 
812 }
813 
830 proshade_double ProSHADE_internal_symmetry::checkForMissingPeak ( ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_double x, proshade_double y, proshade_double z, proshade_double angle, proshade_double heightThres, proshade_double axTol )
831 {
832  //================================================ Initialise variables
833  proshade_double ret = 0.0;
834  proshade_unsign arrIndex = 0;
835  proshade_double* rotMat = new proshade_double [9];
836  ProSHADE_internal_misc::checkMemoryAllocation ( rotMat, __FILE__, __LINE__, __func__ );
837  proshade_double pointHeight, euA, euB, euG, xPk, yPk, zPk, anglPk;
838  proshade_double angTol = std::acos ( 1.0 - axTol );
839 
840  //================================================ Search the self-rotation map
841  for ( proshade_unsign xIt = 0; xIt < ( dataObj->getMaxBand() * 2 ); xIt++ )
842  {
843  for ( proshade_unsign yIt = 0; yIt < ( dataObj->getMaxBand() * 2 ); yIt++ )
844  {
845  for ( proshade_unsign zIt = 0; zIt < ( dataObj->getMaxBand() * 2 ); zIt++ )
846  {
847  //==================================== Get height and check against threshold
848  arrIndex = zIt + ( dataObj->getMaxBand() * 2 ) * ( yIt + ( dataObj->getMaxBand() * 2 ) * xIt );
849  pointHeight = pow( dataObj->getInvSO3Coeffs()[arrIndex][0], 2.0 ) + pow( dataObj->getInvSO3Coeffs()[arrIndex][1], 2.0 );
850  if ( pointHeight < heightThres ) { continue; }
851 
852  //==================================== Get angle-axis values
853  ProSHADE_internal_maths::getEulerZXZFromSOFTPosition ( static_cast< proshade_signed > ( dataObj->getMaxBand() ), static_cast<proshade_signed> ( xIt ),
854  static_cast<proshade_signed> ( yIt ), static_cast<proshade_signed> ( zIt ),
855  &euA, &euB, &euG );
857  ProSHADE_internal_maths::getAxisAngleFromRotationMatrix ( rotMat, &xPk, &yPk, &zPk, &anglPk );
858 
859  //==================================== Check for matching angle
860  if ( ( ( std::abs( anglPk ) - angTol ) < std::abs ( angle ) ) && ( ( std::abs( anglPk ) + angTol ) > std::abs ( angle ) ) )
861  {
862  //================================ Make sure vector direction is the same
863  const FloatingPoint< proshade_double > lhs1 ( std::max( std::abs( xPk ), std::max( std::abs( yPk ), std::abs( zPk ) ) ) );
864  const FloatingPoint< proshade_double > rhs1 ( std::abs( xPk ) );
865  const FloatingPoint< proshade_double > rhs2 ( std::abs( yPk ) );
866  const FloatingPoint< proshade_double > rhs3 ( std::abs( zPk ) );
867  if ( ( lhs1.AlmostEquals ( rhs1 ) && ( xPk < 0.0 ) ) ||
868  ( lhs1.AlmostEquals ( rhs2 ) && ( yPk < 0.0 ) ) ||
869  ( lhs1.AlmostEquals ( rhs3 ) && ( zPk < 0.0 ) ) )
870  {
871  xPk *= -1.0;
872  yPk *= -1.0;
873  zPk *= -1.0;
874  anglPk *= -1.0;
875  }
876 
877  //================================ Compare axis elements
878  if ( ProSHADE_internal_maths::vectorOrientationSimilarity ( xPk, yPk, zPk, x, y, z, axTol ) )
879  {
880  if ( ret < pointHeight ) { ret = pointHeight; }
881  }
882  }
883  }
884  }
885  }
886 
887  //================================================ Done
888  return ( ret );
889 
890 }
891 
902 void ProSHADE_internal_symmetry::saveDetectedCSymmetry ( proshade_unsign fold, std::vector< proshade_unsign >* matchedPeaks, std::vector< std::vector< proshade_unsign > >* ret, proshade_signed verbose )
903 {
904  //================================================ Save fold as first vector value
905  std::vector< proshade_unsign > hlpVec;
907 
908  //================================================ and follow it with indices of all symmetry forming peaks
909  for ( proshade_unsign pIt = 0; pIt < static_cast<proshade_unsign> ( matchedPeaks->size() ); pIt++ )
910  {
911  ProSHADE_internal_misc::addToUnsignVector ( &hlpVec, matchedPeaks->at(pIt) );
912  }
914 
915  //================================================ Report finding symmetry
916  std::stringstream hlpS;
917  hlpS << "Found symmetry C" << fold;
918  ProSHADE_internal_messages::printProgressMessage ( verbose, 5, hlpS.str() );
919 
920  //================================================ Done
921  return ;
922 
923 }
924 
942 bool ProSHADE_internal_symmetry::completeMissingCSymmetry ( ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_unsign fold, std::vector< proshade_unsign >* grp, std::vector< proshade_double* >* peaks, std::vector< proshade_unsign >* missingPeaks, std::vector< proshade_double >* expectedAngles, std::vector< proshade_unsign >* matchedPeaks, proshade_double axErrTolerance, proshade_signed verbose )
943 {
944  //================================================ Initialise variables
945  bool ret = true;
946 
947  //================================================ Report searching for missing peaks
948  std::stringstream hlpSSP;
949  hlpSSP << "Searching for missing peaks for symmetry C" << fold;
950  ProSHADE_internal_messages::printProgressMessage ( verbose, 4, hlpSSP.str() );
951 
952  //================================================ Height threshold for missing peak
953  proshade_double heightThreshold = 0.0;
954  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( grp->size() ); grIt++ ) { heightThreshold += peaks->at(grp->at(grIt))[4]; }
955  heightThreshold /= static_cast<proshade_double> ( grp->size() );
956  heightThreshold *= 0.5;
957 
958  //================================================ For each missing value
959  for ( proshade_unsign misPkIt = 0; misPkIt < static_cast<proshade_unsign> ( missingPeaks->size() ); misPkIt++ )
960  {
961  //============================================ Ignore the extra values in the expected values
962  if ( expectedAngles->at(missingPeaks->at(misPkIt)) > M_PI ) { continue; }
963  if ( expectedAngles->at(missingPeaks->at(misPkIt)) < -M_PI ) { continue; }
964 
965  //============================================ Search for the missing peak
966  proshade_double misHeight = ProSHADE_internal_symmetry::checkForMissingPeak ( dataObj, peaks->at(grp->at(0))[0], peaks->at(grp->at(0))[1], peaks->at(grp->at(0))[2], expectedAngles->at(missingPeaks->at(misPkIt)), heightThreshold, axErrTolerance );
967  if ( misHeight != 0.0 )
968  {
969  //======================================== Missing peak detected - save it to the group
970  proshade_double* hlpP = new proshade_double [5];
971  ProSHADE_internal_misc::checkMemoryAllocation ( hlpP, __FILE__, __LINE__, __func__ );
972  hlpP[0] = peaks->at(grp->at(0))[0];
973  hlpP[1] = peaks->at(grp->at(0))[1];
974  hlpP[2] = peaks->at(grp->at(0))[2];
975  hlpP[3] = expectedAngles->at(missingPeaks->at(misPkIt));
976  hlpP[4] = misHeight;
977  ProSHADE_internal_misc::addToUnsignVector ( matchedPeaks, static_cast<proshade_unsign> ( peaks->size() ) );
979  }
980  else
981  {
982  ret = false;
983  }
984  }
985 
986  //================================================ Done
987  return ( ret );
988 
989 }
990 
1007 void ProSHADE_internal_symmetry::findSymmetryUsingFold ( ProSHADE_internal_data::ProSHADE_data* dataObj, std::vector< proshade_unsign >* angsToTry, std::vector< proshade_unsign >* grp, std::vector< proshade_double* >* peaks, std::vector< std::vector< proshade_unsign > >* ret, std::vector< proshade_unsign >* testedAlready, proshade_double axErrTolerance, bool axErrToleranceDefault, proshade_double missPeakThres, proshade_signed verbose )
1008 {
1009  //================================================ Initialise variables
1010  bool skipFold = false;
1011  std::vector< proshade_unsign > matchedPeaks, missingPeaks;
1012  std::vector< proshade_double > expectedAngles;
1013  proshade_double angTolerance = std::acos ( 1.0 - axErrTolerance );
1014 
1015  //================================================ Testing folds for being supported by peaks
1016  for ( proshade_unsign fIt = 0; fIt < static_cast<proshade_unsign> ( angsToTry->size() ); fIt++ )
1017  {
1018  //============================================ Was this fold already found?
1019  skipFold = false;
1020  for ( proshade_unsign ftIt = 0; ftIt < static_cast<proshade_unsign> ( testedAlready->size() ); ftIt++ ) { if ( testedAlready->at(ftIt) == angsToTry->at(fIt) ) { skipFold = true; } }
1021  if ( skipFold ) { continue; }
1022  else { ProSHADE_internal_misc::addToUnsignVector( testedAlready, angsToTry->at(fIt) ); }
1023 
1024  //============================================ Set axis tolerance based on fold (if required)
1025  if ( axErrToleranceDefault )
1026  {
1027  angTolerance = std::max ( std::min ( angTolerance, ( ( (M_PI * 2.0) / static_cast<double> ( angsToTry->at(fIt) ) ) -
1028  ( (M_PI * 2.0) / static_cast<double> ( angsToTry->at(fIt) + 1 ) ) ) * 2.0 ), 0.02 );
1029  axErrTolerance = std::max ( 1.0 - std::cos ( angTolerance ), 0.0008 );
1030  }
1031 
1032  //============================================ Find expected peak rotation angles
1033  expectedAngles.clear ( );
1034  ProSHADE_internal_symmetry::findExpectedPeakRotations ( angsToTry->at(fIt), &expectedAngles );
1035 
1036  //============================================ Compare group to expected angles
1037  matchedPeaks.clear ( );
1038  missingPeaks.clear ( );
1039  proshade_unsign consecMatches = ProSHADE_internal_symmetry::checkExpectedAgainstFound ( *grp, *peaks, &expectedAngles,
1040  &matchedPeaks, &missingPeaks, angTolerance );
1041 
1042  //============================================ If enough consecutive matches, symmetry was found. Save it
1043  if ( consecMatches >= angsToTry->at(fIt) )
1044  {
1045  ProSHADE_internal_symmetry::saveDetectedCSymmetry ( angsToTry->at(fIt), &matchedPeaks, ret, verbose );
1046  }
1047  else
1048  {
1049  if ( ( static_cast<proshade_double> ( matchedPeaks.size() ) / static_cast<proshade_double> ( angsToTry->at(fIt) ) ) >= ( 1.0 - missPeakThres ) )
1050  {
1051  //==================================== Attempt completing the symmetry using missing peaks
1052  if ( ProSHADE_internal_symmetry::completeMissingCSymmetry ( dataObj, angsToTry->at(fIt), grp, peaks, &missingPeaks,
1053  &expectedAngles, &matchedPeaks, axErrTolerance, verbose ) )
1054  {
1055  ProSHADE_internal_symmetry::saveDetectedCSymmetry ( angsToTry->at(fIt), &matchedPeaks, ret, verbose );
1056  }
1057  }
1058  else
1059  {
1060  //=================================== Symmetry not detected
1061  continue;
1062  }
1063  }
1064  }
1065 
1066  //=============================================== Done
1067  return ;
1068 
1069 }
1070 
1077 void ProSHADE_internal_symmetry::printSymmetryGroup ( std::vector< proshade_unsign > grp, std::vector< proshade_double* > peaks, proshade_signed verbose )
1078 {
1079  //================================================ Detected symmetry table header
1080  std::stringstream ss;
1081  ss << "Detected C" << grp.at(0) << " symmetry with following peaks:";
1082  ProSHADE_internal_messages::printProgressMessage ( verbose, 5, ss.str() );
1083  ProSHADE_internal_messages::printProgressMessage ( verbose, 5, "\tx\t y\t z\tAngle\tPeak height" );
1084 
1085  //================================================ Now print all supporting peaks
1086  for ( proshade_unsign pkIt = 1; pkIt < static_cast<proshade_unsign> ( grp.size() ); pkIt++ )
1087  {
1088  std::stringstream SS;
1089  SS << " " << static_cast<int>( peaks.at(grp.at(pkIt))[0] * 100.0 ) / 100.0 << "\t" << static_cast<int>( peaks.at(grp.at(pkIt))[1] * 100.0 ) / 100.0 << "\t" << static_cast<int>( peaks.at(grp.at(pkIt))[2] * 100.0 ) / 100.0 << "\t" << static_cast<int>( peaks.at(grp.at(pkIt))[3] * 100.0 ) / 100.0 << "\t" << static_cast<int>( peaks.at(grp.at(pkIt))[4] * 100.0 ) / 100.0;
1090  ProSHADE_internal_messages::printProgressMessage ( verbose, 5, SS.str() );
1091  }
1092 
1093  //================================================ Done
1094  return ;
1095 
1096 }
1097 
1103 void ProSHADE_internal_symmetry::printSymmetryCompletion ( proshade_unsign noSyms, proshade_signed verbose )
1104 {
1105  //================================================ Report completion of symmetry detection
1106  std::stringstream ss;
1107  ss << "Detected " << noSyms << " Cyclic symmetries.";
1108  ProSHADE_internal_messages::printProgressMessage ( verbose, 2, ss.str() );
1109 
1110  //================================================ If no symmetries were found, print warning
1111  if ( noSyms < 1 )
1112  {
1113  ProSHADE_internal_messages::printWarningMessage ( verbose, "!!! ProSHADE WARNING !!! Failed to detect any symmetries. If you believe there should be one, you can try decreasing the resolution or checking that the map is centred on the centry of symmetry (or use map centering option in ProSHADE).", "WS00030" );
1114  }
1115 
1116  //================================================ Done
1117  return ;
1118 
1119 }
1120 
1133 void ProSHADE_internal_symmetry::saveAllCSymmetries ( std::vector< std::vector< proshade_unsign > > detected, std::vector< proshade_double* > peaks, std::vector< proshade_double* >* ret, proshade_double axErr )
1134 {
1135  //================================================ Initialise variables
1136  proshade_double sumX, sumY, sumZ, sumH;
1137  proshade_signed matchedPos = -1;
1138 
1139  //================================================ Start saving
1140  for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( detected.size() ); symIt++ )
1141  {
1142  //============================================ Allocate the memory
1143  proshade_double* hlpP = new proshade_double [6];
1144  ProSHADE_internal_misc::checkMemoryAllocation ( hlpP, __FILE__, __LINE__, __func__ );
1145 
1146  //============================================ Set obvious values
1147  hlpP[0] = static_cast<proshade_double> ( detected.at(symIt).at(0) );
1148  hlpP[4] = static_cast<proshade_double> ( ( 2.0 * M_PI ) / hlpP[0] );
1149 
1150  //============================================ Compute peak averages for rest
1151  sumX = 0.0; sumY = 0.0; sumZ = 0.0; sumH = 0.0;
1152  for ( proshade_unsign pkIt = 1; pkIt < static_cast<proshade_unsign> ( detected.at(symIt).size() ); pkIt++ )
1153  {
1154  sumX += peaks.at(detected.at(symIt).at(pkIt))[0];
1155  sumY += peaks.at(detected.at(symIt).at(pkIt))[1];
1156  sumZ += peaks.at(detected.at(symIt).at(pkIt))[2];
1157  sumH += peaks.at(detected.at(symIt).at(pkIt))[4];
1158  }
1159  sumX /= static_cast<proshade_double> ( detected.at(symIt).size() - 1 );
1160  sumY /= static_cast<proshade_double> ( detected.at(symIt).size() - 1 );
1161  sumZ /= static_cast<proshade_double> ( detected.at(symIt).size() - 1 );
1162  sumH /= static_cast<proshade_double> ( detected.at(symIt).size() - 1 );
1163 
1164  //============================================ And add these as well
1165  hlpP[1] = sumX;
1166  hlpP[2] = sumY;
1167  hlpP[3] = sumZ;
1168  hlpP[5] = sumH;
1169 
1170  //============================================ Save the complete symmetry description to the vector, unless already there
1171  if ( !ProSHADE_internal_symmetry::isSymmetrySame ( ret, hlpP, axErr, &matchedPos ) )
1172  {
1174  }
1175  else
1176  {
1177  delete[] hlpP;
1178  }
1179  }
1180 
1181  //================================================ Done
1182  return ;
1183 
1184 }
1185 
1197 bool ProSHADE_internal_symmetry::isSymmetrySame ( std::vector< proshade_double* >* ret, proshade_double* sym, proshade_double simThres, proshade_signed* matchedPos )
1198 {
1199  //================================================ Initialise variables
1200  proshade_double dotProduct = 0.0;
1201  *matchedPos = -1;
1202 
1203  //================================================ Check
1204  for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( ret->size() ); symIt++ )
1205  {
1206  //============================================ Minor speed-up => only test for same folds
1207  const FloatingPoint< proshade_double > lhs ( ret->at(symIt)[0] ), rhs ( sym[0] );
1208  if ( lhs.AlmostEquals ( rhs ) )
1209  {
1210  //======================================== Is axis the same?
1211  dotProduct = ProSHADE_internal_maths::computeDotProduct ( &ret->at(symIt)[1], &ret->at(symIt)[2],
1212  &ret->at(symIt)[3], &sym[1], &sym[2], &sym[3] );
1213  if ( ( ( 1.0 > ( dotProduct - simThres ) ) && ( 1.0 < ( dotProduct + simThres ) ) ) || ( ( -1.0 > ( dotProduct - simThres ) ) && ( -1.0 < ( dotProduct + simThres ) ) ) )
1214  {
1215  //==================================== Matched. Save the index
1216  *matchedPos = static_cast< proshade_signed > ( symIt );
1217 
1218  //==================================== Does the already saved have higher height?
1219  if ( ret->at(symIt)[5] >= sym[5] ) { return ( true ); }
1220 
1221  //==================================== In this case, new is better than old - sort it out
1222  ret->at(symIt)[1] = sym[1];
1223  ret->at(symIt)[2] = sym[2];
1224  ret->at(symIt)[3] = sym[3];
1225  ret->at(symIt)[5] = sym[5];
1226  return ( true );
1227  }
1228  }
1229  }
1230 
1231  //================================================ Done - no matches found
1232  return ( false );
1233 
1234 }
1235 
1248 bool ProSHADE_internal_symmetry::isSymmetrySame ( std::vector< proshade_double* >* ret, proshade_double* sym, proshade_double simThres, proshade_signed* matchedPos, proshade_double fscVal )
1249 {
1250  //================================================ Initialise variables
1251  proshade_double dotProduct = 0.0;
1252  *matchedPos = -1;
1253 
1254  //================================================ Check
1255  for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( ret->size() ); symIt++ )
1256  {
1257  //============================================ Minor speed-up => only test for same folds
1258  const FloatingPoint< proshade_double > lhs ( ret->at(symIt)[0] ), rhs ( sym[0] );
1259  if ( lhs.AlmostEquals ( rhs ) )
1260  {
1261  //======================================== Is axis the same?
1262  dotProduct = ProSHADE_internal_maths::computeDotProduct ( &ret->at(symIt)[1], &ret->at(symIt)[2],
1263  &ret->at(symIt)[3], &sym[1], &sym[2], &sym[3] );
1264  if ( ( ( 1.0 > ( dotProduct - simThres ) ) && ( 1.0 < ( dotProduct + simThres ) ) ) || ( ( -1.0 > ( dotProduct - simThres ) ) && ( -1.0 < ( dotProduct + simThres ) ) ) )
1265  {
1266  //==================================== Matched. Save the index
1267  *matchedPos = static_cast< proshade_signed > ( symIt );
1268 
1269  //==================================== Does the already saved have higher height?
1270  if ( ret->at(symIt)[5] >= sym[5] ) { return ( true ); }
1271 
1272  //==================================== In this case, new is better than old - sort it out
1273  ret->at(symIt)[1] = sym[1];
1274  ret->at(symIt)[2] = sym[2];
1275  ret->at(symIt)[3] = sym[3];
1276  ret->at(symIt)[5] = sym[5];
1277  ret->at(symIt)[6] = fscVal;
1278  return ( true );
1279  }
1280  }
1281  }
1282 
1283  //================================================ Done - no matches found
1284  return ( false );
1285 
1286 }
1287 
1299 std::vector< proshade_double* > ProSHADE_internal_data::ProSHADE_data::getDihedralSymmetriesList ( ProSHADE_settings* settings, std::vector< proshade_double* >* CSymList )
1300 {
1301  //================================================ Initialise variables
1302  std::vector< proshade_double* > ret;
1303  proshade_double dotProduct;
1304 
1305  //================================================ Report progress
1306  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting D symmetry detection." );
1307 
1308  //================================================If not enough axes, just end here
1309  if ( CSymList->size() < 2 ) { return ( ret ); }
1310 
1311  //================================================ For each unique pair of axes
1312  for ( proshade_unsign ax1 = 0; ax1 < static_cast<proshade_unsign> ( CSymList->size() ); ax1++ )
1313  {
1314  //============================================ Ignore small axes
1315  const FloatingPoint< proshade_double > lhs1 ( CSymList->at(ax1)[5] ), rhs1 ( -999.9 );
1316  if ( ( CSymList->at(ax1)[5] < settings->minSymPeak ) && !( lhs1.AlmostEquals ( rhs1 ) ) ) { continue; }
1317 
1318  for ( proshade_unsign ax2 = 1; ax2 < static_cast<proshade_unsign> ( CSymList->size() ); ax2++ )
1319  {
1320  //======================================= Use unique pairs only
1321  if ( ax1 >= ax2 ) { continue; }
1322 
1323  //======================================== Ignore small axes
1324  const FloatingPoint< proshade_double > lhs2 ( CSymList->at(ax2)[5] ), rhs2 ( -999.9 );
1325  if ( ( CSymList->at(ax2)[5] < settings->minSymPeak ) && !( lhs2.AlmostEquals ( rhs2 ) ) ) { continue; }
1326 
1327  //======================================= Compute the dot product
1328  dotProduct = ProSHADE_internal_maths::computeDotProduct ( &CSymList->at(ax1)[1], &CSymList->at(ax1)[2],
1329  &CSymList->at(ax1)[3], &CSymList->at(ax2)[1],
1330  &CSymList->at(ax2)[2], &CSymList->at(ax2)[3] );
1331 
1332  //======================================== If close to zero, these two axes are perpendicular
1333  if ( std::abs( dotProduct ) < settings->axisErrTolerance )
1334  {
1335  //==================================== Save
1336  if ( CSymList->at(ax1)[0] >= CSymList->at(ax2)[0] )
1337  {
1338  ProSHADE_internal_symmetry::saveDSymmetry ( &ret, CSymList, ax1, ax2 );
1339 
1340  std::vector< proshade_unsign > DSymInd;
1344 
1345  }
1346  else
1347  {
1348  ProSHADE_internal_symmetry::saveDSymmetry ( &ret, CSymList, ax2, ax1 );
1349 
1350  std::vector< proshade_unsign > DSymInd;
1354  }
1355  }
1356  }
1357  }
1358 
1359  //================================================ Report progress
1360  std::stringstream hlpSS;
1361  hlpSS << "Detected " << ret.size() << " D symmetries.";
1362  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, hlpSS.str() );
1363 
1364  //================================================ Done
1365  return ( ret );
1366 
1367 }
1368 
1380 void ProSHADE_internal_symmetry::saveDSymmetry ( std::vector< proshade_double* >* ret, std::vector< proshade_double* >* CSymList, proshade_unsign axisOne, proshade_unsign axisTwo )
1381 {
1382  //================================================ Allocate the memory
1383  proshade_double* hlpP = new proshade_double [14];
1384  ProSHADE_internal_misc::checkMemoryAllocation ( hlpP, __FILE__, __LINE__, __func__ );
1385 
1386  //================================================ Set the axis and heights
1387  hlpP[0] = CSymList->at(axisOne)[0];
1388  hlpP[1] = CSymList->at(axisOne)[1];
1389  hlpP[2] = CSymList->at(axisOne)[2];
1390  hlpP[3] = CSymList->at(axisOne)[3];
1391  hlpP[4] = CSymList->at(axisOne)[4];
1392  hlpP[5] = CSymList->at(axisOne)[5];
1393  hlpP[6] = CSymList->at(axisOne)[6];
1394  hlpP[7] = CSymList->at(axisTwo)[0];
1395  hlpP[8] = CSymList->at(axisTwo)[1];
1396  hlpP[9] = CSymList->at(axisTwo)[2];
1397  hlpP[10] = CSymList->at(axisTwo)[3];
1398  hlpP[11] = CSymList->at(axisTwo)[4];
1399  hlpP[12] = CSymList->at(axisTwo)[5];
1400  hlpP[13] = CSymList->at(axisTwo)[6];
1401 
1402  //================================================ Save to ret
1404 
1405  //================================================ Done
1406  return ;
1407 
1408 }
1409 
1419 std::vector< proshade_double* > ProSHADE_internal_data::ProSHADE_data::getTetrahedralSymmetriesList ( ProSHADE_settings* settings, std::vector< proshade_double* >* CSymList )
1420 {
1421  //================================================ Initialise variables
1422  std::vector< proshade_double* > ret;
1423 
1424  //================================================ Report progress
1425  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting T symmetry detection." );
1426 
1427  //================================================ Are the basic requirements for tetrahedral symmetry met?
1429  {
1430  //============================================ Search for all the symmetry axes
1431  ProSHADE_internal_symmetry::findTetra4C3s ( CSymList, &ret, settings->axisErrTolerance, this, settings->verbose, settings->minSymPeak );
1432  if ( ret.size() != 4 ) { ProSHADE_internal_messages::printWarningMessage ( settings->verbose, "!!! ProSHADE WARNING !!! Failed to detect some of the polyhedral symmetries, while detecting the correct dihedral angles.", "WS00031" ); return ( ret ); }
1433 
1434  ProSHADE_internal_symmetry::findTetra3C2s ( CSymList, &ret, settings->axisErrTolerance, this, settings->verbose, settings->minSymPeak );
1435  if ( ret.size() != 7 ) { ProSHADE_internal_messages::printWarningMessage ( settings->verbose, "!!! ProSHADE WARNING !!! Failed to detect some of the polyhedral symmetries, while detecting the correct dihedral angles.", "WS00031" ); return ( ret ); }
1436  else
1437  {
1438  for ( proshade_unsign csIt = 0; csIt < static_cast<proshade_unsign> ( CSymList->size() ); csIt++ )
1439  {
1440  for ( proshade_unsign retIt = 0; retIt < static_cast<proshade_unsign> ( ret.size() ); retIt++ )
1441  {
1442  //======================================== Sort ret by fold
1443  std::sort ( ret.begin(), ret.end(), ProSHADE_internal_misc::sortSymInvFoldHlp );
1444 
1445  //======================================== Save indices
1446  const FloatingPoint< proshade_double > lhs1 ( CSymList->at(csIt)[0] ), rhs1 ( ret.at(retIt)[0] );
1447  const FloatingPoint< proshade_double > lhs2 ( CSymList->at(csIt)[1] ), rhs2 ( ret.at(retIt)[1] );
1448  const FloatingPoint< proshade_double > lhs3 ( CSymList->at(csIt)[2] ), rhs3 ( ret.at(retIt)[2] );
1449  const FloatingPoint< proshade_double > lhs4 ( CSymList->at(csIt)[3] ), rhs4 ( ret.at(retIt)[3] );
1450  const FloatingPoint< proshade_double > lhs5 ( CSymList->at(csIt)[4] ), rhs5 ( ret.at(retIt)[4] );
1451  const FloatingPoint< proshade_double > lhs6 ( CSymList->at(csIt)[5] ), rhs6 ( ret.at(retIt)[5] );
1452  if ( ( lhs1.AlmostEquals ( rhs1 ) ) &&
1453  ( lhs2.AlmostEquals ( rhs2 ) ) &&
1454  ( lhs3.AlmostEquals ( rhs3 ) ) &&
1455  ( lhs4.AlmostEquals ( rhs4 ) ) &&
1456  ( lhs5.AlmostEquals ( rhs5 ) ) &&
1457  ( lhs6.AlmostEquals ( rhs6 ) ) )
1458  {
1460  }
1461  }
1462  }
1463  }
1464  }
1465 
1466  //================================================ Report progress
1467  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "T symmetry detection complete." );
1468 
1469  //================================================ Done
1470  return ( ret );
1471 
1472 }
1473 
1485 bool ProSHADE_internal_symmetry::detectTetrahedralSymmetry ( std::vector< proshade_double* >* CSymList, proshade_double axErr, proshade_double minPeakHeight )
1486 {
1487  //================================================ Initialise variables
1488  std::vector< proshade_unsign > C3List;
1489  proshade_double dotProduct;
1490 
1491  //================================================ Find all C3 symmetries
1492  for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
1493  {
1494  const FloatingPoint< proshade_double > lhs ( CSymList->at(cSym)[0] ), rhs ( 3.0 );
1495  if ( lhs.AlmostEquals ( rhs ) && CSymList->at(cSym)[5] >= minPeakHeight ) { ProSHADE_internal_misc::addToUnsignVector ( &C3List, cSym ); }
1496  }
1497 
1498  //================================================ For each unique pair of C3s
1499  for ( proshade_unsign c31 = 0; c31 < static_cast<proshade_unsign> ( C3List.size() ); c31++ )
1500  {
1501  for ( proshade_unsign c32 = 1; c32 < static_cast<proshade_unsign> ( C3List.size() ); c32++ )
1502  {
1503  //================================ Unique pairs only
1504  if ( c31 >= c32 ) { continue; }
1505 
1506  //======================================== Check the angle between the C3 axes
1507  dotProduct = ProSHADE_internal_maths::computeDotProduct ( &CSymList->at(C3List.at(c31))[1], &CSymList->at(C3List.at(c31))[2], &CSymList->at(C3List.at(c31))[3], &CSymList->at(C3List.at(c32))[1], &CSymList->at(C3List.at(c32))[2], &CSymList->at(C3List.at(c32))[3] );
1508 
1509  //================================ Is the angle approximately the dihedral angle
1510  if ( ( ( 1.0 / 3.0 ) > ( dotProduct - axErr ) ) && ( ( 1.0 / 3.0 ) < ( dotProduct + axErr ) ) )
1511  {
1512  return ( true );
1513  }
1514  }
1515  }
1516 
1517  //================================================ Done
1518  return ( false );
1519 
1520 }
1521 
1535 void ProSHADE_internal_symmetry::findTetra4C3s ( std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_signed verbose, proshade_double minPeakHeight )
1536 {
1537  //================================================ Initialise variables
1538  std::vector< proshade_unsign > C3PossibilitiesHlp;
1539  std::vector< std::vector< proshade_unsign > > C3Possibilities;
1540  bool groupMatched;
1541 
1542  //================================================ Report progress
1543  ProSHADE_internal_messages::printProgressMessage ( verbose, 2, "Starting detection of four C3 axes." );
1544 
1545  //================================================ For all symmetries in the C symmetries list
1546  for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
1547  {
1548  //============================================ Search only using C3s
1549  if ( CSymList->at(cIt)[0] != 3.0 || CSymList->at(cIt)[0] < minPeakHeight ) { continue; }
1550 
1551  //============================================ If this is the first C3, then just save it to the first group of the temporary holder
1552  if ( C3Possibilities.size() == 0 ) { ProSHADE_internal_misc::addToUnsignVector ( &C3PossibilitiesHlp, cIt ); ProSHADE_internal_misc::addToUnsignVectorVector ( &C3Possibilities, C3PossibilitiesHlp ); continue; }
1553 
1554  //============================================ If second or more C3, check if it has the correct angle to all other already found C3s for each group
1555  groupMatched = false;
1556  for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( C3Possibilities.size() ); gIt++ )
1557  {
1558  if ( ProSHADE_internal_symmetry::testGroupAgainstSymmetry ( CSymList, &C3Possibilities.at(gIt), CSymList->at(cIt), axErr, 1.0/3.0, true, cIt ) ) { ProSHADE_internal_misc::addToUnsignVector ( &C3Possibilities.at(gIt), cIt ); groupMatched = true; break; }
1559  }
1560 
1561  //============================================ If no group matched, create a new group
1562  if ( !groupMatched ) { C3PossibilitiesHlp.clear(); ProSHADE_internal_misc::addToUnsignVector ( &C3PossibilitiesHlp, cIt ); ProSHADE_internal_misc::addToUnsignVectorVector ( &C3Possibilities, C3PossibilitiesHlp ); continue; }
1563  }
1564 
1565  //================================================ Test for missing symmetry axes, if need be
1566  ProSHADE_internal_symmetry::findMissingAxes ( &C3Possibilities, CSymList, 4, axErr, 1.0/3.0, 3, dataObj, minPeakHeight );
1567 
1568  //================================================ Any group has 4 entries? If more such groups, take the one with highest average height.
1569  proshade_double maxHeight = 0.0; proshade_unsign maxGrp = 0;
1570  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( C3Possibilities.size() ); iter++ ) { if ( C3Possibilities.at(iter).size() == 4 ) { if ( ( ( CSymList->at(C3Possibilities.at(iter).at(0))[5] + CSymList->at(C3Possibilities.at(iter).at(1))[5] + CSymList->at(C3Possibilities.at(iter).at(2))[5] + CSymList->at(C3Possibilities.at(iter).at(3))[5] ) / 4.0 ) > maxHeight ) { maxHeight = ( ( CSymList->at(C3Possibilities.at(iter).at(0))[5] + CSymList->at(C3Possibilities.at(iter).at(1))[5] + CSymList->at(C3Possibilities.at(iter).at(2))[5] + CSymList->at(C3Possibilities.at(iter).at(3))[5] ) / 4.0 ); maxGrp = iter; } } }
1571 
1572  if ( C3Possibilities.at(maxGrp).size() == 4 )
1573  {
1574  //============================================ Success! Save and exit
1575  for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( C3Possibilities.at(maxGrp).size() ); it++ ) { ProSHADE_internal_misc::addToDblPtrVector ( ret, CSymList->at(C3Possibilities.at(maxGrp).at(it)) ); }
1576 
1577  //============================================ Report progress
1578  ProSHADE_internal_messages::printProgressMessage ( verbose, 3, "Detection of four C3 axes successfull." );
1579 
1580  //============================================ Done
1581  return ;
1582  }
1583 
1584  //================================================ Done
1585  return ;
1586 
1587 }
1588 
1605 bool ProSHADE_internal_symmetry::testGroupAgainstSymmetry ( std::vector< proshade_double* >* CSymList, std::vector< proshade_unsign >* grp, proshade_double* sym, proshade_double axErr, proshade_double angle, bool improve, proshade_unsign pos )
1606 {
1607  //================================================ Initialise variables
1608  bool allAnglesMet = true;
1609  proshade_double dotProduct;
1610 
1611  //================================================ Improve if required
1612  if ( improve )
1613  {
1614  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( grp->size() ); mIt++ )
1615  {
1616  dotProduct = ProSHADE_internal_maths::computeDotProduct ( &CSymList->at(grp->at(mIt))[1], &CSymList->at(grp->at(mIt))[2], &CSymList->at(grp->at(mIt))[3], &sym[1], &sym[2], &sym[3] );
1617 
1618  if ( ( ( 1.0 > ( dotProduct - axErr ) ) && ( 1.0 < ( dotProduct + axErr ) ) ) || ( ( -1.0 > ( dotProduct - axErr ) ) && ( -1.0 < ( dotProduct + axErr ) ) ) )
1619  {
1620  if ( sym[5] > CSymList->at(grp->at(mIt))[5] )
1621  {
1622  grp->at(mIt) = pos;
1623  }
1624  else
1625  {
1626  allAnglesMet = false;
1627  return ( allAnglesMet );
1628  }
1629  }
1630  }
1631  }
1632 
1633  //================================================ For all group members
1634  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( grp->size() ); mIt++ )
1635  {
1636  dotProduct = ProSHADE_internal_maths::computeDotProduct ( &CSymList->at(grp->at(mIt))[1], &CSymList->at(grp->at(mIt))[2], &CSymList->at(grp->at(mIt))[3], &sym[1], &sym[2], &sym[3] );
1637 
1638  if ( ( angle > ( std::abs ( dotProduct ) - axErr ) ) &&
1639  ( angle < ( std::abs ( dotProduct ) + axErr ) ) )
1640  {
1641  //======================================== Matching group memner - try next one
1642  }
1643  else
1644  {
1645  //======================================== Group member not matched - try next group
1646  allAnglesMet = false;
1647  break;
1648  }
1649  }
1650 
1651  //================================================ Done
1652  return ( allAnglesMet );
1653 
1654 }
1655 
1672 bool ProSHADE_internal_symmetry::findMissingAxes ( std::vector< std::vector< proshade_unsign > >* possibilities, std::vector< proshade_double* >* CSymList, proshade_unsign requiredNoAxes, proshade_double axErr, proshade_double angle, proshade_unsign fold, ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_double minPeakHeight )
1673 {
1674  //================================================ Initialise variables
1675  std::vector< proshade_double* > hlpVec;
1676  bool atLeastOne = false;
1677 
1678  //================================================ Proceed only if need be
1679  for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( possibilities->size() ); gIt++ )
1680  {
1681  if ( static_cast<proshade_unsign> ( possibilities->at(gIt).size() ) == requiredNoAxes ) { atLeastOne = true; return ( atLeastOne ); }
1682  }
1683 
1684  //================================================ For each possible group
1685  for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( possibilities->size() ); gIt++ )
1686  {
1687  //============================================ This will not work for less than two axes in group
1688  if ( possibilities->at(gIt).size() < 2 ) { continue; }
1689 
1690  //============================================ Prepare iteration
1691  hlpVec.clear ( );
1692 
1693  //============================================ Search for missing axes
1694  ProSHADE_internal_symmetry::searchMissingSymmetrySpace ( dataObj, CSymList, &possibilities->at(gIt), &hlpVec, axErr, angle, fold, minPeakHeight );
1695 
1696  //============================================ Add missing axes
1697  if ( hlpVec.size() > 0 )
1698  {
1699  //======================================== Start adding by highest first
1700  std::sort ( hlpVec.begin(), hlpVec.end(), ProSHADE_internal_misc::sortSymHlpInv );
1701 
1702  //======================================== For each missing axis
1703  for ( proshade_unsign axIt = 0; axIt < static_cast<proshade_unsign> ( hlpVec.size() ); axIt++ )
1704  {
1705  if ( ProSHADE_internal_symmetry::testGroupAgainstSymmetry ( CSymList, &possibilities->at(gIt), hlpVec.at(axIt), axErr, angle, false ) )
1706  {
1707  //================================ Check for uniqueness
1708  if ( ProSHADE_internal_maths::isAxisUnique ( CSymList, hlpVec.at(axIt), axErr ) )
1709  {
1710  //============================ Add
1711  ProSHADE_internal_misc::addToDblPtrVector ( CSymList, hlpVec.at(axIt) );
1712  ProSHADE_internal_misc::addToUnsignVector ( &possibilities->at(gIt), static_cast<proshade_unsign> ( CSymList->size()-1 ) );
1713  }
1714  }
1715  }
1716  }
1717 
1718  if ( possibilities->at(gIt).size() == requiredNoAxes ) { atLeastOne = true; }
1719  }
1720 
1721  //================================================ Done
1722  return ( atLeastOne );
1723 
1724 }
1725 
1732 bool ProSHADE_internal_symmetry::sortArrVecHlp ( const proshade_double* a, const proshade_double* b )
1733 {
1734  //================================================ Compare
1735  return ( a[0] < b[0] );
1736 
1737 }
1738 
1753 proshade_double ProSHADE_internal_symmetry::missingAxisHeight ( proshade_double xVal, proshade_double yVal, proshade_double zVal, ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_unsign fold, proshade_double axErr )
1754 {
1755  //================================================ Initialise variables
1756  proshade_double ret = 0.0;
1757  proshade_double curSum = 0.0;
1758  proshade_double maxVal = 0.0;
1759  proshade_double angStep = std::acos ( 1.0 - axErr ) / 2;
1760  std::vector< proshade_double* > angVec;
1761 
1762  //================================================ Find map points conforming to the axis
1763  angVec = ProSHADE_internal_symmetry::findMissingAxisPoints ( xVal, yVal, zVal, dataObj, axErr );
1764 
1765  //================================================ Sort points by angle
1766  std::sort ( angVec.begin(), angVec.end(), ProSHADE_internal_symmetry::sortArrVecHlp );
1767 
1768  //================================================ Find the best X peaks with correct distances
1769  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( std::floor ( ( 2.0 * M_PI / angStep ) / static_cast< proshade_double > ( fold ) ) ); iter++ )
1770  {
1771  //============================================ Initialise new ang group iteration
1772  curSum = 0.0;
1773 
1774  //============================================ For each of the fold times
1775  for ( proshade_unsign angCmb = 0; angCmb < static_cast<proshade_unsign> ( fold ); angCmb++ )
1776  {
1777  //======================================== Initialise
1778  maxVal = 0.0;
1779 
1780  //======================================== Search
1781  for ( proshade_unsign angIt = 0; angIt < static_cast<proshade_unsign> ( angVec.size() ); angIt++ )
1782  {
1783  if ( angVec.at(angIt)[0] < ( ( static_cast< proshade_double > ( iter ) * angStep ) +
1784  ( ( 2.0 * M_PI / static_cast< proshade_double > ( fold ) ) * static_cast< proshade_double > ( angCmb ) ) ) ) { continue; }
1785  if ( angVec.at(angIt)[0] > ( ( ( static_cast< proshade_double > ( iter ) + 1.0 ) * angStep ) +
1786  ( ( 2.0 * M_PI / static_cast< proshade_double > ( fold ) ) * static_cast< proshade_double > ( angCmb ) ) ) ) { break; }
1787 
1788  if ( angVec.at(angIt)[1] > maxVal ) { maxVal = angVec.at(angIt)[1]; }
1789  }
1790  curSum += maxVal;
1791  }
1792  curSum /= static_cast<proshade_double> ( fold );
1793  if ( ret < curSum ) { ret = curSum; }
1794  }
1795 
1796  //================================================ Release memory
1797  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( angVec.size() ); iter++ ) { delete[] angVec.at(iter); }
1798 
1799  //================================================ Done
1800  return ( ret );
1801 
1802 }
1803 
1817 std::vector < proshade_double* > ProSHADE_internal_symmetry::findMissingAxisPoints ( proshade_double xVal, proshade_double yVal, proshade_double zVal, ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_double axErr )
1818 {
1819  //================================================ Initialise variables
1820  proshade_double euA, euB, euG, xPk, yPk, zPk, anglPk;
1821  proshade_double* rotMat = new proshade_double [9];
1822  ProSHADE_internal_misc::checkMemoryAllocation ( rotMat, __FILE__, __LINE__, __func__ );
1823  proshade_unsign arrIndex;
1824  std::vector< proshade_double* > angVec;
1825 
1826  //================================================ Search the self-rotation map
1827  for ( proshade_unsign xIt = 0; xIt < ( dataObj->getMaxBand() * 2 ); xIt++ )
1828  {
1829  for ( proshade_unsign yIt = 0; yIt < ( dataObj->getMaxBand() * 2 ); yIt++ )
1830  {
1831  for ( proshade_unsign zIt = 0; zIt < ( dataObj->getMaxBand() * 2 ); zIt++ )
1832  {
1833  //==================================== Get height and check against threshold
1834  arrIndex = zIt + ( dataObj->getMaxBand() * 2 ) * ( yIt + ( dataObj->getMaxBand() * 2 ) * xIt );
1835 
1836  //==================================== Get angle-axis values
1837  ProSHADE_internal_maths::getEulerZXZFromSOFTPosition ( static_cast< proshade_signed > ( dataObj->getMaxBand() ), static_cast< proshade_signed > ( xIt ),
1838  static_cast< proshade_signed > ( yIt ), static_cast< proshade_signed > ( zIt ),
1839  &euA, &euB, &euG );
1841  ProSHADE_internal_maths::getAxisAngleFromRotationMatrix ( rotMat, &xPk, &yPk, &zPk, &anglPk );
1842 
1843  //==================================== Set largest axis element to positive
1844  const FloatingPoint< proshade_double > lhs1 ( std::max ( std::abs ( xPk ), std::max( std::abs ( yPk ), std::abs ( zPk ) ) ) );
1845  const FloatingPoint< proshade_double > rhs1 ( std::abs ( xPk ));
1846  const FloatingPoint< proshade_double > rhs2 ( std::abs ( yPk ) );
1847  const FloatingPoint< proshade_double > rhs3 ( std::abs ( zPk ) );
1848  if ( ( lhs1.AlmostEquals ( rhs1 ) && ( xPk < 0.0 ) ) ||
1849  ( lhs1.AlmostEquals ( rhs2 ) && ( yPk < 0.0 ) ) ||
1850  ( lhs1.AlmostEquals ( rhs3 ) && ( zPk < 0.0 ) ) )
1851  {
1852  xPk *= -1.0;
1853  yPk *= -1.0;
1854  zPk *= -1.0;
1855  anglPk *= -1.0;
1856  }
1857 
1858  //==================================== Does the peak match the required axis?
1859  if ( ProSHADE_internal_maths::vectorOrientationSimilarity ( xPk, yPk, zPk, xVal, yVal, zVal, axErr ) )
1860  {
1861  //================================ Matching map point - save it
1862  proshade_double* hlpArr = new proshade_double [2];
1863  ProSHADE_internal_misc::checkMemoryAllocation ( hlpArr, __FILE__, __LINE__, __func__ );
1864  hlpArr[0] = anglPk + M_PI;
1865  hlpArr[1] = pow( dataObj->getInvSO3Coeffs()[arrIndex][0], 2.0 ) +
1866  pow( dataObj->getInvSO3Coeffs()[arrIndex][1], 2.0 );
1867  ProSHADE_internal_misc::addToDblPtrVector ( &angVec, hlpArr );
1868  }
1869  }
1870  }
1871  }
1872 
1873  //================================================ Release memory
1874  delete[] rotMat;
1875 
1876  //================================================ Done
1877  return ( angVec );
1878 
1879 }
1880 
1895 void ProSHADE_internal_symmetry::saveMissingAxisNewOnly ( std::vector< proshade_double* >* axVec, proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double height, proshade_unsign fold, proshade_double axErr )
1896 {
1897  //================================================ Create symmetry array from the inputs
1898  proshade_double* hlpSym = new proshade_double [6];
1899  ProSHADE_internal_misc::checkMemoryAllocation ( hlpSym, __FILE__, __LINE__, __func__ );
1900 
1901  //================================================ Fill it in
1902  hlpSym[0] = static_cast<proshade_double> ( fold );
1903  hlpSym[1] = axX;
1904  hlpSym[2] = axY;
1905  hlpSym[3] = axZ;
1906  hlpSym[4] = ( 2.0 * M_PI ) / static_cast<proshade_double> ( fold );
1907  hlpSym[5] = height;
1908 
1909  //================================================ Check if similar symmetry does not exist already
1910  for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( axVec->size() ); symIt++ )
1911  {
1912  //============================================ Minor speed-up => only test for same folds
1913  const FloatingPoint< proshade_double > lhs1 ( axVec->at(symIt)[0] ), rhs1 ( hlpSym[0] );
1914  if ( lhs1.AlmostEquals ( rhs1 ) )
1915  {
1916  if ( ProSHADE_internal_maths::vectorOrientationSimilarity ( axVec->at(symIt)[1],
1917  axVec->at(symIt)[2],
1918  axVec->at(symIt)[3],
1919  hlpSym[1],
1920  hlpSym[2],
1921  hlpSym[3],
1922  axErr ) )
1923  {
1924  //==================================== Almost identical entry
1925  if ( axVec->at(symIt)[5] < hlpSym[5] )
1926  {
1927  //================================ If higher, save
1928  delete[] axVec->at(symIt);
1929  axVec->at(symIt) = hlpSym;
1930  return ;
1931  }
1932  else
1933  {
1934  //================================ or just terminate if better is already saved
1935  delete[] hlpSym;
1936  return ;
1937  }
1938  }
1939  }
1940  }
1941 
1942  //================================================ Not matched to anything
1944 
1945  //================================================ Done
1946  return ;
1947 
1948 }
1949 
1964 void ProSHADE_internal_symmetry::searchMissingSymmetrySpace ( ProSHADE_internal_data::ProSHADE_data* dataObj, std::vector< proshade_double* >* CSymList, std::vector< proshade_unsign >* grp, std::vector< proshade_double* >* hlpVec, proshade_double axErr, proshade_double angle, proshade_unsign fold, proshade_double minPeakHeight )
1965 {
1966  //================================================ Sanity check
1967  if ( grp->size() < 2 ) { return; }
1968 
1969  //================================================ Initialise variables
1970  proshade_double axHeight = 0.0;
1971  proshade_double* symHlp = new proshade_double[7];
1972  ProSHADE_internal_misc::checkMemoryAllocation ( symHlp, __FILE__, __LINE__, __func__ );
1973 
1974  //================================================ For each axis pair in the group, find the possible solutions
1975  for ( proshade_unsign fAx = 0; fAx < static_cast<proshade_unsign> ( grp->size() ); fAx++ )
1976  {
1977  for ( proshade_unsign sAx = 1; sAx < static_cast<proshade_unsign> ( grp->size() ); sAx++ )
1978  {
1979  //======================================== Only unique pairs
1980  if ( fAx >= sAx ) { continue; }
1981 
1982  //======================================== Find possible axis having the required angle to this pair ( solution 1 )
1983  std::vector< proshade_double > solVec = ProSHADE_internal_maths::findVectorFromTwoVAndTwoD ( CSymList->at(grp->at(fAx))[1],
1984  CSymList->at(grp->at(fAx))[2],
1985  CSymList->at(grp->at(fAx))[3],
1986  CSymList->at(grp->at(sAx))[1],
1987  CSymList->at(grp->at(sAx))[2],
1988  CSymList->at(grp->at(sAx))[3], angle, angle );
1989 
1990  //======================================== Set largest axis element to positive
1991  const FloatingPoint< proshade_double > lhs1 ( std::max ( std::abs ( solVec.at(0) ), std::max( std::abs ( solVec.at(1) ), std::abs ( solVec.at(2) ) ) ) );
1992  const FloatingPoint< proshade_double > rhs1 ( std::abs ( solVec.at(0) ) );
1993  const FloatingPoint< proshade_double > rhs2 ( std::abs ( solVec.at(1) ) );
1994  const FloatingPoint< proshade_double > rhs3 ( std::abs ( solVec.at(2) ) );
1995  if ( ( lhs1.AlmostEquals ( rhs1 ) && ( solVec.at(0) < 0.0 ) ) ||
1996  ( lhs1.AlmostEquals ( rhs2 ) && ( solVec.at(1) < 0.0 ) ) ||
1997  ( lhs1.AlmostEquals ( rhs3 ) && ( solVec.at(2) < 0.0 ) ) )
1998  {
1999  solVec.at(0) *= -1.0;
2000  solVec.at(1) *= -1.0;
2001  solVec.at(2) *= -1.0;
2002  }
2003 
2004  //======================================== Does the solution fit the whole group?
2005  symHlp[1] = solVec.at(0); symHlp[2] = solVec.at(1); symHlp[3] = solVec.at(2); symHlp[6] = -1.0;
2006  if ( ProSHADE_internal_symmetry::testGroupAgainstSymmetry ( CSymList, grp, symHlp, axErr, angle, true ) )
2007  {
2008  //==================================== Find the height for the axis
2009  axHeight = ProSHADE_internal_symmetry::missingAxisHeight ( solVec.at(0), solVec.at(1), solVec.at(2), dataObj, fold, axErr );
2010 
2011  //================================ Save max height result
2012  if ( axHeight >= minPeakHeight ) { ProSHADE_internal_symmetry::saveMissingAxisNewOnly ( hlpVec, solVec.at(0), solVec.at(1), solVec.at(2), axHeight, fold, axErr ); }
2013  }
2014 
2015  //======================================== Find possible axis having the required angle to this pair ( solution 2 )
2016  solVec = ProSHADE_internal_maths::findVectorFromTwoVAndTwoD ( CSymList->at(grp->at(fAx))[1],
2017  CSymList->at(grp->at(fAx))[2],
2018  CSymList->at(grp->at(fAx))[3],
2019  CSymList->at(grp->at(sAx))[1],
2020  CSymList->at(grp->at(sAx))[2],
2021  CSymList->at(grp->at(sAx))[3], -angle, -angle );
2022 
2023  //======================================== Set largest axis element to positive
2024  const FloatingPoint< proshade_double > lhs2 ( std::max ( std::abs ( solVec.at(0) ), std::max( std::abs ( solVec.at(1) ), std::abs ( solVec.at(2) ) ) ) );
2025  const FloatingPoint< proshade_double > rhs4 ( std::abs ( solVec.at(0) ) );
2026  const FloatingPoint< proshade_double > rhs5 ( std::abs ( solVec.at(1) ) );
2027  const FloatingPoint< proshade_double > rhs6 ( std::abs ( solVec.at(2) ) );
2028  if ( ( lhs2.AlmostEquals ( rhs4 ) && ( solVec.at(0) < 0.0 ) ) ||
2029  ( lhs2.AlmostEquals ( rhs5 ) && ( solVec.at(1) < 0.0 ) ) ||
2030  ( lhs2.AlmostEquals ( rhs6 ) && ( solVec.at(2) < 0.0 ) ) )
2031  {
2032  solVec.at(0) *= -1.0;
2033  solVec.at(1) *= -1.0;
2034  solVec.at(2) *= -1.0;
2035  }
2036 
2037  //======================================== Does the solution fit the whole group?
2038  symHlp[1] = solVec.at(0); symHlp[2] = solVec.at(1); symHlp[3] = solVec.at(2); symHlp[6] = -1.0;
2039  if ( ProSHADE_internal_symmetry::testGroupAgainstSymmetry ( CSymList, grp, symHlp, axErr, angle, true ) )
2040  {
2041  //==================================== Find the height for the axis
2042  axHeight = ProSHADE_internal_symmetry::missingAxisHeight ( solVec.at(0), solVec.at(1), solVec.at(2), dataObj, fold, axErr );
2043 
2044  //================================ Save max height result
2045  if ( axHeight >= minPeakHeight ) { ProSHADE_internal_symmetry::saveMissingAxisNewOnly ( hlpVec, solVec.at(0), solVec.at(1), solVec.at(2), axHeight, fold, axErr ); }
2046  }
2047  }
2048  }
2049 
2050  //================================================ Release memory
2051  delete[] symHlp;
2052 
2053  //================================================ Done
2054  return ;
2055 
2056 }
2057 
2071 void ProSHADE_internal_symmetry::findTetra3C2s ( std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_signed verbose, proshade_double minPeakHeight )
2072 {
2073  //================================================ Initialise variables
2074  std::vector< proshade_unsign > C3s, prospectiveC2s, C2PossibilitiesHlp;
2075  std::vector< std::vector< proshade_unsign > > C2Possibilities;
2076  proshade_double dotProd;
2077  bool groupMatched;
2078  for ( proshade_unsign iter = 0; iter < 4; iter++ ) { ProSHADE_internal_misc::addToUnsignVector ( &C3s, iter ); }
2079 
2080  //================================================ Report progress
2081  ProSHADE_internal_messages::printProgressMessage ( verbose, 2, "Starting detection of three C2 axes." );
2082 
2083  //================================================ For each C3
2084  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( ret->size() ); rIt++ )
2085  {
2086  //============================================ For each C2, check it has angle ( acos(0.5) ) to the tested C3
2087  for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
2088  {
2089  //======================================== Search only using C2s
2090  const FloatingPoint< proshade_double > lhs999 ( CSymList->at(cIt)[5] ), rhs999 ( static_cast< proshade_double > ( -999.9 ) );
2091  if ( CSymList->at(cIt)[0] != 2.0 || ( ( CSymList->at(cIt)[5] < minPeakHeight ) && !( lhs999.AlmostEquals( rhs999 ) ) ) ) { continue; }
2092 
2093  //======================================== Check the C2 axis to the C3 ( acos ( 0.5 ) )
2094  dotProd = ProSHADE_internal_maths::computeDotProduct ( &ret->at(rIt)[1], &ret->at(rIt)[2], &ret->at(rIt)[3],
2095  &CSymList->at(cIt)[1], &CSymList->at(cIt)[2], &CSymList->at(cIt)[3] );
2096 
2097  if ( ( std::abs ( dotProd ) > ( 0.5 - axErr ) ) && ( std::abs ( dotProd ) < ( 0.5 + axErr ) ) ) { ProSHADE_internal_misc::addToUnsignVector ( &prospectiveC2s, cIt ); }
2098  }
2099  }
2100 
2101  //================================================ Group the prospective C2s
2102  C2Possibilities.clear(); C2PossibilitiesHlp.clear();
2103  for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( prospectiveC2s.size() ); cIt++ )
2104  {
2105  //============================================ If second or more C2, check if it can be placed in any group with being perpendicular to all its members
2106  groupMatched = false;
2107  for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( C2Possibilities.size() ); gIt++ )
2108  {
2109  if ( ProSHADE_internal_symmetry::testGroupAgainstSymmetry ( CSymList, &C2Possibilities.at(gIt), CSymList->at(prospectiveC2s.at(cIt)), axErr, 0.0, true, prospectiveC2s.at(cIt) ) ) { ProSHADE_internal_misc::addToUnsignVector ( &C2Possibilities.at(gIt), prospectiveC2s.at(cIt) ); groupMatched = true; break; }
2110  }
2111 
2112  //============================================ If no group matched, create a new group
2113  if ( !groupMatched ) { C2PossibilitiesHlp.clear(); ProSHADE_internal_misc::addToUnsignVector ( &C2PossibilitiesHlp, prospectiveC2s.at(cIt) ); ProSHADE_internal_misc::addToUnsignVectorVector ( &C2Possibilities, C2PossibilitiesHlp ); continue; }
2114  }
2115 
2116  //================================================ Find the best group or return empty
2117  while ( C2Possibilities.size() != 0 )
2118  {
2119  //============================================ Test for missing symmetry axes, if need be
2120  ProSHADE_internal_symmetry::findMissingAxes ( &C2Possibilities, CSymList, 3, axErr, 0.0, 2, dataObj, minPeakHeight );
2121 
2122  //============================================ Found 3 C2s?
2123  if ( C2Possibilities.at(0).size() == 3 )
2124  {
2125  //======================================== Success! Save and exit
2126  for ( proshade_unsign it = 0; it < 3; it++ ) { ProSHADE_internal_misc::addToDblPtrVector ( ret, CSymList->at(C2Possibilities.at(0).at(it)) ); }
2127 
2128  //======================================== Report progress
2129  ProSHADE_internal_messages::printProgressMessage ( verbose, 3, "Detection of three C2 axes successfull." );
2130 
2131  //======================================== Done
2132  return ;
2133  }
2134  else { C2Possibilities.erase ( C2Possibilities.begin() ); }
2135  }
2136 
2137  //================================================ Done
2138  return ;
2139 
2140 }
2141 
2156 bool ProSHADE_internal_symmetry::testGroupAgainstGroup ( std::vector< proshade_double* >* GrList1, std::vector< proshade_unsign >* grp1, std::vector< proshade_double* >* GrList2, std::vector< proshade_unsign >* grp2, proshade_double angle, proshade_double axErr )
2157 {
2158  //================================================ Initialise variables
2159  bool ret = false;
2160  proshade_double dotProduct;
2161 
2162  //================================================ For all pairs of axes
2163  for ( proshade_unsign g1It = 0; g1It < static_cast<proshade_unsign> ( grp1->size() ); g1It++ )
2164  {
2165  for ( proshade_unsign g2It = 0; g2It < static_cast<proshade_unsign> ( grp2->size() ); g2It++ )
2166  {
2167  //======================================== Find the angle
2168  dotProduct = ProSHADE_internal_maths::computeDotProduct ( &GrList1->at(grp1->at(g1It))[1],
2169  &GrList1->at(grp1->at(g1It))[2],
2170  &GrList1->at(grp1->at(g1It))[3],
2171  &GrList2->at(grp2->at(g2It))[1],
2172  &GrList2->at(grp2->at(g2It))[2],
2173  &GrList2->at(grp2->at(g2It))[3] );
2174 
2175  //======================================== Check the angle
2176  if ( ( angle > ( dotProduct - axErr ) ) && ( angle < ( dotProduct + axErr ) ) )
2177  {
2178  ret = true;
2179  return ( ret );
2180  }
2181  }
2182  }
2183 
2184  //================================================ Done
2185  return ( ret );
2186 
2187 }
2188 
2199 std::vector< proshade_double* > ProSHADE_internal_data::ProSHADE_data::getOctahedralSymmetriesList ( ProSHADE_settings* settings, std::vector< proshade_double* >* CSymList )
2200 {
2201  //================================================ Initialise variables
2202  std::vector< proshade_double* > ret;
2203 
2204  //================================================ Report progress
2205  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting O symmetry detection." );
2206 
2207  //================================================ Are the basic requirements for tetrahedral symmetry met?
2208  if ( ProSHADE_internal_symmetry::detectOctahedralSymmetry ( CSymList, settings->axisErrTolerance, settings->minSymPeak ) )
2209  {
2210  //============================================ Search for all the symmetry axes
2211  ProSHADE_internal_symmetry::findOcta3C4s ( CSymList, &ret, settings->axisErrTolerance, this, settings->verbose, settings->minSymPeak );
2212  if ( ret.size() != 3 ) { ProSHADE_internal_messages::printWarningMessage ( settings->verbose, "!!! ProSHADE WARNING !!! Failed to detect some of the polyhedral symmetries, while detecting the correct dihedral angles.", "WS00031" ); return ( ret ); }
2213 
2214  ProSHADE_internal_symmetry::findOcta4C3s ( CSymList, &ret, settings->axisErrTolerance, this, settings->verbose, settings->minSymPeak );
2215  if ( ret.size() != 7 ) { ProSHADE_internal_messages::printWarningMessage ( settings->verbose, "!!! ProSHADE WARNING !!! Failed to detect some of the polyhedral symmetries, while detecting the correct dihedral angles.", "WS00031" ); return ( ret ); }
2216 
2217  ProSHADE_internal_symmetry::findOcta6C2s ( CSymList, &ret, settings->axisErrTolerance, this, settings->verbose, settings->minSymPeak );
2218  if ( ret.size() != 13 ) { ProSHADE_internal_messages::printWarningMessage ( settings->verbose, "!!! ProSHADE WARNING !!! Failed to detect some of the polyhedral symmetries, while detecting the correct dihedral angles.", "WS00031" ); return ( ret ); }
2219  else
2220  {
2221  for ( proshade_unsign csIt = 0; csIt < static_cast<proshade_unsign> ( CSymList->size() ); csIt++ )
2222  {
2223  for ( proshade_unsign retIt = 0; retIt < static_cast<proshade_unsign> ( ret.size() ); retIt++ )
2224  {
2225  //======================================== Sort ret by fold
2226  std::sort ( ret.begin(), ret.end(), ProSHADE_internal_misc::sortSymInvFoldHlp );
2227 
2228  //======================================== Save indices
2229  const FloatingPoint< proshade_double > lhs1 ( CSymList->at(csIt)[0] ), rhs1 ( ret.at(retIt)[0] );
2230  const FloatingPoint< proshade_double > lhs2 ( CSymList->at(csIt)[1] ), rhs2 ( ret.at(retIt)[1] );
2231  const FloatingPoint< proshade_double > lhs3 ( CSymList->at(csIt)[2] ), rhs3 ( ret.at(retIt)[2] );
2232  const FloatingPoint< proshade_double > lhs4 ( CSymList->at(csIt)[3] ), rhs4 ( ret.at(retIt)[3] );
2233  const FloatingPoint< proshade_double > lhs5 ( CSymList->at(csIt)[4] ), rhs5 ( ret.at(retIt)[4] );
2234  const FloatingPoint< proshade_double > lhs6 ( CSymList->at(csIt)[5] ), rhs6 ( ret.at(retIt)[5] );
2235  if ( lhs1.AlmostEquals ( rhs1 ) &&
2236  lhs2.AlmostEquals ( rhs2 ) &&
2237  lhs3.AlmostEquals ( rhs3 ) &&
2238  lhs4.AlmostEquals ( rhs4 ) &&
2239  lhs5.AlmostEquals ( rhs5 ) &&
2240  lhs6.AlmostEquals ( rhs6 ) )
2241  {
2243  }
2244  }
2245  }
2246  }
2247  }
2248 
2249  //================================================ Report progress
2250  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "O symmetry detection complete." );
2251 
2252  //================================================ Done
2253  return ( ret );
2254 
2255 }
2256 
2268 bool ProSHADE_internal_symmetry::detectOctahedralSymmetry ( std::vector< proshade_double* >* CSymList, proshade_double axErr, proshade_double minPeakHeight )
2269 {
2270  //================================================ Initialise variables
2271  std::vector< proshade_unsign > C4List;
2272  proshade_double dotProduct;
2273 
2274  //================================================ Find all C4 symmetries
2275  for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
2276  {
2277  const FloatingPoint< proshade_double > lhs1 ( CSymList->at(cSym)[0] ), rhs1 ( 4.0 );
2278  if ( lhs1.AlmostEquals ( rhs1 ) && CSymList->at(cSym)[5] >= minPeakHeight ) { ProSHADE_internal_misc::addToUnsignVector ( &C4List, cSym ); }
2279  }
2280 
2281  //================================================ For each unique pair of C3s
2282  for ( proshade_unsign c4 = 0; c4 < static_cast<proshade_unsign> ( C4List.size() ); c4++ )
2283  {
2284  for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
2285  {
2286  //======================================== Compare only C3s to the C3List
2287  const FloatingPoint< proshade_double > lhs1 ( CSymList->at(cSym)[0] ), rhs1 ( 3.0 );
2288  if ( !lhs1.AlmostEquals ( rhs1 ) ) { continue; }
2289 
2290  //======================================== Check the angle between the C4 and C3 axes
2291  dotProduct = ProSHADE_internal_maths::computeDotProduct ( &CSymList->at(C4List.at(c4))[1],
2292  &CSymList->at(C4List.at(c4))[2],
2293  &CSymList->at(C4List.at(c4))[3],
2294  &CSymList->at(cSym)[1],
2295  &CSymList->at(cSym)[2],
2296  &CSymList->at(cSym)[3] );
2297 
2298  //======================================== Is the angle approximately the dihedral angle
2299  if ( ( ( 1.0 / sqrt ( 3.0 ) ) > ( dotProduct - axErr ) ) && ( ( 1.0 / sqrt ( 3.0 ) ) < ( dotProduct + axErr ) ) )
2300  {
2301  return ( true );
2302  }
2303  }
2304  }
2305 
2306  //================================================ Done
2307  return ( false );
2308 
2309 }
2310 
2324 void ProSHADE_internal_symmetry::findOcta3C4s ( std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_signed verbose, proshade_double minPeakHeight )
2325 {
2326  //================================================ Initialise variables
2327  std::vector< proshade_unsign > C4PossibilitiesHlp;
2328  std::vector< std::vector< proshade_unsign > > C4Possibilities;
2329  bool groupMatched;
2330 
2331  //================================================ Report progress
2332  ProSHADE_internal_messages::printProgressMessage ( verbose, 2, "Starting detection of three C4 axes." );
2333 
2334  //================================================ For all symmetries in the C symmetries list
2335  for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
2336  {
2337  //============================================ Search only using C4s
2338  if ( CSymList->at(cIt)[0] != 4.0 || CSymList->at(cIt)[5] < minPeakHeight ) { continue; }
2339 
2340  //============================================ If second or more C4, check if it has the correct angle to all other already found C4s for each group
2341  groupMatched = false;
2342  for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( C4Possibilities.size() ); gIt++ )
2343  {
2344  if ( ProSHADE_internal_symmetry::testGroupAgainstSymmetry ( CSymList, &C4Possibilities.at(gIt), CSymList->at(cIt), axErr, 0.0, true, cIt ) ) { ProSHADE_internal_misc::addToUnsignVector ( &C4Possibilities.at(gIt), cIt ); groupMatched = true; break; }
2345  }
2346 
2347  //=========================================== If no group matched, create a new group
2348  if ( !groupMatched ) { C4PossibilitiesHlp.clear(); ProSHADE_internal_misc::addToUnsignVector ( &C4PossibilitiesHlp, cIt ); ProSHADE_internal_misc::addToUnsignVectorVector ( &C4Possibilities, C4PossibilitiesHlp ); continue; }
2349  }
2350 
2351  //================================================ Test for missing symmetry axes, if need be
2352  ProSHADE_internal_symmetry::findMissingAxes ( &C4Possibilities, CSymList, 3, axErr, 0.0, 4, dataObj, minPeakHeight );
2353 
2354  //================================================ Any group has 3 entries? If more such groups, take the one with highest average height.
2355  proshade_double maxHeight = 0.0; proshade_unsign maxGrp = 0;
2356  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( C4Possibilities.size() ); iter++ ) { if ( C4Possibilities.at(iter).size() == 3 ) { if ( ( ( CSymList->at(C4Possibilities.at(iter).at(0))[5] + CSymList->at(C4Possibilities.at(iter).at(1))[5] + CSymList->at(C4Possibilities.at(iter).at(2))[5] ) / 3.0 ) > maxHeight ) { maxHeight = ( ( CSymList->at(C4Possibilities.at(iter).at(0))[5] + CSymList->at(C4Possibilities.at(iter).at(1))[5] + CSymList->at(C4Possibilities.at(iter).at(2))[5] ) / 3.0 ); maxGrp = iter; } } }
2357 
2358  if ( C4Possibilities.at(maxGrp).size() == 3 )
2359  {
2360  //============================================ Success! Save and exit
2361  for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( C4Possibilities.at(maxGrp).size() ); it++ ) { ProSHADE_internal_misc::addToDblPtrVector ( ret, CSymList->at(C4Possibilities.at(maxGrp).at(it)) ); }
2362 
2363  //============================================ Report progress
2364  ProSHADE_internal_messages::printProgressMessage ( verbose, 3, "Detection of three C4 axes successfull." );
2365 
2366  //============================================ Done
2367  return ;
2368  }
2369 
2370  //================================================ Done
2371  return ;
2372 
2373 }
2374 
2390 void ProSHADE_internal_symmetry::findOcta4C3s ( std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_signed verbose, proshade_double minPeakHeight )
2391 {
2392  //================================================ Initialise variables
2393  std::vector< proshade_unsign > C4s, prospectiveC3s, C3PossibilitiesHlp;
2394  std::vector< std::vector< proshade_unsign > > C3Possibilities;
2395  proshade_double dotProd;
2396  bool groupMatched;
2397  for ( proshade_unsign iter = 0; iter < 3; iter++ ) { ProSHADE_internal_misc::addToUnsignVector ( &C4s, iter ); }
2398 
2399  //================================================ Report progress
2400  ProSHADE_internal_messages::printProgressMessage ( verbose, 2, "Starting detection of four C3 axes." );
2401 
2402  //================================================ For each C4
2403  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( ret->size() ); rIt++ )
2404  {
2405  //============================================ For each C3, check it has angle ( acos( 1/sqrt(3) ) ) to the tested C4
2406  for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
2407  {
2408  //======================================== Search only using C3s
2409  if ( CSymList->at(cIt)[0] != 3.0 || CSymList->at(cIt)[5] < minPeakHeight ) { continue; }
2410 
2411  //======================================== Check the C3 axis to the C4 ( acos ( 1/sqrt(3) ) )
2412  dotProd = ProSHADE_internal_maths::computeDotProduct ( &ret->at(rIt)[1], &ret->at(rIt)[2], &ret->at(rIt)[3], &CSymList->at(cIt)[1], &CSymList->at(cIt)[2], &CSymList->at(cIt)[3] );
2413 
2414  if ( ( std::abs ( dotProd ) > ( ( 1.0 / sqrt(3.0) ) - axErr ) ) && ( std::abs ( dotProd ) < ( ( 1.0 / sqrt(3.0) ) + axErr ) ) ) { ProSHADE_internal_misc::addToUnsignVector ( &prospectiveC3s, cIt ); }
2415  }
2416  }
2417 
2418  //================================================ Group the prospective C3s
2419  C3Possibilities.clear(); C3PossibilitiesHlp.clear();
2420  for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( prospectiveC3s.size() ); cIt++ )
2421  {
2422  //============================================ If second or more C3, check if it can be placed in any group with having acos (1/3) to all its members
2423  groupMatched = false;
2424  for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( C3Possibilities.size() ); gIt++ )
2425  {
2426  if ( ProSHADE_internal_symmetry::testGroupAgainstSymmetry ( CSymList, &C3Possibilities.at(gIt), CSymList->at(prospectiveC3s.at(cIt)), axErr, 1.0/3.0, true, prospectiveC3s.at(cIt) ) ) { ProSHADE_internal_misc::addToUnsignVector ( &C3Possibilities.at(gIt), prospectiveC3s.at(cIt) ); groupMatched = true; break; }
2427  }
2428 
2429  //============================================ If no group matched, create a new group
2430  if ( !groupMatched ) { C3PossibilitiesHlp.clear(); ProSHADE_internal_misc::addToUnsignVector ( &C3PossibilitiesHlp, prospectiveC3s.at(cIt) ); ProSHADE_internal_misc::addToUnsignVectorVector ( &C3Possibilities, C3PossibilitiesHlp ); continue; }
2431  }
2432 
2433  //================================================ Find the best group or return empty
2434  while ( C3Possibilities.size() != 0 )
2435  {
2436  //============================================ Test for missing symmetry axes, if need be
2437  ProSHADE_internal_symmetry::findMissingAxes ( &C3Possibilities, CSymList, 4, axErr, 1.0/3.0, 3, dataObj, minPeakHeight );
2438 
2439  //============================================ Found four C3s?
2440  if ( C3Possibilities.at(0).size() == 4 )
2441  {
2442  //======================================== Success! Save and exit
2443  for ( proshade_unsign it = 0; it < 4; it++ ) { ProSHADE_internal_misc::addToDblPtrVector ( ret, CSymList->at(C3Possibilities.at(0).at(it)) ); }
2444 
2445  //======================================== Report progress
2446  ProSHADE_internal_messages::printProgressMessage ( verbose, 3, "Detection of four C3 axes successfull." );
2447 
2448  //======================================== Done
2449  return ;
2450  }
2451  else { C3Possibilities.erase ( C3Possibilities.begin() ); }
2452  }
2453 
2454  //================================================ Done
2455  return ;
2456 
2457 }
2458 
2472 void ProSHADE_internal_symmetry::findOcta6C2s ( std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_signed verbose, proshade_double minPeakHeight )
2473 {
2474  //================================================ Initialise variables
2475  std::vector< proshade_unsign > prospectiveC2s, retGrp;
2476  proshade_double dotProd;
2477  proshade_unsign noPerpendicular, noSqrtTwo;
2478 
2479  //================================================ Report progress
2480  ProSHADE_internal_messages::printProgressMessage ( verbose, 2, "Starting detection of six C2 axes." );
2481 
2482  //================================================ For each C2
2483  for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
2484  {
2485  //============================================ Use only C2s
2486  const FloatingPoint< proshade_double > lhs999 ( CSymList->at(cIt)[5] ), rhs999 ( static_cast< proshade_double > ( -999.9 ) );
2487  if ( CSymList->at(cIt)[0] != 2.0 || ( ( CSymList->at(cIt)[5] < minPeakHeight ) && ! ( lhs999.AlmostEquals( rhs999 ) ) ) ) { continue; }
2488 
2489  //============================================ Check the C2 has acos ( 1/sqrt(2) ) to 2 C4s and acos ( 0.0 ) to the third C4
2490  noPerpendicular = 0; noSqrtTwo = 0;
2491  for ( proshade_unsign rIt = 0; rIt < 3; rIt++ )
2492  {
2493  dotProd = ProSHADE_internal_maths::computeDotProduct ( &ret->at(rIt)[1],
2494  &ret->at(rIt)[2],
2495  &ret->at(rIt)[3],
2496  &CSymList->at(cIt)[1],
2497  &CSymList->at(cIt)[2],
2498  &CSymList->at(cIt)[3] );
2499 
2500  if ( ( std::abs ( dotProd ) > ( ( 1.0 / sqrt(2.0) ) - axErr ) ) && ( std::abs ( dotProd ) < ( ( 1.0 / sqrt(2.0) ) + axErr ) ) ) { noSqrtTwo += 1; continue; }
2501  if ( ( std::abs ( dotProd ) > ( 0.0 - axErr ) ) && ( std::abs ( dotProd ) < ( 0.0 + axErr ) ) ) { noPerpendicular += 1; continue; }
2502  }
2503 
2504  //============================================ If correct angles distribution is found, save the axis
2505  if ( ( noSqrtTwo == 2 ) && ( noPerpendicular == 1 ) )
2506  {
2507  ProSHADE_internal_misc::addToUnsignVector ( &prospectiveC2s, cIt );
2508  }
2509  }
2510 
2511  //================================================ Search for missing axes
2512  for ( proshade_unsign iter = 0; iter < 3; iter++ ) { ProSHADE_internal_misc::addToUnsignVector ( &retGrp, iter ); }
2513  if ( !ProSHADE_internal_symmetry::findMissingAxesDual ( &prospectiveC2s, CSymList, ret, &retGrp, 6, axErr, 1, 0.0, 2, 1/sqrt(2.0), 2, dataObj ) )
2514  {
2515  return ;
2516  }
2517 
2518  //================================================ Found correct number of axes! Now save the
2519  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( prospectiveC2s.size() ); iter++ )
2520  {
2521  ProSHADE_internal_misc::addToDblPtrVector ( ret, CSymList->at(prospectiveC2s.at(iter)) );
2522  }
2523 
2524  //================================================ Report progress
2525  ProSHADE_internal_messages::printProgressMessage ( verbose, 3, "Detection of six C2 axes successfull." );
2526 
2527  //================================================ Done
2528  return ;
2529 
2530 }
2531 
2553 bool ProSHADE_internal_symmetry::findMissingAxesDual ( std::vector< proshade_unsign >* possibilities, std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, std::vector< proshade_unsign >* retGroup, proshade_unsign requiredNoAxes, proshade_double axErr, proshade_unsign noMatchesG1, proshade_double angle1, proshade_unsign noMatchesG2, proshade_double angle2, proshade_unsign fold, ProSHADE_internal_data::ProSHADE_data* dataObj )
2554 {
2555  //================================================ Initialise variables
2556  bool atLeastOne = false;
2557  std::vector< proshade_double* > prosp;
2558  std::vector< proshade_double > sol;
2559 
2560  //================================================ Proceed only if need be
2561  if ( static_cast<proshade_unsign> ( possibilities->size() ) == requiredNoAxes ) { atLeastOne = true; return ( atLeastOne ); }
2562 
2563  //================================================ Copy already found to prospective
2564  for ( proshade_unsign prIt = 0; prIt < static_cast<proshade_unsign> ( possibilities->size() ); prIt++ )
2565  {
2566  ProSHADE_internal_symmetry::addAxisUnlessSame ( static_cast< proshade_unsign > ( CSymList->at(possibilities->at(prIt))[0] ),
2567  CSymList->at(possibilities->at(prIt))[1],
2568  CSymList->at(possibilities->at(prIt))[2],
2569  CSymList->at(possibilities->at(prIt))[3],
2570  CSymList->at(possibilities->at(prIt))[5], &prosp, axErr );
2571  }
2572 
2573  //================================================ Start generating possible solutions
2574  for ( proshade_unsign rgIt1 = 0; rgIt1 < static_cast<proshade_unsign> ( retGroup->size() ); rgIt1++ )
2575  {
2576  for ( proshade_unsign rgIt2 = 0; rgIt2 < static_cast<proshade_unsign> ( retGroup->size() ); rgIt2++ )
2577  {
2578  //======================================== Use unique combinations (order matters here!)
2579  if ( rgIt1 == rgIt2 ) { continue; }
2580 
2581  //======================================== Generate possible solution (1)
2582  sol = ProSHADE_internal_maths::findVectorFromTwoVAndTwoD ( ret->at(rgIt1)[1], ret->at(rgIt1)[2], ret->at(rgIt1)[3],
2583  ret->at(rgIt2)[1], ret->at(rgIt2)[2], ret->at(rgIt2)[3], angle1, angle2 );
2584 
2585  //======================================== Check if solution fits the group completely
2586  ProSHADE_internal_symmetry::checkFittingAxisDualAndSave ( retGroup, ret, fold, sol.at(0), sol.at(1), sol.at(2), &prosp, axErr, noMatchesG1, angle1, noMatchesG2, angle2, dataObj );
2587  if ( prosp.size() == requiredNoAxes ) { break; }
2588 
2589  //======================================== Generate possible solution (2)
2590  sol = ProSHADE_internal_maths::findVectorFromTwoVAndTwoD ( ret->at(rgIt1)[1], ret->at(rgIt1)[2], ret->at(rgIt1)[3],
2591  ret->at(rgIt2)[1], ret->at(rgIt2)[2], ret->at(rgIt2)[3], -angle1, -angle2 );
2592 
2593  //======================================== Check if solution fits the group completely
2594  ProSHADE_internal_symmetry::checkFittingAxisDualAndSave ( retGroup, ret, fold, sol.at(0), sol.at(1), sol.at(2), &prosp, axErr, noMatchesG1, angle1, noMatchesG2, angle2, dataObj );
2595  if ( prosp.size() == requiredNoAxes ) { break; }
2596  }
2597 
2598  if ( prosp.size() == requiredNoAxes ) { break; }
2599  }
2600 
2601  //================================================ Found all required axes!
2602  if ( static_cast<proshade_unsign> ( prosp.size() ) == requiredNoAxes )
2603  {
2604  //============================================ Copy the detected axes
2605  for ( proshade_unsign iter = static_cast<proshade_unsign> ( possibilities->size() ); iter < static_cast<proshade_unsign> ( prosp.size() ); iter++ )
2606  {
2607  if ( ProSHADE_internal_maths::isAxisUnique ( CSymList, prosp.at(iter), axErr ) )
2608  {
2609  //==================================== Add
2610  ProSHADE_internal_misc::addToUnsignVector ( possibilities, static_cast< proshade_unsign > ( CSymList->size() ) );
2611  ProSHADE_internal_misc::addToDblPtrVector ( CSymList, prosp.at(iter) );
2612  }
2613  }
2614 
2615  //============================================ Done
2616  atLeastOne = true;
2617  return ( atLeastOne );
2618  }
2619  else
2620  {
2621  //============================================ Delete the created, but not used axes
2622  for ( proshade_unsign iter = static_cast<proshade_unsign> ( possibilities->size() ); iter < static_cast<proshade_unsign> ( prosp.size() ); iter++ )
2623  {
2624  delete[] prosp.at(iter);
2625  }
2626  }
2627 
2628  //================================================ Done
2629  return ( atLeastOne );
2630 
2631 }
2632 
2648 proshade_signed ProSHADE_internal_symmetry::addAxisUnlessSame ( proshade_unsign fold, proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double axHeight, proshade_double averageFSC, std::vector< proshade_double* >* prosp, proshade_double axErr )
2649 {
2650  //================================================ Initialise variables
2651  proshade_double* symHlp = new proshade_double[7];
2652  ProSHADE_internal_misc::checkMemoryAllocation ( symHlp, __FILE__, __LINE__, __func__ );
2653  proshade_signed ret = -1;
2654 
2655  //================================================ Fill in the prospective axis
2656  symHlp[0] = static_cast<proshade_double> ( fold );
2657  symHlp[1] = axX;
2658  symHlp[2] = axY;
2659  symHlp[3] = axZ;
2660  symHlp[4] = 2.0 * M_PI / symHlp[0];
2661  symHlp[5] = axHeight;
2662  symHlp[6] = averageFSC;
2663 
2664  //================================================ If it is not the same as already saved axes
2665  if ( !ProSHADE_internal_symmetry::isSymmetrySame ( prosp, symHlp, axErr, &ret, averageFSC ) )
2666  {
2668  return ( static_cast< proshade_signed > ( prosp->size() - 1 ) );
2669  }
2670  else
2671  {
2672  delete[] symHlp;
2673  return ( ret );
2674  }
2675 
2676  //================================================ Done
2677 
2678 }
2679 
2694 proshade_signed ProSHADE_internal_symmetry::addAxisUnlessSame ( proshade_unsign fold, proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double axHeight, std::vector< proshade_double* >* prosp, proshade_double axErr )
2695 {
2696  //================================================ Initialise variables
2697  proshade_double* symHlp = new proshade_double[7];
2698  ProSHADE_internal_misc::checkMemoryAllocation ( symHlp, __FILE__, __LINE__, __func__ );
2699  proshade_signed ret = -1;
2700 
2701  //================================================ Fill in the prospective axis
2702  symHlp[0] = static_cast<proshade_double> ( fold );
2703  symHlp[1] = axX;
2704  symHlp[2] = axY;
2705  symHlp[3] = axZ;
2706  symHlp[4] = 2.0 * M_PI / symHlp[0];
2707  symHlp[5] = axHeight;
2708  symHlp[6] = -1.0;
2709 
2710  //================================================ If it is not the same as already saved axes
2711  if ( !ProSHADE_internal_symmetry::isSymmetrySame ( prosp, symHlp, axErr, &ret ) )
2712  {
2714  return ( static_cast< proshade_signed > ( prosp->size() - 1 ) );
2715  }
2716  else
2717  {
2718  delete[] symHlp;
2719  return ( ret );
2720  }
2721 
2722  //================================================ Done
2723 
2724 }
2725 
2747 bool ProSHADE_internal_symmetry::checkFittingAxisDualAndSave ( std::vector< proshade_unsign >* retGroup, std::vector< proshade_double* >* ret, proshade_unsign fold, proshade_double axX, proshade_double axY, proshade_double axZ, std::vector< proshade_double* >* prosp, proshade_double axErr, proshade_unsign noMatchesG1, proshade_double angle1, proshade_unsign noMatchesG2, proshade_double angle2, ProSHADE_internal_data::ProSHADE_data* dataObj )
2748 {
2749  //================================================ Initialise variables
2750  proshade_unsign noG1 = 0;
2751  proshade_unsign noG2 = 0;
2752  proshade_double dotProd = 0.0;
2753  proshade_double axHeight = 0.0;
2754 
2755  //================================================ Find the angle and count dual matching frequencies
2756  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( retGroup->size() ); rIt++ )
2757  {
2758  dotProd = ProSHADE_internal_maths::computeDotProduct ( &ret->at(retGroup->at(rIt))[1],
2759  &ret->at(retGroup->at(rIt))[2],
2760  &ret->at(retGroup->at(rIt))[3],
2761  &axX, &axY, &axZ );
2762 
2763  if ( ( std::abs ( dotProd ) > ( angle1 - axErr ) ) && ( std::abs ( dotProd ) < ( angle1 + axErr ) ) ) { noG1 += 1; continue; }
2764  if ( ( std::abs ( dotProd ) > ( angle2 - axErr ) ) && ( std::abs ( dotProd ) < ( angle2 + axErr ) ) ) { noG2 += 1; continue; }
2765  }
2766 
2767  //================================================ If correct frequencies are matched, check height.
2768  if ( ( noG1 == noMatchesG1 ) && ( noG2 == noMatchesG2 ) )
2769  {
2770  //============================================ Is the height good enough?
2771  axHeight = ProSHADE_internal_symmetry::missingAxisHeight ( axX, axY, axZ, dataObj, fold, axErr );
2772 
2773  //============================================ If so, save
2774  if ( axHeight > 0.1 )
2775  {
2776  proshade_unsign prevProsp = static_cast<proshade_unsign> ( prosp->size() );
2777  ProSHADE_internal_symmetry::addAxisUnlessSame ( fold, axX, axY, axZ, axHeight, prosp, axErr );
2778 
2779  if ( static_cast<proshade_unsign> ( prosp->size() ) > prevProsp ) { return ( true ); }
2780  else { return ( false ); }
2781  }
2782  }
2783 
2784  //================================================ Done
2785  return ( false );
2786 
2787 }
2788 
2799 std::vector< proshade_double* > ProSHADE_internal_data::ProSHADE_data::getIcosahedralSymmetriesList ( ProSHADE_settings* settings, std::vector< proshade_double* >* CSymList )
2800 {
2801  //================================================ Initialise variables
2802  std::vector< proshade_double* > ret;
2803 
2804  //================================================ Report progress
2805  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting I symmetry detection." );
2806 
2807  //================================================ Are the basic requirements for icosahedral symmetry met?
2809  {
2810  //============================================ Search for all the symmetry axes
2811  ProSHADE_internal_symmetry::findIcos6C5s ( CSymList, &ret, settings->axisErrTolerance, this, settings->verbose, settings->minSymPeak );
2812  if ( ret.size() != 6 ) { ProSHADE_internal_messages::printWarningMessage ( settings->verbose, "!!! ProSHADE WARNING !!! Failed to detect some of the polyhedral symmetries, while detecting the correct dihedral angles.", "WS00031" ); return ( ret ); }
2813 
2814  ProSHADE_internal_symmetry::findIcos10C3s ( CSymList, &ret, settings->axisErrTolerance, this, settings->verbose, settings->minSymPeak );
2815  if ( ret.size() != 16 ) { ProSHADE_internal_messages::printWarningMessage ( settings->verbose, "!!! ProSHADE WARNING !!! Failed to detect some of the polyhedral symmetries, while detecting the correct dihedral angles.", "WS00031" ); return ( ret ); }
2816 
2817  ProSHADE_internal_symmetry::findIcos15C2s ( CSymList, &ret, settings->axisErrTolerance, this, settings->verbose, settings->minSymPeak );
2818  if ( ret.size() != 31 ) { ProSHADE_internal_messages::printWarningMessage ( settings->verbose, "!!! ProSHADE WARNING !!! Failed to detect some of the polyhedral symmetries, while detecting the correct dihedral angles.", "WS00031" ); return ( ret ); }
2819  else
2820  {
2821  //======================================== Sort ret by fold
2822  std::sort ( ret.begin(), ret.end(), ProSHADE_internal_misc::sortSymInvFoldHlp );
2823 
2824  //======================================== Save indices
2825  for ( proshade_unsign csIt = 0; csIt < static_cast<proshade_unsign> ( CSymList->size() ); csIt++ )
2826  {
2827  for ( proshade_unsign retIt = 0; retIt < static_cast<proshade_unsign> ( ret.size() ); retIt++ )
2828  {
2829  const FloatingPoint< proshade_double > lhs1 ( CSymList->at(csIt)[0] ), rhs1 ( ret.at(retIt)[0] );
2830  const FloatingPoint< proshade_double > lhs2 ( CSymList->at(csIt)[1] ), rhs2 ( ret.at(retIt)[1] );
2831  const FloatingPoint< proshade_double > lhs3 ( CSymList->at(csIt)[2] ), rhs3 ( ret.at(retIt)[2] );
2832  const FloatingPoint< proshade_double > lhs4 ( CSymList->at(csIt)[3] ), rhs4 ( ret.at(retIt)[3] );
2833  const FloatingPoint< proshade_double > lhs5 ( CSymList->at(csIt)[4] ), rhs5 ( ret.at(retIt)[4] );
2834  const FloatingPoint< proshade_double > lhs6 ( CSymList->at(csIt)[5] ), rhs6 ( ret.at(retIt)[5] );
2835  if ( lhs1.AlmostEquals ( rhs1 ) &&
2836  lhs2.AlmostEquals ( rhs2 ) &&
2837  lhs3.AlmostEquals ( rhs3 ) &&
2838  lhs4.AlmostEquals ( rhs4 ) &&
2839  lhs5.AlmostEquals ( rhs5 ) &&
2840  lhs6.AlmostEquals ( rhs6 ) )
2841  {
2843  }
2844  }
2845  }
2846  }
2847  }
2848 
2849  //================================================ Report progress
2850  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "I symmetry detection complete." );
2851 
2852  //================================================ Done
2853  return ( ret );
2854 
2855 }
2856 
2871 std::vector < std::vector< proshade_double* > > ProSHADE_internal_data::ProSHADE_data::getPredictedIcosahedralSymmetriesList ( ProSHADE_settings* settings, std::vector< proshade_double* >* CSymList )
2872 {
2873  //================================================ Initialise variables
2874  std::vector< std::vector< proshade_double* > > ret;
2875 
2876  //================================================ Report progress
2877  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting I symmetry prediction." );
2878 
2879  //================================================ Are the basic requirements for icosahedral symmetry met?
2881  {
2882  //============================================ Generate the rest of the axes
2883  ProSHADE_internal_symmetry::predictIcosAxes ( CSymList, &ret, settings->axisErrTolerance, settings->minSymPeak );
2884 
2885  //============================================ For each possible axes pair
2886  for ( size_t pIt = 0; pIt < ret.size(); pIt++ )
2887  {
2888  //======================================== Get heights for the predicted axes
2889  ProSHADE_internal_symmetry::findPredictedAxesHeights ( &(ret.at(pIt)), this, settings );
2890  }
2891  }
2892 
2893 
2894  //================================================ Report progress
2895  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "I symmetry prediction complete." );
2896 
2897  //================================================ Done
2898  return ( ret );
2899 
2900 }
2901 
2916 std::vector< proshade_double* > ProSHADE_internal_data::ProSHADE_data::getPredictedOctahedralSymmetriesList ( ProSHADE_settings* settings, std::vector< proshade_double* >* CSymList )
2917 {
2918  //================================================ Initialise variables
2919  std::vector< proshade_double* > ret;
2920 
2921  //================================================ Report progress
2922  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting O symmetry prediction." );
2923 
2924  //================================================ Are the basic requirements for icosahedral symmetry met?
2925  if ( ProSHADE_internal_symmetry::detectOctahedralSymmetry ( CSymList, settings->axisErrTolerance, settings->minSymPeak ) )
2926  {
2927  //============================================ Generate the rest of the axes
2928  ProSHADE_internal_symmetry::predictOctaAxes ( CSymList, &ret, settings->axisErrTolerance, settings->minSymPeak );
2929 
2930  //============================================ Get heights for the predicted axes
2932 
2933  //============================================ Add predicted axes to detected C axes list and also to the settings Icosahedral symmetry list
2934  for ( proshade_unsign retIt = 0; retIt < static_cast < proshade_unsign > ( ret.size() ); retIt++ )
2935  {
2936  proshade_signed matchedPos = ProSHADE_internal_symmetry::addAxisUnlessSame ( static_cast< proshade_unsign > ( ret.at(retIt)[0] ), ret.at(retIt)[1], ret.at(retIt)[2], ret.at(retIt)[3], ret.at(retIt)[5], CSymList, settings->axisErrTolerance );
2937  ProSHADE_internal_misc::addToUnsignVector ( &settings->allDetectedOAxes, static_cast < proshade_unsign > ( matchedPos ) );
2938  }
2939  }
2940 
2941  //================================================ Report progress
2942  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "O symmetry prediction complete." );
2943 
2944  //================================================ Done
2945  return ( ret );
2946 
2947 }
2948 
2960 bool ProSHADE_internal_symmetry::detectIcosahedralSymmetry ( std::vector< proshade_double* >* CSymList, proshade_double axErr, proshade_double minPeakHeight )
2961 {
2962  //================================================ Initialise variables
2963  std::vector< proshade_unsign > C5List;
2964  proshade_double dotProduct;
2965 
2966  //================================================ Find all C5 symmetries
2967  for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
2968  {
2969  const FloatingPoint< proshade_double > lhs1 ( CSymList->at(cSym)[0] ), rhs1 ( 5.0 );
2970  if ( lhs1.AlmostEquals ( rhs1 ) && CSymList->at(cSym)[5] >= minPeakHeight ) { ProSHADE_internal_misc::addToUnsignVector ( &C5List, cSym ); }
2971  }
2972 
2973  //================================================ For each unique pair of C5 and C3
2974  for ( proshade_unsign c5 = 0; c5 < static_cast<proshade_unsign> ( C5List.size() ); c5++ )
2975  {
2976  for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
2977  {
2978  //======================================== Compare only C3s to the C5List
2979  const FloatingPoint< proshade_double > lhs1 ( CSymList->at(cSym)[0] ), rhs1 ( 3.0 );
2980  if ( !lhs1.AlmostEquals ( rhs1 ) ) { continue; }
2981 
2982  //======================================== Check the angle between the C5 and C3 axes
2983  dotProduct = ProSHADE_internal_maths::computeDotProduct ( &CSymList->at(C5List.at(c5))[1],
2984  &CSymList->at(C5List.at(c5))[2],
2985  &CSymList->at(C5List.at(c5))[3],
2986  &CSymList->at(cSym)[1],
2987  &CSymList->at(cSym)[2],
2988  &CSymList->at(cSym)[3] );
2989 
2990  //======================================== Is the angle approximately the dihedral angle
2991  if ( std::abs ( std::abs( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) - std::abs( dotProduct ) ) < axErr )
2992  {
2993  return ( true );
2994  }
2995  }
2996  }
2997 
2998  //================================================ Done
2999  return ( false );
3000 
3001 }
3002 
3020 void ProSHADE_internal_symmetry::findIcos6C5s ( std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_signed verbose, proshade_double minPeakHeight )
3021 {
3022  //================================================ Initialise variables
3023  std::vector< proshade_unsign > C5PossibilitiesHlp;
3024  std::vector< std::vector< proshade_unsign > > C5Possibilities;
3025  bool groupMatched;
3026 
3027  //================================================ Report progress
3028  ProSHADE_internal_messages::printProgressMessage ( verbose, 2, "Starting detection of six C5 axes." );
3029 
3030  //================================================ For all symmetries in the C symmetries list
3031  for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
3032  {
3033  //============================================ Search only using C5s and check peak height
3034  const FloatingPoint< proshade_double > lhs1 ( CSymList->at(cIt)[0] ), rhs1 ( 5.0 );
3035  if ( !lhs1.AlmostEquals ( rhs1 ) || CSymList->at(cIt)[5] < minPeakHeight ) { continue; }
3036 
3037  //============================================ If second or more C5, check if it has the correct angle to all other already found C5s for each group
3038  groupMatched = false;
3039  for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( C5Possibilities.size() ); gIt++ )
3040  {
3041  if ( ProSHADE_internal_symmetry::testGroupAgainstSymmetry ( CSymList, &C5Possibilities.at(gIt), CSymList->at(cIt), axErr, 1.0/2.0, true, cIt ) ) { ProSHADE_internal_misc::addToUnsignVector ( &C5Possibilities.at(gIt), cIt ); groupMatched = true; break; }
3042  }
3043 
3044  //============================================ If no group matched, create a new group
3045  if ( !groupMatched ) { C5PossibilitiesHlp.clear(); ProSHADE_internal_misc::addToUnsignVector ( &C5PossibilitiesHlp, cIt ); ProSHADE_internal_misc::addToUnsignVectorVector ( &C5Possibilities, C5PossibilitiesHlp ); continue; }
3046  }
3047 
3048  //================================================ Test for missing symmetry axes, if need be
3049  ProSHADE_internal_symmetry::findMissingAxes ( &C5Possibilities, CSymList, 6, axErr, 1.0 / 2.0, 5, dataObj, minPeakHeight );
3050 
3051  //=================================================Any group has 6 entries? If more such groups, take the one with highest average height.
3052  proshade_double maxHeight = 0.0; proshade_unsign maxGrp = 0;
3053  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( C5Possibilities.size() ); iter++ ) { if ( C5Possibilities.at(iter).size() == 6 ) { if ( ( ( CSymList->at(C5Possibilities.at(iter).at(0))[5] + CSymList->at(C5Possibilities.at(iter).at(1))[5] + CSymList->at(C5Possibilities.at(iter).at(2))[5] + CSymList->at(C5Possibilities.at(iter).at(3))[5] + CSymList->at(C5Possibilities.at(iter).at(4))[5] + CSymList->at(C5Possibilities.at(iter).at(5))[5] ) / 6.0 ) > maxHeight ) { maxHeight = ( ( CSymList->at(C5Possibilities.at(iter).at(0))[5] + CSymList->at(C5Possibilities.at(iter).at(1))[5] + CSymList->at(C5Possibilities.at(iter).at(2))[5] + CSymList->at(C5Possibilities.at(iter).at(3))[5] + CSymList->at(C5Possibilities.at(iter).at(4))[5] + CSymList->at(C5Possibilities.at(iter).at(5))[5] ) / 6.0 ); maxGrp = iter; } } }
3054 
3055  if ( C5Possibilities.at(maxGrp).size() == 6 )
3056  {
3057  //============================================ Success! Save and exit
3058  for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( C5Possibilities.at(maxGrp).size() ); it++ ) { ProSHADE_internal_misc::addToDblPtrVector ( ret, CSymList->at(C5Possibilities.at(maxGrp).at(it)) ); }
3059 
3060  //============================================ Report progress
3061  ProSHADE_internal_messages::printProgressMessage ( verbose, 3, "Detection of six C5 axes successfull." );
3062 
3063  //============================================ Done
3064  return ;
3065  }
3066 
3067  //================================================ Done
3068  return ;
3069 
3070 }
3071 
3079 std::vector < std::pair< proshade_unsign, proshade_unsign > > findBestIcosDihedralPair ( std::vector< proshade_double* >* CSymList, proshade_double minPeakHeight, proshade_double axErr )
3080 {
3081  //================================================ Initialise variables
3082  std::vector < std::pair< proshade_unsign, proshade_unsign > > ret;
3083  std::vector< proshade_unsign > C5List;
3084  proshade_double dotProduct;
3085 
3086  //================================================ Find all C5 symmetries
3087  for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ ) { const FloatingPoint< proshade_double > lhs1 ( CSymList->at(cSym)[0] ), rhs1 ( 5.0 ); if ( lhs1.AlmostEquals ( rhs1 ) && CSymList->at(cSym)[5] >= minPeakHeight ) { ProSHADE_internal_misc::addToUnsignVector ( &C5List, cSym ); } }
3088 
3089  //================================================ For each unique pair of C5 and C3
3090  for ( proshade_unsign c5 = 0; c5 < static_cast<proshade_unsign> ( C5List.size() ); c5++ )
3091  {
3092  for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
3093  {
3094  //======================================== Compare only C3s to the C5List and only with decent average peak height
3095  const FloatingPoint< proshade_double > lhs1 ( CSymList->at(cSym)[0] ), rhs1 ( 3.0 );
3096  if ( !lhs1.AlmostEquals ( rhs1 ) ) { continue; }
3097  if ( CSymList->at(cSym)[5] < minPeakHeight ) { continue; }
3098 
3099  //======================================== Check the angle between the C5 and C3 axes
3100  dotProduct = ProSHADE_internal_maths::computeDotProduct ( &CSymList->at(C5List.at(c5))[1],
3101  &CSymList->at(C5List.at(c5))[2],
3102  &CSymList->at(C5List.at(c5))[3],
3103  &CSymList->at(cSym)[1],
3104  &CSymList->at(cSym)[2],
3105  &CSymList->at(cSym)[3] );
3106 
3107  //======================================== Is the angle approximately the dihedral angle?
3108  if ( std::abs ( std::abs( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) - std::abs( dotProduct ) ) < axErr )
3109  {
3110  std::pair< proshade_unsign, proshade_unsign > hlp;
3111  hlp.first = C5List.at(c5);
3112  hlp.second = cSym;
3113  ret.emplace_back ( hlp );
3114  }
3115  }
3116  }
3117 
3118  //================================================ Done
3119  return ( ret );
3120 }
3121 
3142 void ProSHADE_internal_symmetry::predictIcosAxes ( std::vector< proshade_double* >* CSymList, std::vector< std::vector< proshade_double* > >* ret, proshade_double axErr, proshade_double minPeakHeight )
3143 {
3144  //================================================ Find the best axis combination with dihedral angle and correct folds
3145  std::vector < std::pair< proshade_unsign, proshade_unsign > > initAxes = findBestIcosDihedralPair ( CSymList, minPeakHeight, axErr );
3146 
3147  //================================================ For each pair of possible axis combinations
3148  for ( size_t pIt = 0; pIt < initAxes.size(); pIt++ )
3149  {
3150  //============================================ Create the tetrahedronAxes object
3152 
3153  //============================================ Find rotation between the detected C5 and the model C5 axes.
3154  proshade_double* rotMat = ProSHADE_internal_maths::findRotMatMatchingVectors ( icoAx->getValue ( 0, 1 ),
3155  icoAx->getValue ( 0, 2 ),
3156  icoAx->getValue ( 0, 3 ),
3157  CSymList->at(initAxes.at(pIt).first)[1],
3158  CSymList->at(initAxes.at(pIt).first)[2],
3159  CSymList->at(initAxes.at(pIt).first)[3] );
3160 
3161  //============================================ Rotate the model C3 to the correct orientation relative to the detected C5 axis.
3162  proshade_double* rotModelC3 = ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( rotMat,
3163  icoAx->getValue ( 6, 1 ),
3164  icoAx->getValue ( 6, 2 ),
3165  icoAx->getValue ( 6, 3 ) );
3166 
3167  //============================================ Find the angle betwen the rotated model C3 and the detected C3 axes along the detected C5 axis.
3168  proshade_double bestAng = 0.0, curAngDist, bestAngDist = 999.9;
3169  proshade_double* rotMatHlp = new proshade_double[9];
3170  ProSHADE_internal_misc::checkMemoryAllocation ( rotMatHlp, __FILE__, __LINE__, __func__ );
3171  for ( proshade_double ang = 0.0; ang < ( M_PI * 2.0 ); ang += 0.002 )
3172  {
3173  //============================================ Compute rotation matrix for this angle value
3174  ProSHADE_internal_maths::getRotationMatrixFromAngleAxis ( rotMatHlp, CSymList->at(initAxes.at(pIt).first)[1], CSymList->at(initAxes.at(pIt).first)[2], CSymList->at(initAxes.at(pIt).first)[3], ang );
3175 
3176  //======================================== Rotate the rotated C2 by the matrix
3177  proshade_double* rotRotModelC3 = ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( rotMatHlp,
3178  rotModelC3[0],
3179  rotModelC3[1],
3180  rotModelC3[2] );
3181 
3182  //======================================== Find distance
3183  curAngDist = std::sqrt ( std::pow ( rotRotModelC3[0] - CSymList->at(initAxes.at(pIt).second)[1], 2.0 ) +
3184  std::pow ( rotRotModelC3[1] - CSymList->at(initAxes.at(pIt).second)[2], 2.0 ) +
3185  std::pow ( rotRotModelC3[2] - CSymList->at(initAxes.at(pIt).second)[3], 2.0 ) );
3186 
3187  //======================================== Save best angle
3188  if ( curAngDist < bestAngDist ) { bestAngDist = curAngDist; bestAng = ang; }
3189 
3190  //======================================== Release memory
3191  delete[] rotRotModelC3;
3192  }
3193 
3194  //============================================ Release memory
3195  delete[] rotMatHlp;
3196 
3197  //============================================ For the rotation matrix along the detected C5 axis with the same anlge as is between the rotated model C3 and the detected C3 axes.
3198  proshade_double* rotMat2 = new proshade_double[9];
3199  ProSHADE_internal_misc::checkMemoryAllocation ( rotMat2, __FILE__, __LINE__, __func__ );
3200  ProSHADE_internal_maths::getRotationMatrixFromAngleAxis ( rotMat2, CSymList->at(initAxes.at(pIt).first)[1], CSymList->at(initAxes.at(pIt).first)[2], CSymList->at(initAxes.at(pIt).first)[3], bestAng );
3201 
3202  //============================================ Combine the two rotation matrices into a single rotation matrix
3203  proshade_double* rotMatFin = ProSHADE_internal_maths::compute3x3MatrixMultiplication ( rotMat2, rotMat );
3204 
3205  //============================================ For each model axis
3206  std::vector< proshade_double* > hlpAxes;
3207  for ( proshade_unsign iter = 0; iter < icoAx->getNoAxes ( ); iter++ )
3208  {
3209  //======================================== Rotate the model axis to fit the detected orientation
3210  proshade_double* rotAxis = ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( rotMatFin,
3211  icoAx->getValue ( iter, 1 ),
3212  icoAx->getValue ( iter, 2 ),
3213  icoAx->getValue ( iter, 3 ) );
3214 
3215  //======================================== Create ProSHADE symmetry axis representation
3216  proshade_double* axis = new proshade_double[7];
3217  ProSHADE_internal_misc::checkMemoryAllocation ( axis, __FILE__, __LINE__, __func__ );
3218 
3219  axis[0] = icoAx->getValue ( iter, 0 );
3220  axis[1] = rotAxis[0];
3221  axis[2] = rotAxis[1];
3222  axis[3] = rotAxis[2];
3223  axis[4] = ( 2.0 * M_PI ) / axis[0];
3224  axis[5] = 0.0;
3225  axis[6] = -1.0;
3226 
3227  //======================================== Save axis to ret
3229 
3230  //======================================== Release memory
3231  delete[] rotAxis;
3232  delete[] axis;
3233  }
3234 
3235  //============================================ Save to ret
3236  ret->emplace_back ( hlpAxes );
3237 
3238  //============================================ Release memory
3239  delete[] rotMat;
3240  delete[] rotMat2;
3241  delete[] rotMatFin;
3242  delete[] rotModelC3;
3243  delete icoAx;
3244  }
3245 
3246  //================================================ Done
3247  return ;
3248 
3249 }
3250 
3258 std::pair< proshade_unsign, proshade_unsign > findBestOctaDihedralPair ( std::vector< proshade_double* >* CSymList, proshade_double minPeakHeight, proshade_double axErr )
3259 {
3260  //================================================ Initialise variables
3261  std::pair< proshade_unsign, proshade_unsign > ret;
3262  std::vector< proshade_unsign > C4List;
3263  proshade_double bestHeightSum = 0.0;
3264  proshade_double dotProduct;
3265 
3266  //================================================ Find all C5 symmetries
3267  for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ ) { const FloatingPoint< proshade_double > lhs1 ( CSymList->at(cSym)[0] ), rhs1 ( 4.0 ); if ( lhs1.AlmostEquals ( rhs1 ) && CSymList->at(cSym)[5] >= minPeakHeight ) { ProSHADE_internal_misc::addToUnsignVector ( &C4List, cSym ); } }
3268 
3269  //================================================ For each unique pair of C5 and C3
3270  for ( proshade_unsign c4 = 0; c4 < static_cast<proshade_unsign> ( C4List.size() ); c4++ )
3271  {
3272  for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
3273  {
3274  //======================================== Compare only C3s to the C5List and only with decent average peak height
3275  const FloatingPoint< proshade_double > lhs1 ( CSymList->at(cSym)[0] ), rhs1 ( 3.0 );
3276  if ( !lhs1.AlmostEquals ( rhs1 ) ) { continue; }
3277  if ( CSymList->at(cSym)[5] < minPeakHeight ) { continue; }
3278 
3279  //======================================== Check the angle between the C5 and C3 axes
3280  dotProduct = ProSHADE_internal_maths::computeDotProduct ( &CSymList->at(C4List.at(c4))[1],
3281  &CSymList->at(C4List.at(c4))[2],
3282  &CSymList->at(C4List.at(c4))[3],
3283  &CSymList->at(cSym)[1],
3284  &CSymList->at(cSym)[2],
3285  &CSymList->at(cSym)[3] );
3286 
3287  //======================================== Is the angle approximately the dihedral angle?
3288  if ( ( ( 1.0 / sqrt ( 3.0 ) ) > ( std::abs( dotProduct ) - axErr ) ) && ( ( 1.0 / sqrt ( 3.0 ) ) < ( std::abs( dotProduct ) + axErr ) ) )
3289  {
3290  if ( bestHeightSum < ( CSymList->at(C4List.at(c4))[5] + CSymList->at(cSym)[5] ) )
3291  {
3292  bestHeightSum = ( CSymList->at(C4List.at(c4))[5] + CSymList->at(cSym)[5] );
3293  ret.first = C4List.at(c4);
3294  ret.second = cSym;
3295  }
3296  }
3297  }
3298  }
3299 
3300  //================================================ Done
3301  return ( ret );
3302 
3303 }
3304 
3325 void ProSHADE_internal_symmetry::predictOctaAxes ( std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, proshade_double axErr, proshade_double minPeakHeight )
3326 {
3327  //================================================ Create the tetrahedronAxes object
3329 
3330  //================================================ Find the best axis combination with dihedral angle and correct folds
3331  std::pair< proshade_unsign, proshade_unsign > initAxes = findBestOctaDihedralPair ( CSymList, minPeakHeight, axErr );
3332 
3333  //================================================ Find rotation between the detected C4 and the model C4 axes.
3334  proshade_double* rotMat = ProSHADE_internal_maths::findRotMatMatchingVectors ( octAx->getValue ( 0, 1 ),
3335  octAx->getValue ( 0, 2 ),
3336  octAx->getValue ( 0, 3 ),
3337  CSymList->at(initAxes.first)[1],
3338  CSymList->at(initAxes.first)[2],
3339  CSymList->at(initAxes.first)[3] );
3340 
3341  //================================================ Rotate the model C3 to the correct orientation relative to the detected C4 axis.
3342  proshade_double* rotModelC3 = ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( rotMat,
3343  octAx->getValue ( 3, 1 ),
3344  octAx->getValue ( 3, 2 ),
3345  octAx->getValue ( 3, 3 ) );
3346 
3347  //================================================ Find the angle betwen the rotated model C3 and the detected C3 axes along the detected C4 axis.
3348  proshade_double bestAng = 0.0, curAngDist, bestAngDist = 999.9;
3349  proshade_double* rotMatHlp = new proshade_double[9];
3350  ProSHADE_internal_misc::checkMemoryAllocation ( rotMatHlp, __FILE__, __LINE__, __func__ );
3351  for ( proshade_double ang = 0.0; ang < ( M_PI * 2.0 ); ang += 0.002 )
3352  {
3353  //============================================ Compute rotation matrix for this angle value
3354  ProSHADE_internal_maths::getRotationMatrixFromAngleAxis ( rotMatHlp, CSymList->at(initAxes.first)[1], CSymList->at(initAxes.first)[2], CSymList->at(initAxes.first)[3], ang );
3355 
3356  //============================================ Rotate the rotated C2 by the matrix
3357  proshade_double* rotRotModelC3 = ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( rotMatHlp,
3358  rotModelC3[0],
3359  rotModelC3[1],
3360  rotModelC3[2] );
3361 
3362  //============================================ Find distance
3363  curAngDist = std::sqrt ( std::pow ( rotRotModelC3[0] - CSymList->at(initAxes.second)[1], 2.0 ) +
3364  std::pow ( rotRotModelC3[1] - CSymList->at(initAxes.second)[2], 2.0 ) +
3365  std::pow ( rotRotModelC3[2] - CSymList->at(initAxes.second)[3], 2.0 ) );
3366 
3367  //============================================ Save best angle
3368  if ( curAngDist < bestAngDist ) { bestAngDist = curAngDist; bestAng = ang; }
3369 
3370  //============================================ Release memory
3371  delete[] rotRotModelC3;
3372  }
3373 
3374  //============================================ Release memory
3375  delete[] rotMatHlp;
3376 
3377  //================================================ For the rotation matrix along the detected C5 axis with the same anlge as is between the rotated model C3 and the detected C3 axes.
3378  proshade_double* rotMat2 = new proshade_double[9];
3379  ProSHADE_internal_misc::checkMemoryAllocation ( rotMat2, __FILE__, __LINE__, __func__ );
3380  ProSHADE_internal_maths::getRotationMatrixFromAngleAxis ( rotMat2, CSymList->at(initAxes.first)[1], CSymList->at(initAxes.first)[2], CSymList->at(initAxes.first)[3], bestAng );
3381 
3382  //================================================ Combine the two rotation matrices into a single rotation matrix
3383  proshade_double* rotMatFin = ProSHADE_internal_maths::compute3x3MatrixMultiplication ( rotMat2, rotMat );
3384 
3385  //================================================ For each model axis
3386  for ( proshade_unsign iter = 0; iter < octAx->getNoAxes ( ); iter++ )
3387  {
3388  //============================================ Rotate the model axis to fit the detected orientation
3389  proshade_double* rotAxis = ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( rotMatFin,
3390  octAx->getValue ( iter, 1 ),
3391  octAx->getValue ( iter, 2 ),
3392  octAx->getValue ( iter, 3 ) );
3393 
3394  //============================================ Create ProSHADE symmetry axis representation
3395  proshade_double* axis = new proshade_double[7];
3396  ProSHADE_internal_misc::checkMemoryAllocation ( axis, __FILE__, __LINE__, __func__ );
3397 
3398  axis[0] = octAx->getValue ( iter, 0 );
3399  axis[1] = rotAxis[0];
3400  axis[2] = rotAxis[1];
3401  axis[3] = rotAxis[2];
3402  axis[4] = ( 2.0 * M_PI ) / axis[0];
3403  axis[5] = 0.0;
3404  axis[6] = -1.0;
3405 
3406  //============================================ Save axis to ret
3408 
3409  //============================================ Release memory
3410  delete[] rotAxis;
3411  }
3412 
3413  //================================================ Release memory
3414  delete[] rotMat;
3415  delete[] rotMat2;
3416  delete[] rotMatFin;
3417  delete[] rotModelC3;
3418  delete octAx;
3419 
3420  //================================================ Done
3421  return ;
3422 
3423 }
3424 
3438 void ProSHADE_internal_symmetry::findIcos10C3s ( std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_signed verbose, proshade_double minPeakHeight )
3439 {
3440  //================================================ Initialise variables
3441  std::vector< proshade_unsign > prospectiveC3s, retGrp;
3442  proshade_double dotProd;
3443  proshade_unsign noClose, noAway;
3444 
3445  //================================================ Report progress
3446  ProSHADE_internal_messages::printProgressMessage ( verbose, 2, "Starting detection of ten C3 axes." );
3447 
3448  //================================================ For each C3
3449  for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
3450  {
3451  //============================================ Use only C3s with hight enough average
3452  if ( CSymList->at(cIt)[0] != 3.0 || CSymList->at(cIt)[0] < minPeakHeight ) { continue; }
3453 
3454  //============================================ Check the C3 has acos ( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) to 3 C5s and acos ( 1.0 - ( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) ) to the other three C5s
3455  noClose = 0; noAway = 0;
3456  for ( proshade_unsign rIt = 0; rIt < 6; rIt++ )
3457  {
3458  dotProd = ProSHADE_internal_maths::computeDotProduct ( &ret->at(rIt)[1],
3459  &ret->at(rIt)[2],
3460  &ret->at(rIt)[3],
3461  &CSymList->at(cIt)[1],
3462  &CSymList->at(cIt)[2],
3463  &CSymList->at(cIt)[3] );
3464 
3465  if ( ( std::abs ( dotProd ) > ( ( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) - axErr ) ) && ( std::abs ( dotProd ) < ( ( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) + axErr ) ) ) { noClose += 1; continue; }
3466  if ( ( std::abs ( dotProd ) > ( 1.0 - ( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) - axErr ) ) && ( std::abs ( dotProd ) < ( 1.0 - ( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) + axErr ) ) ) { noAway += 1; continue; }
3467  }
3468 
3469  //============================================ If correct angles distribution is found, save the axis
3470  if ( ( noClose == 3 ) && ( noAway == 3 ) )
3471  {
3472  ProSHADE_internal_misc::addToUnsignVector ( &prospectiveC3s, cIt );
3473  }
3474  }
3475 
3476  //================================================ Search for missing axes
3477  for ( proshade_unsign iter = 0; iter < 6; iter++ ) { ProSHADE_internal_misc::addToUnsignVector ( &retGrp, iter ); }
3478  if ( !ProSHADE_internal_symmetry::findMissingAxesDual ( &prospectiveC3s, CSymList, ret, &retGrp, 10, axErr, 3, std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ), 3, 1.0 - ( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ), 3, dataObj ) )
3479  {
3480  return ;
3481  }
3482 
3483  //================================================ Found correct number of axes! Now save the
3484  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( prospectiveC3s.size() ); iter++ )
3485  {
3486  ProSHADE_internal_misc::addToDblPtrVector ( ret, CSymList->at(prospectiveC3s.at(iter)) );
3487  }
3488 
3489  //================================================ Report progress
3490  ProSHADE_internal_messages::printProgressMessage ( verbose, 3, "Detection of ten C3 axes successfull." );
3491 
3492  //================================================ Done
3493  return ;
3494 
3495 }
3496 
3510 void ProSHADE_internal_symmetry::findIcos15C2s ( std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data* dataObj, proshade_signed verbose, proshade_double minPeakHeight )
3511 {
3512  //================================================ Initialise variables
3513  std::vector< proshade_unsign > prospectiveC2s, retGrp;
3514  proshade_double dotProd;
3515  proshade_unsign noClose, noMidway, noAway;
3516 
3517  //================================================ Report progress
3518  ProSHADE_internal_messages::printProgressMessage ( verbose, 2, "Starting detection of fifteen C2 axes." );
3519 
3520  //================================================ For each C2
3521  for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
3522  {
3523  //============================================ Use only C2s
3524  const FloatingPoint< proshade_double > lhs999 ( CSymList->at(cIt)[5] ), rhs999 ( static_cast< proshade_double > ( -999.9 ) );
3525  if ( CSymList->at(cIt)[0] != 2.0 || ( ( CSymList->at(cIt)[5] < minPeakHeight ) && !( lhs999.AlmostEquals( rhs999 ) ) ) ) { continue; }
3526 
3527  //============================================ Check the C2 has acos ( 0.0 ) to 2 C5s, acos ( 0.5 ) to another 2 C5s and acos ( sqrt ( 3.0 ) / 2.0 ) to the last two C5s
3528  noClose = 0; noMidway = 0; noAway = 0;
3529  for ( proshade_unsign rIt = 0; rIt < 6; rIt++ )
3530  {
3531  dotProd = ProSHADE_internal_maths::computeDotProduct ( &ret->at(rIt)[1],
3532  &ret->at(rIt)[2],
3533  &ret->at(rIt)[3],
3534  &CSymList->at(cIt)[1],
3535  &CSymList->at(cIt)[2],
3536  &CSymList->at(cIt)[3] );
3537 
3538  if ( ( std::abs ( dotProd ) > ( ( sqrt ( 3.0 ) / 2.0 ) - axErr ) ) && ( std::abs ( dotProd ) < ( ( sqrt ( 3.0 ) / 2.0 ) + axErr ) ) ) { noAway += 1; continue; }
3539  if ( ( std::abs ( dotProd ) > ( ( 1.0 / 2.0 ) - axErr ) ) && ( std::abs ( dotProd ) < ( ( 1.0 / 2.0 ) + axErr ) ) ) { noMidway += 1; continue; }
3540  if ( ( std::abs ( dotProd ) > ( ( 0.0 ) - axErr ) ) && ( std::abs ( dotProd ) < ( ( 0.0 ) + axErr ) ) ) { noClose += 1; continue; }
3541  }
3542 
3543  //============================================ If correct angles distribution is found, save the axis
3544  if ( ( noClose == 2 ) && ( noMidway == 2 ) && ( noAway == 2 ) )
3545  {
3546  ProSHADE_internal_misc::addToUnsignVector ( &prospectiveC2s, cIt );
3547  }
3548  }
3549 
3550  //================================================ Search for missing axes
3551  for ( proshade_unsign iter = 0; iter < 6; iter++ ) { ProSHADE_internal_misc::addToUnsignVector ( &retGrp, iter ); }
3552  if ( !ProSHADE_internal_symmetry::findMissingAxesTriple ( &prospectiveC2s, CSymList, ret, &retGrp, 15, axErr, 2, 0.0, 2, 1.0/2.0, 2, sqrt ( 3.0 ) / 2.0, 2, dataObj ) )
3553  {
3554  return ;
3555  }
3556 
3557  //================================================ Found correct number of axes! Now save the
3558  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( prospectiveC2s.size() ); iter++ )
3559  {
3560  ProSHADE_internal_misc::addToDblPtrVector ( ret, CSymList->at(prospectiveC2s.at(iter)) );
3561  }
3562 
3563  //================================================ Report progress
3564  ProSHADE_internal_messages::printProgressMessage ( verbose, 3, "Detection of fifteen C2 axes successfull." );
3565 
3566  //================================================ Done
3567  return ;
3568 
3569 }
3570 
3593 bool ProSHADE_internal_symmetry::findMissingAxesTriple ( std::vector< proshade_unsign >* possibilities, std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, std::vector< proshade_unsign >* retGroup, proshade_unsign requiredNoAxes, proshade_double axErr, proshade_unsign noMatchesG1, proshade_double angle1, proshade_unsign noMatchesG2, proshade_double angle2, proshade_unsign noMatchesG3, proshade_double angle3, proshade_unsign fold, ProSHADE_internal_data::ProSHADE_data* dataObj )
3594 {
3595  //================================================ Initialise variables
3596  bool atLeastOne = false;
3597  std::vector< proshade_double* > prosp;
3598  std::vector< proshade_double > sol;
3599 
3600  //================================================ Proceed only if need be
3601  if ( static_cast<proshade_unsign> ( possibilities->size() ) == requiredNoAxes ) { atLeastOne = true; return ( atLeastOne ); }
3602 
3603  //================================================ Copy already found to prospective
3604  for ( proshade_unsign prIt = 0; prIt < static_cast<proshade_unsign> ( possibilities->size() ); prIt++ )
3605  {
3606  ProSHADE_internal_symmetry::addAxisUnlessSame ( static_cast< proshade_unsign > ( CSymList->at(possibilities->at(prIt))[0] ),
3607  CSymList->at(possibilities->at(prIt))[1],
3608  CSymList->at(possibilities->at(prIt))[2],
3609  CSymList->at(possibilities->at(prIt))[3],
3610  CSymList->at(possibilities->at(prIt))[5], &prosp, axErr );
3611  }
3612 
3613  //================================================ Start generating possible solutions
3614  for ( proshade_unsign rgIt1 = 0; rgIt1 < static_cast<proshade_unsign> ( retGroup->size() ); rgIt1++ )
3615  {
3616  for ( proshade_unsign rgIt2 = 0; rgIt2 < static_cast<proshade_unsign> ( retGroup->size() ); rgIt2++ )
3617  {
3618  //======================================== Use unique combinations (order matters here!)
3619  if ( rgIt1 == rgIt2 ) { continue; }
3620 
3621  for ( proshade_unsign rgIt3 = 0; rgIt3 < static_cast<proshade_unsign> ( retGroup->size() ); rgIt3++ )
3622  {
3623  //==================================== Use unique combinations (order matters here!)
3624  if ( ( rgIt1 == rgIt3 ) || ( rgIt2 == rgIt3 ) ) { continue; }
3625 
3626  //==================================== Generate possible solution (1)
3627  sol = ProSHADE_internal_maths::findVectorFromThreeVAndThreeD ( ret->at(rgIt1)[1], ret->at(rgIt1)[2], ret->at(rgIt1)[3],
3628  ret->at(rgIt2)[1], ret->at(rgIt2)[2], ret->at(rgIt2)[3],
3629  ret->at(rgIt3)[1], ret->at(rgIt3)[2], ret->at(rgIt3)[3], angle1, angle2, angle3 );
3630 
3631  //==================================== Check if solution fits the group completely
3632  ProSHADE_internal_symmetry::checkFittingAxisTripleAndSave ( retGroup, ret, fold, sol.at(0), sol.at(1), sol.at(2), &prosp, axErr, noMatchesG1, angle1, noMatchesG2, angle2, noMatchesG3, angle3, dataObj );
3633  if ( prosp.size() == requiredNoAxes ) { break; }
3634 
3635  //==================================== Generate possible solution (2)
3636  sol = ProSHADE_internal_maths::findVectorFromThreeVAndThreeD ( ret->at(rgIt1)[1], ret->at(rgIt1)[2], ret->at(rgIt1)[3],
3637  ret->at(rgIt2)[1], ret->at(rgIt2)[2], ret->at(rgIt2)[3],
3638  ret->at(rgIt3)[1], ret->at(rgIt3)[2], ret->at(rgIt3)[3], -angle1, -angle2, -angle3 );
3639 
3640  //==================================== Check if solution fits the group completely
3641  ProSHADE_internal_symmetry::checkFittingAxisTripleAndSave ( retGroup, ret, fold, sol.at(0), sol.at(1), sol.at(2), &prosp, axErr, noMatchesG1, angle1, noMatchesG2, angle2, noMatchesG3, angle3, dataObj );
3642  if ( prosp.size() == requiredNoAxes ) { break; }
3643  }
3644 
3645  if ( prosp.size() == requiredNoAxes ) { break; }
3646  }
3647 
3648  if ( prosp.size() == requiredNoAxes ) { break; }
3649  }
3650 
3651  //================================================ Found all required axes
3652  if ( prosp.size() == requiredNoAxes )
3653  {
3654  //============================================ For each found missing axis
3655  for ( proshade_unsign axIt = static_cast<proshade_unsign> ( possibilities->size() ); axIt < static_cast<proshade_unsign> ( prosp.size() ); axIt++ )
3656  {
3657  if ( ProSHADE_internal_maths::isAxisUnique ( CSymList, prosp.at(axIt), axErr ) )
3658  {
3659  //======================================== Add
3660  ProSHADE_internal_misc::addToDblPtrVector ( CSymList, prosp.at(axIt) );
3661  ProSHADE_internal_misc::addToUnsignVector ( possibilities, static_cast<proshade_unsign> ( CSymList->size()-1 ) );
3662  }
3663  }
3664 
3665  atLeastOne = true;
3666  return ( atLeastOne );
3667  }
3668  else
3669  {
3670  //============================================ Delete all found, but unnecessary axes
3671  for ( proshade_unsign axIt = static_cast<proshade_unsign> ( possibilities->size() ); axIt < static_cast<proshade_unsign> ( prosp.size() ); axIt++ )
3672  {
3673  delete[] prosp.at(axIt);
3674  }
3675  }
3676 
3677  //================================================ Done
3678  return ( atLeastOne );
3679 
3680 }
3681 
3704 void ProSHADE_internal_symmetry::checkFittingAxisTripleAndSave ( std::vector< proshade_unsign >* retGroup, std::vector< proshade_double* >* ret, proshade_unsign fold, proshade_double axX, proshade_double axY, proshade_double axZ, std::vector< proshade_double* >* prosp, proshade_double axErr, proshade_unsign noMatchesG1, proshade_double angle1, proshade_unsign noMatchesG2, proshade_double angle2, proshade_unsign noMatchesG3, proshade_double angle3, ProSHADE_internal_data::ProSHADE_data* dataObj )
3705 {
3706  //================================================ Initialise variables
3707  proshade_unsign noG1 = 0;
3708  proshade_unsign noG2 = 0;
3709  proshade_unsign noG3 = 0;
3710  proshade_double dotProd = 0.0;
3711  proshade_double axHeight = 0.0;
3712 
3713  //================================================ Find the angle and count dual matching frequencies
3714  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( retGroup->size() ); rIt++ )
3715  {
3716  dotProd = ProSHADE_internal_maths::computeDotProduct ( &ret->at(retGroup->at(rIt))[1],
3717  &ret->at(retGroup->at(rIt))[2],
3718  &ret->at(retGroup->at(rIt))[3],
3719  &axX, &axY, &axZ );
3720 
3721  if ( ( std::abs ( dotProd ) > ( angle1 - axErr ) ) && ( std::abs ( dotProd ) < ( angle1 + axErr ) ) ) { noG1 += 1; continue; }
3722  if ( ( std::abs ( dotProd ) > ( angle2 - axErr ) ) && ( std::abs ( dotProd ) < ( angle2 + axErr ) ) ) { noG2 += 1; continue; }
3723  if ( ( std::abs ( dotProd ) > ( angle3 - axErr ) ) && ( std::abs ( dotProd ) < ( angle3 + axErr ) ) ) { noG3 += 1; continue; }
3724  }
3725 
3726  //================================================ If correct frequencies are matched, check height.
3727  if ( ( noG1 == noMatchesG1 ) && ( noG2 == noMatchesG2 ) && ( noG3 == noMatchesG3 ) )
3728  {
3729  //============================================ Is the height good enough?
3730  axHeight = ProSHADE_internal_symmetry::missingAxisHeight ( axX, axY, axZ, dataObj, fold, axErr );
3731 
3732  //============================================ If so, save
3733  if ( axHeight > 0.1 )
3734  {
3735  ProSHADE_internal_symmetry::addAxisUnlessSame ( fold, axX, axY, axZ, axHeight, prosp, axErr );
3736  }
3737  }
3738 
3739  //================================================ Done
3740  return ;
3741 
3742 }
3743 
3757 {
3758  //================================================ Initialise variables
3759  std::vector< proshade_unsign > primes = ProSHADE_internal_maths::findAllPrimes ( settings->maxSymmetryFold );
3760  std::vector< proshade_double* > ret, tmpHolder;
3761  std::vector< proshade_unsign > testedFolds;
3762  proshade_double symThres;
3763  proshade_unsign foldToTest;
3764  bool foldDone, anyNewSyms = true;
3765 
3766  //================================================ For each found prime number in the limit
3767  for ( proshade_unsign prIt = 0; prIt < static_cast< proshade_unsign > ( primes.size() ); prIt++ )
3768  {
3769  //============================================ Report progress
3770  std::stringstream hlpSS;
3771  hlpSS << "Searching for prime fold symmetry C" << primes.at(prIt) << ".";
3772  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 3, hlpSS.str() );
3773 
3774  //============================================ Get all symmetries for this prime fold
3775  std::vector< proshade_double* > prSyms = this->findRequestedCSymmetryFromAngleAxis ( settings, primes.at(prIt), &symThres );
3776 
3777  //============================================ Save the detected C symmetries
3778  for ( size_t axIt = 0; axIt < prSyms.size(); axIt++ )
3779  {
3780  //======================================== Is this symmetry passing the threshold?
3781  if ( prSyms.at(axIt)[5] >= symThres )
3782  {
3783  //==================================== Add this symmetry to final list
3784  if ( ProSHADE_internal_maths::isAxisUnique ( &ret, prSyms.at(axIt), settings->axisErrTolerance, true ) )
3785  {
3786  ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( &ret, prSyms.at(axIt) );
3787  }
3788  }
3789 
3790  //======================================== Release memory
3791  delete[] prSyms.at(axIt);
3792  }
3793  }
3794 
3795  //================================================ Was anything found?
3796  if ( ret.size() < 1 ) { return ( ret ); }
3797 
3798  //================================================ Check for prime symmetry fold multiples
3799  while ( anyNewSyms )
3800  {
3801  //============================================ Initialise new iteration
3802  anyNewSyms = false;
3803 
3804  //============================================ For each passing symmetry, look if there are any combinations of symmetries that would contain it
3805  for ( proshade_unsign axIt1 = 0; axIt1 < static_cast< proshade_unsign > ( ret.size() ); axIt1++ )
3806  {
3807  for ( proshade_unsign axIt2 = 0; axIt2 < static_cast< proshade_unsign > ( ret.size() ); axIt2++ )
3808  {
3809  //==================================== Initialise iteration
3810  foldToTest = static_cast< proshade_unsign > ( ret.at(axIt1)[0] * ret.at(axIt2)[0] );
3811  if ( foldToTest > settings->maxSymmetryFold ) { continue; }
3812 
3813  //==================================== Was this fold tested already?
3814  foldDone = false;
3815  for ( proshade_unsign fIt = 0; fIt < static_cast< proshade_unsign > ( testedFolds.size() ); fIt++ ) { if ( testedFolds.at(fIt) == foldToTest ) { foldDone = true; break; } }
3816  if ( foldDone ) { continue; }
3817  else { ProSHADE_internal_misc::addToUnsignVector ( &testedFolds, foldToTest ); }
3818 
3819  //==================================== Report progress
3820  std::stringstream hlpSS2;
3821  hlpSS2 << "Searching for fold combination of detected folds " << ret.at(axIt1)[0] << " and " << ret.at(axIt2)[0] << ", i.e. C" << foldToTest << ".";
3822  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 3, hlpSS2.str() );
3823 
3824  //==================================== Get all symmetries for this fold
3825  std::vector< proshade_double* > prSyms = this->findRequestedCSymmetryFromAngleAxis ( settings, foldToTest, &symThres );
3826 
3827  //==================================== For each detected group with the required fold
3828  for ( size_t newAxIt = 0; newAxIt < prSyms.size(); newAxIt++ )
3829  {
3830  if ( prSyms.at(newAxIt)[5] >= symThres )
3831  {
3832  //================================ Add to detected axes
3833  if ( ProSHADE_internal_maths::isAxisUnique ( &ret, prSyms.at(newAxIt), settings->axisErrTolerance, true ) )
3834  {
3835  ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( &tmpHolder, prSyms.at(newAxIt) );
3836  }
3837  }
3838 
3839  //==================================== Release memory
3840  delete[] prSyms.at(newAxIt);
3841  }
3842  }
3843  }
3844 
3845  //============================================ Add newly found groups and repeat if need be
3846  if ( tmpHolder.size() > 0 )
3847  {
3848  for ( proshade_unsign tmpIt = 0; tmpIt < static_cast< proshade_unsign > ( tmpHolder.size() ); tmpIt++ )
3849  {
3850  ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( &ret, tmpHolder.at(tmpIt) );
3851  delete[] tmpHolder.at(tmpIt);
3852  }
3853 
3854  anyNewSyms = true;
3855  tmpHolder.clear ( );
3856  }
3857  }
3858 
3859  //================================================ Sort the vector
3860  std::sort ( ret.begin(), ret.end(), ProSHADE_internal_misc::sortSymHlpInv );
3861 
3862  //================================================ Done
3863  return ( ret );
3864 }
3865 
3872 bool sortProSHADESymmetryByPeak ( proshade_double* a, proshade_double* b)
3873 {
3874  //================================================ Done
3875  return ( a[5] > b[5] );
3876 
3877 }
3878 
3895 std::vector < proshade_double* > ProSHADE_internal_data::ProSHADE_data::findRequestedCSymmetryFromAngleAxis ( ProSHADE_settings* settings, proshade_unsign fold, proshade_double* peakThres )
3896 {
3897  //================================================ Initialise variables
3898  proshade_double soughtAngle;
3899  std::vector< proshade_double > allPeakHeights;
3900  std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_spherePeakGroup* > peakGroups;
3901  std::vector< proshade_double* > ret;
3902  bool newPeak;
3903 
3904  //================================================ Make sure we have a clean start
3905  this->sphereMappedRotFun.clear();
3906 
3907  //================================================ Convert rotation function to only the required angle-axis space spheres and find all peaks
3908  for ( proshade_double angIt = 1.0; angIt < static_cast<proshade_double> ( fold ); angIt += 1.0 )
3909  {
3910  //============================================ Figure the angles to form the symmetry
3911  soughtAngle = angIt * ( 2.0 * M_PI / static_cast<proshade_double> ( fold ) );
3912 
3913  //============================================ Create the angle-axis sphere with correct radius (angle)
3914  this->sphereMappedRotFun.emplace_back ( new ProSHADE_internal_spheres::ProSHADE_rotFun_sphere ( soughtAngle,
3915  M_PI / static_cast< proshade_double > ( this->maxShellBand ),
3916  this->maxShellBand * 2,
3917  soughtAngle,
3918  static_cast<proshade_unsign> ( angIt - 1.0 ) ) );
3919 
3920  //=========================================== Interpolate rotation function onto the sphere
3921  this->sphereMappedRotFun.at(static_cast<proshade_unsign> ( angIt - 1.0 ))->interpolateSphereValues ( this->getInvSO3Coeffs ( ) );
3922 
3923  //============================================ Find all peaks for this sphere
3924  this->sphereMappedRotFun.at(static_cast<proshade_unsign> ( angIt - 1.0 ))->findAllPeaks ( static_cast< proshade_signed > ( settings->peakNeighbours ), &allPeakHeights );
3925  }
3926 
3927  //============================================ Report progress
3928  std::stringstream hlpSS;
3929  hlpSS << "Found a total of " << std::pow ( static_cast< proshade_double > ( this->maxShellBand ) * 2.0 * ( static_cast< proshade_double > ( fold ) - 1.0 ), 2.0 ) - static_cast< proshade_double > ( allPeakHeights.size() ) << " non-peaks for thresholding.";
3930  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 4, hlpSS.str() );
3931 
3932  //================================================ Determine the threshold for significant peaks
3933  *peakThres = std::max ( settings->minSymPeak, determinePeakThreshold ( allPeakHeights, settings->noIQRsFromMedianNaivePeak ) );
3934 
3935  //============================================ Report progress
3936  std::stringstream hlpSS2;
3937  hlpSS2 << "Determined peak threshold " << *peakThres << ".";
3938  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 4, hlpSS2.str() );
3939 
3940  //================================================ Remove small peaks
3941  for ( proshade_unsign shIt = 0; shIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.size() ); shIt++ )
3942  {
3943  this->sphereMappedRotFun.at(shIt)->removeSmallPeaks ( *peakThres );
3944  }
3945 
3946  //================================================ Group peaks
3947  for ( proshade_unsign sphIt = 0; sphIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.size() ); sphIt++ )
3948  {
3949  //============================================ For each peak
3950  for ( proshade_unsign pkIt = 0; pkIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.at(sphIt)->getPeaks().size() ); pkIt++ )
3951  {
3952  //======================================== Check if peak belongs to an already detected peak group
3953  newPeak = true;
3954  for ( proshade_unsign pkGrpIt = 0; pkGrpIt < static_cast<proshade_unsign> ( peakGroups.size() ); pkGrpIt++ )
3955  {
3956  if ( peakGroups.at(pkGrpIt)->checkIfPeakBelongs ( static_cast< proshade_double > ( this->sphereMappedRotFun.at(sphIt)->getPeaks().at(pkIt).first ),
3957  static_cast< proshade_double > ( this->sphereMappedRotFun.at(sphIt)->getPeaks().at(pkIt).second ),
3958  sphIt, settings->axisErrTolerance, settings->verbose ) ) { newPeak = false; break; }
3959  }
3960 
3961  //======================================== If already added, go to next one
3962  if ( !newPeak ) { continue; }
3963 
3964  //======================================== If not, create a new group with this peak
3965  peakGroups.emplace_back ( new ProSHADE_internal_spheres::ProSHADE_rotFun_spherePeakGroup ( static_cast< proshade_double > ( this->sphereMappedRotFun.at(sphIt)->getPeaks().at(pkIt).first ),
3966  static_cast< proshade_double > ( this->sphereMappedRotFun.at(sphIt)->getPeaks().at(pkIt).second ),
3967  sphIt,
3968  this->sphereMappedRotFun.at(sphIt)->getAngularDim() ) );
3969  }
3970  }
3971 
3972  //================================================ For each peak group, look for the requested fold
3973  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( peakGroups.size() ); grIt++ )
3974  {
3975  //============================================ Report progress
3976  std::stringstream hlpSS3;
3977  hlpSS3 << "Now considering group with LAT " << peakGroups.at(grIt)->getLatFromIndices() << " - " << peakGroups.at(grIt)->getLatToIndices() << " and LON " << peakGroups.at(grIt)->getLonFromIndices() << " - " << peakGroups.at(grIt)->getLonToIndices() << " spanning spheres ";
3978  for ( proshade_unsign sphIt = 0; sphIt < static_cast<proshade_unsign> ( peakGroups.at(grIt)->getSpherePositions().size() ); sphIt++ ) { hlpSS3 << peakGroups.at(grIt)->getSpherePositions().at(sphIt) << " ; "; }
3979  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 5, hlpSS3.str() );
3980 
3981  //============================================ Find point groups in the peak group
3982  peakGroups.at(grIt)->findCyclicPointGroupsGivenFold ( this->sphereMappedRotFun, &ret, settings->useBiCubicInterpolationOnPeaks, fold, settings->verbose );
3983 
3984  //============================================ Release the memory
3985  delete peakGroups.at(grIt);
3986  }
3987 
3988  //================================================ Sort ret by peak height
3989  std::sort ( ret.begin(), ret.end(), sortProSHADESymmetryByPeak );
3990 
3991  //================================================ Done
3992  return ( ret );
3993 
3994 }
3995 
4009 {
4010  //================================================ Initialise variables
4011  std::vector < proshade_unsign > folds;
4012  bool alreadyFound;
4013  proshade_double lat, lon;
4014  proshade_double latSamlUnit = ( 2.0 * M_PI ) / ( static_cast< proshade_double > ( dataObj->maxShellBand ) * 2.0 );
4015  proshade_double lonSamlUnit = ( 1.0 * M_PI ) / ( static_cast< proshade_double > ( dataObj->maxShellBand ) * 2.0 );
4016 
4017  //================================================ Determine all the folds for which rotation function mapping will be required
4018  for ( proshade_unsign iter = 0; iter < static_cast < proshade_unsign > ( ret->size() ); iter++ )
4019  {
4020  alreadyFound = false;
4021  for ( proshade_unsign it = 0; it < static_cast < proshade_unsign > ( folds.size() ); it++ ) { const FloatingPoint< proshade_double > lhs1 ( static_cast< proshade_double > ( folds.at(it) ) ), rhs1 ( ret->at(iter)[0] ); if ( lhs1.AlmostEquals ( rhs1 ) ) { alreadyFound = true; break; } }
4022 
4023  if ( !alreadyFound ) { ProSHADE_internal_misc::addToUnsignVector ( &folds, static_cast< proshade_unsign > ( ret->at(iter)[0] ) ); }
4024  }
4025 
4026  //================================================ For each fold which needs rotation function mapping
4027  for ( proshade_unsign foldIt = 0; foldIt < static_cast < proshade_unsign > ( folds.size() ); foldIt++ )
4028  {
4029  //============================================ Make sure we have a clean start
4030  dataObj->sphereMappedRotFun.clear();
4031 
4032  //============================================ Convert rotation function to only the required angle-axis space spheres and find all peaks
4033  for ( proshade_double angIt = 1.0; angIt < static_cast<proshade_double> ( folds.at(foldIt) ); angIt += 1.0 )
4034  {
4035  //======================================== Create the angle-axis sphere with correct radius (angle)
4036  dataObj->sphereMappedRotFun.emplace_back ( new ProSHADE_internal_spheres::ProSHADE_rotFun_sphere ( angIt * ( 2.0 * M_PI / static_cast<proshade_double> ( folds.at(foldIt) ) ),
4037  M_PI / static_cast < proshade_double > ( folds.at(foldIt) ),
4038  dataObj->maxShellBand * 2,
4039  angIt * ( 2.0 * M_PI / static_cast<proshade_double> ( folds.at(foldIt) ) ),
4040  static_cast<proshade_unsign> ( angIt - 1.0 ) ) );
4041 
4042  //=========================================== Interpolate rotation function onto the sphere
4043  dataObj->sphereMappedRotFun.at( static_cast < proshade_unsign > ( angIt - 1.0 ))->interpolateSphereValues ( dataObj->getInvSO3Coeffs ( ) );
4044  }
4045 
4046  //============================================ For each ret axis with this fold
4047  for ( proshade_unsign axIt = 0; axIt < static_cast< proshade_unsign > ( ret->size() ); axIt++ )
4048  {
4049  //======================================== Ignore different folds
4050  const FloatingPoint< proshade_double > lhs1 ( ret->at(axIt)[0] ), rhs1 ( static_cast< proshade_double > ( folds.at(foldIt) ) );
4051  if ( !lhs1.AlmostEquals ( rhs1 ) ) { continue; }
4052 
4053  //======================================== Convert XYZ to lat and lon INDICES
4054  lat = std::atan2( ret->at(axIt)[2], ret->at(axIt)[1] ) / latSamlUnit;
4055  lon = std::acos ( ret->at(axIt)[3] ) / lonSamlUnit;
4056 
4057  if ( lat < 0.0 ) { lat += ( static_cast< proshade_double > ( dataObj->maxShellBand ) * 2.0 ); }
4058  if ( lon < 0.0 ) { lon += ( static_cast< proshade_double > ( dataObj->maxShellBand ) * 2.0 ); }
4059 
4060  lat = std::round ( lat );
4061  lon = std::round ( lon );
4062 
4063  //======================================== Initialise the peak group
4065 
4066  //======================================== Construct a peak group with entry from each sphere with the axis as the peak
4067  for ( proshade_unsign sphIt = 0; sphIt < static_cast<proshade_unsign> ( dataObj->sphereMappedRotFun.size() ); sphIt++ )
4068  {
4069  if ( sphIt == 0 )
4070  {
4071  //================================ If first sphere, create the peak group
4072  grp = new ProSHADE_internal_spheres::ProSHADE_rotFun_spherePeakGroup ( lat, lon, sphIt, dataObj->sphereMappedRotFun.at(sphIt)->getAngularDim() );
4073  }
4074  else
4075  {
4076  //================================ Add to the existing object
4077  grp->checkIfPeakBelongs ( lat, lon, sphIt, settings->axisErrTolerance, settings->verbose );
4078  }
4079  }
4080 
4081  //======================================== Find the peak height
4082  std::vector < proshade_double* > detectedAxis;
4083  grp->findCyclicPointGroupsGivenFold ( dataObj->sphereMappedRotFun, &detectedAxis, settings->useBiCubicInterpolationOnPeaks, folds.at(foldIt), settings->verbose );
4084 
4085  //======================================== Save it!
4086  if ( detectedAxis.size() > 0 ) { ret->at(axIt)[5] = detectedAxis.at(0)[5]; }
4087  else { ret->at(axIt)[5] = 0.0; }
4088 
4089  //======================================== Release memory
4090  for ( proshade_unsign i = 0; i < static_cast < proshade_unsign > ( detectedAxis.size() ); i++ ) { delete detectedAxis.at(i); }
4091  delete grp;
4092  }
4093  }
4094 
4095  //================================================ Sort axes by fold
4096  std::sort ( ret->begin(), ret->end(), ProSHADE_internal_misc::sortSymInvFoldHlp );
4097 
4098  //================================================ Done
4099  return ;
4100 
4101 }
4102 
4115 proshade_double ProSHADE_internal_symmetry::findPredictedSingleAxisHeight ( proshade_double* axis, proshade_double fold, ProSHADE_internal_data::ProSHADE_data* dataObj, ProSHADE_settings* settings )
4116 {
4117  //================================================ Initialise variables
4118  proshade_double height = 0.0;
4119  proshade_double lat, lon;
4120  proshade_double latSamlUnit = ( 2.0 * M_PI ) / ( static_cast< proshade_double > ( dataObj->maxShellBand ) * 2.0 );
4121  proshade_double lonSamlUnit = ( 1.0 * M_PI ) / ( static_cast< proshade_double > ( dataObj->maxShellBand ) * 2.0 );
4122 
4123  //================================================ Make sure we have a clean start
4124  dataObj->sphereMappedRotFun.clear ( );
4125 
4126  //================================================ Convert rotation function to only the required angle-axis space spheres and find all peaks
4127  for ( proshade_double angIt = 1.0; angIt < fold; angIt += 1.0 )
4128  {
4129  //============================================ Create the angle-axis sphere with correct radius (angle)
4130  dataObj->sphereMappedRotFun.emplace_back ( new ProSHADE_internal_spheres::ProSHADE_rotFun_sphere ( angIt * ( 2.0 * M_PI / fold ),
4131  M_PI / fold,
4132  dataObj->maxShellBand * 2,
4133  angIt * ( 2.0 * M_PI / fold ),
4134  static_cast<proshade_unsign> ( angIt - 1.0 ) ) );
4135 
4136  //============================================ Interpolate rotation function onto the sphere
4137  dataObj->sphereMappedRotFun.at( static_cast < proshade_unsign > ( angIt - 1.0 ))->interpolateSphereValues ( dataObj->getInvSO3Coeffs ( ) );
4138  }
4139 
4140  //================================================ Convert XYZ to lat and lon INDICES
4141  lat = std::atan2( axis[1], axis[0] ) / latSamlUnit;
4142  lon = std::acos ( axis[2] ) / lonSamlUnit;
4143 
4144  if ( lat < 0.0 ) { lat += ( static_cast< proshade_double > ( dataObj->maxShellBand ) * 2.0 ); }
4145  if ( lon < 0.0 ) { lon += ( static_cast< proshade_double > ( dataObj->maxShellBand ) * 2.0 ); }
4146 
4147  lat = std::round ( lat );
4148  lon = std::round ( lon );
4149 
4150  //================================================ Initialise the peak group
4152 
4153  //================================================ Construct a peak group with entry from each sphere with the axis as the peak
4154  for ( proshade_unsign sphIt = 0; sphIt < static_cast<proshade_unsign> ( dataObj->sphereMappedRotFun.size() ); sphIt++ )
4155  {
4156  if ( sphIt == 0 )
4157  {
4158  //======================================== If first sphere, create the peak group
4159  grp = new ProSHADE_internal_spheres::ProSHADE_rotFun_spherePeakGroup ( lat, lon, sphIt, dataObj->sphereMappedRotFun.at(sphIt)->getAngularDim() );
4160  }
4161  else
4162  {
4163  //======================================== Add to the existing object
4164  grp->checkIfPeakBelongs ( lat, lon, sphIt, settings->axisErrTolerance, settings->verbose );
4165  }
4166  }
4167 
4168  //================================================ Find the peak height
4169  std::vector < proshade_double* > detectedAxis;
4170  grp->findCyclicPointGroupsGivenFold ( dataObj->sphereMappedRotFun, &detectedAxis, settings->useBiCubicInterpolationOnPeaks, static_cast< proshade_unsign > ( fold ), settings->verbose );
4171 
4172  //================================================ Save it!
4173  if ( detectedAxis.size() > 0 ) { height = detectedAxis.at(0)[5]; }
4174  else { height = 0.0; }
4175 
4176  //================================================ Release memory
4177  for ( proshade_unsign i = 0; i < static_cast < proshade_unsign > ( detectedAxis.size() ); i++ ) { delete detectedAxis.at(i); }
4178  delete grp;
4179 
4180  //================================================ Done
4181  return ( height );
4182 
4183 }
4184 
4192 std::pair< proshade_unsign, proshade_unsign > findBestTetraDihedralPair ( std::vector< proshade_double* >* CSymList, proshade_double minPeakHeight, proshade_double axErr )
4193 {
4194  //================================================ Initialise variables
4195  std::pair< proshade_unsign, proshade_unsign > ret;
4196  std::vector< proshade_unsign > C3List;
4197  proshade_double bestHeightSum = 0.0;
4198  proshade_double dotProduct;
4199 
4200  //================================================ Find all C3 symmetries
4201  for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ ) { const FloatingPoint< proshade_double > lhs1 ( CSymList->at(cSym)[0] ), rhs1 ( 3.0 ); if ( lhs1.AlmostEquals ( rhs1 ) && CSymList->at(cSym)[5] >= minPeakHeight ) { ProSHADE_internal_misc::addToUnsignVector ( &C3List, cSym ); } }
4202 
4203  //================================================ For each unique pair of C3 and C2
4204  for ( proshade_unsign c3 = 0; c3 < static_cast<proshade_unsign> ( C3List.size() ); c3++ )
4205  {
4206  for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
4207  {
4208  //======================================== Compare only C2s to the C3List and only with decent average peak height
4209  const FloatingPoint< proshade_double > lhs1 ( CSymList->at(cSym)[0] ), rhs1 ( 2.0 );
4210  if ( !lhs1.AlmostEquals ( rhs1 ) ) { continue; }
4211  if ( CSymList->at(cSym)[5] < minPeakHeight ) { continue; }
4212 
4213  //======================================== Check the angle between the C5 and C3 axes
4214  dotProduct = ProSHADE_internal_maths::computeDotProduct ( &CSymList->at(C3List.at(c3))[1],
4215  &CSymList->at(C3List.at(c3))[2],
4216  &CSymList->at(C3List.at(c3))[3],
4217  &CSymList->at(cSym)[1],
4218  &CSymList->at(cSym)[2],
4219  &CSymList->at(cSym)[3] );
4220 
4221  //======================================== Is the angle approximately the dihedral angle?
4222  if ( ( ( 1.0 / sqrt ( 3.0 ) ) > ( std::abs( dotProduct ) - axErr ) ) && ( ( 1.0 / sqrt ( 3.0 ) ) < ( std::abs( dotProduct ) + axErr ) ) )
4223  {
4224  if ( bestHeightSum < ( CSymList->at(C3List.at(c3))[5] + CSymList->at(cSym)[5] ) )
4225  {
4226  bestHeightSum = CSymList->at(C3List.at(c3))[5] + CSymList->at(cSym)[5];
4227  ret.first = C3List.at(c3);
4228  ret.second = cSym;
4229  }
4230  }
4231  }
4232  }
4233 
4234  //================================================ Done
4235  return ( ret );
4236 
4237 }
4238 
4253 std::vector< proshade_double* > ProSHADE_internal_data::ProSHADE_data::getPredictedTetrahedralSymmetriesList ( ProSHADE_settings* settings, std::vector< proshade_double* >* CSymList )
4254 {
4255  //================================================ Initialise variables
4256  std::vector< proshade_double* > ret;
4257 
4258  //================================================ Report progress
4259  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting T symmetry prediction." );
4260 
4261  //================================================ Are the basic requirements for icosahedral symmetry met?
4263  {
4264  //============================================ Generate the rest of the axes
4265  ProSHADE_internal_symmetry::predictTetraAxes ( CSymList, &ret, settings->axisErrTolerance, settings->minSymPeak );
4266 
4267  //============================================ Get heights for the predicted axes
4269 
4270  //============================================ Add predicted axes to detected C axes list and also to the settings Icosahedral symmetry list
4271  for ( proshade_unsign retIt = 0; retIt < static_cast < proshade_unsign > ( ret.size() ); retIt++ )
4272  {
4273  proshade_signed matchedPos = ProSHADE_internal_symmetry::addAxisUnlessSame ( static_cast< proshade_unsign > ( ret.at(retIt)[0] ), ret.at(retIt)[1], ret.at(retIt)[2], ret.at(retIt)[3], ret.at(retIt)[5], CSymList, settings->axisErrTolerance );
4274  ProSHADE_internal_misc::addToUnsignVector ( &settings->allDetectedTAxes, static_cast < proshade_unsign > ( matchedPos ) );
4275  }
4276  }
4277 
4278  //================================================ Report progress
4279  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "T symmetry prediction complete." );
4280 
4281  //================================================ Done
4282  return ( ret );
4283 
4284 }
4285 
4306 void ProSHADE_internal_symmetry::predictTetraAxes ( std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, proshade_double axErr, proshade_double minPeakHeight )
4307 {
4308  //================================================ Create the tetrahedronAxes object
4310 
4311  //================================================ Find the best axis combination with dihedral angle and correct folds
4312  std::pair< proshade_unsign, proshade_unsign > initAxes = findBestTetraDihedralPair ( CSymList, minPeakHeight, axErr );
4313 
4314  //================================================ Find rotation between the detected C3 and the model C3 axes.
4315  proshade_double* rotMat = ProSHADE_internal_maths::findRotMatMatchingVectors ( tetAx->getValue ( 0, 1 ),
4316  tetAx->getValue ( 0, 2 ),
4317  tetAx->getValue ( 0, 3 ),
4318  CSymList->at(initAxes.first)[1],
4319  CSymList->at(initAxes.first)[2],
4320  CSymList->at(initAxes.first)[3] );
4321 
4322  //================================================ Rotate the model C2 to the correct orientation relative to the detected C3 axis.
4323  proshade_double* rotModelC2 = ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( rotMat,
4324  tetAx->getValue ( 4, 1 ),
4325  tetAx->getValue ( 4, 2 ),
4326  tetAx->getValue ( 4, 3 ) );
4327 
4328  //================================================ Find the angle betwen the rotated model C2 and the detected C2 axes along the detected C3 axis.
4329  proshade_double bestAng = 0.0, curAngDist, bestAngDist = 999.9;
4330  proshade_double* rotMatHlp = new proshade_double[9];
4331  ProSHADE_internal_misc::checkMemoryAllocation ( rotMatHlp, __FILE__, __LINE__, __func__ );
4332  for ( proshade_double ang = 0.0; ang < ( M_PI * 2.0 ); ang += 0.002 )
4333  {
4334  //============================================ Compute rotation matrix for this angle value
4335  ProSHADE_internal_maths::getRotationMatrixFromAngleAxis ( rotMatHlp, CSymList->at(initAxes.first)[1], CSymList->at(initAxes.first)[2], CSymList->at(initAxes.first)[3], ang );
4336 
4337  //============================================ Rotate the rotated C2 by the matrix
4338  proshade_double* rotRotModelC2 = ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( rotMatHlp,
4339  rotModelC2[0],
4340  rotModelC2[1],
4341  rotModelC2[2] );
4342 
4343  //============================================ Find distance
4344  curAngDist = std::sqrt ( std::pow ( rotRotModelC2[0] - CSymList->at(initAxes.second)[1], 2.0 ) +
4345  std::pow ( rotRotModelC2[1] - CSymList->at(initAxes.second)[2], 2.0 ) +
4346  std::pow ( rotRotModelC2[2] - CSymList->at(initAxes.second)[3], 2.0 ) );
4347 
4348  //============================================ Save best angle
4349  if ( curAngDist < bestAngDist ) { bestAngDist = curAngDist; bestAng = ang; }
4350 
4351  //============================================ Release memory
4352  delete[] rotRotModelC2;
4353  }
4354 
4355  //================================================ Release memory
4356  delete[] rotMatHlp;
4357 
4358  //================================================ For the rotation matrix along the detected C5 axis with the same anlge as is between the rotated model C3 and the detected C3 axes.
4359  proshade_double* rotMat2 = new proshade_double[9];
4360  ProSHADE_internal_misc::checkMemoryAllocation ( rotMat2, __FILE__, __LINE__, __func__ );
4361  ProSHADE_internal_maths::getRotationMatrixFromAngleAxis ( rotMat2, CSymList->at(initAxes.first)[1], CSymList->at(initAxes.first)[2], CSymList->at(initAxes.first)[3], bestAng );
4362 
4363  //================================================ Combine the two rotation matrices into a single rotation matrix
4364  proshade_double* rotMatFin = ProSHADE_internal_maths::compute3x3MatrixMultiplication ( rotMat2, rotMat );
4365 
4366  //================================================ For each model axis
4367  for ( proshade_unsign iter = 0; iter < tetAx->getNoAxes( ); iter++ )
4368  {
4369  //============================================ Rotate the model axis to fit the detected orientation
4370  proshade_double* rotAxis = ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( rotMatFin,
4371  tetAx->getValue ( iter, 1 ),
4372  tetAx->getValue ( iter, 2 ),
4373  tetAx->getValue ( iter, 3 ) );
4374 
4375  //============================================ Create ProSHADE symmetry axis representation
4376  proshade_double* axis = new proshade_double[7];
4377  ProSHADE_internal_misc::checkMemoryAllocation ( axis, __FILE__, __LINE__, __func__ );
4378 
4379  axis[0] = tetAx->getValue ( iter, 0 );
4380  axis[1] = rotAxis[0];
4381  axis[2] = rotAxis[1];
4382  axis[3] = rotAxis[2];
4383  axis[4] = ( 2.0 * M_PI ) / axis[0];
4384  axis[5] = 0.0;
4385  axis[6] = -1.0;
4386 
4387  //============================================ Save axis to ret
4389 
4390  //============================================ Release memory
4391  delete[] rotAxis;
4392  }
4393 
4394  //================================================ Release memory
4395  delete[] rotMat;
4396  delete[] rotMat2;
4397  delete[] rotMatFin;
4398  delete[] rotModelC2;
4399  delete tetAx;
4400 
4401  //================================================ Done
4402  return ;
4403 
4404 }
ProSHADE_settings::noIQRsFromMedianNaivePeak
proshade_double noIQRsFromMedianNaivePeak
When doing peak searching, how many IQRs from the median the threshold for peak height should be (in ...
Definition: ProSHADE_settings.hpp:118
ProSHADE_internal_distances::normaliseEMatrices
void normaliseEMatrices(ProSHADE_internal_data::ProSHADE_data *obj1, ProSHADE_internal_data::ProSHADE_data *obj2, ProSHADE_settings *settings)
This function normalises the E matrices.
Definition: ProSHADE_distances.cpp:572
sortProSHADESymmetryByPeak
bool sortProSHADESymmetryByPeak(proshade_double *a, proshade_double *b)
This function allows using std::sort to sort vectors of ProSHADE symmetry format..
Definition: ProSHADE_symmetry.cpp:3872
ProSHADE_internal_symmetry::detectOctahedralSymmetry
bool detectOctahedralSymmetry(std::vector< proshade_double * > *CSymList, proshade_double axErr, proshade_double minPeakHeight)
This function takes the list of C symmetries and decides whether basic requirements for octahhedral s...
Definition: ProSHADE_symmetry.cpp:2268
ProSHADE_internal_maths::findVectorFromThreeVAndThreeD
std::vector< proshade_double > findVectorFromThreeVAndThreeD(proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2, proshade_double x3, proshade_double y3, proshade_double z3, proshade_double dot1, proshade_double dot2, proshade_double dot3)
Function for finding a vector which would have a given three dot products to three other vectors.
Definition: ProSHADE_maths.cpp:2197
ProSHADE_internal_symmetry::findExpectedPeakRotations
void findExpectedPeakRotations(proshade_unsign fold, std::vector< proshade_double > *expAngs)
This function computes the expected peak rotations for given fold.
Definition: ProSHADE_symmetry.cpp:719
ProSHADE_internal_precomputedVals::octahedronAxes
Definition: ProSHADE_precomputedValues.hpp:55
ProSHADE_internal_misc::sortSymInvFoldHlp
bool sortSymInvFoldHlp(const proshade_double *a, const proshade_double *b)
This function compares two arrays of two based on the first number, sorting highest first.
Definition: ProSHADE_misc.cpp:296
ProSHADE_internal_symmetry::findMissingAxes
bool findMissingAxes(std::vector< std::vector< proshade_unsign > > *possibilities, std::vector< proshade_double * > *CSymList, proshade_unsign requiredNoAxes, proshade_double axErr, proshade_double angle, proshade_unsign fold, ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_double minPeakHeight)
This function tries to find an axis which would complete a particular group of axes for polyhedral sy...
Definition: ProSHADE_symmetry.cpp:1672
ProSHADE_internal_symmetry::findIcos15C2s
void findIcos15C2s(std::vector< proshade_double * > *CSymList, std::vector< proshade_double * > *ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_signed verbose, proshade_double minPeakHeight)
This function takes the list of C symmetries and finds the fifteen C3 symmetries with correct angles ...
Definition: ProSHADE_symmetry.cpp:3510
ProSHADE_internal_misc::addToDblPtrVector
void addToDblPtrVector(std::vector< proshade_double * > *vecToAddTo, proshade_double *elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:143
ProSHADE_settings::maxSymmetryFold
proshade_unsign maxSymmetryFold
The highest symmetry fold to search for.
Definition: ProSHADE_settings.hpp:133
ProSHADE_internal_symmetry::groupSameAxes
std::vector< std::vector< proshade_unsign > > groupSameAxes(std::vector< proshade_double * > &peaks, proshade_double errTolerance)
This function groups the peaks by their axes of rotation.
Definition: ProSHADE_symmetry.cpp:426
ProSHADE_internal_data::ProSHADE_data::computeRotationFunction
void computeRotationFunction(ProSHADE_settings *settings)
This function computes the self-rotation function for this structure.
Definition: ProSHADE_symmetry.cpp:41
ProSHADE_internal_symmetry::sortArrVecHlp
bool sortArrVecHlp(const proshade_double *a, const proshade_double *b)
This function compares two arrays of two based on the first number.
Definition: ProSHADE_symmetry.cpp:1732
ProSHADE_internal_symmetry::testGroupAgainstGroup
bool testGroupAgainstGroup(std::vector< proshade_double * > *CSymList, std::vector< proshade_unsign > *grp1, std::vector< proshade_double * > *RetList, std::vector< proshade_unsign > *grp2, proshade_double angle, proshade_double axErr)
This function compares two groups of axes for a single pair having the required angle.
Definition: ProSHADE_symmetry.cpp:2156
ProSHADE_internal_symmetry::predictOctaAxes
void predictOctaAxes(std::vector< proshade_double * > *CSymList, std::vector< proshade_double * > *ret, proshade_double axErr, proshade_double minPeakHeight)
This function predicts all octahedral point group symmetry axes from the cyclic point groups list.
Definition: ProSHADE_symmetry.cpp:3325
ProSHADE_internal_symmetry::smallestDistanceBetweenAngles
bool smallestDistanceBetweenAngles(std::vector< proshade_unsign > grp, std::vector< proshade_double * > peaks, std::vector< proshade_double > *tried, proshade_double *dist)
This function finds the smallest distance between the rotation angles within a group.
Definition: ProSHADE_symmetry.cpp:562
ProSHADE_internal_symmetry::predictTetraAxes
void predictTetraAxes(std::vector< proshade_double * > *CSymList, std::vector< proshade_double * > *ret, proshade_double axErr, proshade_double minPeakHeight)
This function predicts all tetrahedral point group symmetry axes from the cyclic point groups list.
Definition: ProSHADE_symmetry.cpp:4306
ProSHADE_internal_symmetry::addAxisUnlessSame
proshade_signed addAxisUnlessSame(proshade_unsign fold, proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double axHeight, proshade_double averageFSC, std::vector< proshade_double * > *prosp, proshade_double axErr)
This function simply creates a new axis from information in aruments and tests if no such axis alread...
Definition: ProSHADE_symmetry.cpp:2648
ProSHADE_internal_spheres::ProSHADE_rotFun_sphere
This class contains all inputed data for the rotation function angle-axis converted spheres.
Definition: ProSHADE_maths.hpp:54
ProSHADE_internal_precomputedVals::icosahedronAxes::getValue
proshade_double getValue(proshade_unsign axis, proshade_unsign element)
Accessor for the icosahedronAxesVals variable.
Definition: ProSHADE_precomputedValues.cpp:223
ProSHADE_internal_symmetry::saveAllCSymmetries
void saveAllCSymmetries(std::vector< std::vector< proshade_unsign > > detected, std::vector< proshade_double * > peaks, std::vector< proshade_double * > *ret, proshade_double axErr)
This function takes the detected symmetries indices and peaks and saves these in the main cyclic symm...
Definition: ProSHADE_symmetry.cpp:1133
ProSHADE_internal_maths::normalDistributionValue
proshade_double normalDistributionValue(proshade_double mean, proshade_double standardDev, proshade_double value)
Function to the heiht of normal distribution given by mean and standard deviation for a given value.
Definition: ProSHADE_maths.cpp:1756
ProSHADE_internal_symmetry::findPredictedAxesHeights
void findPredictedAxesHeights(std::vector< proshade_double * > *ret, ProSHADE_internal_data::ProSHADE_data *dataObj, ProSHADE_settings *settings)
This function finds the rotation function value for all axes supplied in the ret parameter.
Definition: ProSHADE_symmetry.cpp:4008
ProSHADE_settings::allDetectedDAxes
std::vector< std::vector< proshade_unsign > > allDetectedDAxes
The vector of all detected dihedral symmetry axes indices in allDetectedCAxes.
Definition: ProSHADE_settings.hpp:148
ProSHADE_internal_misc::addToUnsignVectorVector
void addToUnsignVectorVector(std::vector< std::vector< proshade_unsign > > *vecToAddTo, std::vector< proshade_unsign > elementToAdd)
Adds the element to the vector of vectors.
Definition: ProSHADE_misc.cpp:211
ProSHADE_internal_peakSearch::getAllPeaksNaive
std::vector< proshade_double * > getAllPeaksNaive(proshade_complex *map, proshade_unsign dim, proshade_signed peakSize, proshade_double noIQRs)
This function finds peaks in the 3D map using the "naive" approach.
Definition: ProSHADE_peakSearch.cpp:316
ProSHADE_internal_data::ProSHADE_data
This class contains all inputed and derived data for a single structure.
Definition: ProSHADE_data.hpp:49
ProSHADE_internal_symmetry::findOcta4C3s
void findOcta4C3s(std::vector< proshade_double * > *CSymList, std::vector< proshade_double * > *ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_signed verbose, proshade_double minPeakHeight)
This function takes the list of C symmetries and finds the four C3 symmetries with correct angles req...
Definition: ProSHADE_symmetry.cpp:2390
ProSHADE_internal_symmetry::printSymmetryCompletion
void printSymmetryCompletion(proshade_unsign noSyms, proshade_signed verbose)
This function simply prints the summary and warnings for cyclic symmetries detection completion.
Definition: ProSHADE_symmetry.cpp:1103
ProSHADE_internal_symmetry::findPeaksByHeightBoundaries
std::vector< proshade_double > findPeaksByHeightBoundaries(std::vector< proshade_double * > allPeaks, proshade_double smoothing)
This function groups the peaks by height and returns the boundaries between such groups.
Definition: ProSHADE_symmetry.cpp:317
ProSHADE_internal_symmetry::missingAxisHeight
proshade_double missingAxisHeight(proshade_double xVal, proshade_double yVal, proshade_double zVal, ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_unsign fold, proshade_double axErr)
This function searches for the highest peaks average that would produce the required axis and fold.
Definition: ProSHADE_symmetry.cpp:1753
ProSHADE_internal_symmetry::searchMissingSymmetrySpace
void searchMissingSymmetrySpace(ProSHADE_internal_data::ProSHADE_data *dataObj, std::vector< proshade_double * > *CSymList, std::vector< proshade_unsign > *grp, std::vector< proshade_double * > *hlpVec, proshade_double axErr, proshade_double angle, proshade_unsign fold, proshade_double minPeakHeight)
This function tests feasible axes against the missing axis criteria, returning a set of matching axes...
Definition: ProSHADE_symmetry.cpp:1964
ProSHADE_settings::minSymPeak
proshade_double minSymPeak
Minimum average peak for symmetry axis to be considered as "real".
Definition: ProSHADE_settings.hpp:127
ProSHADE_internal_symmetry::findMissingAxesDual
bool findMissingAxesDual(std::vector< proshade_unsign > *possibilities, std::vector< proshade_double * > *CSymList, std::vector< proshade_double * > *ret, std::vector< proshade_unsign > *retGroup, proshade_unsign requiredNoAxes, proshade_double axErr, proshade_unsign noMatchesG1, proshade_double angle1, proshade_unsign noMatchesG2, proshade_double angle2, proshade_unsign fold, ProSHADE_internal_data::ProSHADE_data *dataObj)
This function tries to find a particular symmetry axes which would complete a group of symmetries wit...
Definition: ProSHADE_symmetry.cpp:2553
ProSHADE_internal_data::ProSHADE_data::maxShellBand
proshade_unsign maxShellBand
The maximum band for any shell of the object.
Definition: ProSHADE_data.hpp:123
ProSHADE_settings::symMissPeakThres
proshade_double symMissPeakThres
Percentage of peaks that could be missing that would warrant starting the missing peaks search proced...
Definition: ProSHADE_settings.hpp:124
ProSHADE_internal_maths::getRotationMatrixFromEulerZXZAngles
void getRotationMatrixFromEulerZXZAngles(proshade_double eulerAlpha, proshade_double eulerBeta, proshade_double eulerGamma, proshade_double *matrix)
Function to find the rotation matrix from Euler angles (ZXZ convention).
Definition: ProSHADE_maths.cpp:1007
ProSHADE_internal_data::ProSHADE_data::getDihedralSymmetriesList
std::vector< proshade_double * > getDihedralSymmetriesList(ProSHADE_settings *settings, std::vector< proshade_double * > *CSymList)
This function obtains a list of all D symmetries from already computed C symmetries list.
Definition: ProSHADE_symmetry.cpp:1299
ProSHADE_internal_symmetry::addZeroPeakToGroups
void addZeroPeakToGroups(std::vector< std::vector< proshade_unsign > > &grpsVec, std::vector< proshade_double * > &peaks)
This function takes the peak groups and adds zero peak to each of them.
Definition: ProSHADE_symmetry.cpp:635
ProSHADE_symmetry.hpp
This header file declares all the functions required for symmetry detection and construction.
ProSHADE_settings::allDetectedTAxes
std::vector< proshade_unsign > allDetectedTAxes
The vector of all detected tetrahedral symmetry axes indices in allDetectedCAxes.
Definition: ProSHADE_settings.hpp:149
ProSHADE_internal_messages::printWarningMessage
void printWarningMessage(proshade_signed verbose, std::string message, std::string warnCode)
General stderr message printing (used for warnings).
Definition: ProSHADE_messages.cpp:101
ProSHADE_internal_symmetry::findPeaksCSymmetry
std::vector< std::vector< proshade_unsign > > findPeaksCSymmetry(std::vector< proshade_double * > *peaks, proshade_signed verbose, proshade_unsign band, proshade_double missPeakThres, proshade_double axisErrTolerance, bool axisErrToleranceDef, ProSHADE_internal_data::ProSHADE_data *dataObj)
This function searches the list of peaks for presence of cyclic symmetry.
Definition: ProSHADE_symmetry.cpp:375
ProSHADE_settings::allDetectedOAxes
std::vector< proshade_unsign > allDetectedOAxes
The vector of all detected octahedral symmetry axes indices in allDetectedCAxes.
Definition: ProSHADE_settings.hpp:150
ProSHADE_settings::verbose
proshade_signed verbose
Should the software report on the progress, or just be quiet? Value between -1 (nothing) and 4 (loud)
Definition: ProSHADE_settings.hpp:142
ProSHADE_internal_symmetry::isSymmetrySame
bool isSymmetrySame(std::vector< proshade_double * > *ret, proshade_double *sym, proshade_double simThres, proshade_signed *matchedPos)
This function checks if a very similar symmetry is not already saved.
Definition: ProSHADE_symmetry.cpp:1197
ProSHADE_internal_symmetry::printSymmetryGroup
void printSymmetryGroup(std::vector< proshade_unsign > grp, std::vector< proshade_double * > peaks, proshade_signed verbose)
This function simply prints the detected symmetry and all its supporting peaks.
Definition: ProSHADE_symmetry.cpp:1077
ProSHADE_internal_maths::vectorOrientationSimilarity
bool vectorOrientationSimilarity(proshade_double a1, proshade_double a2, proshade_double a3, proshade_double b1, proshade_double b2, proshade_double b3, proshade_double tolerance=0.1)
This function compares two vectors using cosine distance and decides if they are similar using tolera...
Definition: ProSHADE_maths.cpp:2332
ProSHADE_internal_misc::addToDoubleVector
void addToDoubleVector(std::vector< proshade_double > *vecToAddTo, proshade_double elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:77
ProSHADE_internal_symmetry::checkExpectedAgainstFound
proshade_unsign checkExpectedAgainstFound(std::vector< proshade_unsign > grp, std::vector< proshade_double * > peaks, std::vector< proshade_double > *expAngs, std::vector< proshade_unsign > *matchedAngs, std::vector< proshade_unsign > *missingAngs, proshade_double axisTol)
This function computes the expected peak rotations for given fold.
Definition: ProSHADE_symmetry.cpp:748
ProSHADE_internal_precomputedVals::octahedronAxes::getValue
proshade_double getValue(proshade_unsign axis, proshade_unsign element)
Accessor for the octahedronAxesVals variable.
Definition: ProSHADE_precomputedValues.cpp:136
ProSHADE_internal_precomputedVals::tetrahedronAxes::getValue
proshade_double getValue(proshade_unsign axis, proshade_unsign element)
Accessor for the tetrahedronAxesVals variable.
Definition: ProSHADE_precomputedValues.cpp:68
ProSHADE_settings::useBiCubicInterpolationOnPeaks
bool useBiCubicInterpolationOnPeaks
This variable switch decides whether best symmetry is detected from peak indices, or whether bicubic ...
Definition: ProSHADE_settings.hpp:132
determinePeakThreshold
proshade_double determinePeakThreshold(std::vector< proshade_double > inArr, proshade_double noIQRsFromMedian)
This function takes a vector of values and determines the threshold for removing noise from it.
Definition: ProSHADE_symmetry.cpp:71
ProSHADE_internal_symmetry::saveMissingAxisNewOnly
void saveMissingAxisNewOnly(std::vector< proshade_double * > *axVec, proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double height, proshade_unsign fold, proshade_double axErr)
This function saves the recovered information about missing axis into a full symmetry,...
Definition: ProSHADE_symmetry.cpp:1895
ProSHADE_internal_maths::computeDotProduct
proshade_double computeDotProduct(proshade_double *x1, proshade_double *y1, proshade_double *z1, proshade_double *x2, proshade_double *y2, proshade_double *z2)
Simple 3D vector dot product computation.
Definition: ProSHADE_maths.cpp:1773
ProSHADE_internal_symmetry::findTetra3C2s
void findTetra3C2s(std::vector< proshade_double * > *CSymList, std::vector< proshade_double * > *ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_signed verbose, proshade_double minPeakHeight)
This function takes the list of C symmetries and finds the 3 C2 symmetries with correct angles requir...
Definition: ProSHADE_symmetry.cpp:2071
ProSHADE_internal_misc::deepCopyAxisToDblPtrVector
void deepCopyAxisToDblPtrVector(std::vector< proshade_double * > *dblPtrVec, proshade_double *axis)
Does a deep copy of a double array to a vector of double arrays.
Definition: ProSHADE_misc.cpp:310
ProSHADE_internal_precomputedVals::tetrahedronAxes::getNoAxes
proshade_unsign getNoAxes()
Accessor for the tetrahedronAxesVals variable number of axes.
Definition: ProSHADE_precomputedValues.cpp:79
ProSHADE_internal_symmetry::findOcta6C2s
void findOcta6C2s(std::vector< proshade_double * > *CSymList, std::vector< proshade_double * > *ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_signed verbose, proshade_double minPeakHeight)
This function takes the list of C symmetries and finds the six C2 symmetries with correct angles requ...
Definition: ProSHADE_symmetry.cpp:2472
ProSHADE_internal_maths::getRotationMatrixFromAngleAxis
void getRotationMatrixFromAngleAxis(proshade_double *rotMat, proshade_double x, proshade_double y, proshade_double z, proshade_double ang)
This function converts the axis-angle representation to the rotation matrix representation.
Definition: ProSHADE_maths.cpp:1444
ProSHADE_internal_data::ProSHADE_data::getMaxBand
proshade_unsign getMaxBand(void)
This function returns the maximum band value for the object.
Definition: ProSHADE_data.cpp:3599
ProSHADE_internal_symmetry::detectIcosahedralSymmetry
bool detectIcosahedralSymmetry(std::vector< proshade_double * > *CSymList, proshade_double axErr, proshade_double minPeakHeight)
This function takes the list of C symmetries and decides whether basic requirements for isosahedral s...
Definition: ProSHADE_symmetry.cpp:2960
ProSHADE_internal_spheres::ProSHADE_rotFun_spherePeakGroup::checkIfPeakBelongs
bool checkIfPeakBelongs(proshade_double lat, proshade_double lon, proshade_unsign sphPos, proshade_double cosTol, proshade_signed verbose)
This function takes a new prospective peak and tests if it belongs to this peak group or not.
Definition: ProSHADE_spheres.cpp:1134
ProSHADE_settings
This class stores all the settings and is passed to the executive classes instead of a multitude of p...
Definition: ProSHADE_settings.hpp:37
ProSHADE_internal_data::ProSHADE_data::getCyclicSymmetriesList
std::vector< proshade_double * > getCyclicSymmetriesList(ProSHADE_settings *settings)
This function obtains a list of all C symmetries from already computed self-rotation map.
Definition: ProSHADE_symmetry.cpp:198
ProSHADE_internal_data::ProSHADE_data::getPredictedTetrahedralSymmetriesList
std::vector< proshade_double * > getPredictedTetrahedralSymmetriesList(ProSHADE_settings *settings, std::vector< proshade_double * > *CSymList)
This function predicts a list of all T symmetry axes from the already computed C symmetries list.
Definition: ProSHADE_symmetry.cpp:4253
ProSHADE_internal_distances::computeEMatrices
void computeEMatrices(ProSHADE_internal_data::ProSHADE_data *obj1, ProSHADE_internal_data::ProSHADE_data *obj2, ProSHADE_settings *settings)
This function computes the complete E matrices and their weights between any two objects.
Definition: ProSHADE_distances.cpp:510
ProSHADE_internal_symmetry::giveOppositeAxesSameDirection
void giveOppositeAxesSameDirection(std::vector< proshade_double * > peaks)
This function modifiest the axes so that the highest vector element is always positive.
Definition: ProSHADE_symmetry.cpp:498
ProSHADE_internal_data::ProSHADE_data::convertRotationFunction
void convertRotationFunction(ProSHADE_settings *settings)
This function converts the self-rotation function of this structure to angle-axis representation.
Definition: ProSHADE_symmetry.cpp:120
ProSHADE_internal_symmetry::printSymmetryPeaks
void printSymmetryPeaks(std::vector< proshade_unsign > grp, std::vector< proshade_double * > peaks, proshade_signed verbose, proshade_unsign groupNo)
This function simply prints the symmetry axis group supplied in the first parameter from the second p...
Definition: ProSHADE_symmetry.cpp:529
ProSHADE_internal_distances::computeInverseSOFTTransform
void computeInverseSOFTTransform(ProSHADE_internal_data::ProSHADE_data *obj1, ProSHADE_internal_data::ProSHADE_data *obj2, ProSHADE_settings *settings)
This function computes the inverse SO(3) transform.
Definition: ProSHADE_distances.cpp:854
ProSHADE_internal_spheres::ProSHADE_rotFun_spherePeakGroup
This class contains peak groups detected in the rotation function mapped spheres.
Definition: ProSHADE_spheres.hpp:129
ProSHADE_internal_symmetry::findIcos10C3s
void findIcos10C3s(std::vector< proshade_double * > *CSymList, std::vector< proshade_double * > *ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_signed verbose, proshade_double minPeakHeight)
This function takes the list of C symmetries and finds the ten C3 symmetries with correct angles requ...
Definition: ProSHADE_symmetry.cpp:3438
ProSHADE_internal_maths::isAxisUnique
bool isAxisUnique(std::vector< proshade_double * > *CSymList, proshade_double *axis, proshade_double tolerance=0.1, bool improve=false)
This function checks if new axis is unique, or already detected.
Definition: ProSHADE_maths.cpp:2717
ProSHADE_internal_symmetry::findMissingAxesTriple
bool findMissingAxesTriple(std::vector< proshade_unsign > *possibilities, std::vector< proshade_double * > *CSymList, std::vector< proshade_double * > *ret, std::vector< proshade_unsign > *retGroup, proshade_unsign requiredNoAxes, proshade_double axErr, proshade_unsign noMatchesG1, proshade_double angle1, proshade_unsign noMatchesG2, proshade_double angle2, proshade_unsign noMatchesG3, proshade_double angle3, proshade_unsign fold, ProSHADE_internal_data::ProSHADE_data *dataObj)
This function tries to find a particular symmetry axis which would complete a group of symmetries wit...
Definition: ProSHADE_symmetry.cpp:3593
ProSHADE_internal_maths::getAxisAngleFromRotationMatrix
void getAxisAngleFromRotationMatrix(proshade_double *rotMat, proshade_double *x, proshade_double *y, proshade_double *z, proshade_double *ang)
This function converts rotation matrix to the axis-angle representation.
Definition: ProSHADE_maths.cpp:1124
ProSHADE_internal_maths::findRotMatMatchingVectors
proshade_double * findRotMatMatchingVectors(proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2)
Computation of rotation matrix rotating one vector onto the other.
Definition: ProSHADE_maths.cpp:1975
findBestIcosDihedralPair
std::vector< std::pair< proshade_unsign, proshade_unsign > > findBestIcosDihedralPair(std::vector< proshade_double * > *CSymList, proshade_double minPeakHeight, proshade_double axErr)
This function finds all the pairs of axes conforming to the icosahedron dihedral angle.
Definition: ProSHADE_symmetry.cpp:3079
ProSHADE_settings::peakNeighbours
proshade_unsign peakNeighbours
Number of points in any direction that have to be lower than the considered index in order to conside...
Definition: ProSHADE_settings.hpp:117
ProSHADE_internal_symmetry::completeMissingCSymmetry
bool completeMissingCSymmetry(ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_unsign fold, std::vector< proshade_unsign > *grp, std::vector< proshade_double * > *peaks, std::vector< proshade_unsign > *missingPeaks, std::vector< proshade_double > *expectedAngles, std::vector< proshade_unsign > *matchedPeaks, proshade_double axErrTolerance, proshade_signed verbose)
This function does the complete missing peak searching and filling in the missing peaks.
Definition: ProSHADE_symmetry.cpp:942
ProSHADE_settings::smoothingFactor
proshade_double smoothingFactor
This factor decides how small the group sizes should be - larger factor means more smaller groups.
Definition: ProSHADE_settings.hpp:121
ProSHADE_internal_mapManip::myRound
proshade_signed myRound(proshade_double x)
Calls the appropriate version of round function depending on compiler version.
Definition: ProSHADE_mapManip.cpp:31
ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication
proshade_double * compute3x3MatrixVectorMultiplication(proshade_double *mat, proshade_double x, proshade_double y, proshade_double z)
Function for computing a 3x3 matrix to 3x1 vector multiplication.
Definition: ProSHADE_maths.cpp:1857
ProSHADE_internal_precomputedVals::tetrahedronAxes
Definition: ProSHADE_precomputedValues.hpp:69
findBestTetraDihedralPair
std::pair< proshade_unsign, proshade_unsign > findBestTetraDihedralPair(std::vector< proshade_double * > *CSymList, proshade_double minPeakHeight, proshade_double axErr)
This function finds the best pair of axes conforming to the tetrahedron dihedral angle.
Definition: ProSHADE_symmetry.cpp:4192
ProSHADE_internal_symmetry::testGroupAgainstSymmetry
bool testGroupAgainstSymmetry(std::vector< proshade_double * > *CSymList, std::vector< proshade_unsign > *grp, proshade_double *sym, proshade_double axErr, proshade_double angle, bool improve, proshade_unsign pos=0)
This function tests whether a symmetry has particular angle to all members of a group.
Definition: ProSHADE_symmetry.cpp:1605
ProSHADE_internal_precomputedVals::icosahedronAxes
Definition: ProSHADE_precomputedValues.hpp:41
ProSHADE_internal_misc::checkMemoryAllocation
void checkMemoryAllocation(chVar checkVar, std::string fileP, unsigned int lineP, std::string funcP, std::string infoP="This error may occurs when ProSHADE requests memory to be\n : allocated to it and this operation fails. This could\n : happen when not enough memory is available, either due to\n : other processes using a lot of memory, or when the machine\n : does not have sufficient memory available. Re-run to see\n : if this problem persists.")
Checks if memory was allocated properly.
Definition: ProSHADE_misc.hpp:67
ProSHADE_internal_symmetry::saveDetectedCSymmetry
void saveDetectedCSymmetry(proshade_unsign fold, std::vector< proshade_unsign > *matchedPeaks, std::vector< std::vector< proshade_unsign > > *ret, proshade_signed verbose)
This function saves a detected symmetry for reporting to the user.
Definition: ProSHADE_symmetry.cpp:902
ProSHADE_internal_maths::getEulerZXZFromSOFTPosition
void getEulerZXZFromSOFTPosition(proshade_signed band, proshade_signed x, proshade_signed y, proshade_signed z, proshade_double *eulerAlpha, proshade_double *eulerBeta, proshade_double *eulerGamma)
Function to find Euler angles (ZXZ convention) from index position in the inverse SOFT map.
Definition: ProSHADE_maths.cpp:963
ProSHADE_internal_symmetry::saveDSymmetry
void saveDSymmetry(std::vector< proshade_double * > *ret, std::vector< proshade_double * > *CSymList, proshade_unsign axisOne, proshade_unsign axisTwo)
This function saves a detected dihedral symmetry to the dihedral symmetries list.
Definition: ProSHADE_symmetry.cpp:1380
ProSHADE_internal_maths::vectorMedianAndIQR
void vectorMedianAndIQR(std::vector< proshade_double > *vec, proshade_double *&ret)
Function to get vector median and inter-quartile range.
Definition: ProSHADE_maths.cpp:149
ProSHADE_internal_precomputedVals::octahedronAxes::getNoAxes
proshade_unsign getNoAxes()
Accessor for the octahedronAxesVals variable number of axes.
Definition: ProSHADE_precomputedValues.cpp:147
ProSHADE_internal_data::ProSHADE_data::getPredictedOctahedralSymmetriesList
std::vector< proshade_double * > getPredictedOctahedralSymmetriesList(ProSHADE_settings *settings, std::vector< proshade_double * > *CSymList)
This function predicts a list of all O symmetry axes from the already computed C symmetries list.
Definition: ProSHADE_symmetry.cpp:2916
ProSHADE_internal_symmetry::checkForMissingPeak
proshade_double checkForMissingPeak(ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_double x, proshade_double y, proshade_double z, proshade_double angle, proshade_double heightThres, proshade_double axTol)
This function checks for the high of the correlation for particular rotation angle and axis.
Definition: ProSHADE_symmetry.cpp:830
ProSHADE_internal_data::ProSHADE_data::getInvSO3Coeffs
proshade_complex * getInvSO3Coeffs(void)
This function allows access to the inverse SO(3) coefficients array.
Definition: ProSHADE_data.cpp:3826
ProSHADE_settings::allDetectedIAxes
std::vector< proshade_unsign > allDetectedIAxes
The vector of all detected icosahedral symmetry axes indices in allDetectedCAxes.
Definition: ProSHADE_settings.hpp:151
ProSHADE_internal_symmetry::findPredictedSingleAxisHeight
proshade_double findPredictedSingleAxisHeight(proshade_double *axis, proshade_double fold, ProSHADE_internal_data::ProSHADE_data *dataObj, ProSHADE_settings *settings)
This function finds the rotation function value for a single axis.
Definition: ProSHADE_symmetry.cpp:4115
ProSHADE_internal_misc::addToUnsignVector
void addToUnsignVector(std::vector< proshade_unsign > *vecToAddTo, proshade_unsign elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:99
ProSHADE_internal_symmetry::findOcta3C4s
void findOcta3C4s(std::vector< proshade_double * > *CSymList, std::vector< proshade_double * > *ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_signed verbose, proshade_double minPeakHeight)
This function takes the list of C symmetries and finds the 3 C4 symmetries with perpendicular angles ...
Definition: ProSHADE_symmetry.cpp:2324
ProSHADE_internal_data::ProSHADE_data::getIcosahedralSymmetriesList
std::vector< proshade_double * > getIcosahedralSymmetriesList(ProSHADE_settings *settings, std::vector< proshade_double * > *CSymList)
This function obtains a list of all I symmetry axes from the already computed C symmetries list.
Definition: ProSHADE_symmetry.cpp:2799
ProSHADE_internal_symmetry::findTetra4C3s
void findTetra4C3s(std::vector< proshade_double * > *CSymList, std::vector< proshade_double * > *ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_signed verbose, proshade_double minPeakHeight)
This function takes the list of C symmetries and finds the 4 C3 symmetries with correct angles requir...
Definition: ProSHADE_symmetry.cpp:1535
ProSHADE_internal_data::ProSHADE_data::getPredictedIcosahedralSymmetriesList
std::vector< std::vector< proshade_double * > > getPredictedIcosahedralSymmetriesList(ProSHADE_settings *settings, std::vector< proshade_double * > *CSymList)
This function predicts a list of all I symmetry axes from the already computed C symmetries list.
Definition: ProSHADE_symmetry.cpp:2871
ProSHADE_internal_data::ProSHADE_data::getCyclicSymmetriesListFromAngleAxis
std::vector< proshade_double * > getCyclicSymmetriesListFromAngleAxis(ProSHADE_settings *settings)
This function obtains a list of all C symmetries from the angle-axis space mapped rotation function v...
Definition: ProSHADE_symmetry.cpp:3756
ProSHADE_internal_symmetry::checkFittingAxisTripleAndSave
void checkFittingAxisTripleAndSave(std::vector< proshade_unsign > *retGroup, std::vector< proshade_double * > *ret, proshade_unsign fold, proshade_double axX, proshade_double axY, proshade_double axZ, std::vector< proshade_double * > *prosp, proshade_double axErr, proshade_unsign noMatchesG1, proshade_double angle1, proshade_unsign noMatchesG2, proshade_double angle2, proshade_unsign noMatchesG3, proshade_double angle3, ProSHADE_internal_data::ProSHADE_data *dataObj)
This function takes a newly detected "missing" axis and tests it for belonging to the group,...
Definition: ProSHADE_symmetry.cpp:3704
findBestOctaDihedralPair
std::pair< proshade_unsign, proshade_unsign > findBestOctaDihedralPair(std::vector< proshade_double * > *CSymList, proshade_double minPeakHeight, proshade_double axErr)
This function finds the best pair of axes conforming to the octahedron dihedral angle.
Definition: ProSHADE_symmetry.cpp:3258
ProSHADE_internal_data::ProSHADE_data::getTetrahedralSymmetriesList
std::vector< proshade_double * > getTetrahedralSymmetriesList(ProSHADE_settings *settings, std::vector< proshade_double * > *CSymList)
This function obtains a list of all T symmetry axes from the already computed C symmetries list.
Definition: ProSHADE_symmetry.cpp:1419
ProSHADE_internal_spheres::ProSHADE_rotFun_spherePeakGroup::findCyclicPointGroupsGivenFold
void findCyclicPointGroupsGivenFold(std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > sphereVals, std::vector< proshade_double * > *detectedCs, bool bicubicInterp, proshade_unsign fold, proshade_signed verbose)
Function detecting cyclic point groups with a particular fold in a peak group.
Definition: ProSHADE_spheres.cpp:1362
ProSHADE_internal_maths::findVectorFromTwoVAndTwoD
std::vector< proshade_double > findVectorFromTwoVAndTwoD(proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2, proshade_double dot1, proshade_double dot2)
Function for finding a vector which would have a given two dot products to two other vectors.
Definition: ProSHADE_maths.cpp:2051
ProSHADE_internal_symmetry::getPeaksAngleAxisPositions
std::vector< proshade_double * > getPeaksAngleAxisPositions(std::vector< proshade_double * > allPeaks, proshade_signed verbose)
This function converts peaks ZXZ Euler anles to angle-axis representation for further processing.
Definition: ProSHADE_symmetry.cpp:263
ProSHADE_settings::axisErrTolerance
proshade_double axisErrTolerance
Allowed error on vector axis in in dot product ( acos ( 1 - axErr ) is the allowed difference in radi...
Definition: ProSHADE_settings.hpp:125
ProSHADE_internal_precomputedVals::icosahedronAxes::getNoAxes
proshade_unsign getNoAxes()
Accessor for the octahedronAxesVals variable number of axes.
Definition: ProSHADE_precomputedValues.cpp:234
ProSHADE_internal_symmetry::determineFoldToTry
bool determineFoldToTry(proshade_double dist, proshade_double *divBasis, proshade_double *divRem, proshade_double peakErr, proshade_double *symmErr, std::vector< proshade_unsign > *angsToTry)
This function determines the symmetry fold to be searched for.
Definition: ProSHADE_symmetry.cpp:672
ProSHADE_internal_distances::generateSO3CoeffsFromEMatrices
void generateSO3CoeffsFromEMatrices(ProSHADE_internal_data::ProSHADE_data *obj1, ProSHADE_internal_data::ProSHADE_data *obj2, ProSHADE_settings *settings)
This function converts the E matrices to SO(3) coefficients.
Definition: ProSHADE_distances.cpp:701
ProSHADE_internal_symmetry::findIcos6C5s
void findIcos6C5s(std::vector< proshade_double * > *CSymList, std::vector< proshade_double * > *ret, proshade_double axErr, ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_signed verbose, proshade_double minPeakHeight)
This function takes the list of C symmetries and finds the six C5 symmetries with given angles requir...
Definition: ProSHADE_symmetry.cpp:3020
ProSHADE_internal_data::ProSHADE_data::findRequestedCSymmetryFromAngleAxis
std::vector< proshade_double * > findRequestedCSymmetryFromAngleAxis(ProSHADE_settings *settings, proshade_unsign fold, proshade_double *peakThres)
This function searches the angle-axis representation of the rotation function for a cyclic point grou...
Definition: ProSHADE_symmetry.cpp:3895
ProSHADE_internal_maths::findAllPrimes
std::vector< proshade_unsign > findAllPrimes(proshade_unsign upTo)
This function finds all prime numbers up to the supplied limit.
Definition: ProSHADE_maths.cpp:2798
ProSHADE_internal_data::ProSHADE_data::getOctahedralSymmetriesList
std::vector< proshade_double * > getOctahedralSymmetriesList(ProSHADE_settings *settings, std::vector< proshade_double * > *CSymList)
This function obtains a list of all O symmetry axes from the already computed C symmetries list.
Definition: ProSHADE_symmetry.cpp:2199
ProSHADE_internal_messages::printProgressMessage
void printProgressMessage(proshade_signed verbose, proshade_signed messageLevel, std::string message)
General stdout message printing.
Definition: ProSHADE_messages.cpp:70
ProSHADE_internal_symmetry::findSymmetryUsingFold
void findSymmetryUsingFold(ProSHADE_internal_data::ProSHADE_data *dataObj, std::vector< proshade_unsign > *angsToTry, std::vector< proshade_unsign > *grp, std::vector< proshade_double * > *peaks, std::vector< std::vector< proshade_unsign > > *ret, std::vector< proshade_unsign > *testedAlready, proshade_double axErrTolerance, bool axErrToleranceDefault, proshade_double missPeakThres, proshade_signed verbose)
This function tests all supplied folds for being supported by the peaks (i.e. and being complete pres...
Definition: ProSHADE_symmetry.cpp:1007
ProSHADE_internal_symmetry::predictIcosAxes
void predictIcosAxes(std::vector< proshade_double * > *CSymList, std::vector< std::vector< proshade_double * > > *ret, proshade_double axErr, proshade_double minPeakHeight)
This function predicts all possible icosahedral point groups symmetry axes from the cyclic point grou...
Definition: ProSHADE_symmetry.cpp:3142
ProSHADE_internal_symmetry::detectTetrahedralSymmetry
bool detectTetrahedralSymmetry(std::vector< proshade_double * > *CSymList, proshade_double axErr, proshade_double minPeakHeight)
This function takes the list of C symmetries and decides whether basic requirements for tetrahedral s...
Definition: ProSHADE_symmetry.cpp:1485
ProSHADE_internal_misc::sortSymHlpInv
bool sortSymHlpInv(const proshade_double *a, const proshade_double *b)
This function compares two arrays of two based on the fifth number, sorting highest first.
Definition: ProSHADE_misc.cpp:266
ProSHADE_internal_symmetry::checkFittingAxisDualAndSave
bool checkFittingAxisDualAndSave(std::vector< proshade_unsign > *retGroup, std::vector< proshade_double * > *ret, proshade_unsign fold, proshade_double axX, proshade_double axY, proshade_double axZ, std::vector< proshade_double * > *prosp, proshade_double axErr, proshade_unsign noMatchesG1, proshade_double angle1, proshade_unsign noMatchesG2, proshade_double angle2, ProSHADE_internal_data::ProSHADE_data *dataObj)
This function takes a newly detected "missing" axis and tests it for belonging to the group,...
Definition: ProSHADE_symmetry.cpp:2747
ProSHADE_internal_symmetry::findMissingAxisPoints
std::vector< proshade_double * > findMissingAxisPoints(proshade_double xVal, proshade_double yVal, proshade_double zVal, ProSHADE_internal_data::ProSHADE_data *dataObj, proshade_double axErr)
This function searches for all the self-rotation map points conforming to the axis,...
Definition: ProSHADE_symmetry.cpp:1817
ProSHADE_internal_maths::compute3x3MatrixMultiplication
proshade_double * compute3x3MatrixMultiplication(proshade_double *mat1, proshade_double *mat2)
Function for computing a 3x3 matrix multiplication.
Definition: ProSHADE_maths.cpp:1827