ProSHADE  0.7.5.4 (MAR 2021)
Protein Shape Detection
ProSHADE_sphericalHarmonics.cpp
Go to the documentation of this file.
1 
22 //==================================================== ProSHADE
24 
40 void ProSHADE_internal_sphericalHarmonics::allocateComputationMemory ( proshade_unsign band, proshade_double*& inputReal, proshade_double*& inputImag, proshade_double*& outputReal, proshade_double*& outputImag, proshade_double*& shWeights, double**& tableSpace, proshade_double*& tableSpaceHelper, fftw_complex*& workspace )
41 {
42  //================================================ Initialise local variables
43  proshade_unsign oneDimmension = 2 * band;
44 
45  //================================================ Allocate Input Memory
46  inputReal = new proshade_double [oneDimmension * oneDimmension];
47  inputImag = new proshade_double [oneDimmension * oneDimmension];
48 
49  //================================================ Allocate Output Memory
50  outputReal = new proshade_double [oneDimmension * oneDimmension];
51  outputImag = new proshade_double [oneDimmension * oneDimmension];
52 
53  //================================================ Allocate Working Memory
54  shWeights = new proshade_double [band * 4];
55  workspace = new fftw_complex [( 8 * band * band ) + ( 10 * band )];
56 
57  //================================================ Allocate table
58  tableSpaceHelper = new proshade_double [static_cast<proshade_unsign> ( Reduced_Naive_TableSize ( band, band ) +
59  Reduced_SpharmonicTableSize ( band, band ) )];
60 
61  //================================================ Check memory allocation success
62  ProSHADE_internal_misc::checkMemoryAllocation ( inputReal, __FILE__, __LINE__, __func__ );
63  ProSHADE_internal_misc::checkMemoryAllocation ( inputImag, __FILE__, __LINE__, __func__ );
64  ProSHADE_internal_misc::checkMemoryAllocation ( outputReal, __FILE__, __LINE__, __func__ );
65  ProSHADE_internal_misc::checkMemoryAllocation ( outputImag, __FILE__, __LINE__, __func__ );
66  ProSHADE_internal_misc::checkMemoryAllocation ( shWeights, __FILE__, __LINE__, __func__ );
67  ProSHADE_internal_misc::checkMemoryAllocation ( tableSpaceHelper, __FILE__, __LINE__, __func__ );
68  ProSHADE_internal_misc::checkMemoryAllocation ( workspace, __FILE__, __LINE__, __func__ );
69 
70  //================================================ Done
71  return ;
72 
73 }
74 
87 void ProSHADE_internal_sphericalHarmonics::placeWithinWorkspacePointers ( fftw_complex*& workspace, proshade_unsign oDim, proshade_double*& rres, proshade_double*& ires, proshade_double*& fltres, proshade_double*& scratchpad )
88 {
89  //================================================ Place pointers as required by SOFT2.0
90  rres = reinterpret_cast<proshade_double*> ( workspace );
91  ires = rres + ( oDim * oDim );
92  fltres = ires + ( oDim * oDim );
93  scratchpad = fltres + ( oDim / 2 );
94 
95  //================================================ Done
96  return ;
97 
98 }
99 
114 void ProSHADE_internal_sphericalHarmonics::initialiseFFTWPlans ( proshade_unsign band, fftw_plan& fftPlan, fftw_plan& dctPlan, proshade_double*& inputReal, proshade_double*& inputImag, proshade_double*& rres, proshade_double*& ires, proshade_double*& scratchpad )
115 {
116  //================================================ Initialize fft plan along phi angles
117  fftw_iodim dims[1];
118  fftw_iodim howmany_dims[1];
119 
120  int rank = 1;
121  int howmany_rank = 1;
122 
123  dims[0].n = static_cast<int> ( band * 2 );
124  dims[0].is = 1;
125  dims[0].os = static_cast<int> ( band * 2 );
126 
127  howmany_dims[0].n = static_cast<int> ( band * 2 );
128  howmany_dims[0].is = static_cast<int> ( band * 2 );
129  howmany_dims[0].os = 1;
130 
131  //================================================ Plan fft transform
132  fftPlan = fftw_plan_guru_split_dft ( rank,
133  dims,
134  howmany_rank,
135  howmany_dims,
136  inputReal,
137  inputImag,
138  rres,
139  ires,
140  FFTW_ESTIMATE );
141 
142  //================================================ Initialize dct plan for SHT
143  dctPlan = fftw_plan_r2r_1d ( static_cast<int> ( band * 2 ),
144  scratchpad,
145  scratchpad + static_cast<int> ( band * 2 ),
146  FFTW_REDFT10,
147  FFTW_ESTIMATE ) ;
148 
149  //================================================ Done
150  return ;
151 
152 }
153 
154 
171 void ProSHADE_internal_sphericalHarmonics::releaseSphericalMemory ( proshade_double*& inputReal, proshade_double*& inputImag, proshade_double*& outputReal, proshade_double*& outputImag, double*& tableSpaceHelper, double**& tableSpace, double*& shWeights, fftw_complex*& workspace, fftw_plan& fftPlan, fftw_plan& dctPlan, proshade_unsign band )
172 {
173  //================================================ Release all memory related to SH
174  delete[] inputReal;
175  delete[] inputImag;
176  delete[] outputReal;
177  delete[] outputImag;
178  delete[] tableSpaceHelper;
179  delete[] shWeights;
180  delete[] workspace;
181 
182  //================================================ Set pointers to NULL
183  tableSpaceHelper = NULL;
184  tableSpace = NULL;
185  shWeights = NULL;
186  workspace = NULL;
187 
188  //================================================ Delete fftw plans
189  fftw_destroy_plan ( dctPlan );
190  fftw_destroy_plan ( fftPlan );
191 
192  //================================================ Done
193  return ;
194 
195 }
196 
218 void ProSHADE_internal_sphericalHarmonics::initialiseAllMemory ( proshade_unsign band, proshade_double*& inputReal, proshade_double*& inputImag, proshade_double*& outputReal, proshade_double*& outputImag, double*& shWeights, double**& tableSpace, double*& tableSpaceHelper, fftw_complex*& workspace, proshade_double*& rres, proshade_double*& ires, proshade_double*& fltres, proshade_double*& scratchpad, fftw_plan& fftPlan, fftw_plan& dctPlan )
219 {
220  //================================================ Initialise local variables
221  proshade_unsign oneDim = band * 2;
222 
223  //================================================ Allocate memory for local pointers
224  allocateComputationMemory ( band, inputReal, inputImag, outputReal, outputImag, shWeights, tableSpace, tableSpaceHelper, workspace );
225 
226  //================================================ Within workspace pointers
227  placeWithinWorkspacePointers ( workspace, oneDim, rres, ires, fltres, scratchpad );
228 
229  //================================================ Generate Seminaive and naive tables for Legendre Polynomials
230  tableSpace = SemiNaive_Naive_Pml_Table ( band, band, tableSpaceHelper, reinterpret_cast<double*> ( workspace ) );
231 
232  //================================================ Make weights for spherical transform
233  makeweights ( band, shWeights );
234 
235  //================================================ Initialize FFTW Plans
236  initialiseFFTWPlans ( band, fftPlan, dctPlan, inputReal, inputImag, rres, ires, scratchpad );
237 
238  //================================================ Done
239  return ;
240 
241 }
242 
258 void ProSHADE_internal_sphericalHarmonics::initialSplitDiscreteTransform ( proshade_unsign oneDim, proshade_double*& inputReal, proshade_double*& inputImag, proshade_double*& rres, proshade_double*& ires, proshade_double* mappedData, fftw_plan& fftPlan, proshade_double normCoeff )
259 {
260  //================================================ Load mapped data to decomposition array
261  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( oneDim * oneDim ); iter++ )
262  {
263  inputReal[iter] = mappedData[iter];
264  inputImag[iter] = 0.0;
265  }
266 
267  //================================================ Execute fft plan along phi
268  fftw_execute_split_dft ( fftPlan, inputReal, inputImag, rres, ires ) ;
269 
270  //================================================ Normalize
271  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( oneDim * oneDim ); iter++ )
272  {
273  rres[iter] *= normCoeff;
274  ires[iter] *= normCoeff;
275  }
276 
277  //================================================ Done
278  return ;
279 
280 }
281 
301 void ProSHADE_internal_sphericalHarmonics::computeSphericalTransformCoeffs ( proshade_unsign band, proshade_double*& rdataptr, proshade_double*& idataptr, proshade_double*& outputReal, proshade_double*& outputImag, proshade_double*& rres, proshade_double*& ires, proshade_double*& fltres, proshade_double*& scratchpad, double**& tablePml, double*& shWeights, fftw_plan& dctPlan )
302 {
303  //================================================ Calculate the coefficients for each band
304  rdataptr = outputReal;
305  idataptr = outputImag;
306  for ( proshade_unsign bandIter = 0; bandIter < band; bandIter++ )
307  {
308  //============================================ Real part calculation
309  SemiNaiveReduced ( rres + ( bandIter * ( band * 2 ) ),
310  band,
311  bandIter,
312  fltres,
313  scratchpad,
314  tablePml[bandIter],
315  shWeights,
316  &dctPlan);
317 
318  //============================================ Save the real results to temporary holder
319  memcpy ( rdataptr, fltres, sizeof(proshade_double) * ( band - bandIter ) );
320  rdataptr += band - bandIter;
321 
322  //============================================ Imaginary part calculation
323  SemiNaiveReduced ( ires + ( bandIter * ( band * 2 ) ),
324  band,
325  bandIter,
326  fltres,
327  scratchpad,
328  tablePml[bandIter],
329  shWeights,
330  &dctPlan );
331 
332  //============================================ Save the imaginary results
333  memcpy ( idataptr, fltres, sizeof(proshade_double) * ( band - bandIter ) );
334  idataptr += band - bandIter;
335  }
336 
337  //================================================ DONE
338  return ;
339 
340 }
341 
353 void ProSHADE_internal_sphericalHarmonics::applyCondonShortleyPhase ( proshade_unsign band, proshade_double* outputReal, proshade_double* outputImag, proshade_complex*& shArray )
354 {
355  //================================================ Copy the results into the final holder
356  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( (band * 2) * (band * 2) ); iter++ )
357  {
358  shArray[iter][0] = outputReal[iter];
359  shArray[iter][1] = outputImag[iter];
360  }
361 
362  //================================================ Apply the Condon-Shortley phase sign
363  proshade_double powerOne = 1.0;
364  proshade_unsign hlp1 = 0;
365  proshade_unsign hlp2 = 0;
366  for ( proshade_signed order = 1; order < static_cast<proshade_signed> ( band ); order++)
367  {
368  powerOne *= -1.0;
369  for ( proshade_signed bandIter = order; bandIter < static_cast<proshade_signed> ( band ); bandIter++)
370  {
371  hlp1 = seanindex ( order, bandIter, band );
372  hlp2 = seanindex ( -order, bandIter, band );
373 
374  shArray[hlp2][0] = powerOne * static_cast<proshade_double> ( outputReal[hlp1] );
375  shArray[hlp2][1] = -powerOne * static_cast<proshade_double> ( outputImag[hlp1] );
376  }
377  }
378 
379  //================================================ DONE
380  return ;
381 
382 }
383 
395 void ProSHADE_internal_sphericalHarmonics::computeSphericalHarmonics ( proshade_unsign band, proshade_double* sphereMappedData, proshade_complex*& shArray )
396 {
397  //================================================ Initialise local variables
398  proshade_double *inputReal = NULL, *inputImag = NULL, *outputReal = NULL, *outputImag = NULL;
399  double *shWeights = NULL, *tableSpaceHelper = NULL;
400  double** tablePml = NULL;
401  fftw_complex* workspace = NULL;
402  proshade_unsign oneDim = static_cast<proshade_unsign> ( band * 2 );
403  proshade_double normCoeff = ( 1.0 / ( static_cast<proshade_double> ( band * 2 ) ) ) * sqrt( 2.0 * M_PI );
404 
405  //================================================ Set output to zeroes (so that all unfilled data are not random)
406  for ( proshade_unsign i = 0; i < ( 2 * band * 2 * band); i++ )
407  {
408  shArray[i][0] = 0.0;
409  shArray[i][1] = 0.0;
410  }
411 
412  //================================================ Within workspace pointers
413  proshade_double *rres = NULL, *ires = NULL, *fltres = NULL, *scratchpad = NULL, *rdataptr = NULL, *idataptr = NULL;
414 
415  //================================================ FFTW Plans
416  fftw_plan fftPlan = NULL;
417  fftw_plan dctPlan = NULL;
418 
419  //================================================ Initialise all memory
420  initialiseAllMemory ( band, inputReal, inputImag, outputReal, outputImag, shWeights, tablePml, tableSpaceHelper, workspace,
421  rres, ires, fltres, scratchpad, fftPlan, dctPlan );
422 
423  //================================================ Do the initial discrete split transform
424  initialSplitDiscreteTransform ( oneDim, inputReal, inputImag, rres, ires, sphereMappedData, fftPlan, normCoeff );
425 
426  //================================================ Complete the spherical harmonics transform
427  computeSphericalTransformCoeffs ( band, rdataptr, idataptr, outputReal, outputImag, rres, ires, fltres, scratchpad, tablePml, shWeights, dctPlan );
428 
429  //================================================ Apply the Condon-Shortley phase and save result to the final array
430  applyCondonShortleyPhase ( band, outputReal, outputImag, shArray );
431 
432  //================================================ Free memory
433  releaseSphericalMemory ( inputReal, inputImag, outputReal, outputImag, tableSpaceHelper, tablePml, shWeights, workspace, fftPlan, dctPlan, band );
434 
435  //================================================ Done
436  return ;
437 
438 }
ProSHADE_internal_sphericalHarmonics::initialiseAllMemory
void initialiseAllMemory(proshade_unsign band, proshade_double *&inputReal, proshade_double *&inputImag, proshade_double *&outputReal, proshade_double *&outputImag, double *&shWeights, double **&tableSpace, double *&tableSpaceHelper, fftw_complex *&workspace, proshade_double *&rres, proshade_double *&ires, proshade_double *&fltres, proshade_double *&scratchpad, fftw_plan &fftPlan, fftw_plan &dctPlan)
This function initialises all the memory required for spherical harmonics computation.
Definition: ProSHADE_sphericalHarmonics.cpp:218
ProSHADE_internal_sphericalHarmonics::initialSplitDiscreteTransform
void initialSplitDiscreteTransform(proshade_unsign oneDim, proshade_double *&inputReal, proshade_double *&inputImag, proshade_double *&rres, proshade_double *&ires, proshade_double *mappedData, fftw_plan &fftPlan, proshade_double normCoeff)
This function computes the spherical harmonics of a aingle shell, saving them in supplied pointer.
Definition: ProSHADE_sphericalHarmonics.cpp:258
ProSHADE_internal_sphericalHarmonics::releaseSphericalMemory
void releaseSphericalMemory(proshade_double *&inputReal, proshade_double *&inputImag, proshade_double *&outputReal, proshade_double *&outputImag, double *&tableSpaceHelper, double **&tableSpace, double *&shWeights, fftw_complex *&workspace, fftw_plan &fftPlan, fftw_plan &dctPlan, proshade_unsign band)
This function computes the spherical harmonics of a aingle shell, saving them in supplied pointer.
Definition: ProSHADE_sphericalHarmonics.cpp:171
ProSHADE_internal_sphericalHarmonics::computeSphericalTransformCoeffs
void computeSphericalTransformCoeffs(proshade_unsign band, proshade_double *&rdataptr, proshade_double *&idataptr, proshade_double *&outputReal, proshade_double *&outputImag, proshade_double *&rres, proshade_double *&ires, proshade_double *&fltres, proshade_double *&scratchpad, double **&tablePml, double *&shWeights, fftw_plan &dctPlan)
This function takes the split discrete transform and proceeds to complete the spherical harmonics dec...
Definition: ProSHADE_sphericalHarmonics.cpp:301
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:65
ProSHADE_internal_sphericalHarmonics::placeWithinWorkspacePointers
void placeWithinWorkspacePointers(fftw_complex *&workspace, proshade_unsign oDim, proshade_double *&rres, proshade_double *&ires, proshade_double *&fltres, proshade_double *&scratchpad)
This function takes the workspace pointer and correctly places the other internal pointers.
Definition: ProSHADE_sphericalHarmonics.cpp:87
ProSHADE_internal_sphericalHarmonics::applyCondonShortleyPhase
void applyCondonShortleyPhase(proshade_unsign band, proshade_double *outputReal, proshade_double *outputImag, proshade_complex *&shArray)
This is the final step in computing the full spherical harmonics decomposition of the input data.
Definition: ProSHADE_sphericalHarmonics.cpp:353
ProSHADE_internal_sphericalHarmonics::allocateComputationMemory
void allocateComputationMemory(proshade_unsign band, proshade_double *&inputReal, proshade_double *&inputImag, proshade_double *&outputReal, proshade_double *&outputImag, double *&shWeights, double **&tableSpace, double *&tableSpaceHelper, fftw_complex *&workspace)
This function determines the integration order for the between spheres integration.
Definition: ProSHADE_sphericalHarmonics.cpp:40
ProSHADE_internal_sphericalHarmonics::initialiseFFTWPlans
void initialiseFFTWPlans(proshade_unsign band, fftw_plan &fftPlan, fftw_plan &dctPlan, proshade_double *&inputReal, proshade_double *&inputImag, proshade_double *&rres, proshade_double *&ires, proshade_double *&scratchpad)
This function initialises the FFTW plans.
Definition: ProSHADE_sphericalHarmonics.cpp:114
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:395
ProSHADE_sphericalHarmonics.hpp
This header file declares the functions required to compute the spherical harmonics decomposition for...