ProSHADE  0.6.5 (NOV 2018)
Protein Shape Descriptors and Symmetry Detection
ProSHADE_wigner.cpp
Go to the documentation of this file.
1 
22 //============================================ FFTW3 + SOFT
23 #ifdef __cplusplus
24 extern "C" {
25 #endif
26 #include <fftw3.h>
27 #include <wrap_fftw.h>
28 #include <makeweights.h>
29 #include <s2_primitive.h>
30 #include <s2_cospmls.h>
31 #include <s2_legendreTransforms.h>
32 #include <s2_semi_fly.h>
33 #include <rotate_so3_utils.h>
34 #include <utils_so3.h>
35 #include <soft_fftw.h>
36 #include <rotate_so3_fftw.h>
37 #ifdef __cplusplus
38 }
39 #endif
40 
41 //============================================ RVAPI
42 #include <rvapi_interface.h>
43 
44 //============================================ ProSHADE
45 #include "ProSHADE.h"
46 #include "ProSHADE_internal.h"
47 
58 void ProSHADE_internal::ProSHADE_comparePairwise::generateWignerMatrices ( ProSHADE::ProSHADE_settings* settings )
59 {
60  //======================================== Sanity check
61  if ( !this->_eulerAnglesFound )
62  {
63  std::cerr << "!!! ProSHADE ERROR !!! Error in file comparison !!! Attempted to get Wigner d matrices without computing the Euler angles first. Call the getEulerAngles function BEFORE the generateWignerMatrices function." << std::endl;
64 
65  if ( settings->htmlReport )
66  {
67  std::stringstream hlpSS;
68  hlpSS << "<font color=\"red\">" << "Attempted to get Wigner d matrices without computing the Euler angles first. This looks like an internal bug, please report this case." << "</font>";
69  rvapi_set_text ( hlpSS.str().c_str(),
70  "ProgressSection",
71  settings->htmlReportLineProgress,
72  1,
73  1,
74  1 );
75  settings->htmlReportLineProgress += 1;
76  rvapi_flush ( );
77  }
78 
79  exit ( -1 );
80  }
81 
82  //======================================== Initialise local variables
83  double trigVals[2];
84  unsigned int noOrders = 0;
85  unsigned int arrConvIter = 0;
86  double *expARStart;
87  double *expAIStart;
88  double *expGRStart;
89  double *expGIStart;
90  double Dij;
91  double eARi;
92  double eAIi;
93  double eGRj;
94  double eGIj;
95  double iSign = 1.0;
96  double rSign = 1.0;
97 
98  double angAlpha = this->_eulerAngles[0];
99  double angBeta = this->_eulerAngles[1];
100  double angGamma = this->_eulerAngles[2];
101 
102  if ( angAlpha < 0.0 ) { angAlpha = M_PI + ( angAlpha - M_PI ); }
103  if ( angBeta < 0.0 ) { angBeta = M_PI + ( angBeta - M_PI ); }
104  if ( angGamma < 0.0 ) { angGamma = M_PI + ( angGamma - M_PI ); }
105 
106  //======================================== Initialise internal values
107  this->_wignerMatrices = std::vector< std::vector< std::vector< std::array<double,2> > > > ( this->_bandwidthLimit );
108 
109  //======================================== Allocate memory
110  double *matIn = new double[static_cast<unsigned int> ( 4 * pow( this->_bandwidthLimit, 2.0 ) - 4 * this->_bandwidthLimit + 1 )];
111  double *matOut = new double[static_cast<unsigned int> ( 4 * pow( this->_bandwidthLimit, 2.0 ) - 4 * this->_bandwidthLimit + 1 )];
112  double *sqrts = new double[static_cast<unsigned int> ( 2 * this->_bandwidthLimit )];
113  double *workspace = new double[static_cast<unsigned int> ( 4 * pow( this->_bandwidthLimit, 2.0 ) )];
114  double *exponentAlphaReal = new double[static_cast<unsigned int> ( 2 * this->_bandwidthLimit - 1 )];
115  double *exponentAlphaImag = new double[static_cast<unsigned int> ( 2 * this->_bandwidthLimit - 1 )];
116  double *exponentGammaReal = new double[static_cast<unsigned int> ( 2 * this->_bandwidthLimit - 1 )];
117  double *exponentGammaImag = new double[static_cast<unsigned int> ( 2 * this->_bandwidthLimit - 1 )];
118 
119  //======================================== Compute the square roots
120  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( 2 * this->_bandwidthLimit ); iter++ )
121  {
122  sqrts[iter] = static_cast<double> ( sqrt ( static_cast<double> ( iter ) ) );
123  }
124 
125  //======================================== Compute the trig values
126  trigVals[0] = static_cast<double> ( cos ( 0.5 * -angBeta ) );
127  trigVals[1] = static_cast<double> ( sin ( 0.5 * -angBeta ) );
128 
129  //======================================== Get alpha and gamma exponents
130  genExp ( this->_bandwidthLimit, angAlpha, exponentAlphaReal, exponentAlphaImag );
131  genExp ( this->_bandwidthLimit, angGamma, exponentGammaReal, exponentGammaImag );
132 
133  //======================================== For each band, find the Wigned d matrix
134  for ( unsigned int bandIter = 0; bandIter < this->_bandwidthLimit; bandIter++ )
135  {
136  //==================================== Initialise
137  noOrders = 2 * bandIter + 1;
138  arrConvIter = 0;
139 
140  //==================================== Initialise exponent pointers
141  expARStart = &exponentAlphaReal[ (this->_bandwidthLimit - 1) - bandIter ];
142  expAIStart = &exponentAlphaImag[ (this->_bandwidthLimit - 1) - bandIter ];
143  expGRStart = &exponentGammaReal[ (this->_bandwidthLimit - 1) - bandIter ];
144  expGIStart = &exponentGammaImag[ (this->_bandwidthLimit - 1) - bandIter ];
145  iSign = 1.0;
146  rSign = 1.0;
147 
148  //==================================== Initialise the matrix holding structure
149  this->_wignerMatrices.at(bandIter) = std::vector< std::vector< std::array<double,2> > > ( noOrders );
150  for ( unsigned int iter = 0; iter < noOrders; iter++ )
151  {
152  this->_wignerMatrices.at(bandIter).at(iter) = std::vector< std::array<double,2> > ( noOrders );
153  }
154 
155  //==================================== Get wigner d matrix values using beta angles only
156  wignerdmat ( bandIter, matIn, matOut, trigVals, sqrts, workspace );
157 
158  //==================================== Multiply the wigner d matrix by alpha and gamma values and save the wigner D matrix to output array
159  for ( unsigned int d1Iter = 0; d1Iter < noOrders; d1Iter++ )
160  {
161  eARi = expARStart[d1Iter];
162  eAIi = expAIStart[d1Iter];
163 
164  for ( unsigned int d2Iter = 0; d2Iter < noOrders; d2Iter++ )
165  {
166  Dij = matOut[arrConvIter];
167  eGRj = expGRStart[d2Iter];
168  eGIj = expGIStart[d2Iter];
169 
170  this->_wignerMatrices.at(bandIter).at(d1Iter).at(d2Iter)[0] = ( Dij * eGRj * eARi - Dij * eGIj * eAIi ) * rSign;
171  this->_wignerMatrices.at(bandIter).at(d1Iter).at(d2Iter)[1] = ( Dij * eGRj * eAIi + Dij * eGIj * eARi ) * iSign;
172 
173  arrConvIter += 1;
174  iSign *= -1.0;
175  rSign *= -1.0;
176  }
177  }
178 
179  //==================================== Get ready for next wigner matrix calculation
180  memcpy ( matIn, matOut, sizeof ( double ) * ( noOrders * noOrders ) );
181  }
182 
183  //======================================== Free memory
184  delete[] matIn;
185  delete[] matOut;
186  delete[] sqrts;
187  delete[] workspace;
188  delete[] exponentAlphaReal;
189  delete[] exponentAlphaImag;
190  delete[] exponentGammaReal;
191  delete[] exponentGammaImag;
192 
193  //======================================== Done
194  this->_wignerMatricesComputed = true;
195 }
196 
207 void ProSHADE_internal::ProSHADE_compareOneAgainstAll::generateWignerMatrices ( ProSHADE::ProSHADE_settings* settings )
208 {
209  //======================================== Sanity check
210  if ( !this->_eulerAnglesFound )
211  {
212  std::cerr << "!!! ProSHADE ERROR !!! Error in file comparison !!! Attempted to get Wigner d matrices without computing the Euler angles first. Call the getEulerAngles function BEFORE the generateWignerMatrices function." << std::endl;
213 
214  if ( settings->htmlReport )
215  {
216  std::stringstream hlpSS;
217  hlpSS << "<font color=\"red\">" << "Attempted to get Wigner d matrices without computing the Euler angles first. This looks like an internal bug, please report this case." << "</font>";
218  rvapi_set_text ( hlpSS.str().c_str(),
219  "ProgressSection",
220  settings->htmlReportLineProgress,
221  1,
222  1,
223  1 );
224  settings->htmlReportLineProgress += 1;
225  rvapi_flush ( );
226  }
227 
228  exit ( -1 );
229  }
230 
231  //======================================== Declare internal values
232  double trigVals[2];
233  unsigned int noOrders;
234  unsigned int arrConvIter;
235  double *expARStart;
236  double *expAIStart;
237  double *expGRStart;
238  double *expGIStart;
239  double Dij;
240  double eARi;
241  double eAIi;
242  double eGRj;
243  double eGIj;
244  double iSign;
245  double rSign;
246 
247  double angAlpha;
248  double angBeta;
249  double angGamma;
250 
251  int localComparisonBandLim;
252 
253  std::vector< std::vector< std::vector< std::array<double,2> > > > wignerMatricesHlp;
254  std::vector< std::vector< std::array<double,2> > > wignerMatricesHlpHlp;
255 
256  for ( unsigned int strIt = 0; strIt < static_cast<unsigned int> ( this->all->size() ); strIt++ )
257  {
258  //==================================== Initialise local variables
259  noOrders = 0;
260  arrConvIter = 0;
261 
262  iSign = 1.0;
263  rSign = 1.0;
264 
265  angAlpha = this->_eulerAngles.at(strIt)[0];
266  angBeta = this->_eulerAngles.at(strIt)[1];
267  angGamma = this->_eulerAngles.at(strIt)[2];
268 
269  wignerMatricesHlp.clear ( );
270  wignerMatricesHlpHlp.clear ( );
271 
272  //==================================== Find the appropriate bandwidth
273  localComparisonBandLim = std::min ( this->one->_bandwidthLimit, this->all->at(strIt)->_bandwidthLimit );
274 
275  //==================================== Allocate memory
276  double *matIn = new double[static_cast<unsigned int> ( 4 * pow( localComparisonBandLim, 2.0 ) - 4 * localComparisonBandLim + 1 )];
277  double *matOut = new double[static_cast<unsigned int> ( 4 * pow( localComparisonBandLim, 2.0 ) - 4 * localComparisonBandLim + 1 )];
278  double *sqrts = new double[static_cast<unsigned int> ( 2 * localComparisonBandLim )];
279  double *workspace = new double[static_cast<unsigned int> ( 4 * pow( localComparisonBandLim, 2.0 ) )];
280  double *exponentAlphaReal = new double[static_cast<unsigned int> ( 2 * localComparisonBandLim - 1 )];
281  double *exponentAlphaImag = new double[static_cast<unsigned int> ( 2 * localComparisonBandLim - 1 )];
282  double *exponentGammaReal = new double[static_cast<unsigned int> ( 2 * localComparisonBandLim - 1 )];
283  double *exponentGammaImag = new double[static_cast<unsigned int> ( 2 * localComparisonBandLim - 1 )];
284 
285  //==================================== Compute the square roots
286  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( 2 * localComparisonBandLim ); iter++ )
287  {
288  sqrts[iter] = static_cast<double> ( sqrt ( static_cast<double> ( iter ) ) );
289  }
290 
291  //==================================== Compute the trig values
292  trigVals[0] = static_cast<double> ( cos ( 0.5 * -angBeta ) );
293  trigVals[1] = static_cast<double> ( sin ( 0.5 * -angBeta ) );
294 
295  //==================================== Get alpha and gamma exponents
296  genExp ( localComparisonBandLim, angAlpha, exponentAlphaReal, exponentAlphaImag );
297  genExp ( localComparisonBandLim, angGamma, exponentGammaReal, exponentGammaImag );
298 
299  //==================================== For each band, find the Wigned d matrix
300  for ( int bandIter = 0; bandIter < localComparisonBandLim; bandIter++ )
301  {
302  //================================ Initialise
303  noOrders = 2 * bandIter + 1;
304  arrConvIter = 0;
305 
306  //================================ Initialise exponent pointers
307  expARStart = &exponentAlphaReal[ (localComparisonBandLim - 1) - bandIter ];
308  expAIStart = &exponentAlphaImag[ (localComparisonBandLim - 1) - bandIter ];
309  expGRStart = &exponentGammaReal[ (localComparisonBandLim - 1) - bandIter ];
310  expGIStart = &exponentGammaImag[ (localComparisonBandLim - 1) - bandIter ];
311  iSign = 1.0;
312  rSign = 1.0;
313 
314  //================================ Initialise the matrix holding structure
315  wignerMatricesHlpHlp = std::vector< std::vector< std::array<double,2> > > ( noOrders );
316  for ( unsigned int iter = 0; iter < noOrders; iter++ )
317  {
318  wignerMatricesHlpHlp.at(iter) = std::vector< std::array<double,2> > ( noOrders );
319  }
320 
321  //================================ Get wigner d matrix values using beta angles only
322  wignerdmat ( bandIter, matIn, matOut, trigVals, sqrts, workspace );
323 
324  //================================ Multiply the wigner d matrix by alpha and gamma values and save the wigner D matrix to output array
325  for ( unsigned int d1Iter = 0; d1Iter < noOrders; d1Iter++ )
326  {
327  eARi = expARStart[d1Iter];
328  eAIi = expAIStart[d1Iter];
329 
330  for ( unsigned int d2Iter = 0; d2Iter < noOrders; d2Iter++ )
331  {
332  Dij = matOut[arrConvIter];
333  eGRj = expGRStart[d2Iter];
334  eGIj = expGIStart[d2Iter];
335 
336  wignerMatricesHlpHlp.at(d1Iter).at(d2Iter)[0] = ( Dij * eGRj * eARi - Dij * eGIj * eAIi ) * rSign;
337  wignerMatricesHlpHlp.at(d1Iter).at(d2Iter)[1] = ( Dij * eGRj * eAIi + Dij * eGIj * eARi ) * iSign;
338 
339  arrConvIter += 1;
340  iSign *= -1.0;
341  rSign *= -1.0;
342  }
343  }
344 
345  //================================ Get ready for next wigner matrix calculation
346  memcpy ( matIn, matOut, sizeof ( double ) * ( noOrders * noOrders ) );
347 
348  //================================ Save results
349  wignerMatricesHlp.emplace_back ( wignerMatricesHlpHlp );
350  }
351 
352  //==================================== Free memory
353  delete[] matIn;
354  delete[] matOut;
355  delete[] sqrts;
356  delete[] workspace;
357  delete[] exponentAlphaReal;
358  delete[] exponentAlphaImag;
359  delete[] exponentGammaReal;
360  delete[] exponentGammaImag;
361 
362  //==================================== Save results
363  this->_wignerMatrices.emplace_back ( wignerMatricesHlp );
364  }
365 
366  //======================================== Done
367  this->_wignerMatricesComputed = true;
368 }
369 
ProSHADE_internal::ProSHADE_comparePairwise::generateWignerMatrices
void generateWignerMatrices(ProSHADE::ProSHADE_settings *settings)
This function is responsible for computing the Wigner D matrices for symmetry detection purposes...
Definition: ProSHADE_wigner.cpp:58
ProSHADE_internal::ProSHADE_compareOneAgainstAll::generateWignerMatrices
void generateWignerMatrices(ProSHADE::ProSHADE_settings *settings)
This function is responsible for computing the Wigner D matrices for full rotation function distance ...
Definition: ProSHADE_wigner.cpp:207
ProSHADE::ProSHADE_settings::htmlReport
bool htmlReport
Should HTML report for the run be created?
Definition: ProSHADE.h:186
ProSHADE.h
The main header file containing all declarations the user of the library needs.
ProSHADE_internal.h
The main header file containing all declarations for the innter workings of the library.
ProSHADE::ProSHADE_settings::htmlReportLineProgress
int htmlReportLineProgress
Iterator for current HTML line in the progress bar.
Definition: ProSHADE.h:188
ProSHADE::ProSHADE_settings
This class stores all the settings and is passed to the executive classes instead of multitude of par...
Definition: ProSHADE.h:74