ProSHADE  0.7.6.1 (AUG 2021)
Protein Shape Detection
ProSHADE_data.cpp
Go to the documentation of this file.
1 
24 //==================================================== ProSHADE
25 #include "ProSHADE_data.hpp"
26 
27 //==================================================== Do not use the following flags for the included files - this causes a lot of warnings that have nothing to do with ProSHADE
28 #if defined ( __GNUC__ )
29  #pragma GCC diagnostic push
30  #pragma GCC diagnostic ignored "-Wpedantic"
31  #pragma GCC diagnostic ignored "-Wshadow"
32  #pragma GCC diagnostic ignored "-Wall"
33  #pragma GCC diagnostic ignored "-Wextra"
34  #pragma GCC diagnostic ignored "-Wdouble-promotion"
35  #pragma GCC diagnostic ignored "-Wconversion"
36 #endif
37 
38 //==================================================== Do not use the following flags for the included files - this causes a lot of warnings that have nothing to do with ProSHADE
39 #if defined ( __clang__ )
40  #pragma clang diagnostic push
41  #pragma clang diagnostic ignored "-Wpedantic"
42  #pragma clang diagnostic ignored "-Wshadow"
43  #pragma clang diagnostic ignored "-Wall"
44  #pragma clang diagnostic ignored "-Wextra"
45  #pragma clang diagnostic ignored "-Wdouble-promotion"
46  #pragma clang diagnostic ignored "-Weverything"
47 #endif
48 
49 //==================================================== Remove MSVC C4996 Warnings caused by Gemmi code
50 #if defined ( _MSC_VER )
51  #pragma warning ( disable:4996 )
52 #endif
53 
54 //==================================================== Gemmi PDB output - this cannot be with the rest of includes for some stb_sprintf library related reasons ...
55 #define GEMMI_WRITE_IMPLEMENTATION
56 #include <gemmi/to_pdb.hpp>
57 
58 //==================================================== Enable MSVC C4996 Warnings for the rest of the code
59 #if defined ( _MSC_VER )
60  #pragma warning ( default:4996 )
61 #endif
62 
63 //==================================================== Now the flags can be restored and used as per the CMakeLists.txt file.
64 #if defined ( __GNUC__ )
65  #pragma GCC diagnostic pop
66 #endif
67 
68 //==================================================== Now the flags can be restored and used as per the CMakeLists.txt file.
69 #if defined ( __clang__ )
70  #pragma clang diagnostic pop
71 #endif
72 
73 //==================================================== Forward declarations
75 {
76  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 );
77 }
78 
79 //==================================================== Local functions prototypes
80 void axesToGroupTypeSanityCheck ( proshade_unsign requiredAxes, proshade_unsign obtainedAxes, std::string groupType );
81 bool checkElementAlreadyExists ( std::vector<std::vector< proshade_double > >* elements, std::vector< proshade_double >* elem, proshade_double matrixTolerance );
82 bool checkElementsFormGroup ( std::vector<std::vector< proshade_double > >* elements, proshade_double matrixTolerance );
83 bool sortProSHADESymmetryByFSC ( proshade_double* a, proshade_double* b );
84 
94 {
95  //================================================ Initialise variables
96  // ... Variables regarding input file
97  this->fileName = "";
98  this->fileType = ProSHADE_internal_io::UNKNOWN;
99 
100  // ... Variables regarding map
101  this->internalMap = nullptr;
102 
103  // ... Variables regarding map information
104  this->xDimSize = 0.0;
105  this->yDimSize = 0.0;
106  this->zDimSize = 0.0;
107  this->aAngle = 0.0;
108  this->bAngle = 0.0;
109  this->cAngle = 0.0;
110  this->xDimIndices = 0;
111  this->yDimIndices = 0;
112  this->zDimIndices = 0;
113  this->xGridIndices = 0;
114  this->yGridIndices = 0;
115  this->zGridIndices = 0;
116  this->xAxisOrder = 1;
117  this->yAxisOrder = 2;
118  this->zAxisOrder = 3;
119  this->xAxisOrigin = 0;
120  this->yAxisOrigin = 0;
121  this->zAxisOrigin = 0;
122  this->xCom = 0.0;
123  this->yCom = 0.0;
124  this->zCom = 0.0;
125 
126  // ... Variables regarding original input values (i.e. these do not change with ProSHADE manipulations)
127  this->xDimSizeOriginal = 0.0;
128  this->yDimSizeOriginal = 0.0;
129  this->zDimSizeOriginal = 0.0;
130  this->xDimIndicesOriginal = 0;
131  this->yDimIndicesOriginal = 0;
132  this->zDimIndicesOriginal = 0;
133  this->xAxisOriginOriginal = 0;
134  this->yAxisOriginOriginal = 0;
135  this->zAxisOriginOriginal = 0;
136  this->originalMapXCom = 0.0;
137  this->originalMapYCom = 0.0;
138  this->originalMapZCom = 0.0;
139  this->mapMovFromsChangeX = 0.0;
140  this->mapMovFromsChangeY = 0.0;
141  this->mapMovFromsChangeZ = 0.0;
142  this->mapCOMProcessChangeX = 0.0;
143  this->mapCOMProcessChangeY = 0.0;
144  this->mapCOMProcessChangeZ = 0.0;
145 
146  // ... Variables regarding rotation and translation of original input files
147  this->originalPdbRotCenX = 0.0;
148  this->originalPdbRotCenY = 0.0;
149  this->originalPdbRotCenZ = 0.0;
150  this->originalPdbTransX = 0.0;
151  this->originalPdbTransY = 0.0;
152  this->originalPdbTransZ = 0.0;
153 
154  // ... Variables regarding iterator positions
155  this->xFrom = 0;
156  this->yFrom = 0;
157  this->zFrom = 0;
158  this->xTo = 0;
159  this->yTo = 0;
160  this->zTo = 0;
161 
162  // ... Variables regarding SH mapping spheres
163  this->spherePos = std::vector<proshade_single> ( );
164  this->noSpheres = 0;
165  this->spheres = nullptr;
166  this->sphericalHarmonics = nullptr;
167  this->rotSphericalHarmonics = nullptr;
168  this->maxShellBand = 0;
169 
170  // ... Variables regarding shape distance computations
171  this->rrpMatrices = nullptr;
172  this->eMatrices = nullptr;
173  this->so3Coeffs = nullptr;
174  this->so3CoeffsInverse = nullptr;
175  this->wignerMatrices = nullptr;
176  this->integrationWeight = 0.0;
177  this->maxCompBand = 0;
178  this->translationMap = nullptr;
179 
180 
181  // ... Control variables
182  this->isEmpty = true;
183 
184  //================================================ Done
185 
186 }
187 
214 ProSHADE_internal_data::ProSHADE_data::ProSHADE_data ( std::string strName, double *mapVals, int len, proshade_single xDmSz, proshade_single yDmSz, proshade_single zDmSz, proshade_unsign xDmInd, proshade_unsign yDmInd, proshade_unsign zDmInd, proshade_signed xFr, proshade_signed yFr, proshade_signed zFr, proshade_signed xT, proshade_signed yT, proshade_signed zT, proshade_unsign inputO )
215 {
216  //================================================ Initialise variables
217  // ... Variables regarding input file
218  this->fileName = strName;
219  this->fileType = ProSHADE_internal_io::MAP;
220 
221  // ... Variables regarding map
222  this->internalMap = nullptr;
223 
224  // ... Variables regarding map information
225  this->xDimSize = xDmSz;
226  this->yDimSize = yDmSz;
227  this->zDimSize = zDmSz;
228  this->aAngle = 90.0;
229  this->bAngle = 90.0;
230  this->cAngle = 90.0;
231  this->xDimIndices = xDmInd;
232  this->yDimIndices = yDmInd;
233  this->zDimIndices = zDmInd;
234  this->xGridIndices = xDmInd;
235  this->yGridIndices = yDmInd;
236  this->zGridIndices = zDmInd;
237  this->xAxisOrder = 1;
238  this->yAxisOrder = 2;
239  this->zAxisOrder = 3;
240  this->xAxisOrigin = xFr;
241  this->yAxisOrigin = yFr;
242  this->zAxisOrigin = zFr;
243  this->xCom = 0.0;
244  this->yCom = 0.0;
245  this->zCom = 0.0;
246 
247  // ... Variables regarding original input values (i.e. these do not change with ProSHADE manipulations)
248  this->xDimSizeOriginal = 0.0;
249  this->yDimSizeOriginal = 0.0;
250  this->zDimSizeOriginal = 0.0;
251  this->xDimIndicesOriginal = 0;
252  this->yDimIndicesOriginal = 0;
253  this->zDimIndicesOriginal = 0;
254  this->xAxisOriginOriginal = 0;
255  this->yAxisOriginOriginal = 0;
256  this->zAxisOriginOriginal = 0;
257  this->originalMapXCom = 0.0;
258  this->originalMapYCom = 0.0;
259  this->originalMapZCom = 0.0;
260  this->mapMovFromsChangeX = 0.0;
261  this->mapMovFromsChangeY = 0.0;
262  this->mapMovFromsChangeZ = 0.0;
263  this->mapCOMProcessChangeX = 0.0;
264  this->mapCOMProcessChangeY = 0.0;
265  this->mapCOMProcessChangeZ = 0.0;
266 
267  // ... Variables regarding rotation and translation of original input files
268  this->originalPdbRotCenX = 0.0;
269  this->originalPdbRotCenY = 0.0;
270  this->originalPdbRotCenZ = 0.0;
271  this->originalPdbTransX = 0.0;
272  this->originalPdbTransY = 0.0;
273  this->originalPdbTransZ = 0.0;
274 
275  // ... Variables regarding iterator positions
276  this->xFrom = xFr;
277  this->yFrom = yFr;
278  this->zFrom = zFr;
279  this->xTo = xT;
280  this->yTo = yT;
281  this->zTo = zT;
282 
283  // ... Variables regarding SH mapping spheres
284  this->spherePos = std::vector<proshade_single> ( );
285  this->noSpheres = 0;
286  this->spheres = nullptr;
287  this->sphericalHarmonics = nullptr;
288  this->rotSphericalHarmonics = nullptr;
289  this->maxShellBand = 0;
290 
291  // ... Variables regarding shape distance computations
292  this->rrpMatrices = nullptr;
293  this->eMatrices = nullptr;
294  this->so3Coeffs = nullptr;
295  this->so3CoeffsInverse = nullptr;
296  this->wignerMatrices = nullptr;
297  this->integrationWeight = 0.0;
298  this->maxCompBand = 0;
299  this->translationMap = nullptr;
300 
301  // ... Control variables
302  this->isEmpty = false;
303  this->inputOrder = inputO;
304 
305  //================================================ Sanity checks
306  if ( static_cast<proshade_unsign> ( len ) != ( xDmInd * yDmInd * zDmInd ) )
307  {
308  throw ProSHADE_exception ( "Structure class input map has wrong dimensions.", "EP00044", __FILE__, __LINE__, __func__, "The supplied map array size has different dimensions to\n : the required map dimensions." );
309  }
310 
311  if ( ( static_cast<proshade_signed> ( xT - xFr ) != static_cast<proshade_signed> ( xDmInd - 1 ) ) ||
312  ( static_cast<proshade_signed> ( yT - yFr ) != static_cast<proshade_signed> ( yDmInd - 1 ) ) ||
313  ( static_cast<proshade_signed> ( zT - zFr ) != static_cast<proshade_signed> ( zDmInd - 1 ) ) )
314  {
315  throw ProSHADE_exception ( "Structure class input dimensions not in line with map\n : to/from indices.", "EP00045", __FILE__, __LINE__, __func__, "The supplied map information does not add up. The\n : dimensions are not in line with the indexing start/stop\n : position distances and therefore proper map indexing\n : cannot be done. Please check the input values." );
316  }
317 
318  //================================================ Allocate the map memory
319  this->internalMap = new proshade_double [this->xDimIndices * this->yDimIndices * this->zDimIndices];
320  ProSHADE_internal_misc::checkMemoryAllocation ( this->internalMap, __FILE__, __LINE__, __func__ );
321 
322  //================================================ Copy the values into the map
323  proshade_unsign arrPos = 0;
324  for ( proshade_unsign xIt = 0; xIt < this->xDimIndices; xIt++ )
325  {
326  for ( proshade_unsign yIt = 0; yIt < this->yDimIndices; yIt++ )
327  {
328  for ( proshade_unsign zIt = 0; zIt < this->zDimIndices; zIt++ )
329  {
330  arrPos = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
331  this->internalMap[arrPos] = static_cast<proshade_double> ( mapVals[arrPos] );
332  }
333  }
334  }
335 
336  //================================================ Release memory (it was allocated by the PyBind11 lambda function and needs to be released)
337  delete[] mapVals;
338 
339  //================================================ Done
340 
341 }
342 
350 {
351  //================================================ Release the internal map
352  if ( this->internalMap != nullptr )
353  {
354  delete[] this->internalMap;
355  }
356 
357  //================================================ Release the sphere mapping
358  if ( this->spheres != nullptr )
359  {
360  for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
361  {
362  if ( this->spheres[iter] != nullptr )
363  {
364  delete this->spheres[iter];
365  this->spheres[iter] = nullptr;
366  }
367  }
368  delete[] this->spheres;
369  }
370 
371  //================================================ Release the spherical harmonics
372  if ( this->sphericalHarmonics != nullptr )
373  {
374  for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
375  {
376  if ( this->sphericalHarmonics[iter] != nullptr )
377  {
378  delete[] this->sphericalHarmonics[iter];
379  this->sphericalHarmonics[iter] = nullptr;
380  }
381  }
382  delete[] this->sphericalHarmonics;
383  }
384 
385  //================================================ Release the rotated spherical harmonics
386  if ( this->rotSphericalHarmonics != nullptr )
387  {
388  for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
389  {
390  if ( this->rotSphericalHarmonics[iter] != nullptr )
391  {
392  delete[] this->rotSphericalHarmonics[iter];
393  this->rotSphericalHarmonics[iter] = nullptr;
394  }
395  }
396  delete[] this->rotSphericalHarmonics;
397  }
398 
399  //================================================ Release the RRP matrices (pre-computation for the energy levels descriptor)
400  if ( this->rrpMatrices != nullptr )
401  {
402  for ( proshade_unsign bwIt = 0; bwIt < this->maxShellBand; bwIt++ )
403  {
404  if ( this->rrpMatrices[bwIt] != nullptr )
405  {
406  for ( proshade_unsign shIt = 0; shIt < this->noSpheres; shIt++ )
407  {
408  if ( this->rrpMatrices[bwIt][shIt] != nullptr )
409  {
410  delete[] this->rrpMatrices[bwIt][shIt];
411  }
412  }
413 
414  delete[] this->rrpMatrices[bwIt];
415  }
416  }
417 
418  delete[] this->rrpMatrices;
419  }
420 
421  //================================================ Release the E matrices
422  if ( this->eMatrices != nullptr )
423  {
424  for ( proshade_unsign bandIter = 0; bandIter < this->maxCompBand; bandIter++ )
425  {
426  if ( this->eMatrices[bandIter] != nullptr )
427  {
428  for ( proshade_unsign band2Iter = 0; band2Iter < static_cast<proshade_unsign> ( ( bandIter * 2 ) + 1 ); band2Iter++ )
429  {
430  if ( this->eMatrices[bandIter][band2Iter] != nullptr )
431  {
432  delete[] this->eMatrices[bandIter][band2Iter];
433  }
434  }
435 
436  delete[] this->eMatrices[bandIter];
437  }
438  }
439 
440  delete[] this->eMatrices;
441  }
442 
443  //================================================ Release SOFT and inverse SOFT coefficients
444  if ( this->so3Coeffs != nullptr )
445  {
446  delete[] this->so3Coeffs;
447  }
448  if ( this->so3CoeffsInverse != nullptr )
449  {
450  delete[] this->so3CoeffsInverse;
451  }
452 
453  //================================================ Release Wigner matrices
454  if ( this->wignerMatrices != nullptr )
455  {
456  for ( proshade_unsign bandIter = 1; bandIter < this->maxCompBand; bandIter++ )
457  {
458  if ( this->wignerMatrices[bandIter] != nullptr )
459  {
460  for ( proshade_unsign order1Iter = 0; order1Iter < ( (bandIter * 2) + 1 ); order1Iter++ )
461  {
462  if ( this->wignerMatrices[bandIter][order1Iter] != nullptr )
463  {
464  delete[] this->wignerMatrices[bandIter][order1Iter];
465  }
466  }
467  delete[] this->wignerMatrices[bandIter];
468  }
469  }
470  delete[] wignerMatrices;
471  }
472 
473  //================================================ Release translation map
474  if ( this->translationMap != nullptr )
475  {
476  delete[] this->translationMap;
477  }
478 
479  //================================================ Release the angle-axis space rotation function
480  if ( this->sphereMappedRotFun.size() > 0 )
481  {
482  for ( proshade_unsign spIt = 0; spIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.size() ); spIt++ )
483  {
484  delete this->sphereMappedRotFun.at(spIt);
485  }
486  }
487 
488  //================================================ Done
489 
490 }
491 
509 void ProSHADE_internal_data::ProSHADE_data::readInStructure ( std::string fName, proshade_unsign inputO, ProSHADE_settings* settings, proshade_double* maskArr, proshade_unsign maskXDim, proshade_unsign maskYDim, proshade_unsign maskZDim, proshade_double* weightsArr, proshade_unsign weigXDim, proshade_unsign weigYDim, proshade_unsign weigZDim )
510 {
511  //================================================ Report function start
512  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting to read the structure: " + fName );
513 
514  //================================================ Check if instance is empty
515  if ( !this->isEmpty )
516  {
517  throw ProSHADE_exception ( "Structure data class not empty.", "E000005", __FILE__, __LINE__, __func__, "Attempted to read in structure into a ProSHADE_data\n : object which already does have structure read in\n : i.e. " + this->fileName );
518  }
519 
520  //================================================ Save the filename
521  this->fileName = fName;
522 
523  //================================================ Check what is the input format
524  this->fileType = ProSHADE_internal_io::figureDataType ( this->fileName );
525 
526  //================================================ Save input order
527  this->inputOrder = inputO;
528 
529  //================================================ Decide how to proceed
530  switch ( this->fileType )
531  {
532  case ProSHADE_internal_io::UNKNOWN:
533  throw ProSHADE_exception ( "Unknown file type.", "E000006", __FILE__, __LINE__, __func__, "When attempting to read the file\n : " + this->fileName + "\n : the file extension was determined as unknown. This could\n : mean either that the file does not exist, or that it is\n : not one of the supported extensions." );
534 
535  case ProSHADE_internal_io::GEMMI:
536  throw ProSHADE_exception ( "Unknown file type.", "E000006", __FILE__, __LINE__, __func__, "When attempting to read the file\n : " + this->fileName + "\n : the file extension was determined as unknown. This could\n : mean either that the file does not exist, or that it is\n : not one of the supported extensions." );
537 
538  case ProSHADE_internal_io::PDB:
539  this->readInPDB ( settings );
540  break;
541 
542  case ProSHADE_internal_io::MAP:
543  this->readInMAP ( settings, maskArr, maskXDim, maskYDim, maskZDim, weightsArr, weigXDim, weigYDim, weigZDim );
544  break;
545  }
546 
547  //================================================ This structure is now full
548  this->isEmpty = false;
549 
550  //================================================ Report function completion
551  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Structure read in successfully." );
552 
553  //================================================ Done
554  return ;
555 
556 }
557 
567 void ProSHADE_internal_data::ProSHADE_data::readInStructure ( gemmi::Structure gemmiStruct, proshade_unsign inputO, ProSHADE_settings* settings )
568 {
569  //================================================ Report function start
570  std::stringstream ss;
571  ss << "Starting to load the structure from Gemmi object " << inputO;
573 
574  //================================================ Check if instance is empty
575  if ( !this->isEmpty )
576  {
577  throw ProSHADE_exception ( "Structure data class not empty.", "E000005", __FILE__, __LINE__, __func__, "Attempted to read in structure into a ProSHADE_data\n : object which already does have structure read in\n : i.e. " + this->fileName );
578  }
579 
580  //================================================ Save the filename
581  this->fileName = gemmiStruct.name;
582 
583  //================================================ Check what is the input format
584  this->fileType = ProSHADE_internal_io::GEMMI;
585 
586  //================================================ Save input order
587  this->inputOrder = inputO;
588 
589  //================================================ Decide how to proceed
590  this->readInGemmi ( gemmiStruct, settings );
591 
592  //================================================ This structure is now full
593  this->isEmpty = false;
594 
595  //================================================ Report function completion
596  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Structure read in successfully." );
597 
598  //================================================ Done
599  return ;
600 
601 }
602 
619 void ProSHADE_internal_data::ProSHADE_data::readInMAP ( ProSHADE_settings* settings, proshade_double* maskArr, proshade_unsign maskXDim, proshade_unsign maskYDim, proshade_unsign maskZDim, proshade_double* weightsArr, proshade_unsign weigXDim, proshade_unsign weigYDim, proshade_unsign weigZDim )
620 {
621  //================================================ Open the file
622  gemmi::Ccp4<float> map;
623  map.read_ccp4 ( gemmi::MaybeGzipped ( this->fileName.c_str() ) );
624 
625  //================================================ Convert to XYZ and create complete map, if need be
626  map.setup ( gemmi::GridSetup::ReorderOnly, NAN );
627 
628  //================================================ Read in the rest of the map file header
630  &this->xDimIndices, &this->yDimIndices, &this->zDimIndices,
631  &this->xDimSize, &this->yDimSize, &this->zDimSize,
632  &this->aAngle, &this->bAngle, &this->cAngle,
633  &this->xFrom, &this->yFrom, &this->zFrom,
634  &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin,
635  &this->xAxisOrder, &this->yAxisOrder, &this->zAxisOrder,
636  &this->xGridIndices, &this->yGridIndices, &this->zGridIndices );
637 
638  //================================================ Save the map density to ProSHADE variable
639  ProSHADE_internal_io::readInMapData ( &map, this->internalMap, this->xDimIndices, this->yDimIndices, this->zDimIndices, this->xAxisOrder, this->yAxisOrder, this->zAxisOrder );
640 
641  //================================================ If mask is supplied and the correct task is used
642  ProSHADE_internal_io::applyMask ( this->internalMap, settings->appliedMaskFileName, this->xDimIndices, this->yDimIndices, this->zDimIndices, settings->verbose,
643  maskArr, maskXDim, maskYDim, maskZDim );
644 
645  //================================================ Apply Fourier weights
646  ProSHADE_internal_io::applyWeights ( this->internalMap, settings->fourierWeightsFileName, this->xDimIndices, this->yDimIndices, this->zDimIndices, settings->verbose,
647  weightsArr, weigXDim, weigYDim, weigZDim );
648 
649  //================================================ Remove negative values if so required
650  if ( settings->removeNegativeDensity ) { for ( size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { if ( this->internalMap[iter] < 0.0 ) { this->internalMap[iter] = 0.0; } } }
651 
652  //================================================ Set resolution if need be
653  if ( settings->requestedResolution < 0.0f )
654  {
655  settings->setResolution ( std::min ( static_cast<proshade_single> ( this->xDimSize ) / static_cast<proshade_single> ( this->xDimIndices ),
656  std::min ( static_cast<proshade_single> ( this->yDimSize ) / static_cast<proshade_single> ( this->yDimIndices ),
657  static_cast<proshade_single> ( this->zDimSize ) / static_cast<proshade_single> ( this->zDimIndices ) ) ) * 2.0f );
658  }
659 
660  //================================================ Set iterators from and to
661  this->figureIndexStartStop ( );
662 
663  //================================================ If specific resolution is requested, make sure the map has it
664  if ( settings->changeMapResolution || settings->changeMapResolutionTriLinear )
665  {
666  //============================================ Before re-sampling sampling rate
667  proshade_single xSampRate = this->xDimSize / static_cast< proshade_single > ( this->xTo - this->xFrom );
668  proshade_single ySampRate = this->yDimSize / static_cast< proshade_single > ( this->yTo - this->yFrom );
669  proshade_single zSampRate = this->zDimSize / static_cast< proshade_single > ( this->zTo - this->zFrom );
670 
671  //============================================ Bofore re-sampling first index position
672  proshade_single xStartPosBefore = static_cast< proshade_single > ( this->xFrom ) * xSampRate;
673  proshade_single yStartPosBefore = static_cast< proshade_single > ( this->yFrom ) * ySampRate;
674  proshade_single zStartPosBefore = static_cast< proshade_single > ( this->zFrom ) * zSampRate;
675 
676  //============================================ Find COM before map re-sampling
677  proshade_double xMapCOMPreReSampl = 0.0, yMapCOMPreReSampl = 0.0, zMapCOMPreReSampl = 0.0;
678  ProSHADE_internal_mapManip::findMAPCOMValues ( this->internalMap, &xMapCOMPreReSampl, &yMapCOMPreReSampl, &zMapCOMPreReSampl, this->xDimSize, this->yDimSize, this->zDimSize, this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
679 
680  //============================================ Re-sample map
681  this->reSampleMap ( settings );
682 
683  //============================================ After re-sampling sampling rate
684  xSampRate = this->xDimSize / static_cast< proshade_single > ( this->xTo - this->xFrom );
685  ySampRate = this->yDimSize / static_cast< proshade_single > ( this->yTo - this->yFrom );
686  zSampRate = this->zDimSize / static_cast< proshade_single > ( this->zTo - this->zFrom );
687 
688  //============================================ After re-sampling first index position
689  proshade_single xStartPosAfter = static_cast< proshade_single > ( this->xFrom ) * xSampRate;
690  proshade_single yStartPosAfter = static_cast< proshade_single > ( this->yFrom ) * ySampRate;
691  proshade_single zStartPosAfter = static_cast< proshade_single > ( this->zFrom ) * zSampRate;
692 
693  //============================================ Translate by change in corners to make the boxes as similarly placed as possible
694  proshade_single xMov = static_cast< proshade_single > ( xStartPosAfter - xStartPosBefore );
695  proshade_single yMov = static_cast< proshade_single > ( yStartPosAfter - yStartPosBefore );
696  proshade_single zMov = static_cast< proshade_single > ( zStartPosAfter - zStartPosBefore );
697  ProSHADE_internal_mapManip::moveMapByIndices ( &xMov, &yMov, &zMov, this->xDimSize, this->yDimSize, this->zDimSize,
698  &this->xFrom, &this->xTo, &this->yFrom, &this->yTo, &this->zFrom, &this->zTo,
699  &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin );
700 
701  //============================================ Find COM after map re-sampling and corner move
702  proshade_double xMapCOMPostReSampl = 0.0, yMapCOMPostReSampl = 0.0, zMapCOMPostReSampl = 0.0;
703  ProSHADE_internal_mapManip::findMAPCOMValues ( this->internalMap, &xMapCOMPostReSampl, &yMapCOMPostReSampl, &zMapCOMPostReSampl, this->xDimSize, this->yDimSize, this->zDimSize, this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
704 
705  //============================================ Match the COMs to get as close position of the re-sampled structure to the original as possible
707  static_cast< proshade_single > ( xMapCOMPreReSampl - xMapCOMPostReSampl ),
708  static_cast< proshade_single > ( yMapCOMPreReSampl - yMapCOMPostReSampl ),
709  static_cast< proshade_single > ( zMapCOMPreReSampl - zMapCOMPostReSampl ),
710  this->xDimSize, this->yDimSize, this->zDimSize,
711  static_cast< proshade_signed > ( this->xDimIndices ),
712  static_cast< proshade_signed > ( this->yDimIndices ),
713  static_cast< proshade_signed > ( this->zDimIndices ) );
714  }
715 
716  //================================================ Save the original sizes
717  this->xDimSizeOriginal = this->xDimSize;
718  this->yDimSizeOriginal = this->yDimSize;
719  this->zDimSizeOriginal = this->zDimSize;
720 
721  //================================================ Save the original index counts
722  this->xDimIndicesOriginal = this->xDimIndices;
723  this->yDimIndicesOriginal = this->yDimIndices;
724  this->zDimIndicesOriginal = this->zDimIndices;
725 
726  //================================================ Save the original axis origins
727  this->xAxisOriginOriginal = this->xAxisOrigin;
728  this->yAxisOriginOriginal = this->yAxisOrigin;
729  this->zAxisOriginOriginal = this->zAxisOrigin;
730 
731  //================================================ Compute and save the COM
732  this->findMapCOM ( );
733  this->originalMapXCom = this->xCom;
734  this->originalMapYCom = this->yCom;
735  this->originalMapZCom = this->zCom;
736 
737  //================================================ Done
738 
739 }
740 
751 {
752  //================================================ Set resolution if need be
753  if ( settings->requestedResolution < 0.0f )
754  {
755  settings->setResolution ( 8.0 );
756  }
757 
758  //================================================ Open PDB file for reading
759  gemmi::Structure pdbFile = gemmi::read_structure ( gemmi::MaybeGzipped ( this->fileName ) );
760 
761  //================================================ Once we have Gemmi object, run the Gemmi function
762  this->readInGemmi ( pdbFile, settings );
763 
764  //================================================ Done
765  return;
766 
767 }
768 
785 void ProSHADE_internal_data::ProSHADE_data::readInGemmi ( gemmi::Structure gemmiStruct, ProSHADE_settings* settings )
786 {
787  //================================================ Set resolution if need be
788  if ( settings->requestedResolution < 0.0f )
789  {
790  settings->setResolution ( 8.0 );
791  }
792 
793  //================================================ Change B-factors if need be
794  if ( settings->pdbBFactorNewVal >= 0.0 )
795  {
797  }
798 
799  //================================================ Remove waters if required
800  if ( settings->removeWaters )
801  {
802  ProSHADE_internal_mapManip::removeWaters ( &gemmiStruct, settings->firstModelOnly );
803  }
804 
805  //================================================ Get PDB COM values
806  proshade_double xCOMPdb, yCOMPdb, zCOMPdb;
807  ProSHADE_internal_mapManip::findPDBCOMValues ( gemmiStruct, &xCOMPdb, &yCOMPdb, &zCOMPdb, settings->firstModelOnly );
808 
809  //================================================ Find the ranges
810  proshade_single xF, xT, yF, yT, zF, zT;
811  ProSHADE_internal_mapManip::determinePDBRanges ( gemmiStruct, &xF, &xT, &yF, &yT, &zF, &zT, settings->firstModelOnly );
812 
813  //================================================ Move ranges to have all FROM values 20
814  proshade_single xMov = static_cast< proshade_single > ( 20.0f - xF );
815  proshade_single yMov = static_cast< proshade_single > ( 20.0f - yF );
816  proshade_single zMov = static_cast< proshade_single > ( 20.0f - zF );
817  ProSHADE_internal_mapManip::movePDBForMapCalc ( &gemmiStruct, xMov, yMov, zMov, settings->firstModelOnly );
818 
819  //================================================ Set the angstrom sizes
820  this->xDimSize = static_cast< proshade_single > ( xT - xF + 40.0f );
821  this->yDimSize = static_cast< proshade_single > ( yT - yF + 40.0f );
822  this->zDimSize = static_cast< proshade_single > ( zT - zF + 40.0f );
823 
824  //================================================ Generate map from nicely placed atoms (cell size will be range + 40)
825  ProSHADE_internal_mapManip::generateMapFromPDB ( gemmiStruct, this->internalMap, settings->requestedResolution, this->xDimSize, this->yDimSize, this->zDimSize, &this->xTo, &this->yTo, &this->zTo, settings->forceP1, settings->firstModelOnly );
826 
827  //================================================ Remove negative values if so required
828  if ( settings->removeNegativeDensity ) { for ( size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { if ( this->internalMap[iter] < 0.0 ) { this->internalMap[iter] = 0.0; } } }
829 
830  //================================================ Set the internal variables to correct values
831  this->setPDBMapValues ( );
832 
833  //================================================ Compute reverse movement based on COMs. If there is more than 1 models, simply moving back the xyzMov is not enough.
834  proshade_double xCOMMap, yCOMMap, zCOMMap;
835  ProSHADE_internal_mapManip::findMAPCOMValues ( this->internalMap, &xCOMMap, &yCOMMap, &zCOMMap,
836  this->xDimSize, this->yDimSize, this->zDimSize,
837  this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
838 
839  if ( gemmiStruct.models.size() > 1 )
840  {
841  xMov = static_cast< proshade_single > ( xCOMMap - xCOMPdb );
842  yMov = static_cast< proshade_single > ( yCOMMap - yCOMPdb );
843  zMov = static_cast< proshade_single > ( zCOMMap - zCOMPdb );
844  }
845 
846  //================================================ Move map back to the original PDB location
847  ProSHADE_internal_mapManip::moveMapByIndices ( &xMov, &yMov, &zMov, this->xDimSize, this->yDimSize, this->zDimSize,
848  &this->xFrom, &this->xTo, &this->yFrom, &this->yTo, &this->zFrom, &this->zTo,
849  &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin );
850  ProSHADE_internal_mapManip::moveMapByFourier ( this->internalMap, xMov, yMov, zMov, this->xDimSize, this->yDimSize, this->zDimSize,
851  static_cast< proshade_signed > ( this->xDimIndices ), static_cast< proshade_signed > ( this->yDimIndices ),
852  static_cast< proshade_signed > ( this->zDimIndices ) );
853 
854  //================================================ If specific resolution is requested, make sure the map has it
855  if ( settings->changeMapResolution || settings->changeMapResolutionTriLinear )
856  {
857  //============================================ Before re-sampling sampling rate
858  proshade_single xSampRate = this->xDimSize / static_cast< proshade_single > ( this->xTo - this->xFrom );
859  proshade_single ySampRate = this->yDimSize / static_cast< proshade_single > ( this->yTo - this->yFrom );
860  proshade_single zSampRate = this->zDimSize / static_cast< proshade_single > ( this->zTo - this->zFrom );
861 
862  //============================================ Bofore re-sampling first index position
863  proshade_single xStartPosBefore = static_cast< proshade_single > ( this->xFrom ) * xSampRate;
864  proshade_single yStartPosBefore = static_cast< proshade_single > ( this->yFrom ) * ySampRate;
865  proshade_single zStartPosBefore = static_cast< proshade_single > ( this->zFrom ) * zSampRate;
866 
867  //============================================ Find COM before map re-sampling
868  proshade_double xMapCOMPreReSampl = 0.0, yMapCOMPreReSampl = 0.0, zMapCOMPreReSampl = 0.0;
869  ProSHADE_internal_mapManip::findMAPCOMValues ( this->internalMap, &xMapCOMPreReSampl, &yMapCOMPreReSampl, &zMapCOMPreReSampl, this->xDimSize, this->yDimSize, this->zDimSize, this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
870 
871  //============================================ Re-sample map
872  this->reSampleMap ( settings );
873 
874  //============================================ After re-sampling sampling rate
875  xSampRate = this->xDimSize / static_cast< proshade_single > ( this->xTo - this->xFrom );
876  ySampRate = this->yDimSize / static_cast< proshade_single > ( this->yTo - this->yFrom );
877  zSampRate = this->zDimSize / static_cast< proshade_single > ( this->zTo - this->zFrom );
878 
879  //============================================ After re-sampling first index position
880  proshade_single xStartPosAfter = static_cast< proshade_single > ( this->xFrom ) * xSampRate;
881  proshade_single yStartPosAfter = static_cast< proshade_single > ( this->yFrom ) * ySampRate;
882  proshade_single zStartPosAfter = static_cast< proshade_single > ( this->zFrom ) * zSampRate;
883 
884  //============================================ Translate by change in corners to make the boxes as similarly placed as possible
885  proshade_single xMovHlp = static_cast< proshade_single > ( xStartPosAfter - xStartPosBefore );
886  proshade_single yMovHlp = static_cast< proshade_single > ( yStartPosAfter - yStartPosBefore );
887  proshade_single zMovHlp = static_cast< proshade_single > ( zStartPosAfter - zStartPosBefore );
888  ProSHADE_internal_mapManip::moveMapByIndices ( &xMovHlp, &yMovHlp, &zMovHlp, this->xDimSize, this->yDimSize, this->zDimSize,
889  &this->xFrom, &this->xTo, &this->yFrom, &this->yTo, &this->zFrom, &this->zTo,
890  &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin );
891 
892  //============================================ Find COM after map re-sampling and corner move
893  proshade_double xMapCOMPostReSampl = 0.0, yMapCOMPostReSampl = 0.0, zMapCOMPostReSampl = 0.0;
894  ProSHADE_internal_mapManip::findMAPCOMValues ( this->internalMap, &xMapCOMPostReSampl, &yMapCOMPostReSampl, &zMapCOMPostReSampl, this->xDimSize, this->yDimSize, this->zDimSize, this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
895 
896  //============================================ Match the COMs to get as close position of the re-sampled structure to the original as possible
898  static_cast< proshade_single > ( xMapCOMPreReSampl - xMapCOMPostReSampl ),
899  static_cast< proshade_single > ( yMapCOMPreReSampl - yMapCOMPostReSampl ),
900  static_cast< proshade_single > ( zMapCOMPreReSampl - zMapCOMPostReSampl ),
901  this->xDimSize, this->yDimSize, this->zDimSize,
902  static_cast< proshade_signed > ( this->xDimIndices ),
903  static_cast< proshade_signed > ( this->yDimIndices ),
904  static_cast< proshade_signed > ( this->zDimIndices ) );
905  }
906 
907  //================================================ Save the original sizes
908  this->xDimSizeOriginal = this->xDimSize;
909  this->yDimSizeOriginal = this->yDimSize;
910  this->zDimSizeOriginal = this->zDimSize;
911 
912  //================================================ Save the original index counts
913  this->xDimIndicesOriginal = this->xDimIndices;
914  this->yDimIndicesOriginal = this->yDimIndices;
915  this->zDimIndicesOriginal = this->zDimIndices;
916 
917  //================================================ Save the original axis origins
918  this->xAxisOriginOriginal = this->xAxisOrigin;
919  this->yAxisOriginOriginal = this->yAxisOrigin;
920  this->zAxisOriginOriginal = this->zAxisOrigin;
921 
922  //================================================ Compute and save the COM
923  this->findMapCOM ( );
924  this->originalMapXCom = this->xCom;
925  this->originalMapYCom = this->yCom;
926  this->originalMapZCom = this->zCom;
927 
928  //================================================ Done
929  return;
930 
931 }
932 
938 {
939  //================================================ Set starts to 0
940  this->xFrom = 0;
941  this->yFrom = 0;
942  this->zFrom = 0;
943 
944  //================================================ Set angles to 90 degrees
945  this->aAngle = 90.0;
946  this->bAngle = 90.0;
947  this->cAngle = 90.0;
948 
949  //================================================ Set dimension sizes in indices
950  this->xDimIndices = static_cast< proshade_unsign > ( this->xTo );
951  this->yDimIndices = static_cast< proshade_unsign > ( this->yTo );
952  this->zDimIndices = static_cast< proshade_unsign > ( this->zTo );
953 
954  //================================================ Set the to indices properly
955  this->xTo -= 1;
956  this->yTo -= 1;
957  this->zTo -= 1;
958 
959  //================================================ Set grid indexing to cell indexing
960  this->xGridIndices = this->xDimIndices;
961  this->yGridIndices = this->yDimIndices;
962  this->zGridIndices = this->zDimIndices;
963 
964  //================================================ Set axis order
965  this->xAxisOrder = 1;
966  this->yAxisOrder = 2;
967  this->zAxisOrder = 3;
968 
969  //================================================ Set origin to the first index
970  this->xAxisOrigin = this->xFrom;
971  this->yAxisOrigin = this->yFrom;
972  this->zAxisOrigin = this->zFrom;
973 
974  //================================================ Done
975  return ;
976 
977 }
978 
984 {
985  //================================================ Set ends to origin + size - 1
986  this->xTo = this->xFrom + static_cast< proshade_signed > ( this->xDimIndices ) - 1;
987  this->yTo = this->yFrom + static_cast< proshade_signed > ( this->yDimIndices ) - 1;
988  this->zTo = this->zFrom + static_cast< proshade_signed > ( this->zDimIndices ) - 1;
989 
990  //================================================ Done
991  return ;
992 
993 }
994 
1005 void ProSHADE_internal_data::ProSHADE_data::writeMap ( std::string fName, std::string title, int mode )
1006 {
1007  //================================================ Create and prepare new Grid gemmi object
1008  gemmi::Grid<float> mapData;
1009  mapData.set_unit_cell ( static_cast< double > ( this->xDimSize ), static_cast< double > ( this->yDimSize ), static_cast< double > ( this->zDimSize ), static_cast< double > ( this->aAngle ), static_cast< double > ( this->bAngle ), static_cast< double > ( this->cAngle ) );
1010  mapData.set_size_without_checking ( static_cast< int > ( this->xDimIndices ), static_cast< int > ( this->yDimIndices ), static_cast< int > ( this->zDimIndices ) );
1011  mapData.axis_order = gemmi::AxisOrder::XYZ;
1012  mapData.spacegroup = &gemmi::get_spacegroup_p1();
1013 
1014  //================================================ Create and prepare new Ccp4 gemmi object
1015  gemmi::Ccp4<float> map;
1016  map.grid = mapData;
1017  map.update_ccp4_header ( mode );
1018 
1019  //================================================ Fill in the header
1021  this->xDimIndices, this->yDimIndices, this->zDimIndices,
1022  this->xDimSize, this->yDimSize, this->zDimSize,
1023  this->aAngle, this->bAngle, this->cAngle,
1024  this->xFrom, this->yFrom, this->zFrom,
1025  this->xAxisOrigin, this->yAxisOrigin, this->zAxisOrigin,
1026  this->xAxisOrder, this->yAxisOrder, this->zAxisOrder,
1027  this->xGridIndices, this->yGridIndices, this->zGridIndices,
1028  title, mode );
1029 
1030  //================================================ Copy internal map to grid
1031  proshade_unsign arrPos = 0;
1032  for ( proshade_unsign uIt = 0; uIt < this->xDimIndices; uIt++ )
1033  {
1034  for ( proshade_unsign vIt = 0; vIt < this->yDimIndices; vIt++ )
1035  {
1036  for ( proshade_unsign wIt = 0; wIt < this->zDimIndices; wIt++ )
1037  {
1038  arrPos = wIt + this->zDimIndices * ( vIt + this->yDimIndices * uIt );
1039  map.grid.set_value ( static_cast< int > ( uIt ), static_cast< int > ( vIt ), static_cast< int > ( wIt ), static_cast<float> ( this->internalMap[arrPos] ) );
1040  }
1041  }
1042  }
1043 
1044  //================================================ Update the statistics in the header
1045  map.update_ccp4_header ( mode, true );
1046 
1047  //================================================ Write out the map
1048  map.write_ccp4_map ( fName );
1049 
1050  //================================================ Done
1051  return ;
1052 
1053 }
1054 
1072 void ProSHADE_internal_data::ProSHADE_data::writePdb ( std::string fName, proshade_double euA, proshade_double euB, proshade_double euG, proshade_double trsX, proshade_double trsY, proshade_double trsZ, proshade_double rotX, proshade_double rotY, proshade_double rotZ, bool firstModel )
1073 {
1074  //================================================ Check for co-ordinate origin
1075  if ( !ProSHADE_internal_io::isFilePDB ( this->fileName ) )
1076  {
1077  throw ProSHADE_exception ( "Cannot write co-ordinate file if the input file did not contain co-ordinates.", "EP00047", __FILE__, __LINE__, __func__, "You have called the WritePDB function on structure which\n : was created by reading in a map. This is not allowed as\n : ProSHADE cannot create co-ordinates from map file." );
1078  }
1079 
1080  //================================================ Open PDB file for reading
1081  gemmi::Structure pdbFile = gemmi::read_structure ( gemmi::MaybeGzipped ( this->fileName ) );
1082 
1083  //================================================ Write out using the gemmi::Structure object
1084  this->writeGemmi ( fName, pdbFile, euA, euB, euG, trsX, trsY, trsZ, rotX, rotY, rotZ, firstModel );
1085 
1086  //================================================ Done
1087  return ;
1088 }
1089 
1109 void ProSHADE_internal_data::ProSHADE_data::writeGemmi ( std::string fName, gemmi::Structure gemmiStruct, proshade_double euA, proshade_double euB, proshade_double euG, proshade_double trsX, proshade_double trsY, proshade_double trsZ, proshade_double rotX, proshade_double rotY, proshade_double rotZ, bool firstModel )
1110 {
1111  //================================================ If the map was rotated, do the same for the co-ordinates, making sure we take into account the rotation centre of the map
1112  if ( ( euA != 0.0 ) || ( euB != 0.0 ) || ( euG != 0.0 ) )
1113  {
1114  //============================================ Rotate the co-ordinates
1115  ProSHADE_internal_mapManip::rotatePDBCoordinates ( &gemmiStruct, euA, euB, euG, rotX, rotY, rotZ, firstModel );
1116  }
1117 
1118  //================================================ Translate by required translation and the map centering (if applied)
1119  ProSHADE_internal_mapManip::translatePDBCoordinates ( &gemmiStruct, trsX, trsY, trsZ, firstModel );
1120 
1121  //================================================ Write the PDB file
1122  std::ofstream outCoOrdFile;
1123  outCoOrdFile.open ( fName.c_str() );
1124 
1125  if ( outCoOrdFile.is_open() )
1126  {
1127  gemmi::PdbWriteOptions opt;
1128  write_pdb ( gemmiStruct, outCoOrdFile, opt );
1129  }
1130  else
1131  {
1132  std::stringstream hlpMessage;
1133  hlpMessage << "Failed to open the PDB file " << fName << " for output.";
1134  throw ProSHADE_exception ( hlpMessage.str().c_str(), "EP00048", __FILE__, __LINE__, __func__, "ProSHADE has failed to open the PDB output file. This is\n : likely caused by either not having the write privileges\n : to the required output path, or by making a mistake in\n : the path." );
1135  }
1136 
1137  outCoOrdFile.close ( );
1138 
1139  //================================================ Done
1140  return ;
1141 }
1142 
1151 void ProSHADE_internal_data::ProSHADE_data::writeMask ( std::string fName, proshade_double* mask )
1152 {
1153  //================================================ Allocate the memory
1154  proshade_double* hlpMap = new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1155  ProSHADE_internal_misc::checkMemoryAllocation ( hlpMap, __FILE__, __LINE__, __func__ );
1156 
1157  //================================================ Copy original map and over-write with the mask
1158  for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1159  {
1160  hlpMap[iter] = this->internalMap[iter];
1161  this->internalMap[iter] = mask[iter];
1162  }
1163 
1164  //================================================ Write out the mask
1165  this->writeMap ( fName );
1166 
1167  //================================================ Copy the original map values back
1168  for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1169  {
1170  this->internalMap[iter] = hlpMap[iter];
1171  }
1172 
1173  //================================================ Release memory
1174  delete[] hlpMap;
1175 
1176  //================================================ Done
1177  return ;
1178 }
1179 
1189 {
1190  //================================================ Report function start
1191  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Map inversion." );
1192 
1193  //================================================ Initialise variables
1194  proshade_signed arrayPos, invPos;
1195 
1196  //================================================ Create helper map
1197  proshade_double* hlpMap = new proshade_double [this->xDimIndices * this->yDimIndices * this->zDimIndices];
1198  ProSHADE_internal_misc::checkMemoryAllocation ( hlpMap, __FILE__, __LINE__, __func__ );
1199 
1200  //================================================ Save map values to the helper map
1201  for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1202  {
1203  hlpMap[iter] = this->internalMap[iter];
1204  }
1205 
1206  //================================================ Invert the values
1207  for ( proshade_signed xIt = 0; xIt < static_cast<proshade_signed> ( this->xDimIndices ); xIt++ )
1208  {
1209  for ( proshade_signed yIt = 0; yIt < static_cast<proshade_signed> ( this->yDimIndices ); yIt++ )
1210  {
1211  for ( proshade_signed zIt = 0; zIt < static_cast<proshade_signed> ( this->zDimIndices ); zIt++ )
1212  {
1213  //==================================== Var init
1214  arrayPos = zIt + static_cast< proshade_signed > ( this->zDimIndices ) * ( yIt + static_cast< proshade_signed > ( this->yDimIndices ) * xIt );
1215  invPos = ( static_cast< proshade_signed > ( this->zDimIndices - 1 ) - zIt ) + static_cast< proshade_signed > ( this->zDimIndices ) * ( ( static_cast< proshade_signed > ( this->yDimIndices - 1 ) - yIt ) + static_cast< proshade_signed > ( this->yDimIndices ) * ( static_cast< proshade_signed > ( this->xDimIndices - 1 ) - xIt ) );
1216 
1217  //==================================== And save
1218  this->internalMap[invPos] = hlpMap[arrayPos];
1219  }
1220  }
1221  }
1222 
1223  //================================================ Release memory
1224  delete[] hlpMap;
1225 
1226  //================================================ Report function completion
1227  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Map inversion completed." );
1228 
1229  //================================================ Done
1230  return ;
1231 
1232 }
1233 
1243 {
1244  //================================================ Report function start
1245  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Map normalisation." );
1246 
1247  //================================================ Initialise vector of map values
1248  std::vector<proshade_double> mapVals ( this->xDimIndices * this->yDimIndices * this->zDimIndices, 0.0 );
1249 
1250  //================================================ Get all map values
1251  for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1252  {
1253  mapVals.at(iter) = this->internalMap[iter];
1254  }
1255 
1256  //================================================ Get mean and sd
1257  proshade_double* meanSD = new proshade_double[2];
1258  ProSHADE_internal_maths::vectorMeanAndSD ( &mapVals, meanSD );
1259 
1260  //================================================ Normalise the values
1261  for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1262  {
1263  this->internalMap[iter] = ( this->internalMap[iter] - meanSD[0] ) / meanSD[1];
1264  }
1265 
1266  //================================================ Clear the vector
1267  mapVals.clear ( );
1268 
1269  //================================================ Release memory
1270  delete[] meanSD;
1271 
1272  //================================================ Report function completion
1273  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Map normalisation completed." );
1274 
1275  //================================================ Done
1276  return ;
1277 
1278 }
1279 
1280 
1290 {
1291  //================================================ Report function start
1292  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Computing mask." );
1293 
1294  //================================================ Initialise the blurred map
1295  proshade_double* blurredMap = new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1296  ProSHADE_internal_misc::checkMemoryAllocation ( blurredMap, __FILE__, __LINE__, __func__ );
1297 
1298  //================================================ Compute blurred map
1299  ProSHADE_internal_mapManip::blurSharpenMap ( this->internalMap, blurredMap, this->xDimIndices, this->yDimIndices, this->zDimIndices,
1300  this->xDimSize, this->yDimSize, this->zDimSize, settings->blurFactor );
1301 
1302  //================================================ Compute mask from blurred map and save it into the original map
1303  ProSHADE_internal_mapManip::getMaskFromBlurr ( blurredMap, this->internalMap, this->xDimIndices, this->yDimIndices, this->zDimIndices, settings->maskingThresholdIQRs );
1304 
1305  //================================================ Print the mask if need be
1306  if ( settings->saveMask ) { if ( settings->maskFileName == "" ) { this->writeMask ( "proshade_mask.map", blurredMap ); } else { std::stringstream ss; ss << settings->maskFileName << "_" << this->inputOrder << ".map"; this->writeMask ( ss.str(), blurredMap ); } }
1307 
1308  //================================================ Release memory
1309  delete[] blurredMap;
1310 
1311  //================================================ Report function completion
1312  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Mask computed." );
1313 
1314  //================================================ Done
1315  return ;
1316 
1317 }
1318 
1330 {
1331  //================================================ Report function start
1332  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Finding new boundaries." );
1333 
1334  //================================================ If same bounds as first one are required, test if possible and return these instead
1335  if ( settings->useSameBounds && ( this->inputOrder != 0 ) )
1336  {
1337  for ( proshade_unsign iter = 0; iter < 6; iter++ ) { ret[iter] = settings->forceBounds[iter]; }
1338  }
1339  //================================================ In this case, bounds need to be found de novo
1340  else
1341  {
1342  //============================================ Find the non-zero bounds
1344  static_cast< proshade_signed > ( this->xDimIndices ),
1345  static_cast< proshade_signed > ( this->yDimIndices ),
1346  static_cast< proshade_signed > ( this->zDimIndices ),
1347  ret );
1348 
1349  //============================================ Add the extra space
1350  ProSHADE_internal_mapManip::addExtraBoundSpace ( this->xDimIndices, this->yDimIndices, this->zDimIndices,
1351  this->xDimSize, this->yDimSize, this->zDimSize, ret, settings->boundsExtraSpace );
1352 
1353  //============================================ Beautify boundaries
1354  ProSHADE_internal_mapManip::beautifyBoundaries ( ret, this->xDimIndices, this->yDimIndices, this->zDimIndices, settings->boundsSimilarityThreshold );
1355 
1356  //============================================ Report function results
1357  std::stringstream ssHlp;
1358  ssHlp << "New boundaries are: " << ret[1] - ret[0] + 1 << " x " << ret[3] - ret[2] + 1 << " x " << ret[5] - ret[4] + 1;
1359  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 3, ssHlp.str() );
1360 
1361  //============================================ If need be, save boundaries to be used for all other structure
1362  if ( settings->useSameBounds && ( this->inputOrder == 0 ) )
1363  {
1364  for ( proshade_unsign iter = 0; iter < 6; iter++ ) { settings->forceBounds[iter] = ret[iter]; }
1365  }
1366  }
1367 
1368  //================================================ Report function completion
1369  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "New boundaries determined." );
1370 
1371  //================================================ Done
1372  return ;
1373 
1374 }
1375 
1388 {
1389  //================================================ Report function start
1390  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Creating new structure according to the new bounds." );
1391 
1392  //================================================ Fill in basic info
1393  newStr->fileName = "N/A";
1394  newStr->fileType = ProSHADE_internal_io::MAP;
1395 
1396  //================================================ Fill in new structure values
1397  newStr->xDimIndices = static_cast< proshade_unsign > ( newBounds[1] ) - static_cast< proshade_unsign > ( newBounds[0] ) + 1;
1398  newStr->yDimIndices = static_cast< proshade_unsign > ( newBounds[3] ) - static_cast< proshade_unsign > ( newBounds[2] ) + 1;
1399  newStr->zDimIndices = static_cast< proshade_unsign > ( newBounds[5] ) - static_cast< proshade_unsign > ( newBounds[4] ) + 1;
1400 
1401  newStr->aAngle = this->aAngle;
1402  newStr->bAngle = this->aAngle;
1403  newStr->cAngle = this->aAngle;
1404 
1405  newStr->xDimSize = static_cast<proshade_single> ( newStr->xDimIndices ) * ( this->xDimSize / static_cast<proshade_single> ( this->xDimIndices ) );
1406  newStr->yDimSize = static_cast<proshade_single> ( newStr->yDimIndices ) * ( this->yDimSize / static_cast<proshade_single> ( this->yDimIndices ) );
1407  newStr->zDimSize = static_cast<proshade_single> ( newStr->zDimIndices ) * ( this->zDimSize / static_cast<proshade_single> ( this->zDimIndices ) );
1408 
1409  newStr->xGridIndices = newStr->xDimIndices;
1410  newStr->yGridIndices = newStr->yDimIndices;
1411  newStr->zGridIndices = newStr->zDimIndices;
1412 
1413  newStr->xAxisOrder = this->xAxisOrder;
1414  newStr->yAxisOrder = this->yAxisOrder;
1415  newStr->zAxisOrder = this->zAxisOrder;
1416 
1417  newStr->xAxisOrigin = this->xAxisOrigin + newBounds[0];
1418  newStr->yAxisOrigin = this->yAxisOrigin + newBounds[2];
1419  newStr->zAxisOrigin = this->zAxisOrigin + newBounds[4];
1420 
1421  newStr->xFrom = this->xFrom + newBounds[0];
1422  newStr->yFrom = this->yFrom + newBounds[2];
1423  newStr->zFrom = this->zFrom + newBounds[4];
1424 
1425  newStr->xTo = this->xTo - ( static_cast< proshade_signed > ( this->xDimIndices - 1 ) - newBounds[1] );
1426  newStr->yTo = this->yTo - ( static_cast< proshade_signed > ( this->yDimIndices - 1 ) - newBounds[3] );
1427  newStr->zTo = this->zTo - ( static_cast< proshade_signed > ( this->zDimIndices - 1 ) - newBounds[5] );
1428 
1429  //================================================ Allocate new structure map
1430  newStr->internalMap = new proshade_double[newStr->xDimIndices * newStr->yDimIndices * newStr->zDimIndices];
1431  ProSHADE_internal_misc::checkMemoryAllocation ( newStr->internalMap, __FILE__, __LINE__, __func__ );
1432 
1433  //================================================ Copy the map
1434  ProSHADE_internal_mapManip::copyMapByBounds ( newStr->xFrom, newStr->xTo, newStr->yFrom, newStr->yTo, newStr->zFrom, newStr->zTo,
1435  this->xFrom, this->yFrom, this->zFrom, newStr->yDimIndices, newStr->zDimIndices,
1436  this->xDimIndices, this->yDimIndices, this->zDimIndices, newStr->internalMap, this->internalMap );
1437 
1438  //================================================ Report function completion
1439  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "New structure created." );
1440 
1441  //================================================ Done
1442  return ;
1443 
1444 }
1445 
1454 {
1455  //================================================ Initialise the return variable
1456  proshade_single* changeVals = new proshade_single[6];
1457 
1458  //================================================ Now re-sample the map
1459  if ( settings->changeMapResolution )
1460  {
1461  ProSHADE_internal_mapManip::reSampleMapToResolutionFourier ( this->internalMap, settings->requestedResolution, this->xDimIndices, this->yDimIndices, this->zDimIndices,
1462  this->xDimSize, this->yDimSize, this->zDimSize, changeVals );
1463 
1464  if ( settings->changeMapResolutionTriLinear )
1465  {
1466  ProSHADE_internal_messages::printWarningMessage ( settings->verbose, "!!! ProSHADE WARNING !!! Requested both Fourier-space and real-space map re-sampling. Defaulting to only Fourier space re-samplling.", "WM00049" );
1467  }
1468  }
1469  if ( settings->changeMapResolutionTriLinear && !settings->changeMapResolution )
1470  {
1471  ProSHADE_internal_mapManip::reSampleMapToResolutionTrilinear ( this->internalMap, settings->requestedResolution, this->xDimIndices, this->yDimIndices, this->zDimIndices,
1472  this->xDimSize, this->yDimSize, this->zDimSize, changeVals );
1473 
1474  }
1475 
1476  //================================================ Set the internal values to reflect the new map size
1477  this->xDimIndices += static_cast<proshade_unsign> ( changeVals[0] );
1478  this->yDimIndices += static_cast<proshade_unsign> ( changeVals[1] );
1479  this->zDimIndices += static_cast<proshade_unsign> ( changeVals[2] );
1480 
1481  this->xGridIndices = this->xDimIndices;
1482  this->yGridIndices = this->yDimIndices;
1483  this->zGridIndices = this->zDimIndices;
1484 
1485  this->xTo += static_cast<proshade_unsign> ( changeVals[0] );
1486  this->yTo += static_cast<proshade_unsign> ( changeVals[1] );
1487  this->zTo += static_cast<proshade_unsign> ( changeVals[2] );
1488 
1489  this->xDimSize = changeVals[3];
1490  this->yDimSize = changeVals[4];
1491  this->zDimSize = changeVals[5];
1492 
1493  //================================================ Figure how much the new map moved
1494  proshade_single xMov = -( ( static_cast<proshade_single> ( this->xFrom ) * ( this->xDimSize / static_cast<proshade_single> ( this->xDimIndices ) - changeVals[0] ) ) -
1495  ( static_cast<proshade_single> ( this->xFrom ) * ( this->xDimSize / static_cast<proshade_single> ( this->xDimIndices ) ) ) );
1496  proshade_single yMov = -( ( static_cast<proshade_single> ( this->yFrom ) * ( this->yDimSize / static_cast<proshade_single> ( this->yDimIndices ) - changeVals[1] ) ) -
1497  ( static_cast<proshade_single> ( this->yFrom ) * ( this->yDimSize / static_cast<proshade_single> ( this->yDimIndices ) ) ) );
1498  proshade_single zMov = -( ( static_cast<proshade_single> ( this->zFrom ) * ( this->zDimSize / static_cast<proshade_single> ( this->zDimIndices ) - changeVals[2] ) ) -
1499  ( static_cast<proshade_single> ( this->zFrom ) * ( this->zDimSize / static_cast<proshade_single> ( this->zDimIndices ) ) ) );
1500 
1501  //================================================ Move by indices (this should be sufficient)
1502  ProSHADE_internal_mapManip::moveMapByIndices ( &xMov, &yMov, &zMov, this->xDimSize, this->yDimSize, this->zDimSize, &this->xFrom, &this->xTo,
1503  &this->yFrom, &this->yTo, &this->zFrom, &this->zTo, &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin );
1504 
1505  ProSHADE_internal_mapManip::moveMapByFourier ( this->internalMap, xMov, yMov, zMov, this->xDimSize, this->yDimSize, this->zDimSize,
1506  static_cast< proshade_signed > ( this->xDimIndices ), static_cast< proshade_signed > ( this->yDimIndices ), static_cast< proshade_signed > ( this->zDimIndices ) );
1507 
1508  //================================================ Release memory
1509  delete[] changeVals;
1510 
1511  //================================================ Done
1512  return ;
1513 
1514 }
1515 
1526 {
1527  //================================================ Report function start
1528  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Centering map onto its COM." );
1529 
1530  //================================================ Initialise local variables
1531  proshade_double xCOM = 0.0;
1532  proshade_double yCOM = 0.0;
1533  proshade_double zCOM = 0.0;
1534 
1535  //================================================ Find the COM location
1536  ProSHADE_internal_mapManip::findMAPCOMValues ( this->internalMap, &xCOM, &yCOM, &zCOM,
1537  this->xDimSize, this->yDimSize, this->xDimSize, this->xFrom,
1538  this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
1539 
1540  //================================================ Find the sampling rates
1541  proshade_single xSampRate = static_cast< proshade_single > ( this->xDimSize ) / static_cast< proshade_single > ( this->xTo - this->xFrom );
1542  proshade_single ySampRate = static_cast< proshade_single > ( this->yDimSize ) / static_cast< proshade_single > ( this->yTo - this->yFrom );
1543  proshade_single zSampRate = static_cast< proshade_single > ( this->zDimSize ) / static_cast< proshade_single > ( this->zTo - this->zFrom );
1544 
1545  //================================================ Convert to position in indices starting from 0
1546  xCOM /= static_cast< proshade_double > ( xSampRate );
1547  yCOM /= static_cast< proshade_double > ( ySampRate );
1548  zCOM /= static_cast< proshade_double > ( zSampRate );
1549 
1550  xCOM -= static_cast< proshade_double > ( this->xFrom );
1551  yCOM -= static_cast< proshade_double > ( this->yFrom );
1552  zCOM -= static_cast< proshade_double > ( this->zFrom );
1553 
1554  //================================================ Find distance from COM to map centre in Angstroms
1555  proshade_double xDist = ( ( static_cast<proshade_double> ( this->xDimIndices ) / 2.0 ) - xCOM ) * static_cast<proshade_double> ( this->xDimSize ) / static_cast<proshade_double> ( this->xDimIndices );
1556  proshade_double yDist = ( ( static_cast<proshade_double> ( this->yDimIndices ) / 2.0 ) - yCOM ) * static_cast<proshade_double> ( this->yDimSize ) / static_cast<proshade_double> ( this->yDimIndices );
1557  proshade_double zDist = ( ( static_cast<proshade_double> ( this->zDimIndices ) / 2.0 ) - zCOM ) * static_cast<proshade_double> ( this->zDimSize ) / static_cast<proshade_double> ( this->zDimIndices );
1558 
1559  //================================================ Move the map within the box
1561  static_cast< proshade_single > ( xDist ),
1562  static_cast< proshade_single > ( yDist ),
1563  static_cast< proshade_single > ( zDist ),
1564  this->xDimSize, this->yDimSize, this->zDimSize,
1565  static_cast< proshade_signed > ( this->xDimIndices ),
1566  static_cast< proshade_signed > ( this->yDimIndices ),
1567  static_cast< proshade_signed > ( this->zDimIndices ) );
1568 
1569  //================================================ Note the change due to centering
1570  this->mapCOMProcessChangeX -= xDist;
1571  this->mapCOMProcessChangeY -= yDist;
1572  this->mapCOMProcessChangeZ -= zDist;
1573 
1574  //================================================ Report function completion
1575  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Map centered." );
1576 
1577  //================================================ Done
1578  return ;
1579 
1580 }
1581 
1591 {
1592  //================================================ Report function start
1593  std::stringstream hlpSS;
1594  hlpSS << "Adding extra " << settings->addExtraSpace << " angstroms.";
1595  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, hlpSS.str() );
1596 
1597  //================================================ Figure how much indices need to change
1598  proshade_unsign xAddIndices = static_cast< proshade_unsign > ( ProSHADE_internal_mapManip::myRound ( settings->addExtraSpace / this->xDimSize / static_cast<proshade_single> ( this->xDimIndices ) ) );
1599  proshade_unsign yAddIndices = static_cast< proshade_unsign > ( ProSHADE_internal_mapManip::myRound ( settings->addExtraSpace / this->yDimSize / static_cast<proshade_single> ( this->yDimIndices ) ) );
1600  proshade_unsign zAddIndices = static_cast< proshade_unsign > ( ProSHADE_internal_mapManip::myRound ( settings->addExtraSpace / this->zDimSize / static_cast<proshade_single> ( this->zDimIndices ) ) );
1601 
1602  //================================================ Update internal data variables
1603  this->xDimSize += 2 * static_cast<proshade_single> ( xAddIndices ) * this->xDimSize / static_cast<proshade_single> ( this->xDimIndices );
1604  this->yDimSize += 2 * static_cast<proshade_single> ( yAddIndices ) * this->yDimSize / static_cast<proshade_single> ( this->yDimIndices );
1605  this->zDimSize += 2 * static_cast<proshade_single> ( zAddIndices ) * this->zDimSize / static_cast<proshade_single> ( this->zDimIndices );
1606 
1607  this->xDimIndices += 2 * xAddIndices;
1608  this->yDimIndices += 2 * yAddIndices;
1609  this->zDimIndices += 2 * zAddIndices;
1610 
1611  this->xGridIndices = this->xDimIndices;
1612  this->yGridIndices = this->yDimIndices;
1613  this->zGridIndices = this->zDimIndices;
1614 
1615  this->xAxisOrigin -= xAddIndices;
1616  this->yAxisOrigin -= yAddIndices;
1617  this->zAxisOrigin -= zAddIndices;
1618 
1619  this->xFrom -= xAddIndices;
1620  this->yFrom -= yAddIndices;
1621  this->zFrom -= zAddIndices;
1622 
1623  this->xTo += xAddIndices;
1624  this->yTo += yAddIndices;
1625  this->zTo += zAddIndices;
1626 
1627  //================================================ Allocate new map
1628  proshade_double* newMap = new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1629  ProSHADE_internal_misc::checkMemoryAllocation ( newMap, __FILE__, __LINE__, __func__ );
1630 
1631  //================================================ Set new map to zeroes
1632  for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1633  {
1634  newMap[iter] = 0.0;
1635  }
1636 
1637  //================================================ Update the map
1638  proshade_unsign newMapIndex, oldMapIndex;
1639  for ( proshade_unsign xIt = 0; xIt < (this->xDimIndices - xAddIndices); xIt++ )
1640  {
1641  //============================================ Check if point is applicable
1642  if ( xIt < xAddIndices ) { continue; }
1643 
1644  for ( proshade_unsign yIt = 0; yIt < (this->yDimIndices - yAddIndices); yIt++ )
1645  {
1646  //======================================== Check if point is applicable
1647  if ( yIt < yAddIndices ) { continue; }
1648 
1649  for ( proshade_unsign zIt = 0; zIt < (this->zDimIndices - zAddIndices); zIt++ )
1650  {
1651  //==================================== Check if point is applicable
1652  if ( zIt < zAddIndices ) { continue; }
1653 
1654  //==================================== Var init
1655  newMapIndex = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
1656  oldMapIndex = (zIt - zAddIndices) + (this->zDimIndices - ( 2 * zAddIndices ) ) * ( (yIt - yAddIndices) + (this->yDimIndices - ( 2 * yAddIndices ) ) * (xIt - xAddIndices) );
1657 
1658  newMap[newMapIndex] = this->internalMap[oldMapIndex];
1659  }
1660  }
1661  }
1662 
1663  //================================================ Copy new to old
1664  delete[] this->internalMap;
1665 
1666  this->internalMap = new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1667  ProSHADE_internal_misc::checkMemoryAllocation ( this->internalMap, __FILE__, __LINE__, __func__ );
1668 
1669  for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1670  {
1671  this->internalMap[iter] = newMap[iter];
1672  }
1673 
1674  //================================================ Release memory
1675  delete[] newMap;
1676 
1677  //================================================ Report function completion
1678  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Extra space added." );
1679 
1680  //================================================ Done
1681  return ;
1682 
1683 }
1684 
1700 {
1701  //================================================ Invert map
1702  if ( settings->invertMap ) { this->invertMirrorMap ( settings ); }
1703  else { ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Map inversion (mirror image) not requested." ); }
1704 
1705  //================================================ Normalise map
1706  if ( settings->normaliseMap ) { this->normaliseMap ( settings ); }
1707  else { ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Map normalisation not requested." ); }
1708 
1709  //================================================ Compute mask
1710  if ( settings->maskMap ) { this->maskMap ( settings ); }
1711  else { ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Masking not requested." ); }
1712 
1713  //================================================ Centre map
1714  if ( settings->moveToCOM ) { this->centreMapOnCOM ( settings ); }
1715  else { ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Map centering not requested." ); }
1716 
1717  //================================================ Add extra space
1718  if ( settings->addExtraSpace != 0.0f ) { this->addExtraSpace ( settings ); }
1719  else { ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Extra space not requested." ); }
1720 
1721  //================================================ Remove phase, if required
1722  if ( !settings->usePhase ) { this->removePhaseInormation ( settings ); ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Phase information removed from the data." ); }
1723  else { ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Phase information retained in the data." ); }
1724 
1725  //================================================ Set settings values which were left on AUTO by user and will not be set later
1726  settings->setVariablesLeftOnAuto ( );
1727 
1728  //================================================ Done
1729  return ;
1730 
1731 }
1732 
1743 {
1744  //================================================ Check the current settings value is set to auto
1745  if ( this->spherePos.size() != 0 )
1746  {
1747  std::stringstream hlpSS;
1748  hlpSS << "The sphere distances were determined as " << this->spherePos.at(0);
1749  for ( proshade_unsign iter = 1; iter < static_cast<proshade_unsign> ( this->spherePos.size() ); iter++ ) { hlpSS << "; " << this->spherePos.at(iter); }
1750  hlpSS << " Angstroms.";
1751  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 3, hlpSS.str() );
1752  return ;
1753  }
1754 
1755  //================================================ Find maximum diagonal
1756  proshade_unsign maxDim = static_cast< proshade_unsign > ( std::max ( this->xDimSize, std::max ( this->yDimSize, this->zDimSize ) ) );
1757  proshade_unsign minDim = static_cast< proshade_unsign > ( std::min ( this->xDimSize, std::min ( this->yDimSize, this->zDimSize ) ) );
1758  proshade_unsign midDim = static_cast< proshade_unsign > ( 0 );
1759  if ( ( this->xDimSize < static_cast< proshade_single > ( maxDim ) ) && ( this->xDimSize > static_cast< proshade_single > ( minDim ) ) ) { midDim = static_cast< proshade_unsign > ( this->xDimSize ); }
1760  else if ( ( this->yDimSize < static_cast< proshade_single > ( maxDim ) ) && ( this->yDimSize > static_cast< proshade_single > ( minDim ) ) ) { midDim = static_cast< proshade_unsign > ( this->yDimSize ); }
1761  else { midDim = static_cast< proshade_unsign > ( this->zDimSize ); }
1762 
1763  proshade_single maxDiag = static_cast< proshade_single > ( std::sqrt ( std::pow ( static_cast<proshade_single> ( maxDim ), 2.0 ) +
1764  std::pow ( static_cast<proshade_single> ( midDim ), 2.0 ) ) );
1765 
1766  //================================================ Set between the points
1767  for ( proshade_single iter = 0.5f; ( iter * settings->maxSphereDists ) < ( maxDiag / 2.0f ); iter += 1.0f )
1768  {
1769  ProSHADE_internal_misc::addToSingleVector ( &this->spherePos, ( iter * settings->maxSphereDists ) );
1770  }
1771 
1772  //================================================ Save the number of spheres
1773  this->noSpheres = static_cast<proshade_unsign> ( this->spherePos.size() );
1774 
1775  //================================================ Report progress
1776  std::stringstream hlpSS;
1777  hlpSS << "The sphere distances were determined as " << this->spherePos.at(0);
1778  for ( proshade_unsign iter = 1; iter < static_cast<proshade_unsign> ( this->spherePos.size() ); iter++ ) { hlpSS << "; " << this->spherePos.at(iter); }
1779  hlpSS << " Angstroms.";
1780  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 3, hlpSS.str() );
1781 
1782  //================================================ Done
1783  return ;
1784 
1785 }
1786 
1787 
1800 {
1801  //================================================ Report progress
1802  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting sphere mapping procedure." );
1803 
1804  //================================================ Determine spherical harmonics variables
1805  settings->determineAllSHValues ( this->xDimIndices, this->yDimIndices,
1806  this->xDimSize, this->yDimSize, this->zDimSize );
1807  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Sphere settings determined." );
1808 
1809  //================================================ Find number of spheres supported
1810  this->getSpherePositions ( settings );
1811  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Sphere positions obtained." );
1812 
1813  //================================================ Create sphere objects and map the density
1814  this->spheres = new ProSHADE_internal_spheres::ProSHADE_sphere* [ this->noSpheres ];
1815  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( this->spherePos.size() ); iter++ )
1816  {
1817  std::stringstream ss;
1818  ss << "Now mapping sphere " << iter << " .";
1819  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 4, ss.str() );
1820 
1821  this->spheres[iter] = new ProSHADE_internal_spheres::ProSHADE_sphere ( this->xDimIndices, this->yDimIndices, this->zDimIndices,
1822  this->xDimSize, this->yDimSize, this->zDimSize, iter,
1823  &this->spherePos, settings->progressiveSphereMapping, settings->maxBandwidth,
1824  this->internalMap, &this->maxShellBand );
1825  }
1826 
1827  //================================================ Report completion
1828  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Sphere mapping procedure completed." );
1829 
1830  //================================================ Done
1831  return ;
1832 
1833 }
1834 
1844 {
1845  //================================================ Report progress
1846  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting spherical harmonics decomposition." );
1847 
1848  //================================================ Initialise memory
1849  this->sphericalHarmonics = new proshade_complex* [this->noSpheres];
1850  ProSHADE_internal_misc::checkMemoryAllocation ( this->sphericalHarmonics, __FILE__, __LINE__, __func__ );
1851  for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
1852  {
1853  this->sphericalHarmonics[iter] = new proshade_complex [(this->spheres[iter]->getLocalBandwidth() * 2) * (this->spheres[iter]->getLocalBandwidth() * 2)];
1854  ProSHADE_internal_misc::checkMemoryAllocation ( this->sphericalHarmonics[iter], __FILE__, __LINE__, __func__ );
1855  }
1856 
1857  //================================================ Compute the spherical harmonics
1858  for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
1859  {
1860  //============================================ Report progress
1861  std::stringstream ss;
1862  ss << "Now decomposing sphere " << iter << ". " << "( Band is: " << this->spheres[iter]->getLocalBandwidth() << ").";
1863  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 4, ss.str() );
1864 
1865  //============================================ Compute
1866  ProSHADE_internal_sphericalHarmonics::computeSphericalHarmonics ( this->spheres[iter]->getLocalBandwidth(), this->spheres[iter]->getMappedData(), this->sphericalHarmonics[iter] );
1867  }
1868 
1869  //================================================ Report completion
1870  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Spherical harmonics decomposition complete." );
1871 
1872  //======================================== Done
1873  return ;
1874 
1875 }
1876 
1883 bool sortProSHADESymmetryByFSC ( proshade_double* a, proshade_double* b)
1884 {
1885  //================================================ Done
1886  return ( a[6] > b[6] );
1887 
1888 }
1889 
1903 void ProSHADE_internal_data::ProSHADE_data::detectSymmetryFromAngleAxisSpace ( ProSHADE_settings* settings, std::vector< proshade_double* >* axes, std::vector < std::vector< proshade_double > >* allCs )
1904 {
1905  //================================================ Modify axis tolerance and matrix tolerance by sampling, if required by user
1906  if ( settings->axisErrToleranceDefault )
1907  {
1908  settings->axisErrTolerance = std::min ( std::max ( 0.01, ( 2.0 * M_PI ) / static_cast< proshade_double > ( this->maxShellBand ) ), 0.05 );
1909  }
1910 
1911  //================================================ Prepare FSC computation memory and variables
1912  fftw_complex *mapData, *origCoeffs, *fCoeffs;
1913  fftw_plan planForwardFourier;
1914  proshade_double **bindata;
1915  proshade_signed *binIndexing, *binCounts, noBins;
1916  this->prepareFSCFourierMemory ( mapData, origCoeffs, fCoeffs, binIndexing, &noBins, bindata, binCounts, &planForwardFourier );
1917 
1918  //================================================ If C was requested, we will do it immediately - this allows for a significant speed-up.
1919  if ( settings->requestedSymmetryType == "C" )
1920  {
1921  //============================================ Report progress
1922  std::stringstream hlpSS;
1923  hlpSS << "Starting detection of cyclic point group C" << settings->requestedSymmetryFold;
1924  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, hlpSS.str() );
1925 
1926  //============================================ Do simplified search only in the applicable data
1927  proshade_double symThres = 0.0;
1928  std::vector< proshade_double* > CSyms = this->findRequestedCSymmetryFromAngleAxis ( settings, settings->requestedSymmetryFold, &symThres );
1929 
1930  //============================================ Deal with the rotation function 0 1 0 PI issue
1931  bool possible010PIIssue = false;
1932  for ( size_t iter = 0; iter < CSyms.size(); iter++ ) { if ( ProSHADE_internal_maths::isAxisUnique ( &CSyms, 0.0, 1.0, 0.0, 2.0, 0.1 ) ) { possible010PIIssue = true; } }
1933  if ( possible010PIIssue )
1934  {
1935  proshade_double* addAxis = new proshade_double[7];
1936  addAxis[0] = 2.0; addAxis[1] = 0.0; addAxis[2] = 1.0; addAxis[3] = 0.0; addAxis[4] = M_PI; addAxis[5] = -999.9; addAxis[6] = -1.0;
1938  delete[] addAxis;
1939  }
1940 
1941  //============================================ Compute FSC for all possible axes
1942  for ( size_t cIt = 0; cIt < CSyms.size(); cIt++ ) { const FloatingPoint< proshade_double > lhs ( CSyms.at(cIt)[5] ), rhs ( -999.9 ); if ( ( CSyms.at(cIt)[5] > settings->peakThresholdMin ) || ( lhs.AlmostEquals ( rhs ) ) ) { this->computeFSC ( settings, &CSyms, cIt, mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts ); } }
1943 
1944  //============================================ Sort by FSC
1945  std::sort ( CSyms.begin(), CSyms.end(), sortProSHADESymmetryByFSC );
1946 
1947  //============================================ Save the best axis as the recommended one
1948  if ( settings->detectedSymmetry.size() == 0 ) { if ( CSyms.size() > 0 ) { settings->setDetectedSymmetry ( CSyms.at(0) ); } }
1949  if ( CSyms.size() > 0 )
1950  {
1951  bool passedTests = false;
1952  for ( size_t cIt = 0; cIt < CSyms.size(); cIt++ )
1953  {
1954  if ( CSyms.at(0)[6] > settings->fscThreshold )
1955  {
1956  settings->setRecommendedSymmetry ( "C" );
1957  settings->setRecommendedFold ( settings->requestedSymmetryFold );
1958 
1960  this->saveDetectedSymmetries ( settings, &CSyms, allCs );
1961 
1962  passedTests = true;
1963  break;
1964  }
1965  }
1966 
1967  if ( !passedTests )
1968  {
1969  settings->setRecommendedSymmetry ( "" );
1970  settings->setRecommendedFold ( 0 );
1971  }
1972  }
1973  else
1974  {
1975  settings->setRecommendedSymmetry ( "" );
1976  settings->setRecommendedFold ( 0 );
1977  }
1978 
1979  //============================================ Release memory after FSC computation
1980  delete[] mapData;
1981  delete[] origCoeffs;
1982  delete[] fCoeffs;
1983  fftw_destroy_plan ( planForwardFourier );
1984  delete[] binIndexing;
1985  for (size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ ) { delete[] bindata[binIt]; }
1986  delete[] bindata;
1987  delete[] binCounts;
1988 
1989  //============================================ Done
1990  return ;
1991  }
1992 
1993  //================================================ Report progress
1994  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting C symmetry detection." );
1995 
1996  //================================================ Initialise variables
1997  std::vector< proshade_double* > CSyms = getCyclicSymmetriesListFromAngleAxis ( settings );
1998 
1999  //============================================ Deal with the rotation function 0 1 0 PI issue
2000  bool possible010PIIssue = false;
2001  for ( size_t iter = 0; iter < CSyms.size(); iter++ ) { if ( ProSHADE_internal_maths::isAxisUnique ( &CSyms, 0.0, 1.0, 0.0, 2.0, 0.1 ) ) { possible010PIIssue = true; } }
2002  if ( possible010PIIssue )
2003  {
2004  proshade_double* addAxis = new proshade_double[7];
2005  addAxis[0] = 2.0; addAxis[1] = 0.0; addAxis[2] = 1.0; addAxis[3] = 0.0; addAxis[4] = M_PI; addAxis[5] = -999.9; addAxis[6] = -1.0;
2007  delete[] addAxis;
2008  }
2009 
2010  for ( size_t cIt = 0; cIt < CSyms.size(); cIt++ ) { const FloatingPoint< proshade_double > lhs ( CSyms.at(cIt)[5] ), rhs ( -999.9 ); if ( lhs.AlmostEquals ( rhs ) ) { this->computeFSC ( settings, &CSyms, cIt, mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts ); } }
2011 
2012  //================================================ Report progress
2013  std::stringstream ss;
2014  ss << "Detected " << CSyms.size() << " C symmetries.";
2015  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, ss.str() );
2016 
2017  //================================================ Sanity check - was the rotation function mapped properly?
2018  if ( this->sphereMappedRotFun.size() < 1 )
2019  {
2020  throw ProSHADE_exception ( "Rotation function was not converted into angle-axis space.", "ES00062", __FILE__, __LINE__, __func__, "It seems that the convertRotationFunction() function was\n : not yet called. Therefore, there are no data to detect the\n : symmetry from; please call the convertRotationFunction()\n : function before the detectSymmetryFromAngleAxisSpace()\n : function." );
2021  }
2022 
2023  //================================================ Sanity check - was any symmetry requested?
2024  if ( ( settings->requestedSymmetryType != "" ) && ( settings->requestedSymmetryType != "C" ) && ( settings->requestedSymmetryType != "D" ) && ( settings->requestedSymmetryType != "T" ) && ( settings->requestedSymmetryType != "O" ) && ( settings->requestedSymmetryType != "I" ) )
2025  {
2026  throw ProSHADE_exception ( "Requested symmetry supplied, but not recognised.", "ES00032", __FILE__, __LINE__, __func__, "There are only the following value allowed for the\n : symmetry type request: \"C\", \"D\", \"T\", \"O\" and \"I\". Any\n : other value will result in this error." );
2027  }
2028 
2029  //================================================ Are we doing general search?
2030  if ( settings->requestedSymmetryType == "" )
2031  {
2032  //============================================ Run the symmetry detection functions for C, D, T, O and I symmetries
2033  std::vector< proshade_double* > DSyms = this->getDihedralSymmetriesList ( settings, &CSyms );
2034  std::vector< proshade_double* > ISyms;
2035  std::vector< proshade_double* > OSyms = this->getPredictedOctahedralSymmetriesList ( settings, &CSyms );
2036  std::vector< proshade_double* > TSyms = this->getPredictedTetrahedralSymmetriesList ( settings, &CSyms );
2037 
2038  //============================================ Find which of the I groups is the correct one
2039  proshade_double fscMax = 0.0;
2040  size_t fscMaxInd = 0;
2041  std::vector < std::vector< proshade_double* > > ISymsHlp = this->getPredictedIcosahedralSymmetriesList ( settings, &CSyms );
2042 
2043  for ( size_t icoIt = 0; icoIt < ISymsHlp.size(); icoIt++ )
2044  {
2045  //======================================== Is this a complete icosahedron?
2046  if ( ISymsHlp.at(icoIt).size() != 31 ) { continue; }
2047 
2048  //======================================== Initialise decision vars
2049  proshade_double fscVal = 0.0;
2050  proshade_double fscValAvg = 0.0;
2051 
2052  //======================================== For each axis
2053  for ( size_t aIt = 0; aIt < ISymsHlp.at(icoIt).size(); aIt++ )
2054  {
2055  //==================================== Match to CSyms
2056  proshade_signed matchedPos = ProSHADE_internal_symmetry::addAxisUnlessSame ( static_cast< proshade_unsign > ( ISymsHlp.at(icoIt).at(aIt)[0] ), ISymsHlp.at(icoIt).at(aIt)[1], ISymsHlp.at(icoIt).at(aIt)[2], ISymsHlp.at(icoIt).at(aIt)[3], ISymsHlp.at(icoIt).at(aIt)[5], ISymsHlp.at(icoIt).at(aIt)[6], &CSyms, settings->axisErrTolerance );
2057 
2058  //==================================== Compute FSC
2059  fscVal = this->computeFSC ( settings, &CSyms, static_cast< size_t > ( matchedPos ), mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts );
2060  ISymsHlp.at(fscMaxInd).at(aIt)[6] = fscVal;
2061  fscValAvg += fscVal;
2062  }
2063 
2064  //==================================== Get FSC average over all axes
2065  fscValAvg /= 31.0;
2066 
2067  //======================================== Is this the best
2068  if ( fscValAvg > fscMax )
2069  {
2070  fscMax = fscValAvg;
2071  fscMaxInd = icoIt;
2072  }
2073  }
2074 
2075  //============================================ If C3 and C5 are found and have correct angle (must have if they are both in ISym)
2076  if ( fscMax >= ( settings->fscThreshold * 0.7 ) )
2077  {
2078  //======================================== Add predicted axes to detected C axes list and also to the settings Icosahedral symmetry list
2079  for ( proshade_unsign retIt = 0; retIt < static_cast < proshade_unsign > ( ISymsHlp.at(fscMaxInd).size() ); retIt++ )
2080  {
2081  //==================================== Add the correct index to the settings object
2082  proshade_signed matchedPos = ProSHADE_internal_symmetry::addAxisUnlessSame ( static_cast< proshade_unsign > ( ISymsHlp.at(fscMaxInd).at(retIt)[0] ), ISymsHlp.at(fscMaxInd).at(retIt)[1], ISymsHlp.at(fscMaxInd).at(retIt)[2], ISymsHlp.at(fscMaxInd).at(retIt)[3], ISymsHlp.at(fscMaxInd).at(retIt)[5], ISymsHlp.at(fscMaxInd).at(retIt)[6], &CSyms, settings->axisErrTolerance );
2083  ProSHADE_internal_misc::addToUnsignVector ( &settings->allDetectedIAxes, static_cast < proshade_unsign > ( matchedPos ) );
2084 
2085  //==================================== Set ISyms for saving
2086  ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( &ISyms, ISymsHlp.at(fscMaxInd).at(retIt) );
2087  }
2088  }
2089 
2090  //============================================ Decide on recommended symmetry
2091  this->saveRecommendedSymmetry ( settings, &CSyms, &DSyms, &TSyms, &OSyms, &ISyms, axes, mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts );
2092  }
2093 
2094  if ( settings->requestedSymmetryType == "D" )
2095  {
2096  //============================================ Run only the D symmetry detection and search for requested fold
2097  std::vector< proshade_double* > DSyms = this->getDihedralSymmetriesList ( settings, &CSyms );
2098  this->saveRequestedSymmetryD ( settings, &DSyms, axes, mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts );
2099  }
2100 
2101  if ( settings->requestedSymmetryType == "T" )
2102  {
2103  //============================================ Run only the T symmetry detection
2104  std::vector< proshade_double* > TSyms = this->getPredictedTetrahedralSymmetriesList ( settings, &CSyms );
2105 
2106  if ( TSyms.size() == 7 )
2107  {
2108  //======================================== Initialise decision vars
2109  proshade_double fscVal = 0.0;
2110  proshade_double fscValAvg = 0.0;
2111 
2112  //======================================== Check if axes have high enough FSC and peak height
2113  for ( size_t tIt = 0; tIt < 7; tIt++ ) { if ( CSyms.at(settings->allDetectedTAxes.at(tIt))[5] > settings->peakThresholdMin ) { fscVal = this->computeFSC ( settings, &CSyms, settings->allDetectedTAxes.at(tIt), mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts ); fscValAvg += fscVal; } }
2114  fscValAvg /= 7.0;
2115 
2116  //======================================== If C3 and C5 are found and have correct angle (must have if they are both in ISym)
2117  if ( fscValAvg >= ( settings->fscThreshold * 0.8 ) )
2118  {
2119  //==================================== The decision is T
2120  settings->setRecommendedSymmetry ( "T" );
2121  settings->setRecommendedFold ( 0 );
2122  for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( settings->allDetectedTAxes.size() ); it++ ) { ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( axes, CSyms.at(settings->allDetectedTAxes.at(it)) ); }
2123  if ( settings->detectedSymmetry.size() == 0 ) { for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( settings->allDetectedTAxes.size() ); it++ ) { settings->setDetectedSymmetry ( CSyms.at(settings->allDetectedTAxes.at(it)) ); } }
2124  }
2125  }
2126  }
2127 
2128  if ( settings->requestedSymmetryType == "O" )
2129  {
2130  //============================================ Run only the O symmetry detection
2131  std::vector< proshade_double* > OSyms = this->getPredictedOctahedralSymmetriesList ( settings, &CSyms );
2132 
2133  if ( OSyms.size() == 13 )
2134  {
2135  //======================================== Initialise decision vars
2136  proshade_double fscVal = 0.0;
2137  proshade_double fscValAvg = 0.0;
2138 
2139  //======================================== Check if at least one C5 and one C3 with the correct angle have high FSC and peak height
2140  for ( size_t oIt = 0; oIt < 13; oIt++ ) { if ( CSyms.at(settings->allDetectedOAxes.at(oIt))[5] > settings->peakThresholdMin ) { fscVal = this->computeFSC ( settings, &CSyms, settings->allDetectedOAxes.at(oIt), mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts ); fscValAvg += fscVal; } }
2141  fscValAvg /= 13.0;
2142 
2143  //======================================== If C3 and C5 are found and have correct angle (must have if they are both in ISym)
2144  if ( fscValAvg >= ( settings->fscThreshold * 0.8 ) )
2145  {
2146  //==================================== The decision is O
2147  settings->setRecommendedSymmetry ( "O" );
2148  settings->setRecommendedFold ( 0 );
2149  for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( settings->allDetectedOAxes.size() ); it++ ) { ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( axes, CSyms.at(settings->allDetectedOAxes.at(it)) ); }
2150  if ( settings->detectedSymmetry.size() == 0 ) { for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( settings->allDetectedOAxes.size() ); it++ ) { settings->setDetectedSymmetry ( CSyms.at(settings->allDetectedOAxes.at(it)) ); } }
2151  }
2152  }
2153  }
2154 
2155  if ( settings->requestedSymmetryType == "I" )
2156  {
2157  //============================================ Find which of the I groups is the correct one
2158  proshade_double fscMax = 0.0;
2159  size_t fscMaxInd = 0;
2160  std::vector < std::vector< proshade_double* > > ISymsHlp = this->getPredictedIcosahedralSymmetriesList ( settings, &CSyms );
2161  std::vector< proshade_double* > ISyms;
2162 
2163  for ( size_t icoIt = 0; icoIt < ISymsHlp.size(); icoIt++ )
2164  {
2165  //======================================== Is this a complete icosahedron?
2166  if ( ISymsHlp.at(icoIt).size() != 31 ) { continue; }
2167 
2168  //======================================== Initialise decision vars
2169  proshade_double fscVal = 0.0;
2170  proshade_double fscValAvg = 0.0;
2171 
2172  //======================================== For each axis
2173  for ( size_t aIt = 0; aIt < ISymsHlp.at(icoIt).size(); aIt++ )
2174  {
2175  //==================================== Match to CSyms
2176  proshade_signed matchedPos = ProSHADE_internal_symmetry::addAxisUnlessSame ( static_cast< proshade_unsign > ( ISymsHlp.at(icoIt).at(aIt)[0] ), ISymsHlp.at(icoIt).at(aIt)[1], ISymsHlp.at(icoIt).at(aIt)[2], ISymsHlp.at(icoIt).at(aIt)[3], ISymsHlp.at(icoIt).at(aIt)[5], ISymsHlp.at(icoIt).at(aIt)[6], &CSyms, settings->axisErrTolerance );
2177 
2178  //==================================== Compute FSC
2179  fscVal = this->computeFSC ( settings, &CSyms, static_cast< size_t > ( matchedPos ), mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts );
2180  ISymsHlp.at(fscMaxInd).at(aIt)[6] = fscVal;
2181  fscValAvg += fscVal;
2182  }
2183 
2184  //==================================== Get FSC average over all axes
2185  fscValAvg /= 31.0;
2186 
2187  //======================================== Is this the best
2188  if ( fscValAvg > fscMax )
2189  {
2190  fscMax = fscValAvg;
2191  fscMaxInd = icoIt;
2192  }
2193  }
2194 
2195  //============================================ If C3 and C5 are found and have correct angle (must have if they are both in ISym)
2196  if ( fscMax >= ( settings->fscThreshold * 0.7 ) )
2197  {
2198  //======================================== Add predicted axes to detected C axes list and also to the settings Icosahedral symmetry list
2199  for ( proshade_unsign retIt = 0; retIt < static_cast < proshade_unsign > ( ISymsHlp.at(fscMaxInd).size() ); retIt++ )
2200  {
2201  //==================================== Add the correct index to the settings object
2202  proshade_signed matchedPos = ProSHADE_internal_symmetry::addAxisUnlessSame ( static_cast< proshade_unsign > ( ISymsHlp.at(fscMaxInd).at(retIt)[0] ), ISymsHlp.at(fscMaxInd).at(retIt)[1], ISymsHlp.at(fscMaxInd).at(retIt)[2], ISymsHlp.at(fscMaxInd).at(retIt)[3], ISymsHlp.at(fscMaxInd).at(retIt)[5], ISymsHlp.at(fscMaxInd).at(retIt)[6], &CSyms, settings->axisErrTolerance );
2203  ProSHADE_internal_misc::addToUnsignVector ( &settings->allDetectedIAxes, static_cast < proshade_unsign > ( matchedPos ) );
2204 
2205  //==================================== Set ISyms for saving
2206  ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( &ISyms, ISymsHlp.at(fscMaxInd).at(retIt) );
2207  }
2208  }
2209 
2210  //============================================ Delete all memory from ISymsHlp
2211  for ( size_t gIt = 0; gIt < ISymsHlp.size(); gIt++ ) { for ( size_t aIt = 0; aIt < ISymsHlp.at(gIt).size(); aIt++ ) { delete[] ISymsHlp.at(gIt).at(aIt); } }
2212 
2213  //============================================ If C3 and C5 are found and have correct angle (must have if they are both in ISym)
2214  if ( fscMax >= ( settings->fscThreshold * 0.7 ) )
2215  {
2216  //======================================== The decision is I
2217  settings->setRecommendedSymmetry ( "I" );
2218  settings->setRecommendedFold ( 0 );
2219  for ( size_t it = 0; it < ISyms.size(); it++ ) { ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( axes, ISyms.at(it) ); }
2220  if ( settings->detectedSymmetry.size() == 0 ) { for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( settings->allDetectedIAxes.size() ); it++ ) { settings->setDetectedSymmetry ( CSyms.at(settings->allDetectedIAxes.at(it)) ); } }
2221  }
2222  }
2223 
2224  //================================================ Save C symmetries to argument and if different from settings, to the settings as well
2225  this->saveDetectedSymmetries ( settings, &CSyms, allCs );
2226 
2227  //================================================ Release memory after FSC computation
2228  delete[] mapData;
2229  delete[] origCoeffs;
2230  delete[] fCoeffs;
2231  fftw_destroy_plan ( planForwardFourier );
2232  delete[] binIndexing;
2233  for (size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ ) { delete[] bindata[binIt]; }
2234  delete[] bindata;
2235  delete[] binCounts;
2236 
2237  //================================================ Done
2238  return ;
2239 
2240 }
2241 
2255 void ProSHADE_internal_data::ProSHADE_data::saveDetectedSymmetries ( ProSHADE_settings* settings, std::vector< proshade_double* >* CSyms, std::vector < std::vector< proshade_double > >* allCs )
2256 {
2257  //================================================ Initialise variables
2258  bool isArgSameAsSettings = true;
2259 
2260  //================================================ For each detected point group
2261  for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSyms->size() ); cIt++ )
2262  {
2263  //============================================ Create vector to replace the pointer
2264  std::vector< proshade_double > nextSym;
2265  ProSHADE_internal_misc::addToDoubleVector ( &nextSym, CSyms->at(cIt)[0] );
2266  ProSHADE_internal_misc::addToDoubleVector ( &nextSym, CSyms->at(cIt)[1] );
2267  ProSHADE_internal_misc::addToDoubleVector ( &nextSym, CSyms->at(cIt)[2] );
2268  ProSHADE_internal_misc::addToDoubleVector ( &nextSym, CSyms->at(cIt)[3] );
2269  ProSHADE_internal_misc::addToDoubleVector ( &nextSym, CSyms->at(cIt)[4] );
2270  ProSHADE_internal_misc::addToDoubleVector ( &nextSym, CSyms->at(cIt)[5] );
2271  ProSHADE_internal_misc::addToDoubleVector ( &nextSym, CSyms->at(cIt)[6] );
2273 
2274  //============================================ Copy the vector to output variable and if different, then also to settings object
2275  if ( ( cIt == 0 ) && ( settings->allDetectedCAxes.size() == 0 ) ) { isArgSameAsSettings = false; }
2276  if ( !isArgSameAsSettings ) { ProSHADE_internal_misc::addToDoubleVectorVector ( &settings->allDetectedCAxes, nextSym ); }
2277 
2278  //============================================ Release memory
2279  nextSym.clear ( );
2280  delete[] CSyms->at(cIt);
2281  }
2282 
2283  //================================================ Done
2284  return ;
2285 
2286 }
2287 
2301 proshade_double ProSHADE_internal_data::ProSHADE_data::findTopGroupSmooth ( std::vector< proshade_double* >* CSym, size_t peakPos, proshade_double step, proshade_double sigma, proshade_signed windowSize )
2302 {
2303  //================================================ Initialise local variables
2304  proshade_double threshold = 0.0;
2305  proshade_signed totSize = static_cast< proshade_signed > ( ( 1.0 / step ) + 1 );
2306  std::vector< std::pair < proshade_double, proshade_unsign > > vals;
2307  std::vector< proshade_double > hist ( static_cast< unsigned long int > ( totSize ), 0.0 );
2308  proshade_unsign histPos = 0;
2309 
2310  //================================================ Make sure window size is odd
2311  if ( windowSize % 2 == 0 ) { windowSize += 1; }
2312 
2313  //================================================ Get vector of pairs of peak heights and indices in CSym array
2314  for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( CSym->size() ); symIt++ ) { vals.emplace_back ( std::pair < proshade_double, proshade_unsign > ( CSym->at(symIt)[peakPos], symIt ) ); }
2315 
2316  //================================================ Convert all found heights to histogram from 0.0 to 1.0 by step
2317  for ( proshade_double it = 0.0; it <= 1.0; it = it + step )
2318  {
2319  for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( vals.size() ); symIt++ )
2320  {
2321  //======================================== Is this height in the range?
2322  if ( ( vals.at(symIt).first > it ) && ( vals.at(symIt).first <= ( it + step ) ) ) { hist.at(histPos) += 1.0; }
2323  }
2324 
2325  //============================================ Update counter and continue
2326  histPos += 1;
2327  }
2328 
2329  //================================================ Smoothen the distribution
2330  std::vector< proshade_double > smoothened = ProSHADE_internal_maths::smoothen1D ( step, windowSize, sigma, hist );
2331 
2332  //================================================ Find peaks in smoothened data
2333  std::vector< proshade_signed > peaks = ProSHADE_internal_peakSearch::findPeaks1D ( smoothened );
2334 
2335  //================================================ Take best peaks surroundings and produce a new set of "high" axes
2336  proshade_signed bestHistPos;
2337  if ( peaks.size() > 0 ) { bestHistPos = peaks.at(peaks.size()-1) + ( ( windowSize - 1 ) / 2 ); }
2338  else { bestHistPos = 0.0; }
2339 
2340  threshold = ( static_cast< proshade_double > ( bestHistPos ) * step ) - ( static_cast< proshade_double > ( windowSize ) * step );
2341 
2342  //================================================ Done
2343  return ( threshold );
2344 
2345 }
2346 
2357 void ProSHADE_internal_data::ProSHADE_data::prepareFSCFourierMemory ( fftw_complex*& mapData, fftw_complex*& origCoeffs, fftw_complex*& fCoeffs, proshade_signed*& binIndexing, proshade_signed* noBins, proshade_double**& bindata, proshade_signed*& binCounts, fftw_plan* planForwardFourier )
2358 {
2359  //================================================ Decide number of bins and allocate which reflection belongs to which bin
2360  ProSHADE_internal_maths::binReciprocalSpaceReflections ( this->xDimIndices, this->yDimIndices, this->zDimIndices, noBins, binIndexing );
2361 
2362  //================================================ Allocate memory for FSC sums
2363  bindata = new proshade_double*[*noBins];
2364  binCounts = new proshade_signed [*noBins];
2365 
2366  //================================================ Allcate memory for bin sumation
2367  for ( size_t binIt = 0; binIt < static_cast< size_t > ( *noBins ); binIt++ )
2368  {
2369  bindata[binIt] = new proshade_double[12];
2370  ProSHADE_internal_misc::checkMemoryAllocation ( bindata[binIt], __FILE__, __LINE__, __func__ );
2371  }
2372 
2373  //================================================ Allocate memory for Fourier transform imputs and outputs
2374  mapData = new fftw_complex [this->xDimIndices * this->yDimIndices * this->zDimIndices];
2375  origCoeffs = new fftw_complex [this->xDimIndices * this->yDimIndices * this->zDimIndices];
2376  fCoeffs = new fftw_complex [this->xDimIndices * this->yDimIndices * this->zDimIndices];
2377 
2378  //================================================ Check memory allocation
2379  ProSHADE_internal_misc::checkMemoryAllocation ( mapData, __FILE__, __LINE__, __func__ );
2380  ProSHADE_internal_misc::checkMemoryAllocation ( origCoeffs, __FILE__, __LINE__, __func__ );
2381  ProSHADE_internal_misc::checkMemoryAllocation ( fCoeffs, __FILE__, __LINE__, __func__ );
2382  ProSHADE_internal_misc::checkMemoryAllocation ( bindata, __FILE__, __LINE__, __func__ );
2383  ProSHADE_internal_misc::checkMemoryAllocation ( binCounts, __FILE__, __LINE__, __func__ );
2384 
2385  //================================================ Prepare memory for Fourier transform
2386  *planForwardFourier = fftw_plan_dft_3d ( static_cast< int > ( this->xDimIndices ), static_cast< int > ( this->yDimIndices ), static_cast< int > ( this->zDimIndices ), mapData, fCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
2387 
2388  //================================================ Compute Fourier transform of the original map
2389  for ( size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { mapData[iter][0] = this->internalMap[iter]; mapData[iter][1] = 0.0; }
2390  fftw_execute ( *planForwardFourier );
2391  for ( size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { origCoeffs[iter][0] = fCoeffs[iter][0]; origCoeffs[iter][1] = fCoeffs[iter][1]; }
2392 
2393  //================================================ Done
2394  return ;
2395 
2396 }
2397 
2425 proshade_double ProSHADE_internal_data::ProSHADE_data::computeFSC ( ProSHADE_settings* settings, std::vector< proshade_double* >* CSym, size_t symIndex, fftw_complex* mapData, fftw_complex* fCoeffs, fftw_complex* origCoeffs, fftw_plan* planForwardFourier, proshade_signed noBins, proshade_signed *binIndexing, proshade_double**& bindata, proshade_signed*& binCounts )
2426 {
2427  //================================================ Sanity check
2428  if ( symIndex >= CSym->size() )
2429  {
2430  std::cerr << "The supplied symmetry axes vector does not contain element number " << symIndex << ". Returning FSC 0.0." << std::endl;
2431  return ( -2.0 );
2432  }
2433 
2434  //================================================ Ignore if already computed
2435  const FloatingPoint< proshade_double > lhs1 ( CSym->at(symIndex)[6] ), rhs1 ( -1.0 );
2436  if ( !lhs1.AlmostEquals ( rhs1 ) ) { return ( CSym->at(symIndex)[6] ); }
2437 
2438  //================================================ Report progress
2439  std::stringstream ss2;
2440  ss2 << "Computing FSC for symmetry C" << CSym->at(symIndex)[0] << " ( " << CSym->at(symIndex)[1] << " ; " << CSym->at(symIndex)[2] << " ; " << CSym->at(symIndex)[3] << " ) with peak height " << CSym->at(symIndex)[5];
2441  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 4, ss2.str() );
2442 
2443  //================================================ Initialise local variables
2444  proshade_double *rotMap;
2445 
2446  //================================================ For each rotation along the axis
2447  proshade_double averageFSC = 0.0;
2448  for ( proshade_double rotIter = 1.0; rotIter < CSym->at(symIndex)[0]; rotIter += 1.0 )
2449  {
2450  //============================================ Get rotated map by the smallest fold angle along the symmetry axis
2451  this->rotateMapRealSpace ( CSym->at(symIndex)[1], CSym->at(symIndex)[2], CSym->at(symIndex)[3], ( ( 2.0 * M_PI ) / CSym->at(symIndex)[0] ) * rotIter, rotMap );
2452 
2453  //============================================ Get Fourier for the rotated map
2454  for ( size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { mapData[iter][0] = rotMap[iter]; mapData[iter][1] = 0.0; }
2455  fftw_execute ( *planForwardFourier );
2456 
2457  //============================================ Clean FSC computation memory
2458  for ( size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ ) { for ( size_t valIt = 0; valIt < 12; valIt++ ) { bindata[binIt][valIt] = 0.0; } }
2459  for ( size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ ) { binCounts[binIt] = 0; }
2460 
2461  //============================================ Compute FSC
2462  averageFSC += ProSHADE_internal_maths::computeFSC ( origCoeffs, fCoeffs, this->xDimIndices, this->yDimIndices, this->zDimIndices, noBins, binIndexing, bindata, binCounts );
2463 
2464  //============================================ Release memory
2465  delete[] rotMap;
2466  }
2467 
2468  //================================================ Convert sum to average
2469  averageFSC /= ( CSym->at(symIndex)[0] - 1.0 );
2470 
2471  //================================================ Save result to the axis
2472  CSym->at(symIndex)[6] = averageFSC;
2473 
2474  //================================================ Report progress
2475  std::stringstream ss3;
2476  ss3 << "FSC value is " << averageFSC << " .";
2477  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 5, ss3.str() );
2478 
2479  //================================================ Done
2480  return ( averageFSC );
2481 
2482 }
2483 
2510 proshade_double ProSHADE_internal_data::ProSHADE_data::computeFSC ( ProSHADE_settings* settings, proshade_double* sym, fftw_complex* mapData, fftw_complex* fCoeffs, fftw_complex* origCoeffs, fftw_plan* planForwardFourier, proshade_signed noBins, proshade_signed *binIndexing, proshade_double**& bindata, proshade_signed*& binCounts )
2511 {
2512  //================================================ Ignore if already computed
2513  const FloatingPoint< proshade_double > lhs1 ( sym[6] ), rhs1 ( -1.0 );
2514  if ( !lhs1.AlmostEquals ( rhs1 ) ) { return ( sym[6] ); }
2515 
2516  //================================================ Report progress
2517  std::stringstream ss2;
2518  ss2 << "Computing FSC for symmetry C" << sym[0] << " ( " << sym[1] << " ; " << sym[2] << " ; " << sym[3] << " ) with peak height " << sym[5];
2519  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 4, ss2.str() );
2520 
2521  //================================================ Initialise local variables
2522  proshade_double *rotMap;
2523 
2524  //================================================ For each rotation along the axis
2525  proshade_double averageFSC = 0.0;
2526  for ( proshade_double rotIter = 1.0; rotIter < sym[0]; rotIter += 1.0 )
2527  {
2528  //============================================ Get rotated map by the smallest fold angle along the symmetry axis
2529  this->rotateMapRealSpace ( sym[1], sym[2], sym[3], ( ( 2.0 * M_PI ) / sym[0] ) * rotIter, rotMap );
2530 
2531  //============================================ Get Fourier for the rotated map
2532  for ( size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { mapData[iter][0] = rotMap[iter]; mapData[iter][1] = 0.0; }
2533  fftw_execute ( *planForwardFourier );
2534 
2535  //============================================ Clean FSC computation memory
2536  for ( size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ ) { for ( size_t valIt = 0; valIt < 12; valIt++ ) { bindata[binIt][valIt] = 0.0; } }
2537  for ( size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ ) { binCounts[binIt] = 0; }
2538 
2539  //============================================ Compute FSC
2540  averageFSC += ProSHADE_internal_maths::computeFSC ( origCoeffs, fCoeffs, this->xDimIndices, this->yDimIndices, this->zDimIndices, noBins, binIndexing, bindata, binCounts );
2541 
2542  //============================================ Release memory
2543  delete[] rotMap;
2544  }
2545 
2546  //================================================ Convert sum to average
2547  averageFSC /= ( sym[0] - 1.0 );
2548 
2549  //================================================ Save result to the axis
2550  sym[6] = averageFSC;
2551 
2552  //================================================ Report progress
2553  std::stringstream ss3;
2554  ss3 << "FSC value is " << averageFSC << " .";
2555  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 5, ss3.str() );
2556 
2557  //================================================ Done
2558  return ( averageFSC );
2559 
2560 }
2561 
2596 void ProSHADE_internal_data::ProSHADE_data::saveRecommendedSymmetry ( ProSHADE_settings* settings, std::vector< proshade_double* >* CSym, std::vector< proshade_double* >* DSym, std::vector< proshade_double* >* TSym, std::vector< proshade_double* >* OSym, std::vector< proshade_double* >* ISym, std::vector< proshade_double* >* axes, fftw_complex* mapData, fftw_complex* origCoeffs, fftw_complex* fCoeffs, fftw_plan* planForwardFourier, proshade_signed noBins, proshade_signed* binIndexing, proshade_double** bindata, proshade_signed* binCounts )
2597 {
2598  //================================================ Report progress
2599  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting recommended symmetry decision procedure." );
2600 
2601  //================================================ If no C symmetries, nothing to save...
2602  if ( CSym->size() == 0 )
2603  {
2604  settings->setRecommendedSymmetry ( "" );
2605  settings->setRecommendedFold ( 0 );
2606  return;
2607  }
2608 
2609  //================================================ Find the top group minimum threshold using smoothened histogram
2610  proshade_double step = 0.01;
2611  proshade_double sigma = 0.03;
2612  proshade_signed windowSize = 9;
2613  proshade_double bestHistPeakStart = this->findTopGroupSmooth ( CSym, 5, step, sigma, windowSize );
2614  if ( bestHistPeakStart > settings->peakThresholdMin ) { bestHistPeakStart = settings->peakThresholdMin; }
2615 
2616  //================================================ Report progress
2617  proshade_unsign noPassed = 0; for ( size_t cIt = 0; cIt < CSym->size(); cIt++ ) { if ( CSym->at(cIt)[5] > bestHistPeakStart ) { noPassed += 1; } }
2618  std::stringstream ss;
2619  ss << "Smoothening has resolved in " << noPassed << " C symmetries.";
2620  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, ss.str() );
2621  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Starting FSC computation to confirm the C symmetries existence." );
2622 
2623  //================================================ Decide if I is the answer
2624  bool alreadyDecided = false;
2625  if ( ISym->size() == 31 )
2626  {
2627  //============================================ Initialise decision vars
2628  proshade_double fscVal = 0.0;
2629  proshade_double fscValAvg = 0.0;
2630 
2631  //============================================ Check if at least one C5 and one C3 with the correct angle have high FSC and peak height
2632  for ( size_t iIt = 0; iIt < 31; iIt++ ) { if ( CSym->at(settings->allDetectedIAxes.at(iIt))[5] > bestHistPeakStart ) { fscVal = this->computeFSC ( settings, CSym, settings->allDetectedIAxes.at(iIt), mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts ); fscValAvg += fscVal; } }
2633  fscValAvg /= 31.0;
2634 
2635  //============================================ If C3 and C5 are found and have correct angle (must have if they are both in ISym)
2636  if ( fscValAvg >= ( settings->fscThreshold * 0.60 ) )
2637  {
2638  //======================================== The decision is I
2639  settings->setRecommendedSymmetry ( "I" );
2640  settings->setRecommendedFold ( 0 );
2641  for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( settings->allDetectedIAxes.size() ); it++ ) { ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( axes, CSym->at(settings->allDetectedIAxes.at(it)) ); }
2642  if ( settings->detectedSymmetry.size() == 0 ) { for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( settings->allDetectedIAxes.size() ); it++ ) { settings->setDetectedSymmetry ( CSym->at(settings->allDetectedIAxes.at(it)) ); } }
2643 
2644  //======================================== Done
2645  alreadyDecided = true;
2646  }
2647  }
2648 
2649  //================================================ Decide if O is the answer
2650  if ( ( OSym->size() == 13 ) && !alreadyDecided )
2651  {
2652  //============================================ Initialise decision vars
2653  proshade_double fscVal = 0.0;
2654  proshade_double fscValAvg = 0.0;
2655 
2656  //============================================ Check if at least one C5 and one C3 with the correct angle have high FSC and peak height
2657  for ( size_t oIt = 0; oIt < 13; oIt++ ) { if ( CSym->at(settings->allDetectedOAxes.at(oIt))[5] > bestHistPeakStart ) { fscVal = this->computeFSC ( settings, CSym, settings->allDetectedOAxes.at(oIt), mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts ); fscValAvg += fscVal; } }
2658  fscValAvg /= 13.0;
2659 
2660  //============================================ If C3 and C5 are found and have correct angle (must have if they are both in ISym)
2661  if ( fscValAvg >= ( settings->fscThreshold * 0.70 ) )
2662  {
2663  //======================================== The decision is O
2664  settings->setRecommendedSymmetry ( "O" );
2665  settings->setRecommendedFold ( 0 );
2666  for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( settings->allDetectedOAxes.size() ); it++ ) { ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( axes, CSym->at(settings->allDetectedOAxes.at(it)) ); }
2667  if ( settings->detectedSymmetry.size() == 0 ) { for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( settings->allDetectedOAxes.size() ); it++ ) { settings->setDetectedSymmetry ( CSym->at(settings->allDetectedOAxes.at(it)) ); } }
2668 
2669  //======================================== Done
2670  alreadyDecided = true;
2671  }
2672  }
2673 
2674  //================================================ Decide if T is the answer
2675  if ( ( TSym->size() == 7 ) && !alreadyDecided )
2676  {
2677  //============================================ Initialise decision vars
2678  proshade_double fscVal = 0.0;
2679  proshade_double fscValAvg = 0.0;
2680 
2681  //============================================ Check if at least one C5 and one C3 with the correct angle have high FSC and peak height
2682  for ( size_t tIt = 0; tIt < 7; tIt++ ) { if ( CSym->at(settings->allDetectedTAxes.at(tIt))[5] > bestHistPeakStart ) { fscVal = this->computeFSC ( settings, CSym, settings->allDetectedTAxes.at(tIt), mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts ); fscValAvg += fscVal; } }
2683  fscValAvg /= 7.0;
2684 
2685  //============================================ If C3 and C5 are found and have correct angle (must have if they are both in ISym)
2686  if ( fscValAvg >= ( settings->fscThreshold * 0.80 ) )
2687  {
2688  //======================================== The decision is T
2689  settings->setRecommendedSymmetry ( "T" );
2690  settings->setRecommendedFold ( 0 );
2691  for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( settings->allDetectedTAxes.size() ); it++ ) { ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( axes, CSym->at(settings->allDetectedTAxes.at(it)) ); }
2692  if ( settings->detectedSymmetry.size() == 0 ) { for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( settings->allDetectedTAxes.size() ); it++ ) { settings->setDetectedSymmetry ( CSym->at(settings->allDetectedTAxes.at(it)) ); } }
2693 
2694  //======================================== Done
2695  alreadyDecided = true;
2696  }
2697  }
2698 
2699  //================================================ Decide if D is the answer
2700  if ( ( settings->allDetectedDAxes.size() > 0 ) && ( DSym->size() > 0 ) && !alreadyDecided )
2701  {
2702  //============================================ Initialise decision vars
2703  proshade_signed bestD = -1;
2704  proshade_unsign bestFold = 0;
2705 
2706  //============================================ Find FSCs
2707  for ( size_t dIt = 0; dIt < settings->allDetectedDAxes.size(); dIt++ )
2708  {
2709  //======================================== Do not consider more than top 20, takes time and is unlikely to produce anything...
2710  if ( dIt > 20 ) { continue; }
2711 
2712  //======================================== Check the peak heights
2713  const FloatingPoint< proshade_double > lhs999a ( CSym->at(settings->allDetectedDAxes.at(dIt).at(0))[5] ), lhs999b ( CSym->at(settings->allDetectedDAxes.at(dIt).at(1))[5] ), rhs999 ( static_cast< proshade_double > ( -999.9 ) );
2714  if ( ( CSym->at(settings->allDetectedDAxes.at(dIt).at(0))[5] < bestHistPeakStart ) && !( lhs999a.AlmostEquals( rhs999 ) ) ) { continue; }
2715  if ( ( CSym->at(settings->allDetectedDAxes.at(dIt).at(1))[5] < bestHistPeakStart ) && !( lhs999b.AlmostEquals( rhs999 ) ) ) { continue; }
2716 
2717  //======================================== Find FSCs
2718  this->computeFSC ( settings, CSym, settings->allDetectedDAxes.at(dIt).at(0), mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts );
2719  this->computeFSC ( settings, CSym, settings->allDetectedDAxes.at(dIt).at(1), mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts );
2720  }
2721 
2722  //============================================ Find FSC top group threshold
2723  proshade_double bestHistFSCStart = this->findTopGroupSmooth ( CSym, 5, step, sigma, windowSize );
2724 
2725  //============================================ Check if both C symmetries are reliable
2726  for ( size_t dIt = 0; dIt < settings->allDetectedDAxes.size(); dIt++ )
2727  {
2728  //======================================== Check the peak heights
2729  const FloatingPoint< proshade_double > lhs999a2 ( CSym->at(settings->allDetectedDAxes.at(dIt).at(0))[5] ), lhs999b2 ( CSym->at(settings->allDetectedDAxes.at(dIt).at(1))[5] ), rhs999 ( static_cast< proshade_double > ( -999.9 ) );
2730  if ( ( CSym->at(settings->allDetectedDAxes.at(dIt).at(0))[5] < bestHistPeakStart ) && !( lhs999a2.AlmostEquals( rhs999 ) ) ) { continue; }
2731  if ( ( CSym->at(settings->allDetectedDAxes.at(dIt).at(1))[5] < bestHistPeakStart ) && !( lhs999b2.AlmostEquals( rhs999 ) ) ) { continue; }
2732 
2733  //======================================== Does this improve the best fold?
2734  if ( ( CSym->at(settings->allDetectedDAxes.at(dIt).at(0))[0] > static_cast< proshade_double > ( bestFold ) ) || ( CSym->at(settings->allDetectedDAxes.at(dIt).at(1))[0] > static_cast< proshade_double > ( bestFold ) ) )
2735  {
2736  //==================================== Check the FSC vals
2737  if ( CSym->at(settings->allDetectedDAxes.at(dIt).at(0))[6] < settings->fscThreshold ) { continue; }
2738  if ( CSym->at(settings->allDetectedDAxes.at(dIt).at(1))[6] < settings->fscThreshold ) { continue; }
2739  if ( std::max ( CSym->at(settings->allDetectedDAxes.at(dIt).at(0))[6], CSym->at(settings->allDetectedDAxes.at(dIt).at(1))[6] ) < bestHistFSCStart ) { continue; }
2740 
2741  //==================================== All good!
2742  bestFold = static_cast< proshade_unsign > ( std::max ( CSym->at(settings->allDetectedDAxes.at(dIt).at(0))[0], CSym->at(settings->allDetectedDAxes.at(dIt).at(1))[0] ) );
2743  bestD = static_cast< proshade_signed > ( dIt );
2744  }
2745  }
2746 
2747  //============================================ Anything?
2748  if ( bestD != -1 )
2749  {
2750  //======================================== The decision is D
2751  settings->setRecommendedSymmetry ( "D" );
2752  settings->setRecommendedFold ( static_cast< proshade_unsign > ( std::max ( CSym->at(settings->allDetectedDAxes.at( static_cast< size_t > ( bestD ) ).at(0))[0], CSym->at(settings->allDetectedDAxes.at( static_cast< size_t > ( bestD ) ).at(1))[0] ) ) );
2753  ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( axes, CSym->at(settings->allDetectedDAxes.at( static_cast< size_t > ( bestD ) ).at(0)) );
2754  ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( axes, CSym->at(settings->allDetectedDAxes.at( static_cast< size_t > ( bestD ) ).at(1)) );
2755  if ( settings->detectedSymmetry.size() == 0 )
2756  {
2757  settings->setDetectedSymmetry ( CSym->at(settings->allDetectedDAxes.at( static_cast< size_t > ( bestD ) ).at(0)) );
2758  settings->setDetectedSymmetry ( CSym->at(settings->allDetectedDAxes.at( static_cast< size_t > ( bestD ) ).at(1)) );
2759  }
2760 
2761  //======================================== Done
2762  alreadyDecided = true;
2763  }
2764  }
2765 
2766  //================================================ Decide if C is the answer
2767  if ( ( CSym->size() > 0 ) && !alreadyDecided )
2768  {
2769  //============================================ Initialise decision vars
2770  proshade_signed bestC = -1;
2771  proshade_unsign bestFold = 0;
2772 
2773  //============================================ Find FSCs for C syms
2774  for ( size_t cIt = 0; cIt < CSym->size(); cIt++ )
2775  {
2776  //======================================== Do not consider more than top 20, takes time and is unlikely to produce anything...
2777  if ( cIt > 20 ) { continue; }
2778 
2779  //======================================== Check the peak height
2780  if ( CSym->at(cIt)[5] < bestHistPeakStart ) { continue; }
2781 
2782  //======================================== Compute FSC
2783  this->computeFSC ( settings, CSym, cIt, mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts );
2784  }
2785 
2786  //============================================ Find FSC top group threshold
2787  proshade_double bestHistFSCStart = this->findTopGroupSmooth ( CSym, 6, step, sigma, windowSize );
2788 
2789  //============================================ Find reliable C syms
2790  for ( size_t cIt = 0; cIt < CSym->size(); cIt++ )
2791  {
2792  //======================================== Check if this improves the best already found fold
2793  if ( CSym->at(cIt)[0] > static_cast< proshade_double > ( bestFold ) )
2794  {
2795  //==================================== If FSC passes
2796  if ( ( CSym->at(cIt)[6] > settings->fscThreshold ) && ( CSym->at(cIt)[6] >= bestHistFSCStart ) )
2797  {
2798  bestFold = static_cast< proshade_unsign > ( CSym->at(cIt)[0] );
2799  bestC = static_cast< proshade_signed > ( cIt );
2800  }
2801  }
2802  }
2803 
2804  //============================================ Anything?
2805  if ( bestC != -1 )
2806  {
2807  //======================================== The decision is C
2808  settings->setRecommendedSymmetry ( "C" );
2809  settings->setRecommendedFold ( static_cast< proshade_unsign > ( CSym->at( static_cast< size_t > ( bestC ) )[0] ) );
2810  ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( axes, CSym->at( static_cast< size_t > ( bestC ) ) );
2811  if ( settings->detectedSymmetry.size() == 0 ) { settings->setDetectedSymmetry ( CSym->at( static_cast< size_t > ( bestC ) ) ); }
2812 
2813  //======================================== Done
2814  alreadyDecided = true;
2815  }
2816  }
2817 
2818  //================================================ Done
2819  return ;
2820 
2821 }
2822 
2834 void ProSHADE_internal_data::ProSHADE_data::saveRequestedSymmetryC ( ProSHADE_settings* settings, std::vector< proshade_double* >* CSym, std::vector< proshade_double* >* axes )
2835 {
2836  //================================================ Initialise variables
2837  proshade_unsign bestIndex = 0;
2838  proshade_double highestSym = 0.0;
2839 
2840  //================================================ Search for best fold
2841  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( CSym->size() ); iter++ )
2842  {
2843  //============================================ Check if it is tbe correct fold
2844  const FloatingPoint< proshade_double > lhs1 ( CSym->at(iter)[0] ), rhs1 ( static_cast< proshade_double > ( settings->requestedSymmetryFold ) );
2845  if ( !lhs1.AlmostEquals ( rhs1 ) ) { continue; }
2846 
2847  //============================================ If correct, is it the highest found?
2848  if ( CSym->at(iter)[5] > highestSym )
2849  {
2850  highestSym = CSym->at(iter)[5];
2851  bestIndex = iter;
2852  }
2853  }
2854 
2855  //================================================ Found?
2856  if ( highestSym > 0.0 )
2857  {
2858  settings->setRecommendedSymmetry ( "C" );
2859  settings->setRecommendedFold ( static_cast< proshade_unsign > ( CSym->at(bestIndex)[0] ) );
2860  ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( axes, CSym->at(bestIndex) );
2861 
2862  if ( settings->detectedSymmetry.size() == 0 ) { settings->setDetectedSymmetry ( CSym->at(bestIndex) ); }
2863  }
2864  else
2865  {
2866  settings->setRecommendedSymmetry ( "" );
2867  settings->setRecommendedFold ( 0 );
2868  }
2869 
2870  //================================================ Done
2871  return ;
2872 
2873 }
2874 
2886 void ProSHADE_internal_data::ProSHADE_data::saveRequestedSymmetryD ( ProSHADE_settings* settings, std::vector< proshade_double* >* DSym, std::vector< proshade_double* >* axes, fftw_complex* mapData, fftw_complex* origCoeffs, fftw_complex* fCoeffs, fftw_plan* planForwardFourier, proshade_signed noBins, proshade_signed* binIndexing, proshade_double** bindata, proshade_signed* binCounts )
2887 {
2888  //================================================ Initialise variables
2889  proshade_unsign bestIndex = 0;
2890  proshade_double highestSym = 0.0;
2891 
2892  //================================================ Search for best fold
2893  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( DSym->size() ); iter++ )
2894  {
2895  //============================================ Check if it is tbe correct fold
2896  const FloatingPoint< proshade_double > lhs1 ( std::max ( DSym->at(iter)[0], DSym->at(iter)[7] ) ), rhs1 ( static_cast< proshade_double > ( settings->requestedSymmetryFold ) );
2897  if ( !lhs1.AlmostEquals ( rhs1 ) ) { continue; }
2898 
2899  //============================================ Check if peak height is decent
2900  const FloatingPoint< proshade_double > lhs999a ( DSym->at(iter)[5] ), lhs999b ( DSym->at(iter)[12] ), rhs999 ( static_cast< proshade_double > ( -999.9 ) );
2901  if ( ( DSym->at(iter)[5] < settings->peakThresholdMin ) && !( lhs999a.AlmostEquals( rhs999 ) ) ) { continue; }
2902  if ( ( DSym->at(iter)[12] < settings->peakThresholdMin ) && !( lhs999b.AlmostEquals( rhs999 ) ) ) { continue; }
2903 
2904  //============================================ If correct, compute FSC
2905  this->computeFSC ( settings, &DSym->at(iter)[0], mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts );
2906  this->computeFSC ( settings, &DSym->at(iter)[7], mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts );
2907 
2908  //============================================ If best, store it
2909  if ( ( DSym->at(iter)[6] + DSym->at(iter)[13] ) > highestSym )
2910  {
2911  highestSym = ( DSym->at(iter)[6] + DSym->at(iter)[13] );
2912  bestIndex = iter;
2913  }
2914  }
2915 
2916  //================================================ Found?
2917  if ( highestSym > 0.0 )
2918  {
2919  settings->setRecommendedSymmetry ( "D" );
2920  settings->setRecommendedFold ( static_cast< proshade_unsign > ( std::max ( DSym->at(bestIndex)[0], DSym->at(bestIndex)[7] ) ) );
2921  ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( axes, &DSym->at(bestIndex)[0] );
2922  ProSHADE_internal_misc::deepCopyAxisToDblPtrVector ( axes, &DSym->at(bestIndex)[7] );
2923 
2924  if ( settings->detectedSymmetry.size() == 0 )
2925  {
2926  settings->setDetectedSymmetry ( &DSym->at(bestIndex)[0] );
2927  settings->setDetectedSymmetry ( &DSym->at(bestIndex)[7] );
2928  }
2929  }
2930  else
2931  {
2932  settings->setRecommendedSymmetry ( "" );
2933  settings->setRecommendedFold ( 0 );
2934  }
2935 
2936  //================================================ Done
2937  return ;
2938 
2939 }
2940 
2949 std::vector<std::vector< proshade_double > > ProSHADE_internal_data::computeGroupElementsForGroup ( proshade_double xAx, proshade_double yAx, proshade_double zAx, proshade_signed fold )
2950 {
2951  //================================================ Initialise variables
2952  std::vector< proshade_double > angList;
2953  std::vector<std::vector< proshade_double > > ret;
2954 
2955  //================================================ Allocate memory
2956  proshade_double* rotMat = new proshade_double[9];
2957  ProSHADE_internal_misc::checkMemoryAllocation ( rotMat, __FILE__, __LINE__, __func__ );
2958 
2959 
2960  //================================================ Normalise the axis to have magnitude of 1.0
2961  proshade_double normF = std::sqrt( std::pow ( xAx, 2.0 ) + std::pow ( yAx, 2.0 ) + std::pow ( zAx, 2.0 ) );
2962  xAx /= normF;
2963  yAx /= normF;
2964  zAx /= normF;
2965 
2966  //================================================ Determine the list of angles
2967  if ( fold % 2 == 0 )
2968  {
2969  //============================================ If fold is even, add the negative angles
2970  for ( proshade_double iter = static_cast < proshade_double > ( -( ( fold / 2 ) - 1 ) ); iter <= static_cast < proshade_double > ( fold / 2 ); iter++ )
2971  {
2972  ProSHADE_internal_misc::addToDoubleVector ( &angList, ( ( 2.0 * M_PI ) / static_cast<proshade_double> ( fold ) ) * iter );
2973  }
2974  }
2975  else
2976  {
2977  //============================================ If fold is odd, do the same as for even, but start one index earlier
2978  for ( proshade_double iter = static_cast < proshade_double > ( -fold / 2 ); iter <= static_cast < proshade_double > ( fold / 2 ); iter++ )
2979  {
2980  ProSHADE_internal_misc::addToDoubleVector ( &angList, ( ( 2.0 * M_PI ) / static_cast<proshade_double> ( fold ) ) * iter );
2981  }
2982  }
2983 
2984  //================================================ For each detected angle
2985  for ( proshade_unsign iter = 0; iter < static_cast < proshade_unsign > ( angList.size() ); iter++ )
2986  {
2987  //============================================ Compute the rotation matrix
2988  ProSHADE_internal_maths::getRotationMatrixFromAngleAxis ( rotMat, xAx, yAx, zAx, angList.at(iter) );
2989 
2990  //============================================ Convert to vector of vectors of doubles and save to ret
2991  std::vector < proshade_double > retEl;
2992  for ( proshade_unsign matIt = 0; matIt < 9; matIt++ )
2993  {
2994  ProSHADE_internal_misc::addToDoubleVector ( &retEl, rotMat[matIt] );
2995  }
2997  }
2998 
2999  //================================================ Release memory
3000  delete[] rotMat;
3001 
3002  //================================================ Done
3003  return ( ret );
3004 
3005 }
3006 
3013 void axesToGroupTypeSanityCheck ( proshade_unsign requiredAxes, proshade_unsign obtainedAxes, std::string groupType )
3014 {
3015  //================================================ Sanity check
3016  if ( obtainedAxes != requiredAxes )
3017  {
3018  std::stringstream hlpSS;
3019  hlpSS << "The supplied number of axes for group element\n : detection ( >" << obtainedAxes << "< ) does not match the group type ( >" << groupType << "< ).";
3020  throw ProSHADE_exception ( "Mismatch between supplied number of axes and\n : symmetry type.", "ES00059", __FILE__, __LINE__, __func__, hlpSS.str() );
3021  }
3022 
3023  //================================================ Done
3024  return ;
3025 
3026 }
3027 
3035 bool checkElementAlreadyExists ( std::vector<std::vector< proshade_double > >* elements, std::vector< proshade_double >* elem, proshade_double matrixTolerance )
3036 {
3037  //================================================ Initialise variables
3038  bool elementFound = false;
3039 
3040  //================================================ For each existing element
3041  for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( elements->size() ); elIt++ )
3042  {
3043  if ( ProSHADE_internal_maths::rotationMatrixSimilarity ( &elements->at(elIt), elem, matrixTolerance ) )
3044  {
3045  elementFound = true;
3046  break;
3047  }
3048  }
3049 
3050  //================================================ Done
3051  return ( elementFound );
3052 
3053 }
3054 
3061 bool checkElementsFormGroup ( std::vector<std::vector< proshade_double > >* elements, proshade_double matrixTolerance )
3062 {
3063  //================================================ Initialise variables
3064  bool isGroup = true;
3065 
3066  //================================================ Multiply all group element pairs
3067  for ( proshade_unsign gr1 = 0; gr1 < static_cast<proshade_unsign> ( elements->size() ); gr1++ )
3068  {
3069  for ( proshade_unsign gr2 = 1; gr2 < static_cast<proshade_unsign> ( elements->size() ); gr2++ )
3070  {
3071  //======================================== Use unique pairs only
3072  if ( gr1 >= gr2 ) { continue; }
3073 
3074  //======================================== Multiply the two rotation matrices
3075  std::vector< proshade_double > product = ProSHADE_internal_maths::multiplyGroupElementMatrices ( &elements->at(gr1), &elements->at(gr2) );
3076 
3077  //======================================== Check the group already contains the produces as an element
3078  if ( !checkElementAlreadyExists ( elements, &product, matrixTolerance ) )
3079  {
3080  isGroup = false;
3081  break;
3082  }
3083  }
3084 
3085  //============================================ Stop if problem was found
3086  if ( !isGroup ) { break; }
3087  }
3088 
3089  //================================================ Done
3090  return ( isGroup );
3091 
3092 }
3093 
3102 std::vector<std::vector< proshade_double > > ProSHADE_internal_data::joinElementsFromDifferentGroups ( std::vector<std::vector< proshade_double > >* first, std::vector<std::vector< proshade_double > >* second, proshade_double matrixTolerance, bool combine )
3103 {
3104  //================================================ Initialise variables
3105  std::vector< std::vector< proshade_double > > ret;
3106 
3107  //================================================ Add the first list to ret, checking for uniqueness
3108  for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( first->size() ); elIt++ )
3109  {
3110  if ( !checkElementAlreadyExists( &ret, &first->at(elIt), matrixTolerance ) )
3111  {
3112  ProSHADE_internal_misc::addToDoubleVectorVector ( &ret, first->at(elIt) );
3113  }
3114  }
3115 
3116  //================================================ Add the second list to ret, checking for uniqueness
3117  for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( second->size() ); elIt++ )
3118  {
3119  if ( !checkElementAlreadyExists( &ret, &second->at(elIt), matrixTolerance ) )
3120  {
3121  ProSHADE_internal_misc::addToDoubleVectorVector ( &ret, second->at(elIt) );
3122  }
3123  }
3124 
3125  //================================================ Multiply all combinations of first and second and check for uniqueness
3126  if ( combine )
3127  {
3128  for ( proshade_unsign gr1 = 0; gr1 < static_cast<proshade_unsign> ( first->size() ); gr1++ )
3129  {
3130  for ( proshade_unsign gr2 = 0; gr2 < static_cast<proshade_unsign> ( second->size() ); gr2++ )
3131  {
3132  //==================================== Multiply the two rotation matrices
3133  std::vector< proshade_double > product = ProSHADE_internal_maths::multiplyGroupElementMatrices ( &first->at(gr1), &second->at(gr2) );
3134 
3135  //==================================== Add
3136  if ( !checkElementAlreadyExists( &ret, &product, matrixTolerance ) )
3137  {
3139  }
3140 
3141  }
3142  }
3143  }
3144 
3145  //================================================ Done
3146  return ( ret );
3147 
3148 }
3149 
3168 std::vector<std::vector< proshade_double > > ProSHADE_internal_data::ProSHADE_data::getAllGroupElements ( ProSHADE_settings* settings, std::vector< proshade_unsign > axesList, std::string groupType, proshade_double matrixTolerance )
3169 {
3170  //================================================ Initialise variables
3171  std::vector<std::vector< proshade_double > > ret;
3172 
3173  //================================================ Select which symmetry type are we computing for
3174  if ( groupType == "C" )
3175  {
3176  //============================================ Sanity check
3177  axesToGroupTypeSanityCheck ( 1, static_cast< proshade_unsign > ( axesList.size() ), groupType );
3178 
3179  //============================================ Generate elements
3180  ret = computeGroupElementsForGroup ( settings->allDetectedCAxes.at(axesList.at(0)).at(1),
3181  settings->allDetectedCAxes.at(axesList.at(0)).at(2),
3182  settings->allDetectedCAxes.at(axesList.at(0)).at(3),
3183  static_cast< proshade_signed > ( settings->allDetectedCAxes.at(axesList.at(0)).at(0) ) );
3184 
3185  //============================================ Check the element to form a group
3186  if ( checkElementsFormGroup ( &ret, matrixTolerance ) ) { return ( ret ); }
3187  else
3188  {
3189  throw ProSHADE_exception ( "Computed point group elements do not form a group.", "ES00060", __FILE__, __LINE__, __func__, "The supplied cyclic groups list does not form a group and\n : therefore such group's elements cannot be obtained. Please\n : check the cyclic groups list supplied to the\n : getAllGroupElements() function." );
3190  }
3191  }
3192  else if ( groupType == "D" )
3193  {
3194  //============================================ Sanity check
3195  axesToGroupTypeSanityCheck ( 2, static_cast<proshade_unsign> ( axesList.size() ), groupType );
3196 
3197  //============================================ Generate elements for both axes
3198  std::vector<std::vector< proshade_double > > first = computeGroupElementsForGroup ( settings->allDetectedCAxes.at(axesList.at(0)).at(1),
3199  settings->allDetectedCAxes.at(axesList.at(0)).at(2),
3200  settings->allDetectedCAxes.at(axesList.at(0)).at(3),
3201  static_cast< proshade_signed > ( settings->allDetectedCAxes.at(axesList.at(0)).at(0) ) );
3202  std::vector<std::vector< proshade_double > > second = computeGroupElementsForGroup ( settings->allDetectedCAxes.at(axesList.at(1)).at(1),
3203  settings->allDetectedCAxes.at(axesList.at(1)).at(2),
3204  settings->allDetectedCAxes.at(axesList.at(1)).at(3),
3205  static_cast< proshade_signed > ( settings->allDetectedCAxes.at(axesList.at(1)).at(0) ) );
3206 
3207  //============================================ Join the element lists
3208  ret = joinElementsFromDifferentGroups ( &first, &second, matrixTolerance, true );
3209 
3210  //============================================ Check the element to form a group
3211  if ( checkElementsFormGroup ( &ret, matrixTolerance ) ) { return ( ret ); }
3212  else
3213  {
3214  throw ProSHADE_exception ( "Computed point group elements do not form a group.", "ES00060", __FILE__, __LINE__, __func__, "The supplied cyclic groups list does not form a group and\n : therefore such group's elements cannot be obtained. Please\n : check the cyclic groups list supplied to the\n : getAllGroupElements() function." );
3215  }
3216  }
3217  else if ( groupType == "T" )
3218  {
3219  //============================================ Sanity check
3220  axesToGroupTypeSanityCheck ( 7, static_cast<proshade_unsign> ( axesList.size() ), groupType );
3221 
3222  //============================================ Generate elements for all four C3 axes first
3223  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3224  {
3225  //======================================== If this is a C3 axis
3226  const FloatingPoint< proshade_double > lhs1 ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 3.0 );
3227  if ( lhs1.AlmostEquals ( rhs1 ) )
3228  {
3229  //==================================== Generate the elements
3230  std::vector<std::vector< proshade_double > > els = computeGroupElementsForGroup ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(1),
3231  settings->allDetectedCAxes.at(axesList.at(grIt)).at(2),
3232  settings->allDetectedCAxes.at(axesList.at(grIt)).at(3),
3233  static_cast< proshade_signed > ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3234 
3235  //==================================== Join the elements to any already found
3236  ret = joinElementsFromDifferentGroups ( &els, &ret, matrixTolerance, false );
3237  }
3238  }
3239 
3240  //============================================ Generate elements for all three C2 axes second
3241  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3242  {
3243  //======================================== If this is a C3 axis
3244  const FloatingPoint< proshade_double > lhs1 ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 2.0 );
3245  if ( lhs1.AlmostEquals ( rhs1 ) )
3246  {
3247  //==================================== Generate the elements
3248  std::vector<std::vector< proshade_double > > els = computeGroupElementsForGroup ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(1),
3249  settings->allDetectedCAxes.at(axesList.at(grIt)).at(2),
3250  settings->allDetectedCAxes.at(axesList.at(grIt)).at(3),
3251  static_cast< proshade_signed > ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3252 
3253  //==================================== Join the elements to any already found
3254  ret = joinElementsFromDifferentGroups ( &els, &ret, matrixTolerance, false );
3255  }
3256  }
3257 
3258  //============================================ Check the element to form a group
3259  if ( checkElementsFormGroup ( &ret, matrixTolerance ) ) { return ( ret ); }
3260  else
3261  {
3262  throw ProSHADE_exception ( "Computed point group elements do not form a group.", "ES00060", __FILE__, __LINE__, __func__, "The supplied cyclic groups list does not form a group and\n : therefore such group's elements cannot be obtained. Please\n : check the cyclic groups list supplied to the\n : getAllGroupElements() function." );
3263  }
3264  }
3265  else if ( groupType == "O" )
3266  {
3267  //============================================ Sanity check
3268  axesToGroupTypeSanityCheck ( 13, static_cast<proshade_unsign> ( axesList.size() ), groupType );
3269 
3270  //============================================ Generate elements for all three C4 axes first
3271  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3272  {
3273  //======================================== If this is a C3 axis
3274  const FloatingPoint< proshade_double > lhs1 ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 4.0 );
3275  if ( lhs1.AlmostEquals ( rhs1 ) )
3276  {
3277  //==================================== Generate the elements
3278  std::vector<std::vector< proshade_double > > els = computeGroupElementsForGroup ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(1),
3279  settings->allDetectedCAxes.at(axesList.at(grIt)).at(2),
3280  settings->allDetectedCAxes.at(axesList.at(grIt)).at(3),
3281  static_cast< proshade_signed > ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3282 
3283  //==================================== Join the elements to any already found
3284  ret = joinElementsFromDifferentGroups ( &els, &ret, matrixTolerance, false );
3285  }
3286  }
3287 
3288  //============================================ Generate elements for all four C3 axes first
3289  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3290  {
3291  //======================================== If this is a C3 axis
3292  const FloatingPoint< proshade_double > lhs1 ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 3.0 );
3293  if ( lhs1.AlmostEquals ( rhs1 ) )
3294  {
3295  //==================================== Generate the elements
3296  std::vector<std::vector< proshade_double > > els = computeGroupElementsForGroup ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(1),
3297  settings->allDetectedCAxes.at(axesList.at(grIt)).at(2),
3298  settings->allDetectedCAxes.at(axesList.at(grIt)).at(3),
3299  static_cast< proshade_signed > ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3300 
3301  //==================================== Join the elements to any already found
3302  ret = joinElementsFromDifferentGroups ( &els, &ret, matrixTolerance, false );
3303  }
3304  }
3305 
3306  //============================================ Generate elements for all six C2 axes next
3307  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3308  {
3309  //======================================== If this is a C3 axis
3310  const FloatingPoint< proshade_double > lhs1 ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 2.0 );
3311  if ( lhs1.AlmostEquals ( rhs1 ) )
3312  {
3313  //==================================== Generate the elements
3314  std::vector<std::vector< proshade_double > > els = computeGroupElementsForGroup ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(1),
3315  settings->allDetectedCAxes.at(axesList.at(grIt)).at(2),
3316  settings->allDetectedCAxes.at(axesList.at(grIt)).at(3),
3317  static_cast< proshade_signed > ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3318 
3319  //==================================== Join the elements to any already found
3320  ret = joinElementsFromDifferentGroups ( &els, &ret, matrixTolerance, false );
3321  }
3322  }
3323 
3324  //============================================ Check the element to form a group
3325  if ( checkElementsFormGroup ( &ret, matrixTolerance ) ) { return ( ret ); }
3326  else
3327  {
3328  throw ProSHADE_exception ( "Computed point group elements do not form a group.", "ES00060", __FILE__, __LINE__, __func__, "The supplied cyclic groups list does not form a group and\n : therefore such group's elements cannot be obtained. Please\n : check the cyclic groups list supplied to the\n : getAllGroupElements() function." );
3329  }
3330  }
3331  else if ( groupType == "I" )
3332  {
3333  //============================================ Sanity check
3334  axesToGroupTypeSanityCheck ( 31, static_cast<proshade_unsign> ( axesList.size() ), groupType );
3335 
3336  //============================================ Generate elements for all six C5 axes first
3337  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3338  {
3339  //======================================== If this is a C5 axis
3340  const FloatingPoint< proshade_double > lhs1 ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 5.0 );
3341  if ( lhs1.AlmostEquals ( rhs1 ) )
3342  {
3343  //==================================== Generate the elements
3344  std::vector<std::vector< proshade_double > > els = computeGroupElementsForGroup ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(1),
3345  settings->allDetectedCAxes.at(axesList.at(grIt)).at(2),
3346  settings->allDetectedCAxes.at(axesList.at(grIt)).at(3),
3347  static_cast< proshade_signed > ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3348 
3349  //==================================== Join the elements to any already found
3350  ret = joinElementsFromDifferentGroups ( &els, &ret, matrixTolerance, false );
3351  }
3352  }
3353 
3354  //============================================ Generate elements for all ten C3 axes next
3355  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3356  {
3357  //======================================== If this is a C3 axis
3358  const FloatingPoint< proshade_double > lhs1 ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 3.0 );
3359  if ( lhs1.AlmostEquals ( rhs1 ) )
3360  {
3361  //==================================== Generate the elements
3362  std::vector<std::vector< proshade_double > > els = computeGroupElementsForGroup ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(1),
3363  settings->allDetectedCAxes.at(axesList.at(grIt)).at(2),
3364  settings->allDetectedCAxes.at(axesList.at(grIt)).at(3),
3365  static_cast< proshade_signed > ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3366 
3367  //==================================== Join the elements to any already found
3368  ret = joinElementsFromDifferentGroups ( &els, &ret, matrixTolerance, false );
3369  }
3370  }
3371 
3372  //============================================ Generate elements for all fifteen C2 axes lastly
3373  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3374  {
3375  //======================================== If this is a C3 axis
3376  const FloatingPoint< proshade_double > lhs1 ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 2.0 );
3377  if ( lhs1.AlmostEquals ( rhs1 ) )
3378  {
3379  //==================================== Generate the elements
3380  std::vector<std::vector< proshade_double > > els = computeGroupElementsForGroup ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(1),
3381  settings->allDetectedCAxes.at(axesList.at(grIt)).at(2),
3382  settings->allDetectedCAxes.at(axesList.at(grIt)).at(3),
3383  static_cast< proshade_signed > ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3384 
3385  //==================================== Join the elements to any already found
3386  ret = joinElementsFromDifferentGroups ( &els, &ret, matrixTolerance, false );
3387  }
3388  }
3389 
3390  //============================================ Check the element to form a group
3391  if ( checkElementsFormGroup ( &ret, matrixTolerance ) ) { return ( ret ); }
3392  else
3393  {
3394  throw ProSHADE_exception ( "Computed point group elements do not form a group.", "ES00060", __FILE__, __LINE__, __func__, "The supplied cyclic groups list does not form a group and\n : therefore such group's elements cannot be obtained. Please\n : check the cyclic groups list supplied to the\n : getAllGroupElements() function." );
3395  }
3396  }
3397  else if ( groupType == "X" )
3398  {
3399  //============================================ User forced no checking for unspecified symmetry
3400  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3401  {
3402  //======================================== Compute group elements
3403  std::vector<std::vector< proshade_double > > els = computeGroupElementsForGroup ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(1),
3404  settings->allDetectedCAxes.at(axesList.at(grIt)).at(2),
3405  settings->allDetectedCAxes.at(axesList.at(grIt)).at(3),
3406  static_cast< proshade_signed > ( settings->allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3407 
3408  //======================================== Join the elements to any already found
3409  ret = joinElementsFromDifferentGroups ( &els, &ret, matrixTolerance, true );
3410  }
3411 
3412  //============================================ Check the element to form a group
3413  if ( checkElementsFormGroup ( &ret, matrixTolerance ) ) { return ( ret ); }
3414  else
3415  {
3416  throw ProSHADE_exception ( "Computed point group elements do not form a group.", "ES00060", __FILE__, __LINE__, __func__, "The supplied cyclic groups list does not form a group and\n : therefore such group's elements cannot be obtained. Please\n : check the cyclic groups list supplied to the\n : getAllGroupElements() function." );
3417  }
3418  }
3419  else
3420  {
3421  std::stringstream hlpSS;
3422  hlpSS << "Unknown symmetry type: >" << groupType << "<";
3423  throw ProSHADE_exception ( hlpSS.str().c_str(), "ES00058", __FILE__, __LINE__, __func__, "Function getAllGroupElements was called with symmetry type\n : value outside of the allowed values C, D, T, O, I\n : or empty for using all supplied axes." );
3424  }
3425 
3426 }
3427 
3435 void ProSHADE_internal_data::ProSHADE_data::deepCopyMap ( proshade_double*& saveTo, proshade_signed verbose )
3436 {
3437  //================================================ Sanity check
3438  if ( saveTo != nullptr )
3439  {
3440  ProSHADE_internal_messages::printWarningMessage ( verbose, "!!! ProSHADE WARNING !!! The deep copy pointer is not set to NULL. Cannot proceed and returning unmodified pointer.", "WB00040" );
3441  return ;
3442  }
3443 
3444  //================================================ Allocate the memory
3445  saveTo = new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
3446 
3447  //================================================ Check memory allocation
3448  ProSHADE_internal_misc::checkMemoryAllocation ( saveTo, __FILE__, __LINE__, __func__ );
3449 
3450  //================================================ Copy internal map to the new pointer
3451  for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
3452  {
3453  saveTo[iter] = this->internalMap[iter];
3454  }
3455 
3456  //================================================ Done
3457  return ;
3458 
3459 }
3460 
3468 {
3469  //================================================ Improve this!
3470  if ( settings->recommendedSymmetryType == "" )
3471  {
3472  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 0, "Did not detect any symmetry!" );
3473  }
3474  else
3475  {
3476  std::stringstream ssHlp;
3477  std::vector< proshade_double > comMove = this->getMapCOMProcessChange ( );
3478  ssHlp << std::endl << "Detected " << settings->recommendedSymmetryType << " symmetry with fold " << settings->recommendedSymmetryFold << " about point [" << comMove.at(0) << " , " << comMove.at(1) << " , " << comMove.at(2) << "] away from centre of mass .";
3479  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 0, ssHlp.str() );
3480 
3481  if ( settings->detectedSymmetry.size() > 0 )
3482  {
3483  ssHlp.clear(); ssHlp.str ( "" );
3484  ssHlp << " Fold X Y Z Angle Height Average FSC";
3485  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 0, ssHlp.str() );
3486  }
3487  for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( settings->detectedSymmetry.size() ); symIt++ )
3488  {
3489  ssHlp.clear(); ssHlp.str ( "" );
3490  ssHlp << std::showpos << std::fixed << std::setprecision(0) << " " << settings->detectedSymmetry.at(symIt)[0] << std::setprecision(5) << " " << settings->detectedSymmetry.at(symIt)[1] << " " << settings->detectedSymmetry.at(symIt)[2] << " " << settings->detectedSymmetry.at(symIt)[3] << " " << settings->detectedSymmetry.at(symIt)[4] << " " << settings->detectedSymmetry.at(symIt)[5] << " " << settings->detectedSymmetry.at(symIt)[6];
3491  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 0, ssHlp.str() );
3492  }
3493 
3494  std::stringstream hlpSS3;
3495  ssHlp.clear(); ssHlp.str ( "" );
3496  hlpSS3 << std::endl << "To facilitate manual checking for symmetries, the following is a list of all detected C symmetries:";
3497  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 0, hlpSS3.str() );
3498 
3499  if ( settings->allDetectedCAxes.size() > 0 )
3500  {
3501  ssHlp.clear(); ssHlp.str ( "" );
3502  ssHlp << " Fold X Y Z Angle Height Average FSC";
3503  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 0, ssHlp.str() );
3504  }
3505  for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( settings->allDetectedCAxes.size() ); symIt++ )
3506  {
3507  ssHlp.clear(); ssHlp.str ( "" );
3508  ssHlp << std::showpos << std::fixed << std::setprecision(0) << " " << settings->allDetectedCAxes.at(symIt)[0] << std::setprecision(5) << " " << settings->allDetectedCAxes.at(symIt)[1] << " " << settings->allDetectedCAxes.at(symIt)[2] << " " << settings->allDetectedCAxes.at(symIt)[3] << " " << settings->allDetectedCAxes.at(symIt)[4] << " " << settings->allDetectedCAxes.at(symIt)[5] << " " << settings->allDetectedCAxes.at(symIt)[6];
3509  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 0, ssHlp.str() );
3510  }
3511 
3512  }
3513 
3514  //================================================ Done
3515  return ;
3516 
3517 }
3518 
3524 {
3525  //================================================ Initialise variables
3526  this->xCom = 0.0;
3527  this->yCom = 0.0;
3528  this->zCom = 0.0;
3529  proshade_double totNonZeroPoints = 0.0;
3530  proshade_signed mapIt = 0;
3531 
3532  //================================================ Compute COM from 0 ; 0 ; 0
3533  for ( proshade_signed xIt = 0; xIt < static_cast<proshade_signed> ( this->xDimIndices ); xIt++ )
3534  {
3535  for ( proshade_signed yIt = 0; yIt < static_cast<proshade_signed> ( this->yDimIndices ); yIt++ )
3536  {
3537  for ( proshade_signed zIt = 0; zIt < static_cast<proshade_signed> ( this->zDimIndices ); zIt++ )
3538  {
3539  //==================================== Find map index
3540  mapIt = zIt + static_cast< proshade_signed > ( this->zDimIndices ) * ( yIt + static_cast< proshade_signed > ( this->yDimIndices ) * xIt );
3541 
3542  //==================================== Use only positive density
3543  if ( this->internalMap[mapIt] <= 0.0 ) { continue; }
3544 
3545  //==================================== Compute Index COM
3546  this->xCom += this->internalMap[mapIt] * static_cast<proshade_double> ( xIt + this->xFrom );
3547  this->yCom += this->internalMap[mapIt] * static_cast<proshade_double> ( yIt + this->yFrom );
3548  this->zCom += this->internalMap[mapIt] * static_cast<proshade_double> ( zIt + this->zFrom );
3549  totNonZeroPoints += this->internalMap[mapIt];
3550  }
3551  }
3552  }
3553 
3554  this->xCom /= totNonZeroPoints;
3555  this->yCom /= totNonZeroPoints;
3556  this->zCom /= totNonZeroPoints;
3557 
3558  //================================================ Convert to real world
3559  this->xCom = static_cast< proshade_double > ( ( static_cast< proshade_single > ( this->xFrom ) * ( this->xDimSizeOriginal / static_cast< proshade_single > ( this->xDimIndicesOriginal ) ) ) +
3560  ( ( static_cast< proshade_single > ( this->xCom ) - static_cast< proshade_single > ( this->xFrom ) ) *
3561  ( static_cast< proshade_single > ( this->xDimSizeOriginal ) / static_cast< proshade_single > ( this->xDimIndicesOriginal ) ) ) );
3562  this->yCom = static_cast< proshade_double > ( ( static_cast< proshade_single > ( this->yFrom ) * ( this->yDimSizeOriginal / static_cast< proshade_single > ( this->yDimIndicesOriginal ) ) ) +
3563  ( ( static_cast< proshade_single > ( this->yCom ) - static_cast< proshade_single > ( this->yFrom ) ) *
3564  ( static_cast< proshade_single > ( this->yDimSizeOriginal ) / static_cast< proshade_single > ( this->yDimIndicesOriginal ) ) ) );
3565  this->zCom = static_cast< proshade_double > ( ( static_cast< proshade_single > ( this->zFrom ) * ( this->zDimSizeOriginal / static_cast< proshade_single > ( this->zDimIndicesOriginal ) ) ) +
3566  ( ( static_cast< proshade_single > ( this->zCom ) - static_cast< proshade_single > ( this->zFrom ) ) *
3567  ( static_cast< proshade_single > ( this->zDimSizeOriginal ) / static_cast< proshade_single > ( this->zDimIndicesOriginal ) ) ) );
3568 
3569  //================================================ Done
3570  return ;
3571 
3572 }
3573 
3579 {
3580  //================================================ Return the value
3581  return ( this->noSpheres );
3582 }
3583 
3589 proshade_double ProSHADE_internal_data::ProSHADE_data::getMapValue ( proshade_unsign pos )
3590 {
3591  //================================================ Return the value
3592  return ( this->internalMap[pos] );
3593 }
3594 
3600 {
3601  //================================================ Return the value
3602  return ( this->maxShellBand );
3603 }
3604 
3609 proshade_double ProSHADE_internal_data::ProSHADE_data::getRRPValue ( proshade_unsign band, proshade_unsign sh1, proshade_unsign sh2 )
3610 {
3611  //================================================ Return the value
3612  return ( this->rrpMatrices[band][sh1][sh2] );
3613 }
3614 
3625 bool ProSHADE_internal_data::ProSHADE_data::shellBandExists ( proshade_unsign shell, proshade_unsign bandVal )
3626 {
3627  if ( this->spheres[shell]->getLocalBandwidth( ) >= bandVal )
3628  {
3629  return ( true );
3630  }
3631  else
3632  {
3633  return ( false );
3634  }
3635 }
3636 
3646 {
3647  //================================================ Report function start
3648  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Centering map onto its COM." );
3649 
3650  //================================================ Copy map for processing
3651  fftw_complex* mapCoeffs = new fftw_complex[this->xDimIndices * this->yDimIndices * this->zDimIndices];
3652  fftw_complex* pattersonMap = new fftw_complex[this->xDimIndices * this->yDimIndices * this->zDimIndices];
3653 
3654  //================================================ Check memory allocation
3655  ProSHADE_internal_misc::checkMemoryAllocation ( mapCoeffs, __FILE__, __LINE__, __func__ );
3656  ProSHADE_internal_misc::checkMemoryAllocation ( pattersonMap, __FILE__, __LINE__, __func__ );
3657 
3658  //================================================ Copy data to mask
3659  for ( proshade_unsign iter = 0; iter < (this->xDimIndices * this->yDimIndices * this->zDimIndices); iter++ )
3660  {
3661  pattersonMap[iter][0] = this->internalMap[iter];
3662  pattersonMap[iter][1] = 0.0;
3663  }
3664 
3665  //================================================ Prepare FFTW plans
3666  fftw_plan forward = fftw_plan_dft_3d ( static_cast< int > ( this->xDimIndices ), static_cast< int > ( this->yDimIndices ), static_cast< int > ( this->zDimIndices ),
3667  pattersonMap, mapCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
3668  fftw_plan inverse = fftw_plan_dft_3d ( static_cast< int > ( this->xDimIndices ), static_cast< int > ( this->yDimIndices ), static_cast< int > ( this->zDimIndices ),
3669  mapCoeffs, pattersonMap, FFTW_BACKWARD, FFTW_ESTIMATE );
3670 
3671  //================================================ Run forward Fourier
3672  fftw_execute ( forward );
3673 
3674  //================================================ Remove the phase
3675  ProSHADE_internal_mapManip::removeMapPhase ( mapCoeffs, this->xDimIndices, this->yDimIndices, this->zDimIndices );
3676 
3677  //================================================ Run inverse Fourier
3678  fftw_execute ( inverse );
3679 
3680  //================================================ Save the results
3681  proshade_signed mapIt, patIt, patX, patY, patZ;
3682  for ( proshade_signed xIt = 0; xIt < static_cast<proshade_signed> ( this->xDimIndices ); xIt++ )
3683  {
3684  for ( proshade_signed yIt = 0; yIt < static_cast<proshade_signed> ( this->yDimIndices ); yIt++ )
3685  {
3686  for ( proshade_signed zIt = 0; zIt < static_cast<proshade_signed> ( this->zDimIndices ); zIt++ )
3687  {
3688  //==================================== Centre patterson map
3689  patX = xIt - ( static_cast<proshade_signed> ( this->xDimIndices ) / 2 ); if ( patX < 0 ) { patX += this->xDimIndices; }
3690  patY = yIt - ( static_cast<proshade_signed> ( this->yDimIndices ) / 2 ); if ( patY < 0 ) { patY += this->yDimIndices; }
3691  patZ = zIt - ( static_cast<proshade_signed> ( this->zDimIndices ) / 2 ); if ( patZ < 0 ) { patZ += this->zDimIndices; }
3692 
3693  //==================================== Find indices
3694  mapIt = zIt + static_cast< proshade_signed > ( this->zDimIndices ) * ( yIt + static_cast< proshade_signed > ( this->yDimIndices ) * xIt );
3695  patIt = patZ + static_cast< proshade_signed > ( this->zDimIndices ) * ( patY + static_cast< proshade_signed > ( this->yDimIndices ) * patX );
3696 
3697  //==================================== Copy
3698  this->internalMap[mapIt] = pattersonMap[patIt][0];
3699  }
3700  }
3701  }
3702 
3703  //================================================ Release memory
3704  delete[] pattersonMap;
3705  delete[] mapCoeffs;
3706 
3707  //================================================ Delete FFTW plans
3708  fftw_destroy_plan ( forward );
3709  fftw_destroy_plan ( inverse );
3710 
3711  //================================================ Report function completion
3712  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Phase information removed." );
3713 
3714  //================================================ Done
3715  return ;
3716 
3717 }
3718 
3723 proshade_double* ProSHADE_internal_data::ProSHADE_data::getRealSphHarmValue ( proshade_unsign band, proshade_unsign order, proshade_unsign shell )
3724 {
3725  //================================================ Done
3726  return ( &this->sphericalHarmonics[shell][seanindex ( static_cast< int > ( order ) - static_cast< int > ( band ),
3727  static_cast< int > ( band ),
3728  static_cast< int > ( this->spheres[shell]->getLocalBandwidth() ) )][0] );
3729 
3730 }
3731 
3736 proshade_double* ProSHADE_internal_data::ProSHADE_data::getImagSphHarmValue ( proshade_unsign band, proshade_unsign order, proshade_unsign shell )
3737 {
3738  //================================================ Done
3739  return ( &this->sphericalHarmonics[shell][seanindex ( static_cast< int > ( order ) - static_cast< int > ( band ),
3740  static_cast< int > ( band ),
3741  static_cast< int > ( this->spheres[shell]->getLocalBandwidth() ) )][1] );
3742 
3743 }
3744 
3749 proshade_double ProSHADE_internal_data::ProSHADE_data::getAnySphereRadius ( proshade_unsign shell )
3750 {
3751  //================================================ Done
3752  return ( this->spheres[shell]->getShellRadius() );
3753 
3754 }
3755 
3761 {
3762  //================================================ Done
3763  return ( this->integrationWeight );
3764 
3765 }
3766 
3772 proshade_unsign ProSHADE_internal_data::ProSHADE_data::getShellBandwidth ( proshade_unsign shell )
3773 {
3774  //================================================ Done
3775  return ( this->spheres[shell]->getLocalBandwidth ( ) );
3776 
3777 }
3778 
3784 proshade_single ProSHADE_internal_data::ProSHADE_data::getSpherePosValue ( proshade_unsign shell )
3785 {
3786  //================================================ Done
3787  return ( this->spherePos.at(shell) );
3788 
3789 }
3790 
3796 proshade_complex** ProSHADE_internal_data::ProSHADE_data::getEMatrixByBand ( proshade_unsign band )
3797 {
3798  //================================================ Done
3799  return ( this->eMatrices[band] );
3800 
3801 }
3802 
3811 void ProSHADE_internal_data::ProSHADE_data::getEMatrixValue ( proshade_unsign band, proshade_unsign order1, proshade_unsign order2, proshade_double* valueReal, proshade_double* valueImag )
3812 {
3813  //================================================ Set pointer
3814  *valueReal = this->eMatrices[band][order1][order2][0];
3815  *valueImag = this->eMatrices[band][order1][order2][1];
3816 
3817  //================================================ Done
3818  return ;
3819 
3820 }
3821 
3827 {
3828  //================================================ Done
3829  return ( this->so3CoeffsInverse );
3830 
3831 }
3832 
3838 {
3839  //================================================ Done
3840  return ( this->so3Coeffs );
3841 
3842 }
3843 
3849 {
3850  //================================================ Done
3851  return ( this->maxCompBand );
3852 
3853 }
3854 
3863 void ProSHADE_internal_data::ProSHADE_data::getWignerMatrixValue ( proshade_unsign band, proshade_unsign order1, proshade_unsign order2, proshade_double* valueReal, proshade_double* valueImag )
3864 {
3865  //================================================ Set pointer
3866  *valueReal = this->wignerMatrices[band][order1][order2][0];
3867  *valueImag = this->wignerMatrices[band][order1][order2][1];
3868 
3869  //================================================ Done
3870  return ;
3871 
3872 }
3873 
3879 {
3880  //================================================ Return the requested value
3881  return ( this->xDimSize );
3882 }
3883 
3889 {
3890  //================================================ Return the requested value
3891  return ( this->yDimSize );
3892 }
3893 
3899 {
3900  //================================================ Return the requested value
3901  return ( this->zDimSize );
3902 }
3903 
3909 {
3910  //================================================ Return the requested value
3911  return ( this->xDimIndices );
3912 }
3913 
3919 {
3920  //================================================ Return the requested value
3921  return ( this->yDimIndices );
3922 }
3923 
3929 {
3930  //================================================ Return the requested value
3931  return ( this->zDimIndices );
3932 }
3933 
3939 {
3940  //================================================ Return the requested value
3941  return ( &this->xFrom );
3942 }
3943 
3949 {
3950  //================================================ Return the requested value
3951  return ( &this->yFrom );
3952 }
3953 
3959 {
3960  //================================================ Return the requested value
3961  return ( &this->zFrom );
3962 }
3963 
3969 {
3970  //================================================ Return the requested value
3971  return ( &this->xTo );
3972 }
3973 
3979 {
3980  //================================================ Return the requested value
3981  return ( &this->yTo );
3982 }
3983 
3989 {
3990  //================================================ Return the requested value
3991  return ( &this->zTo );
3992 }
3993 
3999 {
4000  //================================================ Return the requested value
4001  return ( &this->xAxisOrigin );
4002 }
4003 
4009 {
4010  //================================================ Return the requested value
4011  return ( &this->yAxisOrigin );
4012 }
4013 
4019 {
4020  //================================================ Return the requested value
4021  return ( &this->zAxisOrigin );
4022 }
4023 
4029 {
4030  //================================================ Return the requested value
4031  return ( this->internalMap );
4032 }
4033 
4039 {
4040  //================================================ Return the requested value
4041  return ( this->translationMap );
4042 }
4043 
4049 {
4050  //================================================ Initialise local variables
4051  std::vector< proshade_double > ret;
4052 
4053  //================================================ Save the values
4054  ProSHADE_internal_misc::addToDoubleVector ( &ret, this->mapCOMProcessChangeX );
4055  ProSHADE_internal_misc::addToDoubleVector ( &ret, this->mapCOMProcessChangeY );
4056  ProSHADE_internal_misc::addToDoubleVector ( &ret, this->mapCOMProcessChangeZ );
4057 
4058  //================================================ Return the requested value
4059  return ( ret );
4060 }
4061 
4067 {
4068  //================================================ Mutate
4069  this->integrationWeight = intW;
4070 
4071  //================================================ Done
4072  return ;
4073 
4074 }
4075 
4081 {
4082  //================================================ Mutate
4083  this->integrationWeight += intW;
4084 
4085  //================================================ Done
4086  return ;
4087 
4088 }
4089 
4097 void ProSHADE_internal_data::ProSHADE_data::setEMatrixValue ( proshade_unsign band, proshade_unsign order1, proshade_unsign order2, proshade_complex val )
4098 {
4099  //================================================ Mutate
4100  this->eMatrices[band][order1][order2][0] = val[0];
4101  this->eMatrices[band][order1][order2][1] = val[1];
4102 
4103  //================================================ Done
4104  return ;
4105 
4106 }
4107 
4115 void ProSHADE_internal_data::ProSHADE_data::normaliseEMatrixValue ( proshade_unsign band, proshade_unsign order1, proshade_unsign order2, proshade_double normF )
4116 {
4117  //================================================ Mutate
4118  this->eMatrices[band][order1][order2][0] /= normF;
4119  this->eMatrices[band][order1][order2][1] /= normF;
4120 
4121  //================================================ Done
4122  return ;
4123 
4124 }
4125 
4131 void ProSHADE_internal_data::ProSHADE_data::setSO3CoeffValue ( proshade_unsign position, proshade_complex val )
4132 {
4133  //================================================ Mutate
4134  this->so3Coeffs[position][0] = val[0];
4135  this->so3Coeffs[position][1] = val[1];
4136 
4137  //================================================ Done
4138  return ;
4139 
4140 }
4141 
4149 void ProSHADE_internal_data::ProSHADE_data::setWignerMatrixValue ( proshade_complex val, proshade_unsign band, proshade_unsign order1, proshade_unsign order2 )
4150 {
4151  //================================================ Mutate
4152  this->wignerMatrices[band][order1][order2][0] = val[0];
4153  this->wignerMatrices[band][order1][order2][1] = val[1];
4154 
4155  //================================================ Done
4156  return ;
4157 
4158 }
4159 
4167 void ProSHADE_internal_data::ProSHADE_data::getRealEMatrixValuesForLM ( proshade_signed band, proshade_signed order1, double *eMatsLMReal, int len )
4168 {
4169  //================================================ Save the data into the output array
4170  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4171  {
4172  eMatsLMReal[iter] = static_cast<double> ( this->eMatrices[band][order1][iter][0] );
4173  }
4174 
4175  //================================================ Done
4176  return ;
4177 
4178 }
4179 
4187 void ProSHADE_internal_data::ProSHADE_data::getImagEMatrixValuesForLM ( proshade_signed band, proshade_signed order1, double *eMatsLMImag, int len )
4188 {
4189  //================================================ Save the data into the output array
4190  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4191  {
4192  eMatsLMImag[iter] = static_cast<double> ( this->eMatrices[band][order1][iter][1] );
4193  }
4194 
4195  //================================================ Done
4196  return ;
4197 
4198 }
4199 
4205 void ProSHADE_internal_data::ProSHADE_data::getRealSO3Coeffs ( double *so3CoefsReal, int len )
4206 {
4207  //================================================ Save the data into the output array
4208  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4209  {
4210  so3CoefsReal[iter] = static_cast<double> ( this->so3Coeffs[iter][0] );
4211  }
4212 
4213  //================================================ Done
4214  return ;
4215 
4216 }
4217 
4223 void ProSHADE_internal_data::ProSHADE_data::getImagSO3Coeffs ( double *so3CoefsImag, int len )
4224 {
4225  //================================================ Save the data into the output array
4226  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4227  {
4228  so3CoefsImag[iter] = static_cast<double> ( this->so3Coeffs[iter][1] );
4229  }
4230 
4231  //================================================ Done
4232  return ;
4233 
4234 }
4235 
4245 int ProSHADE_internal_data::ProSHADE_data::so3CoeffsArrayIndex ( proshade_signed order1, proshade_signed order2, proshade_signed band )
4246 {
4247  //================================================ Return the value
4248  return ( static_cast<int> ( so3CoefLoc ( static_cast< int > ( order1 ), static_cast< int > ( order2 ), static_cast< int > ( band ), static_cast< int > ( this->getMaxBand() ) ) ) );
4249 }
4250 
4257 {
4258  //================================================ Save the data into the output array
4259  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4260  {
4261  rotFunReal[iter] = static_cast<double> ( this->so3CoeffsInverse[iter][0] );
4262  }
4263 
4264  //================================================ Done
4265  return ;
4266 
4267 }
4268 
4275 {
4276  //================================================ Save the data into the output array
4277  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4278  {
4279  rotFunImag[iter] = static_cast<double> ( this->so3CoeffsInverse[iter][1] );
4280  }
4281 
4282  //================================================ Done
4283  return ;
4284 
4285 }
4286 
4293 {
4294  //================================================ Save the data into the output array
4295  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4296  {
4297  trsFunReal[iter] = static_cast<double> ( this->translationMap[iter][0] );
4298  }
4299 
4300  //================================================ Done
4301  return ;
4302 
4303 }
4304 
4311 {
4312  //================================================ Save the data into the output array
4313  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4314  {
4315  trsFunImag[iter] = static_cast<double> ( this->translationMap[iter][1] );
4316  }
4317 
4318  //================================================ Done
4319  return ;
4320 
4321 }
4322 
4331 void ProSHADE_internal_data::ProSHADE_data::getRotMatrixFromRotFunInds ( proshade_signed aI, proshade_signed bI, proshade_signed gI, double *rotMat, int len )
4332 {
4333  //================================================ Get Euler angles
4334  proshade_double eA, eB, eG;
4335  ProSHADE_internal_maths::getEulerZXZFromSOFTPosition ( static_cast< int > ( this->getMaxBand() ), aI, bI, gI, &eA, &eB, &eG );
4336 
4337  //================================================ Prepare internal rotation matrix memory
4338  proshade_double* rMat = nullptr;
4339  rMat = new proshade_double[9];
4340  ProSHADE_internal_misc::checkMemoryAllocation ( rMat, __FILE__, __LINE__, __func__ );
4341 
4342  //================================================ Convert to rotation matrix
4344 
4345  //================================================ Copy to output
4346  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4347  {
4348  rotMat[iter] = static_cast<double> ( rMat[iter] );
4349  }
4350 
4351  //================================================ Release internal memory
4352  delete[] rMat;
4353 
4354  //================================================ Done
4355  return ;
4356 
4357 }
4358 
4364 {
4365  //================================================ Return the value
4366  return ( settings->recommendedSymmetryType );
4367 
4368 }
4369 
4375 {
4376  //================================================ Return the value
4377  return ( settings->recommendedSymmetryFold );
4378 
4379 }
4380 
4387 {
4388  //================================================ Return the value
4389  return ( static_cast<proshade_unsign> ( settings->detectedSymmetry.size() ) );
4390 }
4391 
4398 std::vector< std::string > ProSHADE_internal_data::ProSHADE_data::getSymmetryAxis ( ProSHADE_settings* settings, proshade_unsign axisNo )
4399 {
4400  //================================================ Sanity checks
4401  if ( static_cast<proshade_unsign> ( settings->detectedSymmetry.size() ) <= axisNo )
4402  {
4403  ProSHADE_internal_messages::printWarningMessage ( settings->verbose, "!!! ProSHADE WARNING !!! Requested symmetry index does not exist. Returning empty vector.", "WS00039" );
4404  return ( std::vector< std::string > ( ) );
4405  }
4406 
4407  //================================================ Initialise local variables
4408  std::vector< std::string > ret;
4409 
4410  //================================================ Input the axis data as strings
4411  std::stringstream ssHlp;
4412  ssHlp << settings->detectedSymmetry.at(axisNo)[0];
4413  ProSHADE_internal_misc::addToStringVector ( &ret, ssHlp.str() );
4414  ssHlp.str ( "" );
4415 
4416  ssHlp << settings->detectedSymmetry.at(axisNo)[1];
4417  ProSHADE_internal_misc::addToStringVector ( &ret, ssHlp.str() );
4418  ssHlp.str ( "" );
4419 
4420  ssHlp << settings->detectedSymmetry.at(axisNo)[2];
4421  ProSHADE_internal_misc::addToStringVector ( &ret, ssHlp.str() );
4422  ssHlp.str ( "" );
4423 
4424  ssHlp << settings->detectedSymmetry.at(axisNo)[3];
4425  ProSHADE_internal_misc::addToStringVector ( &ret, ssHlp.str() );
4426  ssHlp.str ( "" );
4427 
4428  ssHlp << settings->detectedSymmetry.at(axisNo)[4];
4429  ProSHADE_internal_misc::addToStringVector ( &ret, ssHlp.str() );
4430  ssHlp.str ( "" );
4431 
4432  ssHlp << settings->detectedSymmetry.at(axisNo)[5];
4433  ProSHADE_internal_misc::addToStringVector ( &ret, ssHlp.str() );
4434  ssHlp.str ( "" );
4435 
4436  ssHlp << settings->detectedSymmetry.at(axisNo)[6];
4437  ProSHADE_internal_misc::addToStringVector ( &ret, ssHlp.str() );
4438  ssHlp.str ( "" );
4439 
4440  //================================================ Done
4441  return ( ret );
4442 
4443 }
4444 
4457 void ProSHADE_internal_data::ProSHADE_data::writeOutOverlayFiles ( ProSHADE_settings* settings, proshade_double eulA, proshade_double eulB, proshade_double eulG, std::vector< proshade_double >* rotCentre, std::vector< proshade_double >* ultimateTranslation )
4458 {
4459  //================================================ Write out rotated map
4460  std::stringstream fNameHlp;
4461  fNameHlp << settings->overlayStructureName << ".map";
4462  this->writeMap ( fNameHlp.str() );
4463 
4464  //================================================ Write out rotated co-ordinates if possible
4465  if ( ProSHADE_internal_io::isFilePDB ( this->fileName ) )
4466  {
4467  fNameHlp.str("");
4468  fNameHlp << settings->overlayStructureName << ".pdb";
4469  this->writePdb ( fNameHlp.str(), eulA, eulB, eulG, ultimateTranslation->at(0), ultimateTranslation->at(1), ultimateTranslation->at(2), rotCentre->at(0), rotCentre->at(1), rotCentre->at(2), settings->firstModelOnly );
4470  }
4471 
4472  //================================================ Write out the json file with the results
4473  ProSHADE_internal_io::writeRotationTranslationJSON ( rotCentre->at(0), rotCentre->at(1), rotCentre->at(2),
4474  eulA, eulB, eulG,
4475  ultimateTranslation->at(0), ultimateTranslation->at(1), ultimateTranslation->at(2),
4476  settings->rotTrsJSONFile );
4477 
4478  //================================================ Done
4479  return ;
4480 
4481 }
4482 
4491 void ProSHADE_internal_data::ProSHADE_data::reportOverlayResults ( ProSHADE_settings* settings, std::vector < proshade_double >* rotationCentre, std::vector < proshade_double >* eulerAngles, std::vector < proshade_double >* finalTranslation )
4492 {
4493  //================================================ Empty line
4495 
4496  //================================================ Write out rotation centre translation results
4497  std::stringstream rotCen; rotCen << std::setprecision (3) << std::showpos << "The rotation centre to origin translation vector is: " << -rotationCentre->at(0) << " " << -rotationCentre->at(1) << " " << -rotationCentre->at(2);
4498  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 0, rotCen.str() );
4499 
4500  //================================================ Write out rotation matrix about origin
4501  proshade_double* rotMat = new proshade_double[9];
4502  ProSHADE_internal_misc::checkMemoryAllocation ( rotMat, __FILE__, __LINE__, __func__ );
4503  ProSHADE_internal_maths::getRotationMatrixFromEulerZXZAngles ( eulerAngles->at(0), eulerAngles->at(1), eulerAngles->at(2), rotMat );
4504 
4505  std::stringstream rotMatSS;
4506  rotMatSS << std::setprecision (3) << std::showpos << "The rotation matrix about origin is : " << rotMat[0] << " " << rotMat[1] << " " << rotMat[2] << std::endl;
4507  rotMatSS << std::setprecision (3) << std::showpos << " : " << rotMat[3] << " " << rotMat[4] << " " << rotMat[5] << std::endl;
4508  rotMatSS << std::setprecision (3) << std::showpos << " : " << rotMat[6] << " " << rotMat[7] << " " << rotMat[8];
4509  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 0, rotMatSS.str() );
4510 
4511  delete[] rotMat;
4512 
4513  //================================================ Write out origin to overlay translation results
4514  std::stringstream finTrs; finTrs << std::setprecision (3) << std::showpos << "The rotation centre to overlay translation vector is: " << finalTranslation->at(0) << " " << finalTranslation->at(1) << " " << finalTranslation->at(2);
4515  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 0, finTrs.str() );
4516 
4517  //================================================ Done
4518  return ;
4519 
4520 }
ProSHADE_internal_data::ProSHADE_data::findTopGroupSmooth
proshade_double findTopGroupSmooth(std::vector< proshade_double * > *CSym, size_t peakPos, proshade_double step, proshade_double sigma, proshade_signed windowSize)
This function finds the distinct group of axes with highest peak heights.
Definition: ProSHADE_data.cpp:2301
ProSHADE_internal_data::ProSHADE_data::originalMapYCom
proshade_double originalMapYCom
The COM of the first map to be loaded/computed without any furhter changes being reflacted along the ...
Definition: ProSHADE_data.hpp:92
ProSHADE_settings::maxBandwidth
proshade_unsign maxBandwidth
The bandwidth of spherical harmonics decomposition for the largest sphere.
Definition: ProSHADE_settings.hpp:58
ProSHADE_internal_mapManip::addExtraBoundSpace
void addExtraBoundSpace(proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed *&bounds, proshade_single extraSpace)
This function takes a set of bounds and adds a given number of Angstroms to them.
Definition: ProSHADE_mapManip.cpp:1139
ProSHADE_internal_data::ProSHADE_data::computeSphericalHarmonics
void computeSphericalHarmonics(ProSHADE_settings *settings)
This function computes the spherical harmonics decomposition for the whole structure.
Definition: ProSHADE_data.cpp:1843
ProSHADE_internal_data::ProSHADE_data::mapMovFromsChangeX
proshade_double mapMovFromsChangeX
When the map is translated, the xFrom and xTo values are changed. This variable holds how much they h...
Definition: ProSHADE_data.hpp:94
ProSHADE_settings::recommendedSymmetryType
std::string recommendedSymmetryType
The symmetry type that ProSHADE finds the best fitting for the structure. Possible values are "" for ...
Definition: ProSHADE_settings.hpp:128
ProSHADE_internal_data::ProSHADE_data::getEMatrixByBand
proshade_complex ** getEMatrixByBand(proshade_unsign band)
This function allows access to E matrix for a particular band.
Definition: ProSHADE_data.cpp:3796
ProSHADE_internal_mapManip::determinePDBRanges
void determinePDBRanges(gemmi::Structure pdbFile, proshade_single *xFrom, proshade_single *xTo, proshade_single *yFrom, proshade_single *yTo, proshade_single *zFrom, proshade_single *zTo, bool firstModel)
Function for finding the PDB file ranges.
Definition: ProSHADE_mapManip.cpp:72
ProSHADE_internal_data::ProSHADE_data::fileName
std::string fileName
This is the original file from which the data were obtained.
Definition: ProSHADE_data.hpp:52
ProSHADE_internal_data::ProSHADE_data::zFrom
proshade_signed zFrom
This is the starting index along the z axis.
Definition: ProSHADE_data.hpp:112
ProSHADE_settings::rotTrsJSONFile
std::string rotTrsJSONFile
The filename to which the rotation and translation operations are to be saved into.
Definition: ProSHADE_settings.hpp:139
ProSHADE_internal_data::ProSHADE_data::xDimIndices
proshade_unsign xDimIndices
This is the size of the map cell x dimension in indices.
Definition: ProSHADE_data.hpp:65
ProSHADE_internal_data::ProSHADE_data::getZDimSize
proshade_single getZDimSize(void)
This function allows access to the map size in angstroms along the Z axis.
Definition: ProSHADE_data.cpp:3898
ProSHADE_internal_mapManip::changePDBBFactors
void changePDBBFactors(gemmi::Structure *pdbFile, proshade_double newBFactorValue, bool firstModel)
Function for changing the PDB B-factor values to a specific single value.
Definition: ProSHADE_mapManip.cpp:444
ProSHADE_internal_data::ProSHADE_data::readInMAP
void readInMAP(ProSHADE_settings *settings, proshade_double *maskArr=nullptr, proshade_unsign maskXDim=0, proshade_unsign maskYDim=0, proshade_unsign maskZDim=0, proshade_double *weightsArr=nullptr, proshade_unsign weigXDim=0, proshade_unsign weigYDim=0, proshade_unsign weigZDim=0)
Function for reading map data using gemmi library.
Definition: ProSHADE_data.cpp:619
ProSHADE_internal_data::ProSHADE_data::zDimSizeOriginal
proshade_single zDimSizeOriginal
This is the size of the map cell z dimension in Angstroms.
Definition: ProSHADE_data.hpp:84
ProSHADE_internal_maths::vectorMeanAndSD
void vectorMeanAndSD(std::vector< proshade_double > *vec, proshade_double *&ret)
Function to get vector mean and standard deviation.
Definition: ProSHADE_maths.cpp:121
ProSHADE_internal_data::ProSHADE_data::createNewMapFromBounds
void createNewMapFromBounds(ProSHADE_settings *settings, ProSHADE_data *&newStr, proshade_signed *newBounds)
This function creates a new structure from the calling structure and new bounds values.
Definition: ProSHADE_data.cpp:1387
ProSHADE_internal_io::isFilePDB
bool isFilePDB(std::string fName)
Function determining if the input data type is PDB.
Definition: ProSHADE_io.cpp:32
ProSHADE_internal_data::ProSHADE_data::xGridIndices
proshade_unsign xGridIndices
As far as I know, this is identical to the xDimIndices.
Definition: ProSHADE_data.hpp:68
ProSHADE_internal_data::ProSHADE_data::translationMap
proshade_complex * translationMap
This is where the translation map will be held, if at all used.
Definition: ProSHADE_data.hpp:133
ProSHADE_internal_data::joinElementsFromDifferentGroups
std::vector< std::vector< proshade_double > > joinElementsFromDifferentGroups(std::vector< std::vector< proshade_double > > *first, std::vector< std::vector< proshade_double > > *second, proshade_double matrixTolerance, bool combine)
This function joins two group element lists using only unique elements.
Definition: ProSHADE_data.cpp:3102
ProSHADE_internal_data::ProSHADE_data::xAxisOrder
proshade_unsign xAxisOrder
This is the order of the x axis.
Definition: ProSHADE_data.hpp:71
ProSHADE_internal_data::ProSHADE_data::rotSphericalHarmonics
proshade_complex ** rotSphericalHarmonics
A set of rotated spherical harmonics values arrays for each sphere, used only if map rotation is requ...
Definition: ProSHADE_data.hpp:122
ProSHADE_internal_data::ProSHADE_data::getReBoxBoundaries
void getReBoxBoundaries(ProSHADE_settings *settings, proshade_signed *&ret)
This function finds the boundaries enclosing positive map values and adds some extra space.
Definition: ProSHADE_data.cpp:1329
ProSHADE_settings::determineAllSHValues
void determineAllSHValues(proshade_unsign xDim, proshade_unsign yDim, proshade_single xDimAngs, proshade_single yDimAngs, proshade_single zDimAngs)
This function determines all the required values for spherical harmonics computation.
Definition: ProSHADE.cpp:1615
ProSHADE_internal_data::ProSHADE_data::getComparisonBand
proshade_unsign getComparisonBand(void)
This function allows access to the maximum band for the comparison.
Definition: ProSHADE_data.cpp:3848
ProSHADE_internal_data::ProSHADE_data::getRecommendedSymmetryFold
proshade_unsign getRecommendedSymmetryFold(ProSHADE_settings *settings)
This function simply returns the detected recommended symmetry fold.
Definition: ProSHADE_data.cpp:4374
ProSHADE_internal_data::ProSHADE_data::sphericalHarmonics
proshade_complex ** sphericalHarmonics
A set of spherical harmonics values arrays for each sphere.
Definition: ProSHADE_data.hpp:121
ProSHADE_internal_data::ProSHADE_data::setIntegrationWeight
void setIntegrationWeight(proshade_double intW)
This function allows setting the integration weight for the object.
Definition: ProSHADE_data.cpp:4066
ProSHADE_internal_data::ProSHADE_data::saveRecommendedSymmetry
void saveRecommendedSymmetry(ProSHADE_settings *settings, std::vector< proshade_double * > *CSym, std::vector< proshade_double * > *DSym, std::vector< proshade_double * > *TSym, std::vector< proshade_double * > *OSym, std::vector< proshade_double * > *ISym, std::vector< proshade_double * > *axes, fftw_complex *mapData, fftw_complex *origCoeffs, fftw_complex *fCoeffs, fftw_plan *planForwardFourier, proshade_signed noBins, proshade_signed *binIndexing, proshade_double **bindata, proshade_signed *binCounts)
This function takes all the detected symmetry results and decides on which are to be recommended for ...
Definition: ProSHADE_data.cpp:2596
ProSHADE_internal_data::ProSHADE_data::cAngle
proshade_single cAngle
This is the angle c of the map cell in degrees.
Definition: ProSHADE_data.hpp:64
ProSHADE_internal_data::ProSHADE_data::so3CoeffsInverse
proshade_complex * so3CoeffsInverse
The inverse coefficients obtained by inverse SO(3) Fourier Transform (SOFT) - i.e....
Definition: ProSHADE_data.hpp:130
ProSHADE_internal_data::ProSHADE_data::prepareFSCFourierMemory
void prepareFSCFourierMemory(fftw_complex *&mapData, fftw_complex *&origCoeffs, fftw_complex *&fCoeffs, proshade_signed *&binIndexing, proshade_signed *noBins, proshade_double **&bindata, proshade_signed *&binCounts, fftw_plan *planForwardFourier)
This function allocates the memory and makes all preparations required for FSC computation.
Definition: ProSHADE_data.cpp:2357
ProSHADE_internal_data::ProSHADE_data::getXAxisOrigin
proshade_signed * getXAxisOrigin(void)
This function allows access to the map X axis origin value.
Definition: ProSHADE_data.cpp:3998
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_data::ProSHADE_data::normaliseEMatrixValue
void normaliseEMatrixValue(proshade_unsign band, proshade_unsign order1, proshade_unsign order2, proshade_double normF)
This function allows normalising the E matrix value.
Definition: ProSHADE_data.cpp:4115
ProSHADE_internal_data::ProSHADE_data::yCom
proshade_double yCom
The COM of the map after processing along the Y-axis.
Definition: ProSHADE_data.hpp:78
ProSHADE_settings::boundsExtraSpace
proshade_single boundsExtraSpace
The number of extra angstroms to be added to all re-boxing bounds just for safety.
Definition: ProSHADE_settings.hpp:93
ProSHADE_exception
This class is the representation of ProSHADE exception.
Definition: ProSHADE_exceptions.hpp:37
ProSHADE_internal_data::ProSHADE_data::getImagTranslationFunction
void getImagTranslationFunction(double *trsFunImag, int len)
This function fills the input array with the imaginary translation function values.
Definition: ProSHADE_data.cpp:4310
ProSHADE_internal_data::ProSHADE_data::setWignerMatrixValue
void setWignerMatrixValue(proshade_complex val, proshade_unsign band, proshade_unsign order1, proshade_unsign order2)
This function allows setting the Wigner D matrix value by its band, order1 and order2 co-ordinate.
Definition: ProSHADE_data.cpp:4149
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_settings::blurFactor
proshade_single blurFactor
This is the amount by which B-factors should be increased to create the blurred map for masking.
Definition: ProSHADE_settings.hpp:78
ProSHADE_internal_data::ProSHADE_data::saveDetectedSymmetries
void saveDetectedSymmetries(ProSHADE_settings *settings, std::vector< proshade_double * > *CSyms, std::vector< std::vector< proshade_double > > *allCs)
This function takes the results of point group searches and saves then into the output variables.
Definition: ProSHADE_data.cpp:2255
ProSHADE_settings::maskMap
bool maskMap
Should the map be masked from noise?
Definition: ProSHADE_settings.hpp:80
ProSHADE_settings::requestedSymmetryFold
proshade_unsign requestedSymmetryFold
The fold of the requested symmetry (only applicable to C and D symmetry types).
Definition: ProSHADE_settings.hpp:131
ProSHADE_settings::requestedSymmetryType
std::string requestedSymmetryType
The symmetry type requested by the user. Allowed values are C, D, T, O and I.
Definition: ProSHADE_settings.hpp:130
ProSHADE_internal_data::ProSHADE_data::getEMatrixValue
void getEMatrixValue(proshade_unsign band, proshade_unsign order1, proshade_unsign order2, proshade_double *valueReal, proshade_double *valueImag)
This function allows access to E matrix by knowing the band, order1 and order2 indices.
Definition: ProSHADE_data.cpp:3811
ProSHADE_internal_data::ProSHADE_data::getXDimSize
proshade_single getXDimSize(void)
This function allows access to the map size in angstroms along the X axis.
Definition: ProSHADE_data.cpp:3878
ProSHADE_internal_data::ProSHADE_data::getInternalMap
proshade_double *& getInternalMap(void)
This function allows access to the first map array value address.
Definition: ProSHADE_data.cpp:4028
ProSHADE_internal_data::ProSHADE_data
This class contains all inputed and derived data for a single structure.
Definition: ProSHADE_data.hpp:49
ProSHADE_internal_data::ProSHADE_data::noSpheres
proshade_unsign noSpheres
The number of spheres with map projected onto them.
Definition: ProSHADE_data.hpp:119
ProSHADE_settings::removeNegativeDensity
bool removeNegativeDensity
Should the negative density be removed from input files?
Definition: ProSHADE_settings.hpp:47
ProSHADE_internal_data::ProSHADE_data::getYDimSize
proshade_single getYDimSize(void)
This function allows access to the map size in angstroms along the Y axis.
Definition: ProSHADE_data.cpp:3888
ProSHADE_internal_mapManip::getMaskFromBlurr
void getMaskFromBlurr(proshade_double *&blurMap, proshade_double *&outMap, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single noIQRs)
Function for computing mask from blurred map.
Definition: ProSHADE_mapManip.cpp:1032
ProSHADE_internal_data::ProSHADE_data::computeFSC
proshade_double computeFSC(ProSHADE_settings *settings, std::vector< proshade_double * > *CSym, size_t symIndex, fftw_complex *mapData, fftw_complex *fCoeffs, fftw_complex *origCoeffs, fftw_plan *planForwardFourier, proshade_signed noBins, proshade_signed *binIndexing, proshade_double **&bindata, proshade_signed *&binCounts)
This function computes FSC for any given axis in the supplied CSym symmetry axes vector.
Definition: ProSHADE_data.cpp:2425
ProSHADE_internal_data::ProSHADE_data::zDimSize
proshade_single zDimSize
This is the size of the map cell z dimension in Angstroms.
Definition: ProSHADE_data.hpp:61
ProSHADE_internal_data::ProSHADE_data::so3Coeffs
proshade_complex * so3Coeffs
The coefficients obtained by SO(3) Fourier Transform (SOFT), in this case derived from the E matrices...
Definition: ProSHADE_data.hpp:129
ProSHADE_internal_data::ProSHADE_data::getMaxSpheres
proshade_unsign getMaxSpheres(void)
This function returns the number of spheres which contain the whole object.
Definition: ProSHADE_data.cpp:3578
ProSHADE_internal_data::ProSHADE_data::setSO3CoeffValue
void setSO3CoeffValue(proshade_unsign position, proshade_complex val)
This function allows setting the SOFT coefficient values using array position and value.
Definition: ProSHADE_data.cpp:4131
ProSHADE_settings::saveMask
bool saveMask
Should the mask be saved?
Definition: ProSHADE_settings.hpp:84
ProSHADE_settings::requestedResolution
proshade_single requestedResolution
The resolution to which the calculations are to be done.
Definition: ProSHADE_settings.hpp:50
ProSHADE_internal_data::ProSHADE_data::mapCOMProcessChangeX
proshade_double mapCOMProcessChangeX
The change in X axis between the creation of the structure (originalMapXCom) and just before rotation...
Definition: ProSHADE_data.hpp:97
ProSHADE_internal_mapManip::findMAPCOMValues
void findMAPCOMValues(proshade_double *map, proshade_double *xCom, proshade_double *yCom, proshade_double *zCom, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed xFrom, proshade_signed xTo, proshade_signed yFrom, proshade_signed yTo, proshade_signed zFrom, proshade_signed zTo)
This function finds the Centre of Mass for a map.
Definition: ProSHADE_mapManip.cpp:240
ProSHADE_internal_mapManip::blurSharpenMap
void blurSharpenMap(proshade_double *&map, proshade_double *&maskedMap, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single blurringFactor)
Function for blurring/sharpening maps.
Definition: ProSHADE_mapManip.cpp:932
ProSHADE_internal_data::ProSHADE_data::originalMapXCom
proshade_double originalMapXCom
The COM of the first map to be loaded/computed without any furhter changes being reflacted along the ...
Definition: ProSHADE_data.hpp:91
ProSHADE_internal_data::ProSHADE_data::getRotMatrixFromRotFunInds
void getRotMatrixFromRotFunInds(proshade_signed aI, proshade_signed bI, proshade_signed gI, double *rotMat, int len)
This function takes rotation function indices, converts them to Euler angles and these to rotation ma...
Definition: ProSHADE_data.cpp:4331
ProSHADE_settings::setDetectedSymmetry
void setDetectedSymmetry(proshade_double *sym)
Sets the final detected symmetry axes information.
Definition: ProSHADE.cpp:1315
ProSHADE_internal_data::ProSHADE_data::xDimIndicesOriginal
proshade_unsign xDimIndicesOriginal
This is the size of the map cell x dimension in indices.
Definition: ProSHADE_data.hpp:85
ProSHADE_internal_mapManip::moveMapByFourier
void moveMapByFourier(proshade_double *&map, proshade_single xMov, proshade_single yMov, proshade_single zMov, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim)
Function for moving map back to original PDB location by using Fourier transformation.
Definition: ProSHADE_mapManip.cpp:811
ProSHADE_internal_data::ProSHADE_data::getAnySphereRadius
proshade_double getAnySphereRadius(proshade_unsign shell)
This function allows access to the radius of any particular sphere.
Definition: ProSHADE_data.cpp:3749
ProSHADE_settings::progressiveSphereMapping
bool progressiveSphereMapping
If true, each shell will have its own angular resolution dependent on the actual number of map points...
Definition: ProSHADE_settings.hpp:105
ProSHADE_internal_mapManip::movePDBForMapCalc
void movePDBForMapCalc(gemmi::Structure *pdbFile, proshade_single xMov, proshade_single yMov, proshade_single zMov, bool firstModel)
Function for moving co-ordinate atoms to better suit theoretical map computation.
Definition: ProSHADE_mapManip.cpp:573
ProSHADE_internal_data::ProSHADE_data::getSO3Coeffs
proshade_complex * getSO3Coeffs(void)
This function allows access to the SO(3) coefficients array.
Definition: ProSHADE_data.cpp:3837
ProSHADE_settings::maxSphereDists
proshade_single maxSphereDists
The distance between spheres in spherical mapping for the largest sphere.
Definition: ProSHADE_settings.hpp:65
ProSHADE_internal_data::ProSHADE_data::maxShellBand
proshade_unsign maxShellBand
The maximum band for any shell of the object.
Definition: ProSHADE_data.hpp:123
ProSHADE_internal_data::ProSHADE_data::originalPdbRotCenY
proshade_double originalPdbRotCenY
The centre of rotation as it relates to the original PDB positions (and not the ProSHADE internal map...
Definition: ProSHADE_data.hpp:103
ProSHADE_internal_data::ProSHADE_data::getSpherePositions
void getSpherePositions(ProSHADE_settings *settings)
This function determines the sphere positions (radii) for sphere mapping.
Definition: ProSHADE_data.cpp:1742
ProSHADE_internal_data::ProSHADE_data::invertMirrorMap
void invertMirrorMap(ProSHADE_settings *settings)
Function for inverting the map to its mirror image.
Definition: ProSHADE_data.cpp:1188
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
checkElementsFormGroup
bool checkElementsFormGroup(std::vector< std::vector< proshade_double > > *elements, proshade_double matrixTolerance)
This function checks if all group element products produce another group element.
Definition: ProSHADE_data.cpp:3061
ProSHADE_internal_data::ProSHADE_data::getMapCOMProcessChange
std::vector< proshade_double > getMapCOMProcessChange(void)
This function allows access to the translation caused by structure processing.
Definition: ProSHADE_data.cpp:4048
ProSHADE_internal_data::ProSHADE_data::saveRequestedSymmetryD
void saveRequestedSymmetryD(ProSHADE_settings *settings, std::vector< proshade_double * > *DSym, std::vector< proshade_double * > *axes, fftw_complex *mapData, fftw_complex *origCoeffs, fftw_complex *fCoeffs, fftw_plan *planForwardFourier, proshade_signed noBins, proshade_signed *binIndexing, proshade_double **bindata, proshade_signed *binCounts)
This function takes the D symmetries and searched for the requested symmetry.
Definition: ProSHADE_data.cpp:2886
ProSHADE_internal_data::ProSHADE_data::mapCOMProcessChangeZ
proshade_double mapCOMProcessChangeZ
The change in Z axis between the creation of the structure (originalMapZCom) and just before rotation...
Definition: ProSHADE_data.hpp:99
ProSHADE_settings::maskingThresholdIQRs
proshade_single maskingThresholdIQRs
Number of inter-quartile ranges from the median to be used for thresholding the blurred map for maski...
Definition: ProSHADE_settings.hpp:79
ProSHADE_settings::usePhase
bool usePhase
If true, the full data will be used, if false, Patterson maps will be used instead and phased data wi...
Definition: ProSHADE_settings.hpp:62
ProSHADE_internal_data::ProSHADE_data::xDimSizeOriginal
proshade_single xDimSizeOriginal
This is the size of the map cell x dimension in Angstroms.
Definition: ProSHADE_data.hpp:82
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_data::ProSHADE_data::getIntegrationWeight
proshade_double getIntegrationWeight(void)
This function allows access to the integration weight for the object.
Definition: ProSHADE_data.cpp:3760
ProSHADE_internal_io::figureDataType
InputType figureDataType(std::string fName)
Function determining input data type.
Definition: ProSHADE_io.cpp:918
ProSHADE_internal_data::ProSHADE_data::setIntegrationWeightCumul
void setIntegrationWeightCumul(proshade_double intW)
This function allows setting the cumulative integration weight for the object.
Definition: ProSHADE_data.cpp:4080
ProSHADE_internal_mapManip::reSampleMapToResolutionTrilinear
void reSampleMapToResolutionTrilinear(proshade_double *&map, proshade_single resolution, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single *&corrs)
This function re-samples a map to conform to given resolution using tri-linear interpolation.
Definition: ProSHADE_mapManip.cpp:1175
ProSHADE_internal_data::ProSHADE_data::zAxisOriginOriginal
proshade_signed zAxisOriginOriginal
This is the origin position along the z axis.
Definition: ProSHADE_data.hpp:90
ProSHADE_internal_spheres::ProSHADE_sphere
This class contains all inputed and derived data for a single sphere.
Definition: ProSHADE_spheres.hpp:49
ProSHADE_internal_data::ProSHADE_data::getSymmetryAxis
std::vector< std::string > getSymmetryAxis(ProSHADE_settings *settings, proshade_unsign axisNo)
This function returns a single symmetry axis as a vector of strings from the recommended symmetry axe...
Definition: ProSHADE_data.cpp:4398
ProSHADE_settings::maskFileName
std::string maskFileName
The filename to which mask should be saved.
Definition: ProSHADE_settings.hpp:85
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_internal_data::ProSHADE_data::eMatrices
proshade_complex *** eMatrices
The trace sigma and full rotation function c*conj(c) integral tables.
Definition: ProSHADE_data.hpp:127
ProSHADE_internal_data::ProSHADE_data::centreMapOnCOM
void centreMapOnCOM(ProSHADE_settings *settings)
This function shits the map so that its COM is in the centre of the map.
Definition: ProSHADE_data.cpp:1525
ProSHADE_internal_data::ProSHADE_data::getXFromPtr
proshade_signed * getXFromPtr(void)
This function allows access to the map start along the X axis.
Definition: ProSHADE_data.cpp:3938
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_data::ProSHADE_data::yDimIndices
proshade_unsign yDimIndices
This is the size of the map cell y dimension in indices.
Definition: ProSHADE_data.hpp:66
ProSHADE_internal_data::ProSHADE_data::xAxisOriginOriginal
proshade_signed xAxisOriginOriginal
This is the origin position along the x axis.
Definition: ProSHADE_data.hpp:88
ProSHADE_internal_data::ProSHADE_data::deepCopyMap
void deepCopyMap(proshade_double *&saveTo, proshade_signed verbose)
This function copies the internal map into the supplied pointer, which it also allocates.
Definition: ProSHADE_data.cpp:3435
ProSHADE_internal_data::ProSHADE_data::getWignerMatrixValue
void getWignerMatrixValue(proshade_unsign band, proshade_unsign order1, proshade_unsign order2, proshade_double *valueReal, proshade_double *valueImag)
This function allows access to the Wigner D matrix by knowing the band, order1 and order2 indices.
Definition: ProSHADE_data.cpp:3863
ProSHADE_internal_data::ProSHADE_data::getShellBandwidth
proshade_unsign getShellBandwidth(proshade_unsign shell)
This function allows access to the bandwidth of a particular shell.
Definition: ProSHADE_data.cpp:3772
ProSHADE_internal_data::ProSHADE_data::yGridIndices
proshade_unsign yGridIndices
As far as I know, this is identical to the yDimIndices.
Definition: ProSHADE_data.hpp:69
ProSHADE_settings::changeMapResolutionTriLinear
bool changeMapResolutionTriLinear
Should maps be re-sampled to obtain the required resolution?
Definition: ProSHADE_settings.hpp:52
ProSHADE_internal_data::ProSHADE_data::zAxisOrigin
proshade_signed zAxisOrigin
This is the origin position along the z axis.
Definition: ProSHADE_data.hpp:76
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_data::ProSHADE_data::getNoRecommendedSymmetryAxes
proshade_unsign getNoRecommendedSymmetryAxes(ProSHADE_settings *settings)
This function returns the number of detected recommended symmetry axes.
Definition: ProSHADE_data.cpp:4386
ProSHADE_internal_data::ProSHADE_data::setPDBMapValues
void setPDBMapValues(void)
Function for determining iterator start and stop positions.
Definition: ProSHADE_data.cpp:937
ProSHADE_internal_data::ProSHADE_data::findMapCOM
void findMapCOM(void)
This function finds the centre of mass of the internal map representation.
Definition: ProSHADE_data.cpp:3523
ProSHADE_internal_data::ProSHADE_data::getRealRotFunction
void getRealRotFunction(double *rotFunReal, int len)
This function fills the input array with the real rotation function values.
Definition: ProSHADE_data.cpp:4256
ProSHADE_internal_data::ProSHADE_data::yDimIndicesOriginal
proshade_unsign yDimIndicesOriginal
This is the size of the map cell y dimension in indices.
Definition: ProSHADE_data.hpp:86
ProSHADE_internal_data::ProSHADE_data::getXDim
proshade_unsign getXDim(void)
This function allows access to the map size in indices along the X axis.
Definition: ProSHADE_data.cpp:3908
ProSHADE_internal_mapManip::removeWaters
void removeWaters(gemmi::Structure *pdbFile, bool firstModel)
This function removed all waters from PDB input file.
Definition: ProSHADE_mapManip.cpp:504
ProSHADE_internal_data::ProSHADE_data::originalPdbTransX
proshade_double originalPdbTransX
The optimal translation vector as it relates to the original PDB positions (and not the ProSHADE inte...
Definition: ProSHADE_data.hpp:105
ProSHADE_data.hpp
This is the header file containing internal data representation and manipulation structures and funct...
ProSHADE_internal_data::ProSHADE_data::zTo
proshade_signed zTo
This is the final index along the z axis.
Definition: ProSHADE_data.hpp:115
ProSHADE_internal_data::ProSHADE_data::readInStructure
void readInStructure(std::string fName, proshade_unsign inputO, ProSHADE_settings *settings, proshade_double *maskArr=nullptr, proshade_unsign maskXDim=0, proshade_unsign maskYDim=0, proshade_unsign maskZDim=0, proshade_double *weightsArr=nullptr, proshade_unsign weigXDim=0, proshade_unsign weigYDim=0, proshade_unsign weigZDim=0)
This function initialises the basic ProSHADE_data variables and reads in a single structure.
Definition: ProSHADE_data.cpp:509
ProSHADE_internal_data::ProSHADE_data::xTo
proshade_signed xTo
This is the final index along the x axis.
Definition: ProSHADE_data.hpp:113
ProSHADE_internal_data::ProSHADE_data::writeMask
void writeMask(std::string fName, proshade_double *mask)
Function for writing out a mask in MRC MAP format.
Definition: ProSHADE_data.cpp:1151
ProSHADE_internal_data::ProSHADE_data::integrationWeight
proshade_double integrationWeight
The Pearson's c.c. type weighting for the integration.
Definition: ProSHADE_data.hpp:128
ProSHADE_internal_data::ProSHADE_data::saveRequestedSymmetryC
void saveRequestedSymmetryC(ProSHADE_settings *settings, std::vector< proshade_double * > *CSym, std::vector< proshade_double * > *axes)
This function takes the C symmetries and searched for the requested symmetry.
Definition: ProSHADE_data.cpp:2834
ProSHADE_settings::peakThresholdMin
proshade_double peakThresholdMin
The threshold for peak height above which axes are considered possible.
Definition: ProSHADE_settings.hpp:135
ProSHADE_internal_data::ProSHADE_data::mapToSpheres
void mapToSpheres(ProSHADE_settings *settings)
This function converts the internal map onto a set of concentric spheres.
Definition: ProSHADE_data.cpp:1799
ProSHADE_internal_peakSearch::findPeaks1D
std::vector< proshade_signed > findPeaks1D(std::vector< proshade_double > data)
This function simply finds all the peaks in a 1D data array.
Definition: ProSHADE_peakSearch.cpp:1014
ProSHADE_internal_data::ProSHADE_data::xAxisOrigin
proshade_signed xAxisOrigin
This is the origin position along the x axis.
Definition: ProSHADE_data.hpp:74
ProSHADE_internal_data::ProSHADE_data::originalPdbTransY
proshade_double originalPdbTransY
The optimal translation vector as it relates to the original PDB positions (and not the ProSHADE inte...
Definition: ProSHADE_data.hpp:106
ProSHADE_internal_maths::computeFSC
proshade_double computeFSC(fftw_complex *fCoeffs1, fftw_complex *fCoeffs2, proshade_unsign xInds, proshade_unsign yInds, proshade_unsign zInds, proshade_signed noBins, proshade_signed *binIndexing, proshade_double **&binData, proshade_signed *&binCounts)
This function computes the FSC.
Definition: ProSHADE_maths.cpp:3086
ProSHADE_internal_data::ProSHADE_data::rrpMatrices
proshade_double *** rrpMatrices
The energy levels descriptor shell correlation tables.
Definition: ProSHADE_data.hpp:126
ProSHADE_internal_data::ProSHADE_data::yTo
proshade_signed yTo
This is the final index along the y axis.
Definition: ProSHADE_data.hpp:114
ProSHADE_settings::forceP1
bool forceP1
Should the P1 spacegroup be forced on the input PDB files?
Definition: ProSHADE_settings.hpp:44
ProSHADE_internal_data::ProSHADE_data::xDimSize
proshade_single xDimSize
This is the size of the map cell x dimension in Angstroms.
Definition: ProSHADE_data.hpp:59
ProSHADE_internal_mapManip::findPDBCOMValues
void findPDBCOMValues(gemmi::Structure pdbFile, proshade_double *xCom, proshade_double *yCom, proshade_double *zCom, bool firstModel)
This function finds the Centre of Mass for the co-ordinate file.
Definition: ProSHADE_mapManip.cpp:157
ProSHADE_internal_data::ProSHADE_data::mapMovFromsChangeY
proshade_double mapMovFromsChangeY
When the map is translated, the yFrom and yTo values are changed. This variable holds how much they h...
Definition: ProSHADE_data.hpp:95
ProSHADE_internal_data::ProSHADE_data::yAxisOrigin
proshade_signed yAxisOrigin
This is the origin position along the y axis.
Definition: ProSHADE_data.hpp:75
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_data::ProSHADE_data::yAxisOrder
proshade_unsign yAxisOrder
This is the order of the y axis.
Definition: ProSHADE_data.hpp:72
ProSHADE_internal_data::ProSHADE_data::getRealSphHarmValue
proshade_double * getRealSphHarmValue(proshade_unsign band, proshade_unsign order, proshade_unsign shell)
This function allows access to the private internal real spherical harmonics values.
Definition: ProSHADE_data.cpp:3723
ProSHADE_internal_data::ProSHADE_data::internalMap
proshade_double * internalMap
The internal map data representation, which may be amended as the run progresses.
Definition: ProSHADE_data.hpp:56
ProSHADE_internal_data::ProSHADE_data::writePdb
void writePdb(std::string fName, proshade_double euA=0.0, proshade_double euB=0.0, proshade_double euG=0.0, proshade_double trsX=0.0, proshade_double trsY=0.0, proshade_double trsZ=0.0, proshade_double rotX=0.0, proshade_double rotY=0.0, proshade_double rotZ=0.0, bool firstModel=true)
This function writes out the co-ordinates file with ProSHADE type rotation and translation applied.
Definition: ProSHADE_data.cpp:1072
ProSHADE_internal_data::ProSHADE_data::zDimIndicesOriginal
proshade_unsign zDimIndicesOriginal
This is the size of the map cell z dimension in indices.
Definition: ProSHADE_data.hpp:87
ProSHADE_internal_data::ProSHADE_data::reportOverlayResults
void reportOverlayResults(ProSHADE_settings *settings, std::vector< proshade_double > *rotationCentre, std::vector< proshade_double > *eulerAngles, std::vector< proshade_double > *finalTranslation)
This function reports the results of the overlay mode.
Definition: ProSHADE_data.cpp:4491
ProSHADE_internal_symmetry
This namespace contains the symmetry detection related code.
Definition: ProSHADE_data.cpp:75
ProSHADE_internal_data::ProSHADE_data::xFrom
proshade_signed xFrom
This is the starting index along the x axis.
Definition: ProSHADE_data.hpp:110
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::writeMap
void writeMap(std::string fName, std::string title="Created by ProSHADE and written by GEMMI", int mode=2)
Function for writing out the internal structure representation in MRC MAP format.
Definition: ProSHADE_data.cpp:1005
ProSHADE_internal_data::ProSHADE_data::shellBandExists
bool shellBandExists(proshade_unsign shell, proshade_unsign bandVal)
This function checks if particular shell has a particular band.
Definition: ProSHADE_data.cpp:3625
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_settings::appliedMaskFileName
std::string appliedMaskFileName
The filename from which mask data will be read from.
Definition: ProSHADE_settings.hpp:86
ProSHADE_settings::normaliseMap
bool normaliseMap
Should the map be normalised to mean 0 sd 1?
Definition: ProSHADE_settings.hpp:72
ProSHADE_internal_data::ProSHADE_data::writeGemmi
void writeGemmi(std::string fName, gemmi::Structure gemmiStruct, proshade_double euA=0.0, proshade_double euB=0.0, proshade_double euG=0.0, proshade_double trsX=0.0, proshade_double trsY=0.0, proshade_double trsZ=0.0, proshade_double rotX=0.0, proshade_double rotY=0.0, proshade_double rotZ=0.0, bool firstModel=true)
This function writes out the gemmi::Structure object with ProSHADE type rotation and translation appl...
Definition: ProSHADE_data.cpp:1109
ProSHADE_settings::addExtraSpace
proshade_single addExtraSpace
If this value is non-zero, this many angstroms of empty space will be added to the internal map.
Definition: ProSHADE_settings.hpp:102
ProSHADE_internal_data::ProSHADE_data::getYDim
proshade_unsign getYDim(void)
This function allows access to the map size in indices along the Y axis.
Definition: ProSHADE_data.cpp:3918
ProSHADE_internal_data::ProSHADE_data::maxCompBand
proshade_unsign maxCompBand
The largest comparison band - this variable tells how large arrays will be allocated for the comparis...
Definition: ProSHADE_data.hpp:132
ProSHADE_internal_data::ProSHADE_data::getImagEMatrixValuesForLM
void getImagEMatrixValuesForLM(proshade_signed band, proshade_signed order1, double *eMatsLMImag, int len)
This function fills the input array with the imaginary E matrix values for particular band and order1...
Definition: ProSHADE_data.cpp:4187
ProSHADE_internal_data::ProSHADE_data::maskMap
void maskMap(ProSHADE_settings *settings)
Function for computing the map mask using blurring and X IQRs from median.
Definition: ProSHADE_data.cpp:1289
ProSHADE_internal_data::ProSHADE_data::~ProSHADE_data
~ProSHADE_data(void)
Destructor for the ProSHADE_data class.
Definition: ProSHADE_data.cpp:349
ProSHADE_internal_mapManip::getNonZeroBounds
void getNonZeroBounds(proshade_double *map, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, proshade_signed *&ret)
Function for finding the map boundaries enclosing positive only values.
Definition: ProSHADE_mapManip.cpp:1082
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::xCom
proshade_double xCom
The COM of the map after processing along the X-axis.
Definition: ProSHADE_data.hpp:77
ProSHADE_internal_data::ProSHADE_data::getRealTranslationFunction
void getRealTranslationFunction(double *trsFunReal, int len)
This function fills the input array with the real translation function values.
Definition: ProSHADE_data.cpp:4292
ProSHADE_settings::fscThreshold
proshade_double fscThreshold
The threshold for FSC value under which the axis is considered to be likely noise.
Definition: ProSHADE_settings.hpp:134
ProSHADE_internal_data::ProSHADE_data::readInGemmi
void readInGemmi(gemmi::Structure gemmiStruct, ProSHADE_settings *settings)
Function for reading co-ordinate data from Gemmi object.
Definition: ProSHADE_data.cpp:785
ProSHADE_internal_maths::binReciprocalSpaceReflections
void binReciprocalSpaceReflections(proshade_unsign xInds, proshade_unsign yInds, proshade_unsign zInds, proshade_signed *noBin, proshade_signed *&binIndexing)
This function does binning of the reciprocal space reflections.
Definition: ProSHADE_maths.cpp:2939
ProSHADE_internal_data::ProSHADE_data::getImagSphHarmValue
proshade_double * getImagSphHarmValue(proshade_unsign band, proshade_unsign order, proshade_unsign shell)
This function allows access to the private internal imaginary spherical harmonics values.
Definition: ProSHADE_data.cpp:3736
ProSHADE_internal_data::ProSHADE_data::originalPdbTransZ
proshade_double originalPdbTransZ
The optimal translation vector as it relates to the original PDB positions (and not the ProSHADE inte...
Definition: ProSHADE_data.hpp:107
ProSHADE_internal_data::ProSHADE_data::writeOutOverlayFiles
void writeOutOverlayFiles(ProSHADE_settings *settings, proshade_double eulA, proshade_double eulB, proshade_double eulG, std::vector< proshade_double > *rotCentre, std::vector< proshade_double > *ultimateTranslation)
This function writes out the rotated map, co-ordinates and transformation JSON file.
Definition: ProSHADE_data.cpp:4457
ProSHADE_internal_io::applyMask
void applyMask(proshade_double *&map, std::string maskFile, proshade_unsign xDimInds, proshade_unsign yDimInds, proshade_unsign zDimInds, proshade_signed verbose, proshade_double *maskArray=nullptr, proshade_unsign maXInds=0, proshade_unsign maYInds=0, proshade_unsign maZInds=0)
This function reads and applies the mask to the map.
Definition: ProSHADE_io.cpp:294
ProSHADE_internal_data::ProSHADE_data::getImagRotFunction
void getImagRotFunction(double *rotFunImag, int len)
This function fills the input array with the imaginary rotation function values.
Definition: ProSHADE_data.cpp:4274
ProSHADE_internal_data::ProSHADE_data::processInternalMap
void processInternalMap(ProSHADE_settings *settings)
This function simply clusters several other functions which should be called together.
Definition: ProSHADE_data.cpp:1699
ProSHADE_internal_io::readInMapData
void readInMapData(gemmi::Ccp4< float > *gemmiMap, proshade_double *&map, proshade_unsign xDimInds, proshade_unsign yDimInds, proshade_unsign zDimInds, proshade_unsign xAxOrder, proshade_unsign yAxOrder, proshade_unsign zAxOrder)
This function converts the gemmi Ccp4 object data to ProSHADE internal map representation.
Definition: ProSHADE_io.cpp:178
ProSHADE_internal_data::ProSHADE_data::getTranslationFnPointer
proshade_complex * getTranslationFnPointer(void)
This function allows access to the translation function through a pointer.
Definition: ProSHADE_data.cpp:4038
ProSHADE_settings::setRecommendedSymmetry
void setRecommendedSymmetry(std::string val)
Sets the ProSHADE detected symmetry type.
Definition: ProSHADE.cpp:1231
ProSHADE_internal_data::ProSHADE_data::zCom
proshade_double zCom
The COM of the map after processing along the Z-axis.
Definition: ProSHADE_data.hpp:79
ProSHADE_internal_data::ProSHADE_data::spherePos
std::vector< proshade_single > spherePos
Vector of sphere radii from the centre of the map.
Definition: ProSHADE_data.hpp:118
ProSHADE_settings::invertMap
bool invertMap
Should the map be inverted? Only use this if you think you have the wrong hand in your map.
Definition: ProSHADE_settings.hpp:75
ProSHADE_internal_maths::smoothen1D
std::vector< proshade_double > smoothen1D(proshade_double step, proshade_signed windowSize, proshade_double sigma, std::vector< proshade_double > data)
This function takes a 1D vector and computes smoothened version based on the parameters.
Definition: ProSHADE_maths.cpp:2869
ProSHADE_internal_data::ProSHADE_data::getZDim
proshade_unsign getZDim(void)
This function allows access to the map size in indices along the Z axis.
Definition: ProSHADE_data.cpp:3928
ProSHADE_internal_mapManip::removeMapPhase
void removeMapPhase(fftw_complex *&mapCoeffs, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim)
This function removes the phase from reciprocal (frequency) map.
Definition: ProSHADE_mapManip.cpp:1631
ProSHADE_settings::boundsSimilarityThreshold
proshade_signed boundsSimilarityThreshold
Number of indices which can be added just to make sure same size in indices is achieved.
Definition: ProSHADE_settings.hpp:94
ProSHADE_internal_data::ProSHADE_data::getSpherePosValue
proshade_single getSpherePosValue(proshade_unsign shell)
This function allows access to sphere positions.
Definition: ProSHADE_data.cpp:3784
ProSHADE_settings::changeMapResolution
bool changeMapResolution
Should maps be re-sampled to obtain the required resolution?
Definition: ProSHADE_settings.hpp:51
ProSHADE_internal_data::ProSHADE_data::getRRPValue
proshade_double getRRPValue(proshade_unsign band, proshade_unsign sh1, proshade_unsign sh2)
This function allows access to the priva internal RRP matrices.
Definition: ProSHADE_data.cpp:3609
ProSHADE_internal_io::writeRotationTranslationJSON
void writeRotationTranslationJSON(proshade_double trsX1, proshade_double trsY1, proshade_double trsZ1, proshade_double eulA, proshade_double eulB, proshade_double eulG, proshade_double trsX2, proshade_double trsY2, proshade_double trsZ2, std::string fileName)
Function for writing out the optimal rotation and translation into a JSON file.
Definition: ProSHADE_io.cpp:960
ProSHADE_internal_data::ProSHADE_data::mapMovFromsChangeZ
proshade_double mapMovFromsChangeZ
When the map is translated, the zFrom and zTo values are changed. This variable holds how much they h...
Definition: ProSHADE_data.hpp:96
axesToGroupTypeSanityCheck
void axesToGroupTypeSanityCheck(proshade_unsign requiredAxes, proshade_unsign obtainedAxes, std::string groupType)
This function checks that the required and obtained numbers of axes are correct, printing error if th...
Definition: ProSHADE_data.cpp:3013
ProSHADE_settings::overlayStructureName
std::string overlayStructureName
The filename to which the rotated and translated moving structure is to be saved.
Definition: ProSHADE_settings.hpp:138
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_data::ProSHADE_data::getMapValue
proshade_double getMapValue(proshade_unsign pos)
This function returns the internal map representation value of a particular array position.
Definition: ProSHADE_data.cpp:3589
ProSHADE_internal_data::ProSHADE_data::getImagSO3Coeffs
void getImagSO3Coeffs(double *so3CoefsImag, int len)
This function fills the input array with the imaginary SO(3) coefficient values.
Definition: ProSHADE_data.cpp:4223
ProSHADE_internal_data::ProSHADE_data::getYToPtr
proshade_signed * getYToPtr(void)
This function allows access to the map last position along the Y axis.
Definition: ProSHADE_data.cpp:3978
ProSHADE_settings::moveToCOM
bool moveToCOM
Logical value stating whether the structure should be moved to have its Centre Of Mass (COM) in the m...
Definition: ProSHADE_settings.hpp:99
ProSHADE_internal_data::ProSHADE_data::wignerMatrices
proshade_complex *** wignerMatrices
These matrices are computed for a particular rotation to be done in spherical harmonics.
Definition: ProSHADE_data.hpp:131
ProSHADE_internal_data::ProSHADE_data::fileType
ProSHADE_internal_io::InputType fileType
This is the type of the input file.
Definition: ProSHADE_data.hpp:53
ProSHADE_internal_data::ProSHADE_data::getZAxisOrigin
proshade_signed * getZAxisOrigin(void)
This function allows access to the map Z axis origin value.
Definition: ProSHADE_data.cpp:4018
ProSHADE_internal_io::writeOutMapHeader
void writeOutMapHeader(gemmi::Ccp4< float > *map, proshade_unsign xDimInds, proshade_unsign yDimInds, proshade_unsign zDimInds, proshade_single xDim, proshade_single yDim, proshade_single zDim, proshade_single aAng, proshade_single bAng, proshade_single cAng, proshade_signed xFrom, proshade_signed yFrom, proshade_signed zFrom, proshade_signed xAxOrigin, proshade_signed yAxOrigin, proshade_signed zAxOrigin, proshade_unsign xAxOrder, proshade_unsign yAxOrder, proshade_unsign zAxOrder, proshade_unsign xGridInds, proshade_unsign yGridInds, proshade_unsign zGridInds, std::string title, int mode)
This function parses the CCP4 MAP file header as read in by gemmi.
Definition: ProSHADE_io.cpp:868
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_data::ProSHADE_data::yDimSizeOriginal
proshade_single yDimSizeOriginal
This is the size of the map cell y dimension in Angstroms.
Definition: ProSHADE_data.hpp:83
ProSHADE_internal_data::ProSHADE_data::zAxisOrder
proshade_unsign zAxisOrder
This is the order of the z axis.
Definition: ProSHADE_data.hpp:73
ProSHADE_internal_mapManip::beautifyBoundaries
void beautifyBoundaries(proshade_signed *&bounds, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_signed boundsDiffThres)
Function for modifying boundaries to a mathematically more pleasant values.
Definition: ProSHADE_mapManip.cpp:1916
ProSHADE_internal_data::ProSHADE_data::reSampleMap
void reSampleMap(ProSHADE_settings *settings)
This function changes the internal map sampling to conform to particular resolution value.
Definition: ProSHADE_data.cpp:1453
ProSHADE_internal_mapManip::rotatePDBCoordinates
void rotatePDBCoordinates(gemmi::Structure *pdbFile, proshade_double euA, proshade_double euB, proshade_double euG, proshade_double xCom, proshade_double yCom, proshade_double zCom, bool firstModel)
Function for rotating the PDB file co-ordinates by Euler angles.
Definition: ProSHADE_mapManip.cpp:296
ProSHADE_internal_mapManip::generateMapFromPDB
void generateMapFromPDB(gemmi::Structure pdbFile, proshade_double *&map, proshade_single requestedResolution, proshade_single xCell, proshade_single yCell, proshade_single zCell, proshade_signed *xTo, proshade_signed *yTo, proshade_signed *zTo, bool forceP1, bool firstModel)
This function generates a theoretical map from co-ordinate input files.
Definition: ProSHADE_mapManip.cpp:644
ProSHADE_internal_maths::rotationMatrixSimilarity
bool rotationMatrixSimilarity(std::vector< proshade_double > *mat1, std::vector< proshade_double > *mat2, proshade_double tolerance=0.1)
This function compares the distance between two rotation matrices and decides if they are similar usi...
Definition: ProSHADE_maths.cpp:2295
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_data::ProSHADE_data::getXToPtr
proshade_signed * getXToPtr(void)
This function allows access to the map last position along the X axis.
Definition: ProSHADE_data.cpp:3968
ProSHADE_internal_data::ProSHADE_data::getYFromPtr
proshade_signed * getYFromPtr(void)
This function allows access to the map start along the Y axis.
Definition: ProSHADE_data.cpp:3948
ProSHADE_internal_mapManip::translatePDBCoordinates
void translatePDBCoordinates(gemmi::Structure *pdbFile, proshade_double transX, proshade_double transY, proshade_double transZ, bool firstModel)
Function for translating the PDB file co-ordinates by given distances in Angstroms.
Definition: ProSHADE_mapManip.cpp:381
ProSHADE_internal_data::ProSHADE_data::mapCOMProcessChangeY
proshade_double mapCOMProcessChangeY
The change in Y axis between the creation of the structure (originalMapYCom) and just before rotation...
Definition: ProSHADE_data.hpp:98
ProSHADE_internal_data::ProSHADE_data::bAngle
proshade_single bAngle
This is the angle b of the map cell in degrees.
Definition: ProSHADE_data.hpp:63
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_data::ProSHADE_data::yAxisOriginOriginal
proshade_signed yAxisOriginOriginal
This is the origin position along the y axis.
Definition: ProSHADE_data.hpp:89
ProSHADE_internal_data::ProSHADE_data::aAngle
proshade_single aAngle
This is the angle a of the map cell in degrees.
Definition: ProSHADE_data.hpp:62
ProSHADE_internal_misc::addToDoubleVectorVector
void addToDoubleVectorVector(std::vector< std::vector< proshade_double > > *vecToAddTo, std::vector< proshade_double > elementToAdd)
Adds the element to the vector of vectors.
Definition: ProSHADE_misc.cpp:233
ProSHADE_internal_data::ProSHADE_data::reportSymmetryResults
void reportSymmetryResults(ProSHADE_settings *settings)
This function takes prints the report for symmetry detection.
Definition: ProSHADE_data.cpp:3467
ProSHADE_internal_data::ProSHADE_data::figureIndexStartStop
void figureIndexStartStop(void)
Function for determining iterator start and stop positions.
Definition: ProSHADE_data.cpp:983
ProSHADE_internal_data::ProSHADE_data::isEmpty
bool isEmpty
This variable stated whether the class contains any information.
Definition: ProSHADE_data.hpp:139
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_settings::forceBounds
proshade_signed * forceBounds
These will be the boundaries to be forced upon the map.
Definition: ProSHADE_settings.hpp:96
ProSHADE_settings::fourierWeightsFileName
std::string fourierWeightsFileName
The filename from which Fourier weights data will be read from.
Definition: ProSHADE_settings.hpp:89
ProSHADE_settings::detectedSymmetry
std::vector< proshade_double * > detectedSymmetry
The vector of detected symmetry axes.
Definition: ProSHADE_settings.hpp:146
ProSHADE_internal_data::ProSHADE_data::yFrom
proshade_signed yFrom
This is the starting index along the y axis.
Definition: ProSHADE_data.hpp:111
ProSHADE_settings::recommendedSymmetryFold
proshade_unsign recommendedSymmetryFold
The fold of the recommended symmetry C or D type, 0 otherwise.
Definition: ProSHADE_settings.hpp:129
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_mapManip::reSampleMapToResolutionFourier
void reSampleMapToResolutionFourier(proshade_double *&map, proshade_single resolution, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single *&corrs)
This function re-samples a map to conform to given resolution using Fourier.
Definition: ProSHADE_mapManip.cpp:1376
ProSHADE_internal_mapManip::copyMapByBounds
void copyMapByBounds(proshade_signed xFrom, proshade_signed xTo, proshade_signed yFrom, proshade_signed yTo, proshade_signed zFrom, proshade_signed zTo, proshade_signed origXFrom, proshade_signed origYFrom, proshade_signed origZFrom, proshade_unsign yDimIndices, proshade_unsign zDimIndices, proshade_unsign origXDimIndices, proshade_unsign origYDimIndices, proshade_unsign origZDimIndices, proshade_double *&newMap, proshade_double *origMap)
This function copies an old map to a new map with different boundaries.
Definition: ProSHADE_mapManip.cpp:2082
ProSHADE_internal_data::ProSHADE_data::getYAxisOrigin
proshade_signed * getYAxisOrigin(void)
This function allows access to the map Y axis origin value.
Definition: ProSHADE_data.cpp:4008
ProSHADE_internal_data::ProSHADE_data::originalPdbRotCenZ
proshade_double originalPdbRotCenZ
The centre of rotation as it relates to the original PDB positions (and not the ProSHADE internal map...
Definition: ProSHADE_data.hpp:104
ProSHADE_internal_misc::addToSingleVector
void addToSingleVector(std::vector< proshade_single > *vecToAddTo, proshade_single elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:55
ProSHADE_internal_data::ProSHADE_data::ProSHADE_data
ProSHADE_data()
Constructor for getting empty ProSHADE_data class.
Definition: ProSHADE_data.cpp:93
ProSHADE_settings::firstModelOnly
bool firstModelOnly
Shoud only the first PDB model be used, or should all models be used?
Definition: ProSHADE_settings.hpp:46
sortProSHADESymmetryByFSC
bool sortProSHADESymmetryByFSC(proshade_double *a, proshade_double *b)
This function allows using std::sort to sort vectors of ProSHADE symmetry format.
Definition: ProSHADE_data.cpp:1883
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_settings::useSameBounds
bool useSameBounds
Switch to say that the same boundaries as used for the first should be used for all input maps.
Definition: ProSHADE_settings.hpp:95
ProSHADE_internal_data::ProSHADE_data::zDimIndices
proshade_unsign zDimIndices
This is the size of the map cell z dimension in indices.
Definition: ProSHADE_data.hpp:67
ProSHADE_internal_data::computeGroupElementsForGroup
std::vector< std::vector< proshade_double > > computeGroupElementsForGroup(proshade_double xAx, proshade_double yAx, proshade_double zAx, proshade_signed fold)
This function computes the group elements as rotation matrices (except for the identity element) for ...
Definition: ProSHADE_data.cpp:2949
ProSHADE_settings::pdbBFactorNewVal
proshade_double pdbBFactorNewVal
Change all PDB B-factors to this value (for smooth maps).
Definition: ProSHADE_settings.hpp:55
ProSHADE_internal_data::ProSHADE_data::getRecommendedSymmetryType
std::string getRecommendedSymmetryType(ProSHADE_settings *settings)
This function simply returns the detected recommended symmetry type.
Definition: ProSHADE_data.cpp:4363
ProSHADE_internal_data::ProSHADE_data::zGridIndices
proshade_unsign zGridIndices
As far as I know, this is identical to the zDimIndices.
Definition: ProSHADE_data.hpp:70
ProSHADE_settings::setRecommendedFold
void setRecommendedFold(proshade_unsign val)
Sets the ProSHADE detected symmetry fold.
Definition: ProSHADE.cpp:1254
ProSHADE_internal_io::applyWeights
void applyWeights(proshade_double *&map, std::string weightsFile, proshade_unsign xDimInds, proshade_unsign yDimInds, proshade_unsign zDimInds, proshade_signed verbose, proshade_double *weightsArray=nullptr, proshade_unsign waXInds=0, proshade_unsign waYInds=0, proshade_unsign waZInds=0)
This function reads and applies the Fourier weights to the map.
Definition: ProSHADE_io.cpp:563
ProSHADE_settings::setVariablesLeftOnAuto
void setVariablesLeftOnAuto(void)
Function to determine general values that the user left on auto-determination.
Definition: ProSHADE.cpp:337
ProSHADE_internal_data::ProSHADE_data::setEMatrixValue
void setEMatrixValue(proshade_unsign band, proshade_unsign order1, proshade_unsign order2, proshade_complex val)
This function allows setting the E matrix value.
Definition: ProSHADE_data.cpp:4097
ProSHADE_settings::removeWaters
bool removeWaters
Should all waters be removed from input PDB files?
Definition: ProSHADE_settings.hpp:45
ProSHADE_internal_data::ProSHADE_data::spheres
ProSHADE_internal_spheres::ProSHADE_sphere ** spheres
The set of concentric spheres to which the intermal density map has been projected.
Definition: ProSHADE_data.hpp:120
ProSHADE_internal_data::ProSHADE_data::normaliseMap
void normaliseMap(ProSHADE_settings *settings)
Function for normalising the map values to mean 0 and sd 1..
Definition: ProSHADE_data.cpp:1242
ProSHADE_internal_data::ProSHADE_data::detectSymmetryFromAngleAxisSpace
void detectSymmetryFromAngleAxisSpace(ProSHADE_settings *settings, std::vector< proshade_double * > *axes, std::vector< std::vector< proshade_double > > *allCs)
This function runs the symmetry detection algorithms on this structure using the angle-axis space and...
Definition: ProSHADE_data.cpp:1903
ProSHADE_internal_data::ProSHADE_data::originalPdbRotCenX
proshade_double originalPdbRotCenX
The centre of rotation as it relates to the original PDB positions (and not the ProSHADE internal map...
Definition: ProSHADE_data.hpp:102
ProSHADE_internal_data::ProSHADE_data::originalMapZCom
proshade_double originalMapZCom
The COM of the first map to be loaded/computed without any furhter changes being reflacted along the ...
Definition: ProSHADE_data.hpp:93
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_data::ProSHADE_data::getRealSO3Coeffs
void getRealSO3Coeffs(double *so3CoefsReal, int len)
This function fills the input array with the real SO(3) coefficient values.
Definition: ProSHADE_data.cpp:4205
ProSHADE_internal_data::ProSHADE_data::yDimSize
proshade_single yDimSize
This is the size of the map cell y dimension in Angstroms.
Definition: ProSHADE_data.hpp:60
ProSHADE_settings::setResolution
void setResolution(proshade_single resolution)
Sets the requested resolution in the appropriate variable.
Definition: ProSHADE.cpp:383
ProSHADE_settings::allDetectedCAxes
std::vector< std::vector< proshade_double > > allDetectedCAxes
The vector of all detected cyclic symmetry axes.
Definition: ProSHADE_settings.hpp:147
ProSHADE_internal_data::ProSHADE_data::removePhaseInormation
void removePhaseInormation(ProSHADE_settings *settings)
This function removes phase from the map, effectively converting it to Patterson map.
Definition: ProSHADE_data.cpp:3645
ProSHADE_internal_mapManip::moveMapByIndices
void moveMapByIndices(proshade_single *xMov, proshade_single *yMov, proshade_single *zMov, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed *xFrom, proshade_signed *xTo, proshade_signed *yFrom, proshade_signed *yTo, proshade_signed *zFrom, proshade_signed *zTo, proshade_signed *xOrigin, proshade_signed *yOrigin, proshade_signed *zOrigin)
Function for moving map back to original PDB location by changing the indices.
Definition: ProSHADE_mapManip.cpp:763
ProSHADE_internal_io::readInMapHeader
void readInMapHeader(gemmi::Ccp4< float > *map, proshade_unsign *xDimInds, proshade_unsign *yDimInds, proshade_unsign *zDimInds, proshade_single *xDim, proshade_single *yDim, proshade_single *zDim, proshade_single *aAng, proshade_single *bAng, proshade_single *cAng, proshade_signed *xFrom, proshade_signed *yFrom, proshade_signed *zFrom, proshade_signed *xAxOrigin, proshade_signed *yAxOrigin, proshade_signed *zAxOrigin, proshade_unsign *xAxOrder, proshade_unsign *yAxOrder, proshade_unsign *zAxOrder, proshade_unsign *xGridInds, proshade_unsign *yGridInds, proshade_unsign *zGridInds)
This function parses the CCP4 MAP file header as read in by gemmi.
Definition: ProSHADE_io.cpp:109
ProSHADE_internal_maths::multiplyGroupElementMatrices
std::vector< proshade_double > multiplyGroupElementMatrices(std::vector< proshade_double > *el1, std::vector< proshade_double > *el2)
This function computes matrix multiplication using the ProSHADE group element matrix format as input ...
Definition: ProSHADE_maths.cpp:2241
ProSHADE_internal_data::ProSHADE_data::getAllGroupElements
std::vector< std::vector< proshade_double > > getAllGroupElements(ProSHADE_settings *settings, std::vector< proshade_unsign > axesList, std::string groupType="", proshade_double matrixTolerance=0.05)
This function returns the group elements as rotation matrices of any defined point group.
Definition: ProSHADE_data.cpp:3168
ProSHADE_internal_data::ProSHADE_data::addExtraSpace
void addExtraSpace(ProSHADE_settings *settings)
This function increases the size of the map so that it can add empty space around it.
Definition: ProSHADE_data.cpp:1590
ProSHADE_internal_data::ProSHADE_data::getRealEMatrixValuesForLM
void getRealEMatrixValuesForLM(proshade_signed band, proshade_signed order1, double *eMatsLMReal, int len)
This function fills the input array with the real E matrix values for particular band and order1 (l a...
Definition: ProSHADE_data.cpp:4167
ProSHADE_internal_misc::addToStringVector
void addToStringVector(std::vector< std::string > *vecToAddTo, std::string elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:33
ProSHADE_internal_sphericalHarmonics::computeSphericalHarmonics
void computeSphericalHarmonics(proshade_unsign band, proshade_double *sphereMappedData, proshade_complex *&shArray)
This function computes the spherical harmonics of a aingle shell, saving them in supplied pointer.
Definition: ProSHADE_sphericalHarmonics.cpp:394
checkElementAlreadyExists
bool checkElementAlreadyExists(std::vector< std::vector< proshade_double > > *elements, std::vector< proshade_double > *elem, proshade_double matrixTolerance)
This function checks if the element list already contains a given matrix.
Definition: ProSHADE_data.cpp:3035
ProSHADE_internal_data::ProSHADE_data::readInPDB
void readInPDB(ProSHADE_settings *settings)
Function for reading pdb data.
Definition: ProSHADE_data.cpp:750
ProSHADE_internal_data::ProSHADE_data::getZFromPtr
proshade_signed * getZFromPtr(void)
This function allows access to the map start along the Z axis.
Definition: ProSHADE_data.cpp:3958
ProSHADE_internal_data::ProSHADE_data::getZToPtr
proshade_signed * getZToPtr(void)
This function allows access to the map last position along the Z axis.
Definition: ProSHADE_data.cpp:3988
ProSHADE_internal_data::ProSHADE_data::so3CoeffsArrayIndex
int so3CoeffsArrayIndex(proshade_signed order1, proshade_signed order2, proshade_signed band)
This function gets the SO(3) coefficients array index for a particular so(3) band,...
Definition: ProSHADE_data.cpp:4245