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