ProSHADE  0.7.5.4 (MAR 2021)
Protein Shape Detection
ProSHADE_maths.cpp
Go to the documentation of this file.
1 
23 //==================================================== ProSHADE
24 #include "ProSHADE_maths.hpp"
25 
38 void ProSHADE_internal_maths::complexMultiplication ( proshade_double* r1, proshade_double* i1, proshade_double* r2, proshade_double* i2, proshade_double* retReal, proshade_double* retImag )
39 {
40  //================================================ Multiplication
41  *retReal = (*r1)*(*r2) - (*i1)*(*i2);
42  *retImag = (*r1)*(*i2) + (*i1)*(*r2);
43 
44  //================================================ Return
45  return ;
46 
47 }
48 
62 void ProSHADE_internal_maths::complexMultiplicationConjug ( proshade_double* r1, proshade_double* i1, proshade_double* r2, proshade_double* i2, proshade_double* retReal, proshade_double* retImag )
63 {
64  //================================================ Multiplication
65  *retReal = (*r1)*(*r2) + (*i1)*(*i2);
66  *retImag = -(*r1)*(*i2) + (*i1)*(*r2);
67 
68  //================================================ Return
69  return ;
70 
71 }
72 
83 proshade_double ProSHADE_internal_maths::complexMultiplicationRealOnly ( proshade_double* r1, proshade_double* i1, proshade_double* r2, proshade_double* i2 )
84 {
85  //================================================ Multiplication
86  proshade_double ret = (*r1)*(*r2) - (*i1)*(*i2);
87 
88  //================================================ Return
89  return ( ret );
90 
91 }
92 
103 proshade_double ProSHADE_internal_maths::complexMultiplicationConjugRealOnly ( proshade_double* r1, proshade_double* i1, proshade_double* r2, proshade_double* i2 )
104 {
105  //================================================ Multiplication
106  proshade_double ret = (*r1)*(*r2) + (*i1)*(*i2);
107 
108  //================================================ Return
109  return ( ret );
110 
111 }
112 
121 void ProSHADE_internal_maths::vectorMeanAndSD ( std::vector<proshade_double>* vec, proshade_double*& ret )
122 {
123  //================================================ Get mean
124  ret[0] = std::accumulate ( vec->begin(), vec->end(), 0.0 ) / static_cast<proshade_double> ( vec->size() );
125 
126  //================================================ Get standard deviation
127  proshade_double squaredSum = std::inner_product ( vec->begin(), vec->end(), vec->begin(), 0.0 );
128  ret[1] = std::sqrt ( ( squaredSum / static_cast<proshade_double> ( vec->size() ) ) - std::pow ( ret[0], 2.0 ) );
129 
130  //================================================ Check for NaN's
131  if ( ret[0] != ret[0] ) { ret[0] = 0.0; }
132  if ( ret[1] != ret[1] ) { ret[1] = 0.0; }
133 
134  //================================================ Return
135  return ;
136 
137 }
138 
147 void ProSHADE_internal_maths::vectorMedianAndIQR ( std::vector<proshade_double>* vec, proshade_double*& ret )
148 {
149  //================================================ Sanity check
150  if ( vec->size() < 3 ) { ret[0] = 0.0; ret[1] = 0.0; return; }
151 
152  //================================================ Sort the vector
153  std::sort ( vec->begin(), vec->end() );
154 
155  //================================================ Get median
156  if ( static_cast<proshade_unsign> ( vec->size() ) % 2 == 0)
157  {
158  ret[0] = ( vec->at( ( static_cast<proshade_unsign> ( vec->size() ) / 2 ) - 1 ) +
159  vec->at( static_cast<proshade_unsign> ( vec->size() ) / 2 ) ) / 2.0;
160  }
161  else
162  {
163  ret[0] = vec->at( static_cast<proshade_unsign> ( vec->size() ) / 2 );
164  }
165 
166  //================================================ Get first and third quartile
167  proshade_double Q1, Q3;
168  if ( static_cast<proshade_unsign> ( vec->size() ) % 2 == 0)
169  {
170  Q1 = ( vec->at( ( static_cast<proshade_unsign> ( vec->size() ) / 4 ) - 1 ) +
171  vec->at( static_cast<proshade_unsign> ( vec->size() ) / 4 ) ) / 2.0;
172  Q3 = ( vec->at( ( ( static_cast<proshade_unsign> ( vec->size() ) / 4 ) * 3 ) - 1 ) +
173  vec->at( ( static_cast<proshade_unsign> ( vec->size() ) / 4 ) * 3 ) ) / 2.0;
174  }
175  else
176  {
177  Q1 = vec->at( static_cast<proshade_unsign> ( vec->size() ) / 4 );
178  Q3 = vec->at( ( static_cast<proshade_unsign> ( vec->size() ) / 4 ) * 3 );
179  }
180 
181  //================================================ And now save the IQR
182  ret[1] = Q3 - Q1;
183 
184  //================================================ Return
185  return ;
186 
187 }
188 
198 void ProSHADE_internal_maths::arrayMedianAndIQR ( proshade_double* vec, proshade_unsign vecSize, proshade_double*& ret )
199 {
200  //================================================ Sort the vector
201  std::sort ( vec, vec + vecSize );
202 
203  //================================================ Get median
204  if ( vecSize % 2 == 0)
205  {
206  ret[0] = ( vec[ ( vecSize / 2 ) - 1 ] + vec[ vecSize / 2 ] ) / 2.0;
207  }
208  else
209  {
210  ret[0] = vec[ vecSize / 2 ];
211  }
212 
213  //================================================ Get first and third quartile
214  proshade_double Q1, Q3;
215  if ( vecSize % 2 == 0)
216  {
217  Q1 = ( vec[ ( vecSize / 4 ) - 1 ] + vec[ vecSize / 4 ] ) / 2.0;
218  Q3 = ( vec[ ( ( vecSize / 4 ) * 3 ) - 1 ] + vec[ ( vecSize / 4 ) * 3 ] ) / 2.0;
219  }
220  else
221  {
222  Q1 = vec[ vecSize / 4 ];
223  Q3 = vec[ ( vecSize / 4 ) * 3 ];
224  }
225 
226  //================================================ And now save the IQR
227  ret[1] = Q3 - Q1;
228 
229  //================================================ Return
230  return ;
231 
232 }
233 
244 proshade_double ProSHADE_internal_maths::pearsonCorrCoeff ( proshade_double* valSet1, proshade_double* valSet2, proshade_unsign length )
245 {
246  //================================================ Find vector means
247  proshade_double xMean = 0.0;
248  proshade_double yMean = 0.0;
249  proshade_double zeroCount = 0.0;
250  for ( proshade_unsign iter = 0; iter < length; iter++ )
251  {
252  xMean += valSet1[iter];
253  yMean += valSet2[iter];
254  }
255  xMean /= static_cast<proshade_double> ( length - zeroCount );
256  yMean /= static_cast<proshade_double> ( length - zeroCount );
257 
258  //================================================ Get Pearson's correlation coefficient
259  proshade_double xmmymm = 0.0;
260  proshade_double xmmsq = 0.0;
261  proshade_double ymmsq = 0.0;
262  for ( proshade_unsign iter = 0; iter < length; iter++ )
263  {
264  xmmymm += ( valSet1[iter] - xMean ) * ( valSet2[iter] - yMean );
265  xmmsq += pow( valSet1[iter] - xMean, 2.0 );
266  ymmsq += pow( valSet2[iter] - yMean, 2.0 );
267  }
268 
269  proshade_double ret = xmmymm / ( sqrt(xmmsq) * sqrt(ymmsq) );
270 
271  //================================================ Done
272  if ( std::isnan ( ret ) ) { return ( 0.0 ); }
273  return ( ret );
274 
275 }
276 
287 void ProSHADE_internal_maths::getLegendreAbscAndWeights ( proshade_unsign order, proshade_double* abscissas, proshade_double* weights, proshade_unsign taylorSeriesCap )
288 {
289  //================================================ Sanity check
290  if ( order < 2 )
291  {
292  throw ProSHADE_exception ( "The integration order is too low.", "EI00019", __FILE__, __LINE__, __func__, "The Gauss-Legendre integration order is less than 2. This\n : seems very low; if you have a very small structure or very\n : low resolution, please manually increase the integration\n : order. Otherwise, please report this as a bug." );
293  }
294 
295  //================================================ Initialise
296  proshade_double polyValue = 0.0;
297  proshade_double deriValue = 0.0;
298  proshade_double weightSum = 0.0;
299 
300  //================================================ Find the polynomial and derivative values at 0
301  getGLPolyAtZero ( order,
302  &polyValue,
303  &deriValue );
304 
305  //================================================ If the order is odd, then 0 is a root ...
306  if ( order % 2 == 1 )
307  {
308  abscissas[((order-1)/2)] = polyValue;
309  weights[((order-1)/2)] = deriValue;
310  }
311  else
312  {
313  // ... and if order is even, find the first root
314  getGLFirstEvenRoot ( polyValue, order, &abscissas[(order/2)], &weights[(order/2)], taylorSeriesCap );
315  }
316 
317  //================================================ Now, having computed the first roots, complete the series
318  completeLegendreSeries ( order, abscissas, weights, taylorSeriesCap );
319 
320  //================================================ Correct weights by anscissa values
321  for ( proshade_unsign iter = 0; iter < order; iter++ )
322  {
323  weights[iter] = 2.0 / ( 1.0 - abscissas[iter] ) / ( 1.0 + abscissas[iter] ) / weights[iter] / weights[iter];
324  weightSum = weightSum + weights[iter];
325  }
326 
327  //================================================ Normalise weights
328  for ( proshade_unsign iter = 0; iter < order; iter++ )
329  {
330  weights[iter] = 2.0 * weights[iter] / weightSum;
331  }
332 
333  //================================================ Done
334  return ;
335 }
336 
347 void ProSHADE_internal_maths::getGLPolyAtZero ( proshade_unsign order, proshade_double *polyValue, proshade_double *deriValue )
348 {
349  //================================================ Initialise
350  proshade_double hlpVal = 0.0;
351  proshade_double prevPoly = 1.0;
352  proshade_double prevPrevPoly = 0.0;
353  proshade_double prevDeri = 0.0;
354  proshade_double prevPrevDeri = 0.0;
355 
356  for ( proshade_unsign ordIt = 0; ordIt < order; ordIt++ )
357  {
358  hlpVal = static_cast<proshade_double> ( ordIt );
359  *polyValue = -hlpVal * prevPrevPoly / ( hlpVal + 1.0 );
360  *deriValue = ( ( 2.0 * hlpVal + 1.0 ) * prevPoly - hlpVal * prevPrevDeri ) / ( hlpVal + 1.0 );
361  prevPrevPoly = prevPoly;
362  prevPoly = *polyValue;
363  prevPrevDeri = prevDeri;
364  prevDeri = *deriValue;
365  }
366 
367  //================================================ Done
368  return ;
369 
370 }
371 
384 void ProSHADE_internal_maths::getGLFirstEvenRoot ( proshade_double polyAtZero, proshade_unsign order, proshade_double *abscAtZero, proshade_double *weighAtZero, proshade_unsign taylorSeriesCap )
385 {
386  //================================================ Sanity check
387  if ( taylorSeriesCap < 2 )
388  {
389  throw ProSHADE_exception ( "The Taylor series cap is too low.", "EI00020", __FILE__, __LINE__, __func__, "The Taylor series expansion limit is less than 2. This\n : seems very low; if you have a very small structure or very\n : low resolution, please manually increase the integration\n : order. Otherwise, please report this as a bug." );
390  }
391 
392  //================================================ Initialise variables
393  *abscAtZero = advanceGLPolyValue ( 0.0, -M_PI / 2.0, 0.0, order, taylorSeriesCap );
394  proshade_double hlp = 0.0;
395  proshade_double hlpVal = static_cast<proshade_double> ( order );
396  proshade_double *abscSteps;
397  proshade_double *weightSteps;
398 
399  //================================================ Allocate memory
400  abscSteps = new proshade_double [taylorSeriesCap+2];
401  weightSteps = new proshade_double [taylorSeriesCap+1];
402 
403  //================================================ Pre-set values
404  abscSteps[0] = 0.0;
405  abscSteps[1] = polyAtZero;
406  weightSteps[0] = 0.0;
407 
408  //================================================ Fill in abscissa and weight steps
409  for ( proshade_unsign iter = 0; iter <= taylorSeriesCap - 2; iter = iter + 2 )
410  {
411  hlp = static_cast<proshade_double> ( iter );
412 
413  abscSteps[iter+2] = 0.0;
414  abscSteps[iter+3] = ( hlp * ( hlp + 1.0 ) - hlpVal * ( hlpVal + 1.0 ) ) * abscSteps[iter+1] / (hlp + 1.0) / (hlp + 2.0 );
415 
416  weightSteps[iter+1] = 0.0;
417  weightSteps[iter+2] = ( hlp + 2.0 ) * abscSteps[iter+3];
418  }
419 
420  //================================================ Find abscissa and weights
421  for ( proshade_double iter = 0; iter < 5; iter++ )
422  {
423  *abscAtZero = *abscAtZero - evaluateGLSeries ( abscSteps, *abscAtZero, taylorSeriesCap ) / evaluateGLSeries ( weightSteps, *abscAtZero, taylorSeriesCap-1 );
424  }
425  *weighAtZero = evaluateGLSeries ( weightSteps, *abscAtZero, taylorSeriesCap-1 );
426 
427  //================================================ Free memory
428  delete abscSteps;
429  delete weightSteps;
430 
431  //================================================ Done
432  return ;
433 
434 }
435 
447 proshade_double ProSHADE_internal_maths::evaluateGLSeries ( proshade_double *series, proshade_double target, proshade_unsign terms )
448 {
449  //================================================ Initalise
450  proshade_double factorialValue = 1.0;
451  proshade_double value = 0.0;
452 
453  //================================================ Compute
454  for ( proshade_unsign iter = 1; iter <= terms; iter++ )
455  {
456  value = value + series[iter] * factorialValue;
457  factorialValue = factorialValue * target;
458  }
459 
460  //================================================ Done
461  return ( value );
462 
463 }
464 
477 proshade_double ProSHADE_internal_maths::advanceGLPolyValue ( proshade_double from, proshade_double to, proshade_double valAtFrom, proshade_unsign noSteps, proshade_unsign taylorSeriesCap )
478 {
479  //================================================ Initialise variables
480  proshade_double hlpVal = 0.0;
481  proshade_double stepSize = 0.0;
482  proshade_double valChange = 0.0;
483  proshade_double valSecChange = 0.0;
484  proshade_double squareSteps = 0.0;
485  proshade_double curVal = 0.0;
486 
487  //================================================ Set initial values
488  stepSize = ( to - from ) / static_cast<proshade_double> ( taylorSeriesCap );
489  squareSteps = sqrt ( static_cast<proshade_double> ( noSteps * ( noSteps + 1 ) ) );
490  curVal = from;
491 
492  //================================================ Go through the series and iteratively improve the estimate
493  for ( proshade_unsign iter = 0; iter < taylorSeriesCap; iter++ )
494  {
495  hlpVal = ( 1.0 - valAtFrom ) * ( 1.0 + valAtFrom );
496  valChange = - stepSize * hlpVal / ( squareSteps * sqrt ( hlpVal ) - 0.5 * valAtFrom * sin ( 2.0 * curVal ) );
497  valAtFrom = valAtFrom + valChange;
498 
499  curVal = curVal + stepSize;
500 
501  hlpVal = ( 1.0 - valAtFrom ) * ( 1.0 + valAtFrom );
502  valSecChange = - stepSize * hlpVal / ( squareSteps * sqrt ( hlpVal ) - 0.5 * valAtFrom * sin ( 2.0 * curVal ) );
503  valAtFrom = valAtFrom + 0.5 * ( valSecChange - valChange );
504  }
505 
506  //================================================ Done
507  return valAtFrom;
508 
509 }
510 
521 void ProSHADE_internal_maths::completeLegendreSeries ( proshade_unsign order, proshade_double* abscissas, proshade_double* weights, proshade_unsign taylorSeriesCap )
522 {
523  //================================================ Initialise internal variables
524  proshade_double hlpTaylorVal = 0.0;
525  proshade_double hlpOrderVal = static_cast<proshade_double> ( order );
526  proshade_double abscValueChange = 0.0;
527  proshade_double prevAbsc = 0.0;
528  proshade_double *hlpAbscSeries;
529  proshade_double *hlpWeightSeries;
530  proshade_unsign noSeriesElems = 0;
531  proshade_unsign oddEvenSwitch = 0;
532 
533  //================================================ Pre-set internal values
534  if ( order % 2 == 1 )
535  {
536  noSeriesElems = ( order - 1 ) / 2 - 1;
537  oddEvenSwitch = 1;
538  }
539  else
540  {
541  noSeriesElems = order / 2 - 1;
542  oddEvenSwitch = 0;
543  }
544 
545  //================================================ Allocate memory
546  hlpAbscSeries = new proshade_double[taylorSeriesCap+2];
547  hlpWeightSeries = new proshade_double[taylorSeriesCap+1];
548 
549  //================================================ For each series element
550  for ( proshade_unsign serIt = noSeriesElems + 1; serIt < order - 1; serIt++ )
551  {
552  //============================================ Init loop
553  prevAbsc = abscissas[serIt];
554  abscValueChange = advanceGLPolyValue ( M_PI/2.0, -M_PI/2.0, prevAbsc, order, taylorSeriesCap ) - prevAbsc;
555 
556  //============================================ Init abscissas
557  hlpAbscSeries[0] = 0.0;
558  hlpAbscSeries[1] = 0.0;
559  hlpAbscSeries[2] = weights[serIt];
560 
561  //============================================ Init weights
562  hlpWeightSeries[0] = 0.0;
563  hlpWeightSeries[1] = hlpAbscSeries[2];
564 
565  //============================================ Taylor expansion
566  for ( proshade_unsign tayIt = 0; tayIt <= taylorSeriesCap - 2; tayIt++ )
567  {
568  hlpTaylorVal = static_cast<proshade_double> ( tayIt );
569 
570  hlpAbscSeries[tayIt+3] = ( 2.0 * prevAbsc * ( hlpTaylorVal + 1.0 ) * hlpAbscSeries[tayIt+2] + ( hlpTaylorVal * ( hlpTaylorVal + 1.0 ) - hlpOrderVal *
571  ( hlpOrderVal + 1.0 ) ) * hlpAbscSeries[tayIt+1] / ( hlpTaylorVal + 1.0 ) ) / ( 1.0 - prevAbsc ) / ( 1.0 + prevAbsc ) /
572  ( hlpTaylorVal + 2.0 );
573 
574  hlpWeightSeries[tayIt+2] = ( hlpTaylorVal + 2.0 ) * hlpAbscSeries[tayIt+3];
575  }
576 
577  //============================================ Sum over results
578  for ( proshade_unsign iter = 0; iter < 5; iter++ )
579  {
580  abscValueChange = abscValueChange - evaluateGLSeries ( hlpAbscSeries, abscValueChange, taylorSeriesCap ) /
581  evaluateGLSeries ( hlpWeightSeries, abscValueChange, taylorSeriesCap-1 );
582  }
583 
584  //============================================ Save results
585  abscissas[serIt+1] = prevAbsc + abscValueChange;
586  weights[serIt+1] = evaluateGLSeries ( hlpWeightSeries, abscValueChange, taylorSeriesCap - 1 );
587  }
588 
589  for ( proshade_unsign serIt = 0; serIt <= noSeriesElems + oddEvenSwitch; serIt++ )
590  {
591  abscissas[serIt] = -abscissas[order-serIt-1];
592  weights[serIt] = weights[order-serIt-1];
593  }
594 
595  //================================================ Free memory
596  delete hlpAbscSeries;
597  delete hlpWeightSeries;
598 
599  //================================================ Done
600  return ;
601 
602 }
603 
619 proshade_double ProSHADE_internal_maths::gaussLegendreIntegrationReal ( proshade_double* vals, proshade_unsign valsSize, proshade_unsign order, proshade_double* abscissas, proshade_double* weights, proshade_double integralOverRange, proshade_double maxSphereDists )
620 {
621  //================================================ Initialise local variables
622  proshade_double ret = 0.0;
623  proshade_complex* intData = new proshade_complex[order];
624  ProSHADE_internal_misc::checkMemoryAllocation ( intData, __FILE__, __LINE__, __func__ );
625  proshade_complex posVals;
626  proshade_unsign lesserPos = 0;
627  proshade_unsign upperPos = 0;
628  proshade_double lesserWeight = 0.0;
629  proshade_double upperWeight = 0.0;
630 
631  //================================================ Rescale to <order> points
632  for ( proshade_unsign absIter = 0; absIter < order; absIter++ )
633  {
634  //============================================ Init loop
635  posVals[0] = 0.0;
636  posVals[1] = 0.0;
637 
638  //============================================ Find real position of abscissas
639  posVals[0] = ( ( abscissas[absIter] + 1.0 ) / 2.0 ) * integralOverRange;
640 
641 
642  //============================================ Find lesser and upper bounds
643  for ( proshade_unsign valIt = 0; valIt < valsSize; valIt++ )
644  {
645  if ( ( ( valIt * maxSphereDists ) <= posVals[0] ) && ( ( ( valIt + 1 ) * maxSphereDists ) > posVals[0] ) )
646  {
647  lesserPos = static_cast<proshade_unsign> ( valIt );
648  upperPos = static_cast<proshade_unsign> ( valIt + 1 );
649  break;
650  }
651  }
652 
653  //============================================ Linear Interpolation
654  lesserWeight = 0.0;
655  upperWeight = 0.0;
656  if ( lesserPos != 0 )
657  {
658  //======================================== Here we realise that the lesser and upper bounds were determined on scale 1 ... N, while our values are on scale 0 ... N-1 and therefore after determining the linear interpolation weights, we subtract 1 from both lesserPos and upperPos; however ...
659  lesserWeight = upperPos - ( posVals[0] / maxSphereDists );
660  upperWeight = 1.0 - lesserWeight;
661 
662  posVals[1] = ( lesserWeight * vals[lesserPos-1] ) + ( upperWeight * vals[upperPos-1] );
663  }
664  else
665  {
666  //======================================== ... this then means that we would require position -1 for when the integration value is between 0 and the first shell. To resolve this, we assume that the values are 0 below the first shell and proceed as follows:
667  upperWeight = 1.0 - ( upperPos - ( posVals[0] / maxSphereDists ) );
668 
669  posVals[1] = ( upperWeight * vals[upperPos-1] );
670  }
671 
672  intData[absIter][0] = posVals[0];
673  intData[absIter][1] = posVals[1];
674  }
675 
676  //================================================ Integrate
677  for ( proshade_unsign absPoint = 0; absPoint < order; absPoint++ )
678  {
679  ret += ( weights[absPoint] * intData[absPoint][1] );
680  }
681 
682  //================================================ Normalise
683  ret *= ( integralOverRange / 2.0 );
684 
685  //================================================ Release memory
686  delete[] intData;
687 
688  //================================================ Done
689  return ( ret );
690 
691 }
692 
709 void ProSHADE_internal_maths::gaussLegendreIntegration ( proshade_complex* vals, proshade_unsign valsSize, proshade_unsign order, proshade_double* abscissas, proshade_double* weights, proshade_double integralOverRange, proshade_double maxSphereDists, proshade_double* retReal, proshade_double* retImag )
710 {
711  //================================================ Initialise local variables
712  proshade_triplet* intData = new proshade_triplet [order];
713  ProSHADE_internal_misc::checkMemoryAllocation ( intData, __FILE__, __LINE__, __func__ );
714  proshade_triplet posVals;
715  proshade_unsign lesserPos = 0;
716  proshade_unsign upperPos = 0;
717  proshade_double lesserWeight = 0.0;
718  proshade_double upperWeight = 0.0;
719 
720  //================================================ Rescale to <order> points
721  for ( proshade_unsign absIter = 0; absIter < order; absIter++ )
722  {
723  //============================================ Init loop
724  posVals[0] = 0.0;
725  posVals[1] = 0.0;
726  posVals[2] = 0.0;
727 
728  //============================================ Find real position of abscissas
729  posVals[0] = ( ( abscissas[absIter] + 1.0 ) / 2.0 ) * integralOverRange;
730 
731 
732  //============================================ Find lesser and upper bounds
733  for ( proshade_unsign valIt = 0; valIt < valsSize; valIt++ )
734  {
735  if ( ( ( valIt * maxSphereDists ) <= posVals[0] ) && ( ( ( valIt + 1 ) * maxSphereDists ) > posVals[0] ) )
736  {
737  lesserPos = static_cast<proshade_unsign> ( valIt );
738  upperPos = static_cast<proshade_unsign> ( valIt + 1 );
739  break;
740  }
741  }
742 
743  //============================================ Linear Interpolation
744  lesserWeight = 0.0;
745  upperWeight = 0.0;
746  if ( lesserPos != 0 )
747  {
748  //======================================== Here we realise that the lesser and upper bounds were determined on scale 1 ... N, while our values are on scale 0 ... N-1 and therefore after determining the linear interpolation weights, we subtract 1 from both lesserPos and upperPos; however ...
749  lesserWeight = upperPos - ( posVals[0] / maxSphereDists );
750  upperWeight = 1.0 - lesserWeight;
751 
752  posVals[1] = ( lesserWeight * vals[lesserPos-1][0] ) + ( upperWeight * vals[upperPos-1][0] );
753  posVals[2] = ( lesserWeight * vals[lesserPos-1][1] ) + ( upperWeight * vals[upperPos-1][1] );
754  }
755  else
756  {
757  //======================================== ... this then means that we would require position -1 for when the integration value is between 0 and the first shell. To resolve this, we assume that the values are 0 below the first shell and proceed as follows:
758  upperWeight = 1.0 - ( upperPos - ( posVals[0] / maxSphereDists ) );
759 
760  posVals[1] = ( upperWeight * vals[upperPos-1][0] );
761  posVals[2] = ( upperWeight * vals[upperPos-1][1] );
762  }
763 
764  intData[absIter][0] = posVals[0];
765  intData[absIter][1] = posVals[1];
766  intData[absIter][2] = posVals[2];
767  }
768 
769  //================================================ Integrate
770  *retReal = 0.0;
771  *retImag = 0.0;
772  for ( proshade_unsign absPoint = 0; absPoint < order; absPoint++ )
773  {
774  *retReal += ( weights[absPoint] * intData[absPoint][1] );
775  *retImag += ( weights[absPoint] * intData[absPoint][2] );
776  }
777 
778  //================================================ Normalise
779  *retReal *= ( integralOverRange / 2.0 );
780  *retImag *= ( integralOverRange / 2.0 );
781 
782  //================================================ Release memory
783  delete[] intData;
784 
785  //================================================ Done
786  return ;
787 
788 }
789 
802 void ProSHADE_internal_maths::complexMatrixSVDSigmasOnly ( proshade_complex** mat, int dim, double*& singularValues )
803 {
804  //================================================ Initialise local variables
805  char job = 'N'; // Save computation of parts of U and V matrices, they are not needed here
806  std::complex<double> *rotMatU = new std::complex<double> [dim*dim]; // The U matrix space
807  std::complex<double> *rotMatV = new std::complex<double> [dim*dim]; // The V^T matrix space
808  std::complex<double> *work = new std::complex<double> [( 4 * dim)]; // Workspace, minimum required is 3*dim, using more for performance
809  int workDim = ( 4 * dim); // Formalism stating just that
810  double* rwork = new double[(7 * dim)]; // Required by LAPACK, from 3.7 requires 7 * dim
811  int* iwork = new int[(8 * dim)]; // Required by LAPACK
812  int returnValue = 0; // This will tell if operation succeeded
813 
814  //================================================ Check memory allocation
815  ProSHADE_internal_misc::checkMemoryAllocation ( rotMatU, __FILE__, __LINE__, __func__ );
816  ProSHADE_internal_misc::checkMemoryAllocation ( rotMatV, __FILE__, __LINE__, __func__ );
817  ProSHADE_internal_misc::checkMemoryAllocation ( work, __FILE__, __LINE__, __func__ );
818  ProSHADE_internal_misc::checkMemoryAllocation ( rwork, __FILE__, __LINE__, __func__ );
819  ProSHADE_internal_misc::checkMemoryAllocation ( iwork, __FILE__, __LINE__, __func__ );
820 
821  //================================================ Load input data into array in column-major order
822  std::complex<double> *matrixToDecompose = new std::complex<double>[dim*dim];
823  ProSHADE_internal_misc::checkMemoryAllocation ( matrixToDecompose, __FILE__, __LINE__, __func__ );
824  for ( int rowIt = 0; rowIt < dim; rowIt++ )
825  {
826  for ( int colIt = 0; colIt < dim; colIt++ )
827  {
828  matrixToDecompose[(colIt*dim)+rowIt] = std::complex<double> ( mat[rowIt][colIt][0], mat[rowIt][colIt][1] );
829  }
830  }
831 
832  //================================================ Run LAPACK ZGESDD
833  zgesdd_ ( &job, &dim, &dim, matrixToDecompose, &dim, singularValues, rotMatU, &dim, rotMatV, &dim,
834  work, &workDim, rwork, iwork, &returnValue );
835 
836  //================================================ Free memory
837  delete[] rotMatU;
838  delete[] rotMatV;
839  delete[] work;
840  delete[] rwork;
841  delete[] iwork;
842  delete[] matrixToDecompose;
843 
844  //================================================ Check result
845  if ( returnValue != 0 )
846  {
847  throw ProSHADE_exception ( "The LAPACK complex SVD algorithm did not converge!", "EL00021", __FILE__, __LINE__, __func__, "LAPACK algorithm for computing the singular value\n : decomposition of complex matrices did not converge and\n : therefore it was not possible to combine SH coefficients\n : from multiple shells. Changing the resolution may help,\n : contact me if this error persists." );
848  }
849 
850  //================================================ Done
851  return ;
852 
853 }
854 
869 void ProSHADE_internal_maths::complexMatrixSVDUandVOnly ( proshade_double* mat, int dim, proshade_double* uAndV, bool fail )
870 {
871  //================================================ Initialise local variables
872  char job = 'A'; // Save computation of parts of U and V matrices, they are not needed here
873  double* singularValues = new double[dim]; // The array of singular values
874  std::complex<double> *rotMatU = new std::complex<double> [dim*dim]; // The U matrix space
875  std::complex<double> *rotMatV = new std::complex<double> [dim*dim]; // The V^T matrix space
876  std::complex<double> *work = new std::complex<double> [static_cast<proshade_unsign>( ( 3 * dim) + pow( dim, 2 ) * dim)]; // Workspace, minimum required is 3*dim, using more for performance
877  int workDim = ( 3 * dim) + pow( dim, 2 ); // Formalism stating just that
878  double* rwork = new double[static_cast<proshade_unsign>((5 * dim) + 5 * pow(dim,2))]; // Required by LAPACK
879  int* iwork = new int[(8 * dim)]; // Required by LAPACK
880  int returnValue = 0; // This will tell if operation succeeded
881  ProSHADE_internal_misc::checkMemoryAllocation ( singularValues, __FILE__, __LINE__, __func__ );
882  ProSHADE_internal_misc::checkMemoryAllocation ( rotMatU, __FILE__, __LINE__, __func__ );
883  ProSHADE_internal_misc::checkMemoryAllocation ( rotMatV, __FILE__, __LINE__, __func__ );
884  ProSHADE_internal_misc::checkMemoryAllocation ( work, __FILE__, __LINE__, __func__ );
885  ProSHADE_internal_misc::checkMemoryAllocation ( rwork, __FILE__, __LINE__, __func__ );
886  ProSHADE_internal_misc::checkMemoryAllocation ( iwork, __FILE__, __LINE__, __func__ );
887 
888  //================================================ Load input data into array in column-major order
889  std::complex<double> *matrixToDecompose = new std::complex<double>[dim*dim];
890  ProSHADE_internal_misc::checkMemoryAllocation ( matrixToDecompose, __FILE__, __LINE__, __func__ );
891  for ( int rowIt = 0; rowIt < dim; rowIt++ )
892  {
893  for ( int colIt = 0; colIt < dim; colIt++ )
894  {
895  matrixToDecompose[(colIt*dim)+rowIt] = std::complex<double> ( mat[(rowIt*dim)+colIt], 0.0 );
896  }
897  }
898 
899  //================================================ Run LAPACK ZGESDD
900  zgesdd_ ( &job, &dim, &dim, matrixToDecompose, &dim, singularValues, rotMatU, &dim, rotMatV, &dim,
901  work, &workDim, rwork, iwork, &returnValue );
902 
903  //================================================ Free memory
904  delete[] work;
905  delete[] rwork;
906  delete[] iwork;
907  delete[] matrixToDecompose;
908  delete[] singularValues;
909 
910  //================================================ Check result
911  if ( ( returnValue != 0 ) && ( fail ) )
912  {
913  throw ProSHADE_exception ( "The LAPACK complex SVD algorithm did not converge!", "EL00022", __FILE__, __LINE__, __func__, "LAPACK algorithm for computing the singular value\n : decomposition of complex matrices did not converge and\n : therefore it was not possible to optimise the peak\n : positions in the (self-)rotation function. Changing the\n : resolution may help, contact me if this error persists." );
914  }
915  if ( ( returnValue != 0 ) && ( !fail ) )
916  {
917  uAndV[0] = -777.7;
918  return ;
919  }
920 
921  //================================================ Save U
922  for ( proshade_signed rowIt = 0; rowIt < dim; rowIt++ )
923  {
924  for ( proshade_signed colIt = 0; colIt < dim; colIt++ )
925  {
926  uAndV[(rowIt*3)+colIt] = rotMatU[( rowIt * 3 ) + colIt].real();
927  }
928  }
929 
930  //================================================ Save V
931  for ( proshade_signed rowIt = 0; rowIt < dim; rowIt++ )
932  {
933  for ( proshade_signed colIt = 0; colIt < dim; colIt++ )
934  {
935  uAndV[(rowIt*3)+colIt+9] = rotMatV[( rowIt * 3 ) + colIt].real();
936  }
937  }
938 
939  //================================================ Release the rest of the memory
940  delete[] rotMatU;
941  delete[] rotMatV;
942 
943  //================================================ Done
944  return ;
945 
946 }
947 
961 void ProSHADE_internal_maths::getEulerZXZFromSOFTPosition ( proshade_signed band, proshade_signed x, proshade_signed y, proshade_signed z, proshade_double* eulerAlpha, proshade_double* eulerBeta, proshade_double* eulerGamma )
962 {
963  //================================================ Convert index to Euler angles
964  *eulerGamma = ( M_PI * y / static_cast<proshade_double> ( 1.0 * band ) );
965  *eulerBeta = ( M_PI * x / static_cast<proshade_double> ( 2.0 * band ) );
966  *eulerAlpha = ( M_PI * z / static_cast<proshade_double> ( 1.0 * band ) );
967 
968  //================================================ Done
969  return ;
970 
971 }
972 
986 void ProSHADE_internal_maths::getSOFTPositionFromEulerZXZ ( proshade_signed band, proshade_double eulerAlpha, proshade_double eulerBeta, proshade_double eulerGamma, proshade_double* x, proshade_double* y, proshade_double* z )
987 {
988  //================================================ Convert Euler angles to indices
989  *x = static_cast<proshade_double> ( ( eulerBeta * static_cast<proshade_double> ( 2.0 * band ) ) / M_PI );
990  *y = static_cast<proshade_double> ( ( eulerGamma * static_cast<proshade_double> ( band ) ) / M_PI );
991  *z = static_cast<proshade_double> ( ( eulerAlpha * static_cast<proshade_double> ( band ) ) / M_PI );
992 
993  //================================================ Done
994  return ;
995 
996 }
997 
1005 void ProSHADE_internal_maths::getRotationMatrixFromEulerZXZAngles ( proshade_double eulerAlpha, proshade_double eulerBeta, proshade_double eulerGamma, proshade_double* matrix )
1006 {
1007  //================================================ First row
1008  matrix[0] = cos ( eulerAlpha ) * cos ( eulerBeta ) * cos ( eulerGamma ) - sin ( eulerAlpha ) * sin ( eulerGamma );
1009  matrix[1] = sin ( eulerAlpha ) * cos ( eulerBeta ) * cos ( eulerGamma ) + cos ( eulerAlpha ) * sin ( eulerGamma );
1010  matrix[2] = -sin ( eulerBeta ) * cos ( eulerGamma );
1011 
1012  //================================================ Second row
1013  matrix[3] = -cos ( eulerAlpha ) * cos ( eulerBeta ) * sin ( eulerGamma ) - sin ( eulerAlpha ) * cos ( eulerGamma );
1014  matrix[4] = -sin ( eulerAlpha ) * cos ( eulerBeta ) * sin ( eulerGamma ) + cos ( eulerAlpha ) * cos ( eulerGamma );
1015  matrix[5] = sin ( eulerBeta ) * sin ( eulerGamma );
1016 
1017  //================================================ Third row
1018  matrix[6] = cos ( eulerAlpha ) * sin ( eulerBeta );
1019  matrix[7] = sin ( eulerAlpha ) * sin ( eulerBeta );
1020  matrix[8] = cos ( eulerBeta );
1021 
1022  //================================================ Done
1023  return ;
1024 
1025 }
1026 
1039  void ProSHADE_internal_maths::getAxisAngleFromRotationMatrix ( proshade_double* rotMat, proshade_double* x, proshade_double* y, proshade_double* z, proshade_double* ang )
1040 {
1041  //================================================ Initialise
1042  proshade_double singAtPiCheck = 0.01;
1043  proshade_double singAtIdentity = 0.05;
1044 
1045  //================================================ Check input for singularities
1046  if ( ( std::abs ( rotMat[1] - rotMat[3] ) < singAtPiCheck ) &&
1047  ( std::abs ( rotMat[2] - rotMat[6] ) < singAtPiCheck ) &&
1048  ( std::abs ( rotMat[5] - rotMat[7] ) < singAtPiCheck ) )
1049  {
1050  //============================================ Singularity in input! Check for identity matrix
1051  if ( ( std::abs ( rotMat[1] + rotMat[3] ) < singAtIdentity ) &&
1052  ( std::abs ( rotMat[2] + rotMat[6] ) < singAtIdentity ) &&
1053  ( std::abs ( rotMat[5] + rotMat[7] ) < singAtIdentity ) &&
1054  ( std::abs ( rotMat[0] + rotMat[4] + rotMat[8] - 3.0 ) < singAtIdentity ) )
1055  {
1056  //======================================== Identity matrix. Return 0 angle.
1057  *x = 1.0;
1058  *y = 0.0;
1059  *z = 0.0;
1060  *ang = 0.0;
1061 
1062  //======================================== Done
1063  return ;
1064  }
1065 
1066  //============================================ If we got here, this is the 180deg (pi rad) singularity. Find which axis should the rotation be done along
1067  *ang = M_PI;
1068 
1069  proshade_double xx = ( rotMat[0] + 1.0 ) / 2.0;
1070  proshade_double yy = ( rotMat[4] + 1.0 ) / 2.0;
1071  proshade_double zz = ( rotMat[8] + 1.0 ) / 2.0;
1072  proshade_double xy = ( rotMat[1] + rotMat[3] ) / 4.0;
1073  proshade_double xz = ( rotMat[2] + rotMat[6] ) / 4.0;
1074  proshade_double yz = ( rotMat[5] + rotMat[7] ) / 4.0;
1075 
1076  if ( ( xx > yy ) && ( xx > zz ) ) // XX is the largest diagonal
1077  {
1078  if ( xx < singAtPiCheck ) // and is still 0
1079  {
1080  *x = 0.0;
1081  *y = 1.0 / sqrt(2);
1082  *z = 1.0 / sqrt(2);
1083  }
1084  else
1085  {
1086  *x = sqrt ( xx );
1087  *y = xy / sqrt ( xx );
1088  *z = xz / sqrt ( xx );
1089  }
1090  }
1091 
1092  else if ( yy > zz ) // YY is the largest diagonal
1093  {
1094  if ( yy < singAtPiCheck ) // and is still 0
1095  {
1096  *x = 1.0 / sqrt(2);
1097  *y = 0.0;
1098  *z = 1.0 / sqrt(2);
1099  }
1100  else
1101  {
1102  *y = sqrt ( yy );
1103  *x = xy / sqrt ( yy );
1104  *z = yz / sqrt ( yy );
1105  }
1106  }
1107 
1108  else // ZZ is the largest diagonal
1109  {
1110  if ( zz < singAtPiCheck ) // and is still 0
1111  {
1112  *x = 1.0 / sqrt(2);
1113  *y = 1.0 / sqrt(2);
1114  *z = 0.0;
1115  }
1116  else
1117  {
1118  *z = sqrt ( zz );
1119  *x = xz / sqrt ( zz );
1120  *y = yz / sqrt ( zz );
1121  }
1122  }
1123 
1124  //============================================ Make sure largest axis is positive and so is the angle
1125  if ( ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *x ) ) && ( *x < 0.0 ) ) ||
1126  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *y ) ) && ( *y < 0.0 ) ) ||
1127  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *z ) ) && ( *z < 0.0 ) ) )
1128  {
1129  *x *= -1.0;
1130  *y *= -1.0;
1131  *z *= -1.0;
1132  *ang *= -1.0;
1133  }
1134  if ( *ang < 0.0 ) { *ang = ( 2.0 * M_PI ) + *ang; }
1135 
1136  //============================================ Done
1137  return ;
1138  }
1139 
1140  //================================================ No singularities! Now get angle
1141  *ang = std::acos ( ( std::max ( -1.0, std::min ( 3.0, rotMat[0] + rotMat[4] + rotMat[8] ) ) - 1.0 ) / 2.0 );
1142 
1143  //================================================ Init return values
1144  *x = 1.0;
1145  *y = 0.0;
1146  *z = 0.0;
1147 
1148  //================================================ Is angle 0? This should not happen, but will
1149  if ( std::abs ( *ang ) < singAtPiCheck )
1150  {
1151  *ang = 0.0;
1152  return ;
1153  }
1154 
1155  //================================================ Axis
1156  *x = rotMat[7] - rotMat[5];
1157  *y = rotMat[2] - rotMat[6];
1158  *z = rotMat[3] - rotMat[1];
1159 
1160  proshade_double normFactor = std::sqrt ( pow ( *x, 2.0 ) + pow ( *y, 2.0 ) + pow ( *z, 2.0 ) );
1161  *x /= normFactor;
1162  *y /= normFactor;
1163  *z /= normFactor;
1164 
1165  //================================================ Make sure largest axis is positive and so is the angle
1166  if ( ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *x ) ) && ( *x < 0.0 ) ) ||
1167  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *y ) ) && ( *y < 0.0 ) ) ||
1168  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *z ) ) && ( *z < 0.0 ) ) )
1169  {
1170  *x *= -1.0;
1171  *y *= -1.0;
1172  *z *= -1.0;
1173  *ang *= -1.0;
1174  }
1175  if ( *ang < 0.0 ) { *ang = ( 2.0 * M_PI ) + *ang; }
1176 
1177  //================================================ Done
1178  return ;
1179 
1180 }
1181 
1194  void ProSHADE_internal_maths::getAxisAngleFromRotationMatrix ( std::vector< proshade_double >* rotMat, proshade_double* x, proshade_double* y, proshade_double* z, proshade_double* ang )
1195 {
1196  //================================================ Initialise
1197  proshade_double singAtPiCheck = 0.01;
1198  proshade_double singAtIdentity = 0.05;
1199 
1200  //================================================ Check input for singularities
1201  if ( ( std::abs ( rotMat->at(1) - rotMat->at(3) ) < singAtPiCheck ) &&
1202  ( std::abs ( rotMat->at(2) - rotMat->at(6) ) < singAtPiCheck ) &&
1203  ( std::abs ( rotMat->at(5) - rotMat->at(7) ) < singAtPiCheck ) )
1204  {
1205  //============================================ Singularity in input! Check for identity matrix
1206  if ( ( std::abs ( rotMat->at(1) + rotMat->at(3) ) < singAtIdentity ) &&
1207  ( std::abs ( rotMat->at(2) + rotMat->at(6) ) < singAtIdentity ) &&
1208  ( std::abs ( rotMat->at(5) + rotMat->at(7) ) < singAtIdentity ) &&
1209  ( std::abs ( rotMat->at(0) + rotMat->at(4) + rotMat->at(8) - 3.0 ) < singAtIdentity ) )
1210  {
1211  //======================================== Identity matrix. Return 0 angle.
1212  *x = 1.0;
1213  *y = 0.0;
1214  *z = 0.0;
1215  *ang = 0.0;
1216 
1217  //======================================== Done
1218  return ;
1219  }
1220 
1221  //============================================ If we got here, this is the 180deg (pi rad) singularity. Find which axis should the rotation be done along
1222  *ang = M_PI;
1223 
1224  proshade_double xx = ( rotMat->at(0) + 1.0 ) / 2.0;
1225  proshade_double yy = ( rotMat->at(4) + 1.0 ) / 2.0;
1226  proshade_double zz = ( rotMat->at(8) + 1.0 ) / 2.0;
1227  proshade_double xy = ( rotMat->at(1) + rotMat->at(3) ) / 4.0;
1228  proshade_double xz = ( rotMat->at(2) + rotMat->at(6) ) / 4.0;
1229  proshade_double yz = ( rotMat->at(5) + rotMat->at(7) ) / 4.0;
1230 
1231  if ( ( xx > yy ) && ( xx > zz ) ) // XX is the largest diagonal
1232  {
1233  if ( xx < singAtPiCheck ) // and is still 0
1234  {
1235  *x = 0.0;
1236  *y = 1.0 / sqrt(2);
1237  *z = 1.0 / sqrt(2);
1238  }
1239  else
1240  {
1241  *x = sqrt ( xx );
1242  *y = xy / sqrt ( xx );
1243  *z = xz / sqrt ( xx );
1244  }
1245  }
1246 
1247  else if ( yy > zz ) // YY is the largest diagonal
1248  {
1249  if ( yy < singAtPiCheck ) // and is still 0
1250  {
1251  *x = 1.0 / sqrt(2);
1252  *y = 0.0;
1253  *z = 1.0 / sqrt(2);
1254  }
1255  else
1256  {
1257  *y = sqrt ( yy );
1258  *x = xy / sqrt ( yy );
1259  *z = yz / sqrt ( yy );
1260  }
1261  }
1262 
1263  else // ZZ is the largest diagonal
1264  {
1265  if ( zz < singAtPiCheck ) // and is still 0
1266  {
1267  *x = 1.0 / sqrt(2);
1268  *y = 1.0 / sqrt(2);
1269  *z = 0.0;
1270  }
1271  else
1272  {
1273  *z = sqrt ( zz );
1274  *x = xz / sqrt ( zz );
1275  *y = yz / sqrt ( zz );
1276  }
1277  }
1278 
1279  //============================================ Make sure largest axis is positive and so is the angle
1280  if ( ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *x ) ) && ( *x < 0.0 ) ) ||
1281  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *y ) ) && ( *y < 0.0 ) ) ||
1282  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *z ) ) && ( *z < 0.0 ) ) )
1283  {
1284  *x *= -1.0;
1285  *y *= -1.0;
1286  *z *= -1.0;
1287  *ang *= -1.0;
1288  }
1289 
1290  //============================================ Done
1291  return ;
1292  }
1293 
1294  //================================================ No singularities! Now get angle
1295  *ang = std::acos ( ( std::max ( -1.0, std::min ( 3.0, rotMat->at(0) + rotMat->at(4) + rotMat->at(8) ) ) - 1.0 ) / 2.0 );
1296 
1297  //================================================ Init return values
1298  *x = 1.0;
1299  *y = 0.0;
1300  *z = 0.0;
1301 
1302  //================================================ Is angle 0? This should not happen, but will
1303  if ( std::abs ( *ang ) < singAtPiCheck )
1304  {
1305  *ang = 0.0;
1306  return ;
1307  }
1308 
1309  //================================================ Axis
1310  *x = rotMat->at(7) - rotMat->at(5);
1311  *y = rotMat->at(2) - rotMat->at(6);
1312  *z = rotMat->at(3) - rotMat->at(1);
1313 
1314  proshade_double normFactor = std::sqrt ( pow ( *x, 2.0 ) + pow ( *y, 2.0 ) + pow ( *z, 2.0 ) );
1315  *x /= normFactor;
1316  *y /= normFactor;
1317  *z /= normFactor;
1318 
1319  //================================================ Make sure largest axis is positive and so is the angle
1320  if ( ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *x ) ) && ( *x < 0.0 ) ) ||
1321  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *y ) ) && ( *y < 0.0 ) ) ||
1322  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *z ) ) && ( *z < 0.0 ) ) )
1323  {
1324  *x *= -1.0;
1325  *y *= -1.0;
1326  *z *= -1.0;
1327  *ang *= -1.0;
1328  }
1329 
1330  //================================================ Done
1331  return ;
1332 
1333 }
1334 
1343 void ProSHADE_internal_maths::getRotationMatrixFromAngleAxis ( proshade_double* rotMat, proshade_double x, proshade_double y, proshade_double z, proshade_double ang )
1344 {
1345  //================================================ If angle is 0 or infinity (anything divided by 0), return identity matrix
1346  if ( ( ang == 0.0 ) || ( std::isinf ( ang ) ) )
1347  {
1348  //============================================ Create identity
1349  for ( proshade_unsign i = 0; i < 9; i++ ) { rotMat[i] = 0.0; }
1350  rotMat[0] = 1.0;
1351  rotMat[4] = 1.0;
1352  rotMat[8] = 1.0;
1353 
1354  //============================================ Done
1355  return ;
1356  }
1357 
1358  //================================================ Compute the matrix
1359  proshade_double cAng = cos ( ang );
1360  proshade_double sAng = sin ( ang );
1361  proshade_double tAng = 1.0 - cAng;
1362 
1363  rotMat[0] = cAng + x * x * tAng;
1364  rotMat[4] = cAng + y * y * tAng;
1365  rotMat[8] = cAng + z * z * tAng;
1366 
1367  proshade_double tmp1 = x * y * tAng;
1368  proshade_double tmp2 = z * sAng;
1369  rotMat[3] = tmp1 + tmp2;
1370  rotMat[1] = tmp1 - tmp2;
1371 
1372  tmp1 = x * z * tAng;
1373  tmp2 = y * sAng;
1374  rotMat[6] = tmp1 - tmp2;
1375  rotMat[2] = tmp1 + tmp2;
1376 
1377  tmp1 = y * z * tAng;
1378  tmp2 = x * sAng;
1379  rotMat[7] = tmp1 + tmp2;
1380  rotMat[5] = tmp1 - tmp2;
1381 
1382  //================================================ Done
1383  return ;
1384 
1385 }
1386 
1394 void ProSHADE_internal_maths::getEulerZXZFromRotMatrix ( proshade_double* rotMat, proshade_double* eA, proshade_double* eB, proshade_double* eG )
1395 {
1396  //================================================ Get ZXZ Euler from matrix
1397  *eA = atan2 ( rotMat[7], rotMat[6] );
1398  *eB = acos ( rotMat[8] );
1399  *eG = atan2 ( rotMat[5], -rotMat[2] );
1400 
1401  //================================================ Solve undefined 0,0 inputs (i.e. identity matrix)
1402  proshade_double errLimit = 0.001;
1403  if ( ( ( rotMat[7] < errLimit ) && ( rotMat[7] > -errLimit ) ) && ( ( rotMat[6] < errLimit ) && ( rotMat[6] > -errLimit ) ) )
1404  {
1405  //============================================ atan2 (0,0) is undefined, we want 0.0 here
1406  *eA = 0.0;
1407  }
1408 
1409  if ( ( ( rotMat[5] < errLimit ) && ( rotMat[5] > -errLimit ) ) && ( ( rotMat[2] < errLimit ) && ( rotMat[2] > -errLimit ) ) )
1410  {
1411  //============================================ atan2 (0,0) is undefined, we want 0.0 here
1412  *eG = 0.0;
1413  }
1414 
1415  //================================================ Get the angles to proper range
1416  if ( *eA < 0.0 ) { *eA = 2.0 * M_PI + *eA; }
1417  if ( *eB < 0.0 ) { *eB = M_PI + *eB; }
1418  if ( *eG < 0.0 ) { *eG = 2.0 * M_PI + *eG; }
1419 
1420  //================================================ Done
1421  return ;
1422 
1423 }
1424 
1437 void ProSHADE_internal_maths::getEulerZXZFromAngleAxis ( proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double axAng, proshade_double* eA, proshade_double* eB, proshade_double* eG, proshade_unsign angDim )
1438 {
1439  //================================================ If angle is 0 or infinity (anything divided by 0), return no rotation
1440  if ( ( axAng == 0.0 ) || ( std::isinf ( axAng ) ) )
1441  {
1442  //============================================ Return 0 ; 0 ; 0 for no angle
1443  *eA = 0.0;
1444  *eB = 0.0;
1445  *eG = 0.0;
1446 
1447  //============================================ Done
1448  return ;
1449  }
1450 
1451  //================================================ Compute required rotation matrix elements
1452  proshade_double cAng = std::cos ( axAng );
1453  proshade_double sAng = std::sin ( axAng );
1454  proshade_double tAng = 1.0 - cAng;
1455 
1456  proshade_double element22 = cAng + axZ * axZ * tAng;
1457 
1458  proshade_double tmp1 = axX * axZ * tAng;
1459  proshade_double tmp2 = axY * sAng;
1460  proshade_double element20 = tmp1 - tmp2;
1461  proshade_double element02 = tmp1 + tmp2;
1462 
1463  tmp1 = axY * axZ * tAng;
1464  tmp2 = axX * sAng;
1465  proshade_double element21 = tmp1 + tmp2;
1466  proshade_double element12 = tmp1 - tmp2;
1467 
1468  //================================================ Convert to Eulers
1469  if ( std::abs( element22 ) <= 0.99999 )
1470  {
1471  //============================================ This case occurs when there is no singularity in the rotation matrix (i.e. it does not have 0 or 180 degrees angle)
1472  *eA = std::atan2 ( element21, element20 );
1473  *eB = std::acos ( element22 );
1474  *eG = std::atan2 ( element12, -element02 );
1475  }
1476  else
1477  {
1478  //============================================ Compute some extra rotation matrix elements
1479  tmp1 = axX * axY * tAng;
1480  tmp2 = axZ * sAng;
1481  proshade_double element10 = tmp1 + tmp2;
1482  proshade_double element00 = cAng + axX * axX * tAng;
1483 
1484  //============================================ This case occurs when there is either 0 or 180 degrees rotation angle in the rotation matrix and therefore when beta is zero.
1485  if ( element22 >= 0.99999 )
1486  {
1487  //======================================== In this case, beta = 0 and alpha and gamma are only defined in terms of their sum. So we arbitrarily set gamma to 0 and solve alpha.
1488  *eA = std::atan2 ( element10, element00 );
1489  *eB = 0.0;
1490  *eG = 0.0;
1491  }
1492  if ( element22 <= -0.99999 )
1493  {
1494  //======================================== In this case, beta = 0 and alpha and gamma are only defined in terms of their difference. So we arbitrarily set gamma to 0 and solve alpha.
1495  *eA = std::atan2 ( element10, element00 );
1496  *eB = M_PI / 2.0;
1497  *eG = 0.0;
1498  }
1499  }
1500 
1501  //================================================ Get the angles to proper range
1502  if ( *eA < 0.0 ) { *eA = 2.0 * M_PI + *eA; }
1503  if ( *eB < 0.0 ) { *eB = M_PI + *eB; }
1504  if ( *eG < 0.0 ) { *eG = 2.0 * M_PI + *eG; }
1505 
1506  //================================================ Done
1507  return ;
1508 
1509 }
1510 
1526 void ProSHADE_internal_maths::getEulerZXZFromAngleAxisFullSearch ( proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double axAng, proshade_double* eA, proshade_double* eB, proshade_double* eG, proshade_signed angDim )
1527 {
1528  //================================================ Initialise variables
1529  proshade_double bestDist = 999.9;
1530  proshade_double eAHlp, eBHlp, eGHlp, axXHlp, axYHlp, axZHlp, axAngHlp, axDist;
1531 
1532  //================================================ Allocate memory
1533  proshade_double* rMat = new proshade_double[9];
1534  ProSHADE_internal_misc::checkMemoryAllocation ( rMat, __FILE__, __LINE__, __func__ );
1535 
1536  //================================================ For each rotation function index (i.e. existing Euler angles ZXZ combination)
1537  for ( proshade_signed xIt = 0; xIt < angDim; xIt++ )
1538  {
1539  for ( proshade_signed yIt = 0; yIt < angDim; yIt++ )
1540  {
1541  for ( proshade_signed zIt = 0; zIt < angDim; zIt++ )
1542  {
1543  //==================================== Speed up
1544  if ( bestDist < 0.001 ) { break; }
1545 
1546  //==================================== Get Euler ZXZ from the indices
1547  getEulerZXZFromSOFTPosition ( angDim/2, xIt, yIt, zIt, &eAHlp, &eBHlp, &eGHlp );
1548  getRotationMatrixFromEulerZXZAngles ( eAHlp, eBHlp, eGHlp, rMat );
1549  getAxisAngleFromRotationMatrix ( rMat, &axXHlp, &axYHlp, &axZHlp, &axAngHlp );
1550 
1551  //==================================== If angle is larger than 180 degrees
1552  if ( axAng > M_PI )
1553  {
1554  axAng = ( 2.0 * M_PI ) - axAng;
1555  axAng *= -1.0;
1556  }
1557 
1558  //==================================== Make sure vector direction is the same
1559  if ( ( ( std::max( std::abs( axXHlp ), std::max( std::abs( axYHlp ), std::abs( axZHlp ) ) ) == std::abs( axXHlp ) ) && ( axXHlp < 0.0 ) ) ||
1560  ( ( std::max( std::abs( axXHlp ), std::max( std::abs( axYHlp ), std::abs( axZHlp ) ) ) == std::abs( axYHlp ) ) && ( axYHlp < 0.0 ) ) ||
1561  ( ( std::max( std::abs( axXHlp ), std::max( std::abs( axYHlp ), std::abs( axZHlp ) ) ) == std::abs( axZHlp ) ) && ( axZHlp < 0.0 ) ) )
1562  {
1563  axXHlp *= -1.0;
1564  axYHlp *= -1.0;
1565  axZHlp *= -1.0;
1566  axAngHlp *= -1.0;
1567  }
1568 
1569  if ( ( ( std::max( std::abs( axX ), std::max( std::abs( axY ), std::abs( axZ ) ) ) == std::abs( axX ) ) && ( axX < 0.0 ) ) ||
1570  ( ( std::max( std::abs( axX ), std::max( std::abs( axY ), std::abs( axZ ) ) ) == std::abs( axY ) ) && ( axY < 0.0 ) ) ||
1571  ( ( std::max( std::abs( axX ), std::max( std::abs( axY ), std::abs( axZ ) ) ) == std::abs( axZ ) ) && ( axZ < 0.0 ) ) )
1572  {
1573  axX *= -1.0;
1574  axY *= -1.0;
1575  axZ *= -1.0;
1576  axAng *= -1.0;
1577  }
1578 
1579  //==================================== Compute distance to the requested angle-axis values
1580  axDist = std::abs( axAng - axAngHlp ) + ( 1.0 - std::abs ( ( ( axX * axXHlp ) + ( axY * axYHlp ) + ( axZ * axZHlp ) ) /
1581  ( sqrt( pow( axX, 2.0 ) + pow( axY, 2.0 ) + pow( axZ, 2.0 ) ) * sqrt( pow( axXHlp, 2.0 ) + pow( axYHlp, 2.0 ) + pow( axZHlp, 2.0 ) ) ) ) );
1582 
1583  //==================================== Is this point an improvement
1584  if ( std::abs ( axDist ) < bestDist )
1585  {
1586  //================================ If so, note it
1587  bestDist = std::abs ( axDist );
1588  *eA = eAHlp;
1589  *eB = eBHlp;
1590  *eG = eGHlp;
1591  }
1592  }
1593  }
1594  }
1595 
1596  //================================================ Release memory
1597  delete[] rMat;
1598 
1599  //================================================ Done
1600  return ;
1601 
1602 }
1603 
1613 void ProSHADE_internal_maths::multiplyTwoSquareMatrices ( proshade_double* A, proshade_double* B, proshade_double* res, proshade_unsign dim )
1614 {
1615  //================================================ Set res to 0.0s
1616  for ( proshade_unsign iter = 0; iter < 9; iter++ ) { res[iter] = 0.0; }
1617 
1618  //================================================ Compute the matrix multiplication
1619  for ( proshade_unsign row = 0; row < dim; row++ )
1620  {
1621  for ( proshade_unsign col = 0; col < dim; col++ )
1622  {
1623  for ( proshade_unsign inner = 0; inner < dim; inner++ )
1624  {
1625  res[(row*dim)+col] += A[(inner*dim)+row] * B[(col*dim)+inner];
1626  }
1627  }
1628  }
1629 
1630  //================================================ Done
1631  return ;
1632 
1633 }
1634 
1639 std::vector < proshade_signed > ProSHADE_internal_maths::primeFactorsDecomp ( proshade_signed number )
1640 {
1641  //================================================ Initialise variables
1642  std::vector < proshade_signed > ret;
1643 
1644  //================================================ Deal with negative numbers
1645  bool changeSign = false;
1646  if ( number < 0 ) { changeSign = true; number = -number; }
1647 
1648  //================================================ Deal with zero and one
1649  if ( number == 0 ) { ProSHADE_internal_misc::addToSignedVector ( &ret, 0 ); return ( ret ); }
1650  if ( number == 1 ) { ProSHADE_internal_misc::addToSignedVector ( &ret, 1 ); return ( ret ); }
1651 
1652  //================================================ Divide by 2 as long as you can
1653  while ( number % 2 == 0 )
1654  {
1656  number = number / 2;
1657  }
1658 
1659  //================================================ Check all odd numbers up to the square root
1660  for ( proshade_unsign posDiv = 3; posDiv <= sqrt ( number ); posDiv += 2)
1661  {
1662  // If posDiv is a divisor of the number, save the result
1663  while ( number % posDiv == 0 )
1664  {
1666  number = number / posDiv;
1667  }
1668  }
1669 
1670  //================================================ If the number was a large prime number, save it as it is
1671  if ( number > 2 ) { ProSHADE_internal_misc::addToSignedVector ( &ret, number ); }
1672 
1673  //================================================ Finish dealing with negative numbers
1674  if ( changeSign ) { ret.at(0) = -ret.at(0); }
1675 
1676  //================================================ Done
1677  return ( ret );
1678 
1679 }
1680 
1688 proshade_double ProSHADE_internal_maths::normalDistributionValue ( proshade_double mean, proshade_double standardDev, proshade_double value )
1689 {
1690  //================================================ Compute and return
1691  return ( ( 1.0 / sqrt ( 2.0 * M_PI * pow(standardDev,2.0) ) ) * std::exp ( - pow( value - mean, 2.0 ) / 2.0 * pow(standardDev,2.0) ) );
1692 
1693 }
1694 
1705 proshade_double ProSHADE_internal_maths::computeDotProduct ( proshade_double* x1, proshade_double* y1, proshade_double* z1, proshade_double* x2, proshade_double* y2, proshade_double* z2 )
1706 {
1707  //================================================ Compute and return
1708  return ( (*x1 * *x2) + (*y1 * *y2) + (*z1 * *z2) );
1709 }
1710 
1721 proshade_double ProSHADE_internal_maths::computeDotProduct ( proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2 )
1722 {
1723  //================================================ Compute and return
1724  return ( (x1 * x2) + (y1 * y2) + (z1 * z2) );
1725 }
1726 
1737 proshade_double* ProSHADE_internal_maths::computeCrossProduct ( proshade_double* x1, proshade_double* y1, proshade_double* z1, proshade_double* x2, proshade_double* y2, proshade_double* z2 )
1738 {
1739  //================================================ Allocate memory
1740  proshade_double* crossProd = new proshade_double[3];
1741  ProSHADE_internal_misc::checkMemoryAllocation ( crossProd, __FILE__, __LINE__, __func__ );
1742 
1743  //================================================ Compute
1744  crossProd[0] = ( (*y1) * (*z2) ) - ( (*z1) * (*y2) );
1745  crossProd[1] = ( (*z1) * (*x2) ) - ( (*x1) * (*z2) );
1746  crossProd[2] = ( (*x1) * (*y2) ) - ( (*y1) * (*x2) );
1747 
1748  //================================================ Done
1749  return ( crossProd );
1750 
1751 }
1752 
1759 proshade_double* ProSHADE_internal_maths::compute3x3MatrixMultiplication ( proshade_double* mat1, proshade_double* mat2 )
1760 {
1761  //================================================ Allocate memory
1762  proshade_double* ret = new proshade_double[9];
1763  ProSHADE_internal_misc::checkMemoryAllocation ( ret, __FILE__, __LINE__, __func__ );
1764 
1765  //================================================ Multiply
1766  ret[0] = ( mat1[0] * mat2[0] ) + ( mat1[1] * mat2[3] ) + ( mat1[2] * mat2[6] );
1767  ret[1] = ( mat1[0] * mat2[1] ) + ( mat1[1] * mat2[4] ) + ( mat1[2] * mat2[7] );
1768  ret[2] = ( mat1[0] * mat2[2] ) + ( mat1[1] * mat2[5] ) + ( mat1[2] * mat2[8] );
1769  ret[3] = ( mat1[3] * mat2[0] ) + ( mat1[4] * mat2[3] ) + ( mat1[5] * mat2[6] );
1770  ret[4] = ( mat1[3] * mat2[1] ) + ( mat1[4] * mat2[4] ) + ( mat1[5] * mat2[7] );
1771  ret[5] = ( mat1[3] * mat2[2] ) + ( mat1[4] * mat2[5] ) + ( mat1[5] * mat2[8] );
1772  ret[6] = ( mat1[6] * mat2[0] ) + ( mat1[7] * mat2[3] ) + ( mat1[8] * mat2[6] );
1773  ret[7] = ( mat1[6] * mat2[1] ) + ( mat1[7] * mat2[4] ) + ( mat1[8] * mat2[7] );
1774  ret[8] = ( mat1[6] * mat2[2] ) + ( mat1[7] * mat2[5] ) + ( mat1[8] * mat2[8] );
1775 
1776  //================================================ Done
1777  return ( ret );
1778 
1779 }
1780 
1789 proshade_double* ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( proshade_double* mat, proshade_double x, proshade_double y, proshade_double z )
1790 {
1791  //================================================ Allocate memory
1792  proshade_double* ret = new proshade_double[3];
1793  ProSHADE_internal_misc::checkMemoryAllocation ( ret, __FILE__, __LINE__, __func__ );
1794 
1795  //================================================ Compute the multiplication
1796  ret[0] = ( x * mat[0] ) + ( y * mat[1] ) + ( z * mat[2] );
1797  ret[1] = ( x * mat[3] ) + ( y * mat[4] ) + ( z * mat[5] );
1798  ret[2] = ( x * mat[6] ) + ( y * mat[7] ) + ( z * mat[8] );
1799 
1800  //================================================ Done
1801  return ( ret );
1802 
1803 }
1804 
1810 proshade_double* ProSHADE_internal_maths::compute3x3MatrixInverse ( proshade_double* mat )
1811 {
1812  //================================================ Allocate memory
1813  proshade_double* inverse = new proshade_double[9];
1814  ProSHADE_internal_misc::checkMemoryAllocation ( inverse, __FILE__, __LINE__, __func__ );
1815 
1816  //================================================ Compute determinant
1817  proshade_double matDet = ( mat[0] * mat[4] * mat[8] ) +
1818  ( mat[1] * mat[5] * mat[6] ) +
1819  ( mat[2] * mat[3] * mat[7] ) -
1820  ( mat[0] * mat[5] * mat[7] ) -
1821  ( mat[1] * mat[3] * mat[8] ) -
1822  ( mat[2] * mat[4] * mat[6] );
1823 
1824  //================================================ Compute inverse matrix
1825  inverse[0] = ( mat[4] * mat[8] - mat[5] * mat[7] ) / matDet;
1826  inverse[1] = ( mat[2] * mat[7] - mat[1] * mat[8] ) / matDet;
1827  inverse[2] = ( mat[1] * mat[5] - mat[2] * mat[4] ) / matDet;
1828  inverse[3] = ( mat[5] * mat[6] - mat[3] * mat[8] ) / matDet;
1829  inverse[4] = ( mat[0] * mat[8] - mat[2] * mat[6] ) / matDet;
1830  inverse[5] = ( mat[2] * mat[3] - mat[0] * mat[5] ) / matDet;
1831  inverse[6] = ( mat[3] * mat[7] - mat[4] * mat[6] ) / matDet;
1832  inverse[7] = ( mat[1] * mat[6] - mat[0] * mat[7] ) / matDet;
1833  inverse[8] = ( mat[0] * mat[4] - mat[1] * mat[3] ) / matDet;
1834 
1835  //================================================ Done
1836  return ( inverse );
1837 }
1838 
1844 {
1845  //================================================ Initialise variables
1846  proshade_double tmp;
1847 
1848  //================================================ Transpose the non-diagonal values
1849  tmp = mat[1];
1850  mat[1] = mat[3];
1851  mat[3] = tmp;
1852 
1853  tmp = mat[2];
1854  mat[2] = mat[6];
1855  mat[6] = tmp;
1856 
1857  tmp = mat[5];
1858  mat[5] = mat[7];
1859  mat[7] = tmp;
1860 
1861  //================================================ Done
1862  return ;
1863 
1864 }
1865 
1883 proshade_double* ProSHADE_internal_maths::findRotMatMatchingVectors ( proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2 )
1884 {
1885  //================================================ Allocate required memory
1886  proshade_double* inPlaneRotation = new proshade_double[9];
1887  proshade_double* basisChangeMat = new proshade_double[9];
1888  ProSHADE_internal_misc::checkMemoryAllocation ( inPlaneRotation, __FILE__, __LINE__, __func__ );
1889  ProSHADE_internal_misc::checkMemoryAllocation ( basisChangeMat, __FILE__, __LINE__, __func__ );
1890 
1891  //================================================ Normalise inputs
1892  proshade_double normF = std::sqrt( std::pow( x1, 2.0 ) + std::pow ( y1, 2.0 ) + std::pow ( z1, 2.0 ) );
1893  x1 /= normF; y1 /= normF; z1 /= normF;
1894 
1895  normF = std::sqrt( std::pow( x2, 2.0 ) + std::pow ( y2, 2.0 ) + std::pow ( z2, 2.0 ) );
1896  x2 /= normF; y2 /= normF; z2 /= normF;
1897 
1898  //================================================ Compute cross product's magnitude
1899  proshade_double* crossProd = ProSHADE_internal_maths::computeCrossProduct( &x1, &y1, &z1, &x2, &y2, &z2 );
1900  proshade_double crossProdMag = std::sqrt( std::pow( crossProd[0], 2.0 ) + std::pow ( crossProd[1], 2.0 ) + std::pow ( crossProd[2], 2.0 ) );
1901  delete[] crossProd;
1902 
1903  //================================================ Compute dot product
1904  proshade_double dotProd = ProSHADE_internal_maths::computeDotProduct ( &x1, &y1, &z1, &x2, &y2, &z2 );
1905 
1906  //================================================ Construct the in-plane rotation matrix
1907  inPlaneRotation[0] = dotProd; inPlaneRotation[1] = -crossProdMag; inPlaneRotation[2] = 0.0;
1908  inPlaneRotation[3] = crossProdMag; inPlaneRotation[4] = dotProd; inPlaneRotation[5] = 0.0;
1909  inPlaneRotation[6] = 0.0; inPlaneRotation[7] = 0.0; inPlaneRotation[8] = 1.0;
1910 
1911  //================================================ Construct change of basis matrix
1912  crossProd = ProSHADE_internal_maths::computeCrossProduct( &x2, &y2, &z2, &x1, &y1, &z1 );
1913  normF = std::sqrt ( std::pow ( x2 - ( dotProd * x1 ), 2.0 ) + std::pow ( y2 - ( dotProd * y1 ), 2.0 ) + std::pow ( z2 - ( dotProd * z1 ), 2.0 ) );
1914 
1915  basisChangeMat[0] = x1; basisChangeMat[1] = ( x2 - ( dotProd * x1 ) ) / normF; basisChangeMat[2] = crossProd[0];
1916  basisChangeMat[3] = y1; basisChangeMat[4] = ( y2 - ( dotProd * y1 ) ) / normF; basisChangeMat[5] = crossProd[1];
1917  basisChangeMat[6] = z1; basisChangeMat[7] = ( z2 - ( dotProd * z1 ) ) / normF; basisChangeMat[8] = crossProd[2];
1918 
1919  //================================================ Invert the change of basis matrix
1920  proshade_double* basisChangeMatInverse = ProSHADE_internal_maths::compute3x3MatrixInverse ( basisChangeMat );
1921 
1922  //================================================ Multiply inverse of change of basis matrix with the in plane rotation matrix, then multiply the result with the inverse
1923  proshade_double* tmpMat = ProSHADE_internal_maths::compute3x3MatrixMultiplication ( basisChangeMat, inPlaneRotation );
1924  proshade_double* rotMat = ProSHADE_internal_maths::compute3x3MatrixMultiplication ( tmpMat, basisChangeMatInverse );
1925 
1926  //================================================ Release memory
1927  delete[] crossProd;
1928  delete[] inPlaneRotation;
1929  delete[] basisChangeMat;
1930  delete[] basisChangeMatInverse;
1931  delete[] tmpMat;
1932 
1933  //================================================ Done
1934  return ( rotMat );
1935 
1936 }
1937 
1959 std::vector < proshade_double > ProSHADE_internal_maths::findVectorFromTwoVAndTwoD ( proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2, proshade_double dot1, proshade_double dot2 )
1960 {
1961  //================================================ Initialise variables
1962  std::vector < proshade_double > ret;
1963 
1964  //================================================ Solution
1965  proshade_double solX = ( -sqrt ( pow ( 2.0 * x1 * y1 * dot2 * y2 + 2.0 * x1 * z1 * dot2 * z2 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( y1, 2.0 ) * dot2 * x2 + 2.0 * y1 * dot1 * x2 * y2 - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
1966  4.0 * ( pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * x1 * y1 * x2 * y2 - 2.0 * x1 * z1 * x2 * z2 + pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) ) *
1967  ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 - 2.0 * y1 * dot1 * dot2 * y2 + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) - 2.0 * z1 * dot1 * dot2 * z2 + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) ) ) -
1968  2.0 * x1 * y1 * dot2 * y2 - 2.0 * x1 * z1 * dot2 * z2 + 2.0 * x1 * dot1 * pow ( y2, 2.0 ) + 2.0 * x1 * dot1 * pow ( z2, 2.0 ) + 2.0 * pow ( y1, 2.0 ) * dot2 * x2 - 2.0 * y1 * dot1 * x2 * y2 + 2.0 * pow ( z1, 2.0 ) * dot2 * x2 - 2.0 * z1 * dot1 * x2 * z2 ) /
1969  ( 2.0 * ( pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * x1 * y1 * x2 * y2 - 2.0 * x1 * z1 * x2 * z2 + pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) ) );
1970  proshade_double solY = ( ( dot2 * pow ( x2, 2.0 ) * pow ( z1, 3.0 ) ) /
1971  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1972  ( dot1 * pow ( x2, 2.0 ) * z2 * pow ( z1, 2.0 ) ) /
1973  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1974  ( 2.0 * x1 * dot2 * x2 * z2 * pow ( z1, 2.0 ) ) /
1975  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) - dot2 * z1 -
1976  ( x2 * sqrt ( pow ( - 2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
1977  4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
1978  ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) * z1 ) /
1979  ( 2.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
1980  ( pow ( y1, 2.0 ) * dot2 * pow ( x2, 2.0 ) * z1 ) /
1981  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
1982  ( x1 * dot1 * x2 * pow ( y2, 2.0 ) * z1 ) /
1983  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
1984  ( pow ( x1, 2.0 ) * dot2 * pow ( z2, 2.0 ) * z1 ) /
1985  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
1986  ( 2.0 * x1 * dot1 * x2 * pow ( z2, 2.0 ) * z1 ) /
1987  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1988  ( y1 * dot1 * pow ( x2, 2.0 ) * y2 * z1 ) /
1989  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1990  ( x1 * y1 * dot2 * x2 * y2 * z1 ) /
1991  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) + dot1 * z2 +
1992  ( x1 * z2 * sqrt ( pow ( -2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
1993  4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
1994  ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) ) /
1995  ( 2.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
1996  ( pow ( x1, 2.0 ) * dot1 * pow ( z2, 3.0 ) ) /
1997  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1998  ( pow ( x1, 2.0 ) * dot1 * pow ( y2, 2.0 ) * z2 ) /
1999  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
2000  ( x1 * pow ( y1, 2.0 ) * dot2 * x2 * z2 ) /
2001  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2002  ( pow ( x1, 2.0 ) * y1 * dot2 * y2 * z2 ) /
2003  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2004  ( x1 * y1 * dot1 * x2 * y2 * z2 ) /
2005  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) / ( y1 * z2 - z1 * y2 );
2006  proshade_double solZ = ( - ( dot2 * pow ( x2, 2.0 ) * y2 * pow ( z1, 3.0 ) ) /
2007  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2008  ( dot2 * pow ( x2, 2.0 ) * pow ( z1, 2.0 ) ) /
2009  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2010  ( dot1 * pow ( x2, 2.0 ) * y2 * z2 * pow ( z1, 2.0 ) ) /
2011  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2012  ( 2.0 * x1 * dot2 * x2 * y2 * z2 * pow ( z1, 2.0 ) ) /
2013  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2014  ( x2 * y2 * sqrt ( pow ( -2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
2015  4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
2016  ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) * z1 ) /
2017  ( 2.0 * ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2018  ( dot2 * y2 * z1 ) / ( y1 * z2 - z1 * y2 ) +
2019  ( dot1 * pow ( x2, 2.0 ) * z2 * z1 ) /
2020  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2021  ( x1 * dot2 * x2 * z2 * z1 ) /
2022  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
2023  ( x1 * dot1 * x2 * pow ( y2, 3.0 ) * z1 ) /
2024  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2025  ( y1 * dot1 * pow ( x2, 2.0 ) * pow ( y2, 2.0 ) * z1 ) /
2026  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2027  ( x1 * y1 * dot2 * x2 * pow ( y2, 2.0 ) * z1 ) /
2028  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2029  ( pow ( x1, 2.0 ) * dot2 * y2 * pow ( z2, 2.0 ) * z1 ) /
2030  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2031  ( 2.0 * x1 * dot1 * x2 * y2 * pow ( z2, 2.0 ) * z1 ) /
2032  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2033  ( pow ( y1, 2.0 ) * dot2 * pow ( x2, 2.0 ) * y2 * z1 ) /
2034  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) + dot2 +
2035  ( x2 * sqrt ( pow ( - 2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
2036  4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
2037  ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) ) /
2038  ( 2.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2039  ( x1 * y2 * z2 * sqrt ( pow ( - 2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
2040  4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
2041  ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) ) /
2042  ( 2.0 * ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2043  ( dot1 * y2 * z2 ) / ( y1 * z2 - z1 * y2 ) -
2044  ( pow ( y1, 2.0 ) * dot2 * pow ( x2, 2.0 ) ) /
2045  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
2046  ( x1 * dot1 * x2 * pow ( y2, 2.0 ) ) /
2047  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
2048  ( x1 * dot1 * x2 * pow ( z2, 2.0 ) ) /
2049  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2050  ( y1 * dot1 * pow ( x2, 2.0 ) * y2 ) /
2051  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2052  ( x1 * y1 * dot2 * x2 * y2 ) /
2053  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2054  ( pow ( x1, 2.0 ) * dot1 * y2 * pow ( z2, 3.0 ) ) /
2055  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2056  ( pow ( x1, 2.0 ) * dot1 * pow ( y2, 3.0 ) * z2 ) /
2057  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2058  ( pow ( x1, 2.0 ) * y1 * dot2 * pow ( y2, 2.0 ) * z2 ) /
2059  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2060  ( x1 * y1 * dot1 * x2 * pow ( y2, 2.0 ) * z2 ) /
2061  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2062  ( x1 * pow ( y1, 2.0 ) * dot2 * x2 * y2 * z2 ) /
2063  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) ) / z2;
2064 
2065  //================================================ Set largest axis element to positive (ProSHADE standard)
2066  if ( ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solX ) ) && ( solX < 0.0 ) ) ||
2067  ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solY ) ) && ( solY < 0.0 ) ) ||
2068  ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solZ ) ) && ( solZ < 0.0 ) ) ) { solX *= -1.0; solY *= -1.0; solZ *= -1.0; }
2069 
2070  //================================================ Save solutions
2074 
2075  //================================================ Done
2076  return ( ret );
2077 
2078 }
2079 
2101 std::vector < proshade_double > ProSHADE_internal_maths::findVectorFromThreeVAndThreeD ( proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2, proshade_double x3, proshade_double y3, proshade_double z3, proshade_double dot1, proshade_double dot2, proshade_double dot3 )
2102 {
2103  //================================================ Initialise variables
2104  std::vector < proshade_double > ret;
2105 
2106  //================================================ Solution
2107  proshade_double solX = - ( y1 * dot2 * z3 - y1 * dot3 * z2 - z1 * dot2 * y3 + z1 * dot3 * y2 + dot1 * y3 * z2 - dot1 * z3 * y2 ) /
2108  ( -x1 * y3 * z2 + x1 * z3 * y2 + y1 * x3 * z2 - y1 * z3 * x2 - z1 * x3 * y2 + z1 * y3 * x2 );
2109  proshade_double solY = - ( x1 * dot2 * z3 - x1 * dot3 * z2 - z1 * dot2 * x3 + z1 * dot3 * x2 + dot1 * x3 * z2 - dot1 * z3 * x2 ) /
2110  ( x1 * y3 * z2 - x1 * z3 * y2 - y1 * x3 * z2 + y1 * z3 * x2 + z1 * x3 * y2 - z1 * y3 * x2 );
2111  proshade_double solZ = - ( x1 * dot2 * y3 - x1 * dot3 * y2 - y1 * dot2 * x3 + y1 * dot3 * x2 + dot1 * x3 * y2 - dot1 * y3 * x2 ) /
2112  ( -x1 * y3 * z2 + x1 * z3 * y2 + y1 * x3 * z2 - y1 * z3 * x2 - z1 * x3 * y2 + z1 * y3 * x2 );
2113 
2114  //================================================ Normalise the axis to magnitude 1
2115  proshade_double normFactor = sqrt ( pow ( solX, 2.0 ) + pow ( solY, 2.0 ) + pow ( solZ, 2.0 ) );
2116  solX /= normFactor;
2117  solY /= normFactor;
2118  solZ /= normFactor;
2119 
2120  //================================================ Set largest axis element to positive (ProSHADE standard)
2121  if ( ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solX ) ) && ( solX < 0.0 ) ) ||
2122  ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solY ) ) && ( solY < 0.0 ) ) ||
2123  ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solZ ) ) && ( solZ < 0.0 ) ) ) { solX *= -1.0; solY *= -1.0; solZ *= -1.0; }
2124 
2125  //================================================ Save solutions
2129 
2130  //================================================ Done
2131  return ( ret );
2132 
2133 }
2134 
2141 std::vector< proshade_double > ProSHADE_internal_maths::multiplyGroupElementMatrices ( std::vector< proshade_double >* el1, std::vector< proshade_double >* el2 )
2142 {
2143  //================================================ Initialise variables
2144  std::vector< proshade_double > ret;
2145 
2146  //================================================ Compute
2147  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(0) * el2->at(0) ) +
2148  ( el1->at(1) * el2->at(3) ) +
2149  ( el1->at(2) * el2->at(6) ) );
2150  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(0) * el2->at(1) ) +
2151  ( el1->at(1) * el2->at(4) ) +
2152  ( el1->at(2) * el2->at(7) ) );
2153  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(0) * el2->at(2) ) +
2154  ( el1->at(1) * el2->at(5) ) +
2155  ( el1->at(2) * el2->at(8) ) );
2156 
2157  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(3) * el2->at(0) ) +
2158  ( el1->at(4) * el2->at(3) ) +
2159  ( el1->at(5) * el2->at(6) ) );
2160  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(3) * el2->at(1) ) +
2161  ( el1->at(4) * el2->at(4) ) +
2162  ( el1->at(5) * el2->at(7) ) );
2163  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(3) * el2->at(2) ) +
2164  ( el1->at(4) * el2->at(5) ) +
2165  ( el1->at(5) * el2->at(8) ) );
2166 
2167  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(6) * el2->at(0) ) +
2168  ( el1->at(7) * el2->at(3) ) +
2169  ( el1->at(8) * el2->at(6) ) );
2170  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(6) * el2->at(1) ) +
2171  ( el1->at(7) * el2->at(4) ) +
2172  ( el1->at(8) * el2->at(7) ) );
2173  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(6) * el2->at(2) ) +
2174  ( el1->at(7) * el2->at(5) ) +
2175  ( el1->at(8) * el2->at(8) ) );
2176 
2177 
2178  //================================================ Done
2179  return ( ret );
2180 
2181 }
2182 
2195 bool ProSHADE_internal_maths::rotationMatrixSimilarity ( std::vector< proshade_double >* mat1, std::vector< proshade_double >* mat2, proshade_double tolerance )
2196 {
2197  //================================================ Initialise variables
2198  bool ret = false;
2199 
2200  //================================================ Compute trace of mat1 * mat2^T
2201  proshade_double trace = ( mat1->at(0) * mat2->at(0) ) + ( mat1->at(1) * mat2->at(1) ) + ( mat1->at(2) * mat2->at(2) );
2202  trace += ( mat1->at(3) * mat2->at(3) ) + ( mat1->at(4) * mat2->at(4) ) + ( mat1->at(5) * mat2->at(5) );
2203  trace += ( mat1->at(6) * mat2->at(6) ) + ( mat1->at(7) * mat2->at(7) ) + ( mat1->at(8) * mat2->at(8) );
2204 
2205  //================================================ Subtract 3 (so that we would have 0 in case of idenity matrix)
2206  trace -= 3.0;
2207 
2208  //================================================ Compare to tolerance
2209  if ( tolerance > std::abs ( trace ) ) { ret = true; }
2210 
2211  //================================================ Done
2212  return ( ret );
2213 
2214 }
2215 
2232 bool ProSHADE_internal_maths::vectorOrientationSimilarity ( proshade_double a1, proshade_double a2, proshade_double a3, proshade_double b1, proshade_double b2, proshade_double b3, proshade_double tolerance )
2233 {
2234  //================================================ Initialise variables
2235  bool ret = false;
2236 
2237  //================================================ Cosine distance
2238  proshade_double cosDist = ( ( a1 * b1 ) + ( a2 * b2 ) + ( a3 * b3 ) ) /
2239  ( sqrt( pow( a1, 2.0 ) + pow( a2, 2.0 ) + pow( a3, 2.0 ) ) *
2240  sqrt( pow( b1, 2.0 ) + pow( b2, 2.0 ) + pow( b3, 2.0 ) ) );
2241 
2242  //================================================ Compare the absolute value of distance to 1.0 - tolerance
2243  if ( std::abs( cosDist ) > ( 1.0 - tolerance ) ) { ret = true; }
2244 
2245  //================================================ Done
2246  return ( ret );
2247 
2248 }
2249 
2266 bool ProSHADE_internal_maths::vectorOrientationSimilaritySameDirection ( proshade_double a1, proshade_double a2, proshade_double a3, proshade_double b1, proshade_double b2, proshade_double b3, proshade_double tolerance )
2267 {
2268  //================================================ Initialise variables
2269  bool ret = false;
2270 
2271  //================================================ Cosine distance
2272  proshade_double cosDist = ( ( a1 * b1 ) + ( a2 * b2 ) + ( a3 * b3 ) ) /
2273  ( sqrt( pow( a1, 2.0 ) + pow( a2, 2.0 ) + pow( a3, 2.0 ) ) *
2274  sqrt( pow( b1, 2.0 ) + pow( b2, 2.0 ) + pow( b3, 2.0 ) ) );
2275 
2276  //================================================ Compare the absolute value of distance to 1.0 - tolerance
2277  if ( cosDist > ( 1.0 - tolerance ) ) { ret = true; }
2278 
2279  //================================================ Done
2280  return ( ret );
2281 
2282 }
2283 
2296 void ProSHADE_internal_maths::optimiseAxisBiCubicInterpolation ( proshade_double* bestLattitude, proshade_double* bestLongitude, proshade_double* bestSum, std::vector<proshade_unsign>* sphereList, std::vector<ProSHADE_internal_spheres::ProSHADE_rotFun_sphere*>* sphereMappedRotFun, proshade_double step )
2297 {
2298  //================================================ Initialise variables
2299  proshade_double lonM, lonP, latM, latP, movSum;
2300  std::vector<proshade_double> latVals ( 3 );
2301  std::vector<proshade_double> lonVals ( 3 );
2302  proshade_double learningRate = 0.1;
2303  proshade_double prevVal = *bestSum;
2304  proshade_double valChange = 999.9;
2305  proshade_double origBestLat = std::round ( *bestLattitude );
2306  proshade_double origBestLon = std::round ( *bestLongitude );
2307  proshade_double tmpVal;
2308 
2309  //================================================ Initialise interpolators in all directions around the point of interest
2310  std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsMinusMinus;
2311  std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsMinusPlus;
2312  std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsPlusMinus;
2313  std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsPlusPlus;
2314  prepareBiCubicInterpolatorsMinusMinus ( std::round ( *bestLattitude ), std::round ( *bestLongitude ), sphereList, &interpolsMinusMinus, sphereMappedRotFun );
2315  prepareBiCubicInterpolatorsMinusPlus ( std::round ( *bestLattitude ), std::round ( *bestLongitude ), sphereList, &interpolsMinusPlus, sphereMappedRotFun );
2316  prepareBiCubicInterpolatorsPlusMinus ( std::round ( *bestLattitude ), std::round ( *bestLongitude ), sphereList, &interpolsPlusMinus, sphereMappedRotFun );
2317  prepareBiCubicInterpolatorsPlusPlus ( std::round ( *bestLattitude ), std::round ( *bestLongitude ), sphereList, &interpolsPlusPlus, sphereMappedRotFun );
2318 
2319  //================================================ Start the pseudo gradient ascent (while there is some change)
2320  while ( valChange > 0.0001 )
2321  {
2322  //============================================ Find the surrounding points to the currently best position
2323  lonM = *bestLongitude - step;
2324  lonP = *bestLongitude + step;
2325  latM = *bestLattitude - step;
2326  latP = *bestLattitude + step;
2327 
2328  //============================================ Deal with optimising outside of prepared range - recursion
2329  if ( latM < ( origBestLat - 1.0 ) ) { tmpVal = *bestLattitude; *bestLattitude = origBestLat - 1.0; optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step ); if ( *bestLattitude == origBestLat - 1.0 ) { *bestLattitude = tmpVal; } break; }
2330  if ( latP > ( origBestLat + 1.0 ) ) { tmpVal = *bestLattitude; *bestLattitude = origBestLat + 1.0; optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step ); if ( *bestLattitude == origBestLat + 1.0 ) { *bestLattitude = tmpVal; } break; }
2331  if ( lonM < ( origBestLon - 1.0 ) ) { tmpVal = *bestLongitude; *bestLongitude = origBestLon - 1.0; optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step ); if ( *bestLongitude == origBestLon - 1.0 ) { *bestLongitude = tmpVal; } break; }
2332  if ( lonP > ( origBestLon + 1.0 ) ) { tmpVal = *bestLongitude; *bestLongitude = origBestLon + 1.0; optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step ); if ( *bestLongitude == origBestLon + 1.0 ) { *bestLongitude = tmpVal; } break; }
2333 
2334  //============================================ Prepare vectors of tested positions
2335  latVals.at(0) = latM; latVals.at(1) = *bestLattitude; latVals.at(2) = latP;
2336  lonVals.at(0) = lonM; lonVals.at(1) = *bestLongitude; lonVals.at(2) = lonP;
2337 
2338  //============================================ Find the best change
2339  for ( proshade_unsign laIt = 0; laIt < static_cast<proshade_unsign> ( latVals.size() ); laIt++ )
2340  {
2341  for ( proshade_unsign loIt = 0; loIt < static_cast<proshade_unsign> ( lonVals.size() ); loIt++ )
2342  {
2343  //==================================== For this combination of lat and lon, find sum over spheres
2344  movSum = 1.0;
2345  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( sphereList->size() ); iter++ )
2346  {
2347  //================================ Interpolate using correct interpolators
2348  if ( ( latVals.at(laIt) <= origBestLat ) && ( lonVals.at(loIt) <= origBestLon ) ) { movSum += interpolsMinusMinus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2349  if ( ( latVals.at(laIt) <= origBestLat ) && ( lonVals.at(loIt) > origBestLon ) ) { movSum += interpolsMinusPlus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2350  if ( ( latVals.at(laIt) > origBestLat ) && ( lonVals.at(loIt) <= origBestLon ) ) { movSum += interpolsPlusMinus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2351  if ( ( latVals.at(laIt) > origBestLat ) && ( lonVals.at(loIt) > origBestLon ) ) { movSum += interpolsPlusPlus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2352  }
2353 
2354  //==================================== If position has improved, save it
2355  if ( *bestSum < movSum )
2356  {
2357  *bestSum = movSum;
2358  *bestLongitude = lonVals.at(loIt);
2359  *bestLattitude = latVals.at(laIt);
2360  }
2361  }
2362  }
2363 
2364  //============================================ Prepare for next iteration
2365  valChange = std::floor ( 100000.0 * ( *bestSum - prevVal ) ) / 100000.0;
2366  prevVal = std::floor ( 100000.0 * ( *bestSum ) ) / 100000.0;
2367  step = std::max ( ( valChange / step ) * learningRate, 0.01 );
2368  if ( learningRate >= 0.02 ) { learningRate -= 0.01; }
2369 
2370  }
2371 
2372  //================================================ Release interpolators memory
2373  for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsMinusMinus.size() ); intIt++ ) { delete interpolsMinusMinus.at(intIt); }
2374  for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsMinusPlus.size() ); intIt++ ) { delete interpolsMinusPlus.at(intIt); }
2375  for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsPlusMinus.size() ); intIt++ ) { delete interpolsPlusMinus.at(intIt); }
2376  for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsPlusPlus.size() ); intIt++ ) { delete interpolsPlusPlus.at(intIt); }
2377 
2378  //================================================ Done
2379  return ;
2380 }
2381 
2393 void ProSHADE_internal_maths::prepareBiCubicInterpolatorsMinusMinus ( proshade_double bestLattitude, proshade_double bestLongitude, std::vector<proshade_unsign>* sphereList, std::vector<ProSHADE_internal_maths::BicubicInterpolator*>* interpols, std::vector<ProSHADE_internal_spheres::ProSHADE_rotFun_sphere*>* sphereMappedRotFun )
2394 {
2395  //================================================ Initialise local variables
2396  proshade_signed latHlp, lonHlp;
2397  proshade_signed angDim = sphereMappedRotFun->at(0)->getAngularDim();
2398 
2399  //================================================ Prepare the interpolator objects for interpolation around the position
2400  for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2401  {
2402  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along first dimension)
2403  proshade_double** interpGrid = new proshade_double*[4];
2404  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid, __FILE__, __LINE__, __func__ );
2405 
2406  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along second dimension)
2407  for ( proshade_unsign iter = 0; iter < 4; iter++ )
2408  {
2409  interpGrid[iter] = new proshade_double[4];
2410  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid[iter], __FILE__, __LINE__, __func__ );
2411  }
2412 
2413  //============================================ Fill in the value grid on which the interpolation is to be done
2414  for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2415  {
2416  for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2417  {
2418  latHlp = bestLattitude - 2 + latIt; if ( latHlp < 0.0 ) { latHlp += angDim; } if ( latHlp >= angDim ) { latHlp -= angDim; }
2419  lonHlp = bestLongitude - 2 + lonIt; if ( lonHlp < 0.0 ) { lonHlp += angDim; } if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2420  interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( latHlp, lonHlp );
2421  }
2422  }
2423 
2424  //============================================ Create the interpolators
2425  ProSHADE_internal_maths::BicubicInterpolator* biCubInterp = new ProSHADE_internal_maths::BicubicInterpolator ( interpGrid, bestLattitude - 1.0, bestLongitude - 1.0 );
2426  interpols->emplace_back ( biCubInterp );
2427 
2428  //============================================ Release memory
2429  for ( proshade_unsign iter = 0; iter < 4; iter++ ) { delete[] interpGrid[iter]; }
2430  delete[] interpGrid;
2431  }
2432 
2433  //================================================ Done
2434  return ;
2435 }
2436 
2448 void ProSHADE_internal_maths::prepareBiCubicInterpolatorsMinusPlus ( proshade_double bestLattitude, proshade_double bestLongitude, std::vector<proshade_unsign>* sphereList, std::vector<ProSHADE_internal_maths::BicubicInterpolator*>* interpols, std::vector<ProSHADE_internal_spheres::ProSHADE_rotFun_sphere*>* sphereMappedRotFun )
2449 {
2450  //================================================ Initialise local variables
2451  proshade_signed latHlp, lonHlp;
2452  proshade_signed angDim = sphereMappedRotFun->at(0)->getAngularDim();
2453 
2454  //================================================ Prepare the interpolator objects for interpolation around the position
2455  for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2456  {
2457  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along first dimension)
2458  proshade_double** interpGrid = new proshade_double*[4];
2459  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid, __FILE__, __LINE__, __func__ );
2460 
2461  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along second dimension)
2462  for ( proshade_unsign iter = 0; iter < 4; iter++ )
2463  {
2464  interpGrid[iter] = new proshade_double[4];
2465  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid[iter], __FILE__, __LINE__, __func__ );
2466  }
2467 
2468  //============================================ Fill in the value grid on which the interpolation is to be done
2469  for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2470  {
2471  for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2472  {
2473  latHlp = bestLattitude - 2 + latIt; if ( latHlp < 0.0 ) { latHlp += angDim; } if ( latHlp >= angDim ) { latHlp -= angDim; }
2474  lonHlp = bestLongitude - 1 + lonIt; if ( lonHlp < 0.0 ) { lonHlp += angDim; } if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2475  interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( latHlp, lonHlp );
2476  }
2477  }
2478 
2479  //============================================ Create the interpolators
2480  ProSHADE_internal_maths::BicubicInterpolator* biCubInterp = new ProSHADE_internal_maths::BicubicInterpolator ( interpGrid, bestLattitude - 1.0, bestLongitude );
2481  interpols->emplace_back ( biCubInterp );
2482 
2483  //============================================ Release memory
2484  for ( proshade_unsign iter = 0; iter < 4; iter++ ) { delete[] interpGrid[iter]; }
2485  delete[] interpGrid;
2486  }
2487 
2488  //================================================ Done
2489  return ;
2490 }
2491 
2503 void ProSHADE_internal_maths::prepareBiCubicInterpolatorsPlusMinus ( proshade_double bestLattitude, proshade_double bestLongitude, std::vector<proshade_unsign>* sphereList, std::vector<ProSHADE_internal_maths::BicubicInterpolator*>* interpols, std::vector<ProSHADE_internal_spheres::ProSHADE_rotFun_sphere*>* sphereMappedRotFun )
2504 {
2505  //================================================ Initialise local variables
2506  proshade_signed latHlp, lonHlp;
2507  proshade_signed angDim = sphereMappedRotFun->at(0)->getAngularDim();
2508 
2509  //================================================ Prepare the interpolator objects for interpolation around the position
2510  for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2511  {
2512  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along first dimension)
2513  proshade_double** interpGrid = new proshade_double*[4];
2514  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid, __FILE__, __LINE__, __func__ );
2515 
2516  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along second dimension)
2517  for ( proshade_unsign iter = 0; iter < 4; iter++ )
2518  {
2519  interpGrid[iter] = new proshade_double[4];
2520  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid[iter], __FILE__, __LINE__, __func__ );
2521  }
2522 
2523  //============================================ Fill in the value grid on which the interpolation is to be done
2524  for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2525  {
2526  for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2527  {
2528  latHlp = bestLattitude - 1 + latIt; if ( latHlp < 0.0 ) { latHlp += angDim; } if ( latHlp >= angDim ) { latHlp -= angDim; }
2529  lonHlp = bestLongitude - 2 + lonIt; if ( lonHlp < 0.0 ) { lonHlp += angDim; } if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2530  interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( latHlp, lonHlp );
2531  }
2532  }
2533 
2534  //============================================ Create the interpolators
2535  ProSHADE_internal_maths::BicubicInterpolator* biCubInterp = new ProSHADE_internal_maths::BicubicInterpolator ( interpGrid, bestLattitude, bestLongitude - 1.0 );
2536  interpols->emplace_back ( biCubInterp );
2537 
2538  //============================================ Release memory
2539  for ( proshade_unsign iter = 0; iter < 4; iter++ ) { delete[] interpGrid[iter]; }
2540  delete[] interpGrid;
2541  }
2542 
2543  //================================================ Done
2544  return ;
2545 }
2546 
2558 void ProSHADE_internal_maths::prepareBiCubicInterpolatorsPlusPlus ( proshade_double bestLattitude, proshade_double bestLongitude, std::vector<proshade_unsign>* sphereList, std::vector<ProSHADE_internal_maths::BicubicInterpolator*>* interpols, std::vector<ProSHADE_internal_spheres::ProSHADE_rotFun_sphere*>* sphereMappedRotFun )
2559 {
2560  //================================================ Initialise local variables
2561  proshade_signed latHlp, lonHlp;
2562  proshade_signed angDim = sphereMappedRotFun->at(0)->getAngularDim();
2563 
2564  //================================================ Prepare the interpolator objects for interpolation around the position
2565  for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2566  {
2567  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along first dimension)
2568  proshade_double** interpGrid = new proshade_double*[4];
2569  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid, __FILE__, __LINE__, __func__ );
2570 
2571  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along second dimension)
2572  for ( proshade_unsign iter = 0; iter < 4; iter++ )
2573  {
2574  interpGrid[iter] = new proshade_double[4];
2575  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid[iter], __FILE__, __LINE__, __func__ );
2576  }
2577 
2578  //============================================ Fill in the value grid on which the interpolation is to be done
2579  for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2580  {
2581  for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2582  {
2583  latHlp = bestLattitude - 1 + latIt; if ( latHlp < 0.0 ) { latHlp += angDim; } if ( latHlp >= angDim ) { latHlp -= angDim; }
2584  lonHlp = bestLongitude - 1 + lonIt; if ( lonHlp < 0.0 ) { lonHlp += angDim; } if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2585  interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( latHlp, lonHlp );
2586  }
2587  }
2588 
2589  //============================================ Create the interpolators
2590  ProSHADE_internal_maths::BicubicInterpolator* biCubInterp = new ProSHADE_internal_maths::BicubicInterpolator ( interpGrid, bestLattitude, bestLongitude );
2591  interpols->emplace_back ( biCubInterp );
2592 
2593  //============================================ Release memory
2594  for ( proshade_unsign iter = 0; iter < 4; iter++ ) { delete[] interpGrid[iter]; }
2595  delete[] interpGrid;
2596  }
2597 
2598  //================================================ Done
2599  return ;
2600 }
2601 
2613 bool ProSHADE_internal_maths::isAxisUnique ( std::vector< proshade_double* >* CSymList, proshade_double* axis, proshade_double tolerance, bool improve )
2614 {
2615  //================================================ Initialise variables
2616  bool ret = true;
2617  proshade_unsign whichImprove;
2618 
2619  //================================================ For each already detected member
2620  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( CSymList->size() ); grIt++ )
2621  {
2622  //============================================ Is fold the same?
2623  if ( CSymList->at(grIt)[0] == axis[0] )
2624  {
2625  if ( ProSHADE_internal_maths::vectorOrientationSimilarity ( CSymList->at(grIt)[1], CSymList->at(grIt)[2], CSymList->at(grIt)[3], axis[1], axis[2], axis[3], tolerance ) )
2626  {
2627  ret = false;
2628  whichImprove = grIt;
2629  break;
2630  }
2631  }
2632  }
2633 
2634  //================================================ Improve, if required
2635  if ( improve && !ret )
2636  {
2637  CSymList->at(whichImprove)[1] = axis[1];
2638  CSymList->at(whichImprove)[2] = axis[2];
2639  CSymList->at(whichImprove)[3] = axis[3];
2640  CSymList->at(whichImprove)[4] = axis[4];
2641  CSymList->at(whichImprove)[5] = axis[5];
2642  }
2643 
2644  //================================================ Done
2645  return ( ret );
2646 
2647 }
2648 
2661 bool ProSHADE_internal_maths::isAxisUnique ( std::vector< proshade_double* >* CSymList, proshade_double X, proshade_double Y, proshade_double Z, proshade_double fold, proshade_double tolerance )
2662 {
2663  //================================================ Initialise variables
2664  bool ret = true;
2665 
2666  //================================================ For each already detected member
2667  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( CSymList->size() ); grIt++ )
2668  {
2669  if ( fold == CSymList->at(grIt)[0] )
2670  {
2671  if ( ProSHADE_internal_maths::vectorOrientationSimilarity ( CSymList->at(grIt)[1], CSymList->at(grIt)[2], CSymList->at(grIt)[3], X, Y, Z, tolerance ) )
2672  {
2673  ret = false;
2674  break;
2675  }
2676  }
2677  }
2678  std::cout << std::endl;
2679 
2680  //================================================ Done
2681  return ( ret );
2682 
2683 }
2684 
2693 std::vector< proshade_unsign > ProSHADE_internal_maths::findAllPrimes ( proshade_unsign upTo )
2694 {
2695  //================================================ Initialise variables
2696  std::vector< proshade_unsign > ret;
2697  std::vector< std::pair< proshade_unsign, bool > > sieveOfEratosthenesArray;
2698 
2699  //================================================ Sanity check
2700  if ( upTo < 2 ) { return ( ret ); }
2701 
2702  //================================================ Initialise the sieve array up to the required number
2703  for ( proshade_unsign iter = 2; iter <= upTo; iter++ ) { sieveOfEratosthenesArray.emplace_back ( std::pair< proshade_unsign, bool > ( iter, true ) ); }
2704 
2705  //================================================ For each entry in the array
2706  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( sieveOfEratosthenesArray.size() ); iter++ )
2707  {
2708  //============================================ If this entry is still true
2709  if ( sieveOfEratosthenesArray.at(iter).second )
2710  {
2711  //======================================== Set all entries with the position x * [this entry value] to false
2712  for ( proshade_unsign it = iter + sieveOfEratosthenesArray.at(iter).first; it < static_cast<proshade_unsign> ( sieveOfEratosthenesArray.size() ); it += sieveOfEratosthenesArray.at(iter).first )
2713  {
2714  sieveOfEratosthenesArray.at(it).second = false;
2715  }
2716  }
2717  }
2718 
2719  //================================================ Copy passing results to return vector
2720  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( sieveOfEratosthenesArray.size() ); iter++ )
2721  {
2722  if ( sieveOfEratosthenesArray.at(iter).second ) { ProSHADE_internal_misc::addToUnsignVector ( &ret, sieveOfEratosthenesArray.at(iter).first ); }
2723  }
2724 
2725  //================================================ Done
2726  return ( ret );
2727 
2728 }
ProSHADE_internal_maths::gaussLegendreIntegration
void gaussLegendreIntegration(proshade_complex *vals, proshade_unsign valsSize, proshade_unsign order, proshade_double *abscissas, proshade_double *weights, proshade_double integralOverRange, proshade_double maxSphereDists, proshade_double *retReal, proshade_double *retImag)
Function to compute the complete complex Gauss-Legendre integration over spherical harmonic values in...
Definition: ProSHADE_maths.cpp:709
ProSHADE_internal_maths::findVectorFromThreeVAndThreeD
std::vector< proshade_double > findVectorFromThreeVAndThreeD(proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2, proshade_double x3, proshade_double y3, proshade_double z3, proshade_double dot1, proshade_double dot2, proshade_double dot3)
Function for finding a vector which would have a given three dot products to three other vectors.
Definition: ProSHADE_maths.cpp:2101
ProSHADE_internal_maths::gaussLegendreIntegrationReal
proshade_double gaussLegendreIntegrationReal(proshade_double *vals, proshade_unsign valsSize, proshade_unsign order, proshade_double *abscissas, proshade_double *weights, proshade_double integralOverRange, proshade_double maxSphereDists)
Function to compute real part of the Gauss-Legendre integration over spherical harmonic values in dif...
Definition: ProSHADE_maths.cpp:619
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_maths::prepareBiCubicInterpolatorsMinusPlus
void prepareBiCubicInterpolatorsMinusPlus(proshade_double bestLattitude, proshade_double bestLongitude, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *interpols, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun)
This function prepares the interpolation objects for the bi-cubic interpolation.
Definition: ProSHADE_maths.cpp:2448
ProSHADE_internal_maths::computeCrossProduct
proshade_double * computeCrossProduct(proshade_double *x1, proshade_double *y1, proshade_double *z1, proshade_double *x2, proshade_double *y2, proshade_double *z2)
Simple 3D vector cross product computation.
Definition: ProSHADE_maths.cpp:1737
ProSHADE_internal_maths::getEulerZXZFromRotMatrix
void getEulerZXZFromRotMatrix(proshade_double *rotMat, proshade_double *eA, proshade_double *eB, proshade_double *eG)
This function converts rotation matrix to the Euler ZXZ angles representation.
Definition: ProSHADE_maths.cpp:1394
ProSHADE_exception
This class is the representation of ProSHADE exception.
Definition: ProSHADE_exceptions.hpp:37
ProSHADE_internal_maths::normalDistributionValue
proshade_double normalDistributionValue(proshade_double mean, proshade_double standardDev, proshade_double value)
Function to the heiht of normal distribution given by mean and standard deviation for a given value.
Definition: ProSHADE_maths.cpp:1688
ProSHADE_internal_maths::prepareBiCubicInterpolatorsPlusMinus
void prepareBiCubicInterpolatorsPlusMinus(proshade_double bestLattitude, proshade_double bestLongitude, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *interpols, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun)
This function prepares the interpolation objects for the bi-cubic interpolation.
Definition: ProSHADE_maths.cpp:2503
ProSHADE_internal_maths::evaluateGLSeries
proshade_double evaluateGLSeries(proshade_double *series, proshade_double target, proshade_unsign terms)
This function evaluates the Taylor expansion.
Definition: ProSHADE_maths.cpp:447
ProSHADE_internal_maths::getSOFTPositionFromEulerZXZ
void getSOFTPositionFromEulerZXZ(proshade_signed band, proshade_double eulerAlpha, proshade_double eulerBeta, proshade_double eulerGamma, proshade_double *x, proshade_double *y, proshade_double *z)
Function to find the index position in the inverse SOFT map from given Euler angles (ZXZ convention).
Definition: ProSHADE_maths.cpp:986
ProSHADE_maths.hpp
This header file declares all the functions required for computing various information from the ProSH...
ProSHADE_internal_maths::getEulerZXZFromAngleAxis
void getEulerZXZFromAngleAxis(proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double axAng, proshade_double *eA, proshade_double *eB, proshade_double *eG, proshade_unsign angDim)
This function converts angle-axis representation to the Euler ZXZ angles representation.
Definition: ProSHADE_maths.cpp:1437
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:1005
ProSHADE_internal_maths::complexMultiplication
void complexMultiplication(proshade_double *r1, proshade_double *i1, proshade_double *r2, proshade_double *i2, proshade_double *retReal, proshade_double *retImag)
Function to multiply two complex numbers.
Definition: ProSHADE_maths.cpp:38
ProSHADE_internal_maths::getGLFirstEvenRoot
void getGLFirstEvenRoot(proshade_double polyAtZero, proshade_unsign order, proshade_double *abscAtZero, proshade_double *weighAtZero, proshade_unsign taylorSeriesCap)
This function finds the first root for Legendre polynomials of odd order.
Definition: ProSHADE_maths.cpp:384
ProSHADE_internal_maths::optimiseAxisBiCubicInterpolation
void optimiseAxisBiCubicInterpolation(proshade_double *bestLattitude, proshade_double *bestLongitude, proshade_double *bestSum, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun, proshade_double step=0.05)
This function provides axis optimisation given starting lattitude and longitude indices.
Definition: ProSHADE_maths.cpp:2296
ProSHADE_internal_maths::vectorOrientationSimilarity
bool vectorOrientationSimilarity(proshade_double a1, proshade_double a2, proshade_double a3, proshade_double b1, proshade_double b2, proshade_double b3, proshade_double tolerance=0.1)
This function compares two vectors using cosine distance and decides if they are similar using tolera...
Definition: ProSHADE_maths.cpp:2232
ProSHADE_internal_maths::complexMatrixSVDSigmasOnly
void complexMatrixSVDSigmasOnly(proshade_complex **mat, int dim, double *&singularValues)
Function to compute the complete complex matrix SVD and return only the sigmas.
Definition: ProSHADE_maths.cpp:802
ProSHADE_internal_maths::prepareBiCubicInterpolatorsMinusMinus
void prepareBiCubicInterpolatorsMinusMinus(proshade_double bestLattitude, proshade_double bestLongitude, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *interpols, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun)
This function prepares the interpolation objects for the bi-cubic interpolation.
Definition: ProSHADE_maths.cpp:2393
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_maths::vectorOrientationSimilaritySameDirection
bool vectorOrientationSimilaritySameDirection(proshade_double a1, proshade_double a2, proshade_double a3, proshade_double b1, proshade_double b2, proshade_double b3, proshade_double tolerance=0.1)
This function compares two vectors using cosine distance and decides if they are similar using tolera...
Definition: ProSHADE_maths.cpp:2266
ProSHADE_internal_maths::pearsonCorrCoeff
proshade_double pearsonCorrCoeff(proshade_double *valSet1, proshade_double *valSet2, proshade_unsign length)
Function for computing the Pearson's correlation coefficient.
Definition: ProSHADE_maths.cpp:244
ProSHADE_internal_maths::computeDotProduct
proshade_double computeDotProduct(proshade_double *x1, proshade_double *y1, proshade_double *z1, proshade_double *x2, proshade_double *y2, proshade_double *z2)
Simple 3D vector dot product computation.
Definition: ProSHADE_maths.cpp:1705
ProSHADE_internal_misc::addToSignedVector
void addToSignedVector(std::vector< proshade_signed > *vecToAddTo, proshade_signed elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:121
ProSHADE_internal_maths::complexMultiplicationConjugRealOnly
proshade_double complexMultiplicationConjugRealOnly(proshade_double *r1, proshade_double *i1, proshade_double *r2, proshade_double *i2)
Function to conjuggate multiply two complex numbers and return the real part only.
Definition: ProSHADE_maths.cpp:103
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:1343
ProSHADE_internal_maths::BicubicInterpolator
Definition: ProSHADE_maths.hpp:102
ProSHADE_internal_maths::primeFactorsDecomp
std::vector< proshade_signed > primeFactorsDecomp(proshade_signed number)
Function to find prime factors of an integer.
Definition: ProSHADE_maths.cpp:1639
ProSHADE_internal_maths::getGLPolyAtZero
void getGLPolyAtZero(proshade_unsign order, proshade_double *polyValue, proshade_double *deriValue)
This function obtains the Legendre polynomial values and its derivative at zero for any positive inte...
Definition: ProSHADE_maths.cpp:347
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:2613
ProSHADE_internal_maths::getLegendreAbscAndWeights
void getLegendreAbscAndWeights(proshade_unsign order, proshade_double *abscissas, proshade_double *weights, proshade_unsign taylorSeriesCap)
Function to prepare abscissas and weights for Gauss-Legendre integration.
Definition: ProSHADE_maths.cpp:287
ProSHADE_internal_maths::getAxisAngleFromRotationMatrix
void getAxisAngleFromRotationMatrix(proshade_double *rotMat, proshade_double *x, proshade_double *y, proshade_double *z, proshade_double *ang)
This function converts rotation matrix to the axis-angle representation.
Definition: ProSHADE_maths.cpp:1039
ProSHADE_internal_maths::advanceGLPolyValue
proshade_double advanceGLPolyValue(proshade_double from, proshade_double to, proshade_double valAtFrom, proshade_unsign noSteps, proshade_unsign taylorSeriesCap)
This function finds the next value of the polynomial.
Definition: ProSHADE_maths.cpp:477
ProSHADE_internal_maths::findRotMatMatchingVectors
proshade_double * findRotMatMatchingVectors(proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2)
Computation of rotation matrix rotating one vector onto the other.
Definition: ProSHADE_maths.cpp:1883
ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication
proshade_double * compute3x3MatrixVectorMultiplication(proshade_double *mat, proshade_double x, proshade_double y, proshade_double z)
Function for computing a 3x3 matrix to 3x1 vector multiplication.
Definition: ProSHADE_maths.cpp:1789
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:2195
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_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:961
ProSHADE_internal_maths::vectorMedianAndIQR
void vectorMedianAndIQR(std::vector< proshade_double > *vec, proshade_double *&ret)
Function to get vector median and inter-quartile range.
Definition: ProSHADE_maths.cpp:147
ProSHADE_internal_maths::getEulerZXZFromAngleAxisFullSearch
void getEulerZXZFromAngleAxisFullSearch(proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double axAng, proshade_double *eA, proshade_double *eB, proshade_double *eG, proshade_signed angDim)
This function converts angle-axis representation to the Euler ZXZ angles representation using full se...
Definition: ProSHADE_maths.cpp:1526
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_maths::complexMultiplicationConjug
void complexMultiplicationConjug(proshade_double *r1, proshade_double *i1, proshade_double *r2, proshade_double *i2, proshade_double *retReal, proshade_double *retImag)
Function to multiply two complex numbers by using the second number's conjugate.
Definition: ProSHADE_maths.cpp:62
ProSHADE_internal_maths::transpose3x3MatrixInPlace
void transpose3x3MatrixInPlace(proshade_double *mat)
Transposes 3x3 matrix in place.
Definition: ProSHADE_maths.cpp:1843
ProSHADE_internal_maths::completeLegendreSeries
void completeLegendreSeries(proshade_unsign order, proshade_double *abscissa, proshade_double *weights, proshade_unsign taylorSeriesCap)
This function completes the Legendre polynomial series assuming you have obtained the first values.
Definition: ProSHADE_maths.cpp:521
ProSHADE_internal_maths::findVectorFromTwoVAndTwoD
std::vector< proshade_double > findVectorFromTwoVAndTwoD(proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2, proshade_double dot1, proshade_double dot2)
Function for finding a vector which would have a given two dot products to two other vectors.
Definition: ProSHADE_maths.cpp:1959
ProSHADE_internal_maths::prepareBiCubicInterpolatorsPlusPlus
void prepareBiCubicInterpolatorsPlusPlus(proshade_double bestLattitude, proshade_double bestLongitude, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *interpols, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun)
This function prepares the interpolation objects for the bi-cubic interpolation.
Definition: ProSHADE_maths.cpp:2558
ProSHADE_internal_maths::multiplyTwoSquareMatrices
void multiplyTwoSquareMatrices(proshade_double *A, proshade_double *B, proshade_double *res, proshade_unsign dim=3)
Function to compute matrix multiplication.
Definition: ProSHADE_maths.cpp:1613
ProSHADE_internal_maths::findAllPrimes
std::vector< proshade_unsign > findAllPrimes(proshade_unsign upTo)
This function finds all prime numbers up to the supplied limit.
Definition: ProSHADE_maths.cpp:2693
ProSHADE_internal_maths::compute3x3MatrixInverse
proshade_double * compute3x3MatrixInverse(proshade_double *mat)
Function for computing a 3x3 matrix inverse.
Definition: ProSHADE_maths.cpp:1810
ProSHADE_internal_maths::complexMatrixSVDUandVOnly
void complexMatrixSVDUandVOnly(proshade_double *mat, int dim, proshade_double *uAndV, bool fail=true)
Function to compute the real matrix SVD and return the U and V matrices.
Definition: ProSHADE_maths.cpp:869
ProSHADE_internal_maths::arrayMedianAndIQR
void arrayMedianAndIQR(proshade_double *vec, proshade_unsign vecSize, proshade_double *&ret)
Function to get array median and inter-quartile range.
Definition: ProSHADE_maths.cpp:198
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:2141
ProSHADE_internal_maths::complexMultiplicationRealOnly
proshade_double complexMultiplicationRealOnly(proshade_double *r1, proshade_double *i1, proshade_double *r2, proshade_double *i2)
Function to multiply two complex numbers and return the real part only.
Definition: ProSHADE_maths.cpp:83
ProSHADE_internal_maths::compute3x3MatrixMultiplication
proshade_double * compute3x3MatrixMultiplication(proshade_double *mat1, proshade_double *mat2)
Function for computing a 3x3 matrix multiplication.
Definition: ProSHADE_maths.cpp:1759