Xmipp  v3.23.11-Nereus
numerical_tools.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  *
3  * Authors: Carlos Oscar S. Sorzano (coss@cnb.csic.es)
4  *
5  * Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC
6  *
7  * This program is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
20  * 02111-1307 USA
21  *
22  * All comments concerning this program package may be sent to the
23  * e-mail address 'xmipp@cnb.csic.es'
24  ***************************************************************************/
25 #include <vector>
26 #include "numerical_tools.h"
27 #include "core/matrix2d.h"
28 #include "core/numerical_recipes.h"
29 #include "core/xmipp_funcs.h"
30 #include "core/xmipp_memory.h"
31 
32 /* Random permutation ------------------------------------------------------ */
33 void randomPermutation(int N, MultidimArray<int>& result)
34 {
35  MultidimArray<double> aux;
36  aux.resize(N);
37  aux.initRandom(0,1);
38 
39  aux.indexSort(result);
40  result-=1;
41 }
42 
43 /* Powell's optimizer ------------------------------------------------------ */
44 void powellOptimizer(Matrix1D<double> &p, int i0, int n,
45  double(*f)(double *x, void *), void * prm,
46  double ftol, double &fret,
47  int &iter, const Matrix1D<double> &steps, bool show)
48 {
49  // Adapt indexes of p
50  double *pptr = p.adaptForNumericalRecipes();
51  double *auxpptr = pptr + (i0 - 1);
52 
53  // Form direction matrix
54  std::vector<double> buffer(n*n);
55  auto *xi= buffer.data()-1;
56  for (int i = 1, ptr = 1; i <= n; i++)
57  for (int j = 1; j <= n; j++, ptr++)
58  xi[ptr] = (i == j) ? steps(i - 1) : 0;
59 
60  // Optimize
61  powell(auxpptr, xi -n, n, ftol, iter, fret, f, prm, show); // xi - n because NR works with matrices starting at [1,1]
62 }
63 
64 /* Gaussian interpolator -------------------------------------------------- */
65 void GaussianInterpolator::initialize(double _xmax, int N, bool normalize)
66 {
67  xmax=_xmax;
68  xstep=xmax/N;
69  ixstep=1.0/xstep;
70  v.initZeros(N);
71  double inorm=1.0/sqrt(2*PI);
72  FOR_ALL_ELEMENTS_IN_ARRAY1D(v)
73  {
74  double x=i*xstep;
75  v(i)=exp(-x*x/2);
76  if (normalize)
77  v(i)*=inorm;
78  }
79 }
80 
81 /* Solve Cx=d, nonnegative x */
82 double solveNonNegative(const Matrix2D<double> &C, const Matrix1D<double> &d,
83  Matrix1D<double> &result)
84 {
85  if (C.Xdim() == 0)
86  REPORT_ERROR(ERR_MATRIX_EMPTY, "Solve_nonneg: Matrix is empty");
87  if (C.Ydim() != d.size())
88  REPORT_ERROR(ERR_MATRIX_SIZE, "Solve_nonneg: Different sizes of Matrix and Vector");
89  if (d.isRow())
90  REPORT_ERROR(ERR_MATRIX_DIM, "Solve_nonneg: Not correct vector shape");
91 
92  Matrix2D<double> Ct;
93  Ct=C.transpose();
94 
95  result.initZeros(Ct.Ydim());
96  double rnorm;
97 
98  // Watch out that matrix Ct is transformed.
99  int success = nnls(MATRIX2D_ARRAY(Ct), Ct.Xdim(), Ct.Ydim(),
100  MATRIX1D_ARRAY(d),
101  MATRIX1D_ARRAY(result),
102  &rnorm, nullptr, nullptr, nullptr);
103  if (success == 1)
104  std::cerr << "Warning, too many iterations in nnls\n";
105  else if (success == 2)
106  REPORT_ERROR(ERR_MEM_NOTENOUGH, "Solve_nonneg: Not enough memory");
107  return rnorm;
108 }
109 
110 /* Solve Ax=b, A definite positive and symmetric --------------------------- */
111 void solveViaCholesky(const Matrix2D<double> &A, const Matrix1D<double> &b,
112  Matrix1D<double> &result)
113 {
114  Matrix2D<double> Ap = A;
115  Matrix1D<double> p(A.Xdim());
116  result.resize(A.Xdim());
117  choldc(Ap.adaptForNumericalRecipes2(), A.Xdim(),
118  p.adaptForNumericalRecipes());
119  cholsl(Ap.adaptForNumericalRecipes2(), A.Xdim(),
120  p.adaptForNumericalRecipes(), b.adaptForNumericalRecipes(),
121  result.adaptForNumericalRecipes());
122 }
123 
124 
125 /* Quadratic form ---------------------------------------------------------- */
126 void evaluateQuadratic(const Matrix1D<double> &x, const Matrix1D<double> &c,
127  const Matrix2D<double> &H, double &val, Matrix1D<double> &grad)
128 {
129  if (x.size() != c.size())
130  REPORT_ERROR(ERR_MATRIX_SIZE, "Eval_quadratic: Not compatible sizes in x and c");
131  if (H.Xdim() != x.size())
132  REPORT_ERROR(ERR_MATRIX_SIZE, "Eval_quadratic: Not compatible sizes in x and H");
133 
134  // H*x, store in grad
135  grad.initZeros(x.size());
136  for (size_t i = 0; i < H.Ydim(); i++)
137  for (size_t j = 0; j < x.size(); j++)
138  grad(i) += H(i, j) * x(j);
139 
140  // Now, compute c^t*x+1/2*x^t*H*x
141  // Add c to the gradient
142  double quad = 0;
143  val = 0;
144  for (size_t j = 0; j < x.size(); j++)
145  {
146  quad += grad(j) * grad(j); // quad=x^t*H^t*H*x
147  val += c(j) * x(j); // val=c^t*x
148 
149  grad(j) += c(j); // grad+=c
150  }
151  val += 0.5 * quad;
152 }
153 
154 /* Quadprog and Lsqlin ----------------------------------------------------- */
155 /* Structure to pass the objective function and constraints to cfsqp*/
156 typedef struct
157 {
158  Matrix2D<double> C;
159  Matrix2D<double> D;
160  Matrix2D<double> A;
161  Matrix2D<double> B;
162 }
163 CDAB;
164 
165 /*/////////////////////////////////////////////////////////////////////////
166  Internal functions used by the quadraticProgramming function
168 /* To calculate the value of the objective function */
169 void quadraticProgramming_obj32(int nparam, int j, double* x, double* fj, void* cd)
170 {
171  auto* in = (CDAB *)cd;
172  Matrix2D<double> X(nparam,1);
173  for (int i=0; i<nparam; ++i)
174  X(i,0)=x[i];
175  Matrix2D<double> result;
176  result = 0.5 * X.transpose() * in->C * X + in->D.transpose() * X;
177 
178  *fj = result(0, 0);
179 }
180 
181 /* To calculate the value of the jth constraint */
182 void quadraticProgramming_cntr32(int nparam, int j, double* x, double* gj, void* cd)
183 {
184  auto* in = (CDAB *)cd;
185  *gj = 0;
186  for (int k = 0; k < nparam; k++)
187  *gj += in->A(j - 1, k) * x[k];
188  *gj -= in->B(j - 1, 0);
189 }
190 
191 /* To calculate the value of the derivative of objective function */
192 void quadraticProgramming_grob32(int nparam, double* x, double* gradfj, void *cd)
193 {
194  auto* in = (CDAB *)cd;
195  Matrix2D<double> X(1,nparam);
196  for (int i=0; i<nparam; ++i)
197  X(0,i)=x[i];
198 
199  Matrix2D<double> gradient;
200  gradient = in->C * X + in->D;
201  for (int k = 0; k < nparam; k++)
202  gradfj[k] = gradient(k, 0);
203 }
204 
205 /* To calculate the value of the derivative of jth constraint */
206 void quadraticProgramming_grcn32(int nparam, int j, double *gradgj, void *cd)
207 {
208  auto* in = (CDAB *)cd;
209  for (int k = 0; k < nparam; k++)
210  gradgj[k] = in->A(j - 1, k);
211 }
212 
213 /**************************************************************************
214 
215  Solves Quadratic programming subproblem.
216 
217  min 0.5*x'Cx + d'x subject to: A*x <= b
218  x Aeq*x=beq
219  bl<=x<=bu
220 
221 **************************************************************************/
222 void quadraticProgramming(const Matrix2D<double> &C, const Matrix1D<double> &d,
223  const Matrix2D<double> &A, const Matrix1D<double> &b,
224  const Matrix2D<double> &Aeq, const Matrix1D<double> &beq,
225  Matrix1D<double> &bl, Matrix1D<double> &bu,
226  Matrix1D<double> &x)
227 {
228  CDAB prm;
229  prm.C = C;
230  prm.D.fromVector(d);
231  prm.A.initZeros(A.Ydim() + Aeq.Ydim(), A.Xdim());
232  prm.B.initZeros(prm.A.Ydim(), 1);
233 
234 
235  // Copy Inequalities
236  for (size_t i = 0; i < A.Ydim(); i++)
237  {
238  for (size_t j = 0; j < A.Xdim(); j++)
239  prm.A(i, j) = A(i, j);
240  prm.B(i, 0) = b(i);
241  }
242 
243  // Copy Equalities
244  for (size_t i = 0; i < Aeq.Ydim(); i++)
245  {
246  for (size_t j = 0; j < Aeq.Xdim(); j++)
247  prm.A(i + A.Ydim(), j) = Aeq(i, j);
248  prm.B(i + A.Ydim(), 0) = beq(i);
249  }
250 
251  double bigbnd = 1e30;
252  // Bounds
253  if (bl.size() == 0)
254  {
255  bl.resize(C.Xdim());
256  bl.initConstant(-bigbnd);
257  }
258  if (bu.size() == 0)
259  {
260  bu.resize(C.Xdim());
261  bu.initConstant(bigbnd);
262  }
263 
264  // Define intermediate variables
265  int mode = 100; // CFSQP mode
266  int iprint = 0; // Debugging
267  int miter = 1000; // Maximum number of iterations
268  double eps = 1e-4; // Epsilon
269  double epsneq = 1e-4; // Epsilon for equalities
270  double udelta = 0.e0; // Finite difference approximation
271  // of the gradients. Not used in this function
272  int nparam = C.Xdim(); // Number of variables
273  int nf = 1; // Number of objective functions
274  int neqn = Aeq.Ydim(); // Number of nonlinear equations
275  int nineqn = A.Ydim(); // Number of nonlinear inequations
276  int nineq = A.Ydim(); // Number of linear inequations
277  int neq = Aeq.Ydim(); // Number of linear equations
278  int inform;
279  int ncsrl = 0;
280  int ncsrn = 0;
281  int nfsr = 0;
282  int mesh_pts[] = {0};
283 
284  if (x.size() == 0)
285  x.initZeros(nparam);
286  Matrix1D<double> f(nf);
287  Matrix1D<double> g(nineq + neq);
288  Matrix1D<double> lambda(nineq + neq + nf + nparam);
289 
290  // Call the minimization routine
291  cfsqp(nparam, nf, nfsr, nineqn, nineq, neqn, neq, ncsrl, ncsrn, mesh_pts,
292  mode, iprint, miter, &inform, bigbnd, eps, epsneq, udelta,
293  MATRIX1D_ARRAY(bl), MATRIX1D_ARRAY(bu),
294  MATRIX1D_ARRAY(x),
295  MATRIX1D_ARRAY(f), MATRIX1D_ARRAY(g),
296  MATRIX1D_ARRAY(lambda),
297  // quadprg_obj32,quadprog_cntr32,quadprog_grob32,quadprog_grcn32,
298  quadraticProgramming_obj32, quadraticProgramming_cntr32, grobfd, grcnfd,
299  (void*)&prm);
300 
301 #ifdef DEBUG
302 
303  if (inform == 0)
304  std::cout << "SUCCESSFUL RETURN. \n";
305  if (inform == 1 || inform == 2)
306  std::cout << "\nINITIAL GUESS INFEASIBLE.\n";
307  if (inform == 3)
308  printf("\n MAXIMUM NUMBER OF ITERATIONS REACHED.\n");
309  if (inform > 3)
310  printf("\ninform=%d\n", inform);
311 #endif
312 }
313 
314 /**************************************************************************
315 
316  Solves the least square problem
317 
318  min 0.5*(Norm(C*x-d)) subject to: A*x <= b
319  x Aeq*x=beq
320  bl<=x<=bu
321 **************************************************************************/
322 void leastSquare(const Matrix2D<double> &C, const Matrix1D<double> &d,
323  const Matrix2D<double> &A, const Matrix1D<double> &b,
324  const Matrix2D<double> &Aeq, const Matrix1D<double> &beq,
325  Matrix1D<double> &bl, Matrix1D<double> &bu,
326  Matrix1D<double> &x)
327 {
328  // Convert d to Matrix2D for multiplication
329  Matrix2D<double> P;
330  P.fromVector(d);
331  P = -2.0 * (P.transpose() * C);
332  P = P.transpose();
333 
334  //Convert back to vector for passing it to quadraticProgramming
335  Matrix1D<double> newd;
336  P.toVector(newd);
337 
338  quadraticProgramming(C.transpose()*C, newd, A, b, Aeq, beq, bl, bu, x);
339 }
340 
341 /* Regularized least squares ----------------------------------------------- */
342 void regularizedLeastSquare(const Matrix2D< double >& A,
343  const Matrix1D< double >& d, double lambda,
344  const Matrix2D< double >& G, Matrix1D< double >& x)
345 {
346  int Nx=A.Xdim(); // Number of variables
347 
348  Matrix2D<double> X(Nx,Nx); // X=(A^t * A +lambda *G^t G)
349  // Compute A^t*A
350  FOR_ALL_ELEMENTS_IN_MATRIX2D(X)
351  // Compute the dot product of the i-th and j-th columns of A
352  for (size_t k=0; k<A.Ydim(); k++)
353  X(i,j) += A(k,i) * A(k,j);
354 
355  // Compute lambda*G^t*G
356  if (G.Xdim()==0)
357  for (int i=0; i<Nx; i++)
358  X(i,i)+=lambda;
359  else
360  FOR_ALL_ELEMENTS_IN_MATRIX2D(X)
361  // Compute the dot product of the i-th and j-th columns of G
362  for (size_t k=0; k<G.Ydim(); k++)
363  X(i,j) += G(k,i) * G(k,j);
364 
365  // Compute A^t*d
366  Matrix1D<double> Atd(Nx);
367  FOR_ALL_ELEMENTS_IN_MATRIX1D(Atd)
368  // Compute the dot product of the i-th column of A and d
369  for (size_t k=0; k<A.Ydim(); k++)
370  Atd(i) += A(k,i) * d(k);
371 
372  // Compute the inverse of X
373  Matrix2D<double> Xinv;
374  X.inv(Xinv);
375 
376  // Now multiply Xinv * A^t * d
377  matrixOperation_Ax(Xinv, Atd, x);
378 }
379 
380 
382 
383 #define Element(a,b,c) a[b*nDim+c]
384 #define RowVector(a,b) (&a[b*nDim])
385 #define CopyVector(a,b) memcpy((a),(b),nDim*sizeof(double))
386 
387 DESolver::DESolver(int dim, int popSize) :
388 nDim(dim), nPop(popSize),
389 generations(0), strategy(stRand1Exp),
390 scale(0.7), probability(0.5), bestEnergy(0.0),
391 trialSolution(0), bestSolution(0),
392 popEnergy(0), population(0)
393 {
394  trialSolution = new double[nDim];
395  bestSolution = new double[nDim];
396  popEnergy = new double[nPop];
397  population = new double[nPop * nDim];
398 
399  randomize_random_generator();
400 }
401 
402 DESolver::~DESolver(void)
403 {
404  if (trialSolution)
405  delete [] trialSolution;
406  if (bestSolution)
407  delete [] bestSolution;
408  if (popEnergy)
409  delete [] popEnergy;
410  if (population)
411  delete [] population;
412 
413  trialSolution = bestSolution = popEnergy = population = nullptr;
414  return;
415 }
416 
417 void DESolver::Setup(double min[], double max[],
418  int deStrategy, double diffScale, double crossoverProb)
419 {
420  int i;
421  int j;
422 
423  strategy = deStrategy;
424  scale = diffScale;
425  probability = crossoverProb;
426  bestEnergy = 1.0E20;
427 
428  for (i = 0; i < nPop; i++)
429  {
430  for (j = 0; j < nDim; j++)
431  {
432  population[i*nDim+j] = rnd_unif(min[j], max[j]);
433  //Element(population, i, j) = rnd_unif(min[j], max[j]);
434  }
435 
436  popEnergy[i] = 1.0E20;
437  }
438 
439  for (i = 0; i < nDim; i++)
440  bestSolution[i] = (min[i] + max[i]) / 2.0;
441 
442  switch (strategy)
443  {
444  case stBest1Exp:
445  calcTrialSolution = &DESolver::Best1Exp;
446  break;
447 
448  case stRand1Exp:
449  calcTrialSolution = &DESolver::Rand1Exp;
450  break;
451 
452  case stRandToBest1Exp:
453  calcTrialSolution = &DESolver::RandToBest1Exp;
454  break;
455 
456  case stBest2Exp:
457  calcTrialSolution = &DESolver::Best2Exp;
458  break;
459 
460  case stRand2Exp:
461  calcTrialSolution = &DESolver::Rand2Exp;
462  break;
463 
464  case stBest1Bin:
465  calcTrialSolution = &DESolver::Best1Bin;
466  break;
467 
468  case stRand1Bin:
469  calcTrialSolution = &DESolver::Rand1Bin;
470  break;
471 
472  case stRandToBest1Bin:
473  calcTrialSolution = &DESolver::RandToBest1Bin;
474  break;
475 
476  case stBest2Bin:
477  calcTrialSolution = &DESolver::Best2Bin;
478  break;
479 
480  case stRand2Bin:
481  calcTrialSolution = &DESolver::Rand2Bin;
482  break;
483  }
484 
485  return;
486 }
487 
488 bool DESolver::Solve(int maxGenerations)
489 {
490  bool bAtSolution = false;
491  int generation;
492 
493  for (generation = 0;(generation < maxGenerations) && !bAtSolution;generation++)
494  for (int candidate = 0; candidate < nPop; candidate++)
495  {
496  (this->*calcTrialSolution)(candidate);
497  trialEnergy = EnergyFunction(trialSolution, bAtSolution);
498 
499  if (trialEnergy < popEnergy[candidate])
500  {
501  // New low for this candidate
502  popEnergy[candidate] = trialEnergy;
503  CopyVector(RowVector(population, candidate), trialSolution);
504 
505  // Check if all-time low
506  if (trialEnergy < bestEnergy)
507  {
508  bestEnergy = trialEnergy;
509  CopyVector(bestSolution, trialSolution);
510  }
511  }
512  }
513 
514  generations = generation;
515  return bAtSolution;
516 }
517 
518 void DESolver::Best1Exp(int candidate)
519 {
520  int r1;
521  int r2;
522  int n;
523 
524  SelectSamples(candidate, &r1, &r2);
525  n = (int)rnd_unif(0.0, (double)nDim);
526 
527  CopyVector(trialSolution, RowVector(population, candidate));
528  for (int i = 0; (rnd_unif(0.0, 1.0) < probability) && (i < nDim); i++)
529  {
530  trialSolution[n] = bestSolution[n]
531  + scale * (Element(population, r1, n)
532  - Element(population, r2, n));
533  n = (n + 1) % nDim;
534  }
535 
536  return;
537 }
538 
539 void DESolver::Rand1Exp(int candidate)
540 {
541  int r1;
542  int r2;
543  int r3;
544  int n;
545 
546  SelectSamples(candidate, &r1, &r2, &r3);
547  n = (int)rnd_unif(0.0, (double)nDim);
548 
549  CopyVector(trialSolution, RowVector(population, candidate));
550  for (int i = 0; (rnd_unif(0.0, 1.0) < probability) && (i < nDim); i++)
551  {
552  trialSolution[n] = Element(population, r1, n)
553  + scale * (Element(population, r2, n)
554  - Element(population, r3, n));
555  n = (n + 1) % nDim;
556  }
557 
558  return;
559 }
560 
561 void DESolver::RandToBest1Exp(int candidate)
562 {
563  int r1;
564  int r2;
565  int n;
566 
567  SelectSamples(candidate, &r1, &r2);
568  n = (int)rnd_unif(0.0, (double)nDim);
569 
570  CopyVector(trialSolution, RowVector(population, candidate));
571  for (int i = 0; (rnd_unif(0.0, 1.0) < probability) && (i < nDim); i++)
572  {
573  trialSolution[n] += scale * (bestSolution[n] - trialSolution[n])
574  + scale * (Element(population, r1, n)
575  - Element(population, r2, n));
576  n = (n + 1) % nDim;
577  }
578 
579  return;
580 }
581 
582 void DESolver::Best2Exp(int candidate)
583 {
584  int r1;
585  int r2;
586  int r3;
587  int r4;
588  int n;
589 
590  SelectSamples(candidate, &r1, &r2, &r3, &r4);
591  n = (int)rnd_unif(0.0, (double)nDim);
592 
593  CopyVector(trialSolution, RowVector(population, candidate));
594  for (int i = 0; (rnd_unif(0.0, 1.0) < probability) && (i < nDim); i++)
595  {
596  trialSolution[n] = bestSolution[n] +
597  scale * (Element(population, r1, n)
598  + Element(population, r2, n)
599  - Element(population, r3, n)
600  - Element(population, r4, n));
601  n = (n + 1) % nDim;
602  }
603 
604  return;
605 }
606 
607 void DESolver::Rand2Exp(int candidate)
608 {
609  int r1;
610  int r2;
611  int r3;
612  int r4;
613  int r5;
614  int n;
615 
616  SelectSamples(candidate, &r1, &r2, &r3, &r4, &r5);
617  n = (int)rnd_unif(0.0, (double)nDim);
618 
619  CopyVector(trialSolution, RowVector(population, candidate));
620  for (int i = 0; (rnd_unif(0.0, 1.0) < probability) && (i < nDim); i++)
621  {
622  trialSolution[n] = Element(population, r1, n)
623  + scale * (Element(population, r2, n)
624  + Element(population, r3, n)
625  - Element(population, r4, n)
626  - Element(population, r5, n));
627  n = (n + 1) % nDim;
628  }
629 
630  return;
631 }
632 
633 void DESolver::Best1Bin(int candidate)
634 {
635  int r1;
636  int r2;
637  int n;
638 
639  SelectSamples(candidate, &r1, &r2);
640  n = (int)rnd_unif(0.0, (double)nDim);
641 
642  CopyVector(trialSolution, RowVector(population, candidate));
643  for (int i = 0; i < nDim; i++)
644  {
645  if ((rnd_unif(0.0, 1.0) < probability) || (i == (nDim - 1)))
646  trialSolution[n] = bestSolution[n]
647  + scale * (Element(population, r1, n)
648  - Element(population, r2, n));
649  n = (n + 1) % nDim;
650  }
651 
652  return;
653 }
654 
655 void DESolver::Rand1Bin(int candidate)
656 {
657  int r1;
658  int r2;
659  int r3;
660  int n;
661 
662  SelectSamples(candidate, &r1, &r2, &r3);
663  n = (int)rnd_unif(0.0, (double)nDim);
664 
665  CopyVector(trialSolution, RowVector(population, candidate));
666  for (int i = 0; i < nDim; i++)
667  {
668  if ((rnd_unif(0.0, 1.0) < probability) || (i == (nDim - 1)))
669  trialSolution[n] = Element(population, r1, n)
670  + scale * (Element(population, r2, n)
671  - Element(population, r3, n));
672  n = (n + 1) % nDim;
673  }
674 
675  return;
676 }
677 
678 void DESolver::RandToBest1Bin(int candidate)
679 {
680  int r1;
681  int r2;
682  int n;
683 
684  SelectSamples(candidate, &r1, &r2);
685  n = (int)rnd_unif(0.0, (double)nDim);
686 
687  CopyVector(trialSolution, RowVector(population, candidate));
688  for (int i = 0; i < nDim; i++)
689  {
690  if ((rnd_unif(0.0, 1.0) < probability) || (i == (nDim - 1)))
691  trialSolution[n] += scale * (bestSolution[n] - trialSolution[n])
692  + scale * (Element(population, r1, n)
693  - Element(population, r2, n));
694  n = (n + 1) % nDim;
695  }
696 
697  return;
698 }
699 
700 void DESolver::Best2Bin(int candidate)
701 {
702  int r1;
703  int r2;
704  int r3;
705  int r4;
706  int n;
707 
708  SelectSamples(candidate, &r1, &r2, &r3, &r4);
709  n = (int)rnd_unif(0.0, (double)nDim);
710 
711  CopyVector(trialSolution, RowVector(population, candidate));
712  for (int i = 0; i < nDim; i++)
713  {
714  if ((rnd_unif(0.0, 1.0) < probability) || (i == (nDim - 1)))
715  trialSolution[n] = bestSolution[n]
716  + scale * (Element(population, r1, n)
717  + Element(population, r2, n)
718  - Element(population, r3, n)
719  - Element(population, r4, n));
720  n = (n + 1) % nDim;
721  }
722 
723  return;
724 }
725 
726 void DESolver::Rand2Bin(int candidate)
727 {
728  int r1, r2, r3, r4, r5;
729  int n;
730 
731  SelectSamples(candidate, &r1, &r2, &r3, &r4, &r5);
732  n = (int)rnd_unif(0.0, (double)nDim);
733 
734  CopyVector(trialSolution, RowVector(population, candidate));
735  for (int i = 0; i < nDim; i++)
736  {
737  if ((rnd_unif(0.0, 1.0) < probability) || (i == (nDim - 1)))
738  trialSolution[n] = Element(population, r1, n)
739  + scale * (Element(population, r2, n)
740  + Element(population, r3, n)
741  - Element(population, r4, n)
742  - Element(population, r5, n));
743  n = (n + 1) % nDim;
744  }
745 
746  return;
747 }
748 
749 void DESolver::SelectSamples(int candidate, int *r1, int *r2,
750  int *r3, int *r4, int *r5)
751 {
752  if (r1)
753  {
754  do
755  {
756  *r1 = (int)rnd_unif(0.0, (double)nPop);
757  }
758  while (*r1 == candidate);
759  } else
760  return;
761 
762  if (r2)
763  {
764  do
765  {
766  *r2 = (int)rnd_unif(0.0, (double)nPop);
767  }
768  while ((*r2 == candidate) || (*r2 == *r1));
769  } else
770  return;
771 
772  if (r3)
773  {
774  do
775  {
776  *r3 = (int)rnd_unif(0.0, (double)nPop);
777  }
778  while ((*r3 == candidate) || (*r3 == *r2) || (*r3 == *r1));
779  } else
780  return;
781 
782  if (r4)
783  {
784  do
785  {
786  *r4 = (int)rnd_unif(0.0, (double)nPop);
787  }
788  while ((*r4 == candidate) || (*r4 == *r3) || (*r4 == *r2) || (*r4 == *r1));
789  } else
790  return;
791 
792  if (r5)
793  {
794  do
795  {
796  *r5 = (int)rnd_unif(0.0, (double)nPop);
797  }
798  while ((*r5 == candidate) || (*r5 == *r4) || (*r5 == *r3)
799  || (*r5 == *r2) || (*r5 == *r1));
800  }
801 }
802 
803 /* Check Randomness ------------------------------------------------------- */
804 /* See http://home.ubalt.edu/ntsbarsh/Business-stat/opre504.htm for the formulas */
805 /* This is called runs test */
806 double checkRandomness(const std::string &sequence)
807 {
808  int imax=sequence.size();
809  if (imax<=1)
810  return 0;
811  double n0=1;
812  double n1=0;
813  double R=1;
814  int current=0;
815  for (int i=1; i<imax; ++i)
816  {
817  if (sequence[i]!=sequence[i-1])
818  {
819  current=(current+1)%2;
820  R++;
821  if (current==0)
822  n0++;
823  else
824  n1++;
825  }
826  else
827  {
828  if (current==0)
829  n0++;
830  else
831  n1++;
832  }
833  }
834  double m=1+2*n0*n1/(n0+n1);
835  double s=sqrt(2*n0*n1*(2*n0*n1-n0-n1)/((n0+n1)*(n0+n1)*(n0+n1-1)));
836  double z=(R-m)/s;
837  return ABS(z);
838 }
839 
840 #ifdef NEVERDEFINED
841 double ZernikeSphericalHarmonics(int l1, int n, int l2, int m, double xr, double yr, double zr, double r)
842 {
843  // General variables
844  double r2=r*r,xr2=xr*xr,yr2=yr*yr,zr2=zr*zr;
845 
846  //Variables needed for l>=5
847  double tht=0.0,phi=0.0,costh=0.0,sinth=0.0,costh2=0.0,sinth2=0.0;
848  if (l2>=5)
849  {
850  tht = atan2(yr,xr);
851  phi = atan2(zr,sqrt(xr2 + yr2));
852  sinth = sin(phi); costh = cos(tht);
853  sinth2 = sinth*sinth; costh2 = costh*costh;
854  }
855 
856  // Zernike polynomial
857  double R=0.0;
858 
859  switch (l1)
860  {
861  case 0:
862  R = 1.0;
863  break;
864  case 1:
865  R = r;
866  break;
867  case 2:
868  switch (n)
869  {
870  case 0:
871  R = -1+2*r2;
872  break;
873  case 2:
874  R = r2;
875  break;
876  } break;
877  case 3:
878  switch (n)
879  {
880  case 1:
881  R = 3*r2*r-2*r;
882  break;
883  case 3:
884  R = r2*r;
885  } break;
886  case 4:
887  switch (n)
888  {
889  case 0:
890  R = 6*r2*r2-6*r2+1;
891  break;
892  case 2:
893  R = 4*r2*r2-3*r2;
894  break;
895  case 4:
896  R = r2*r2;
897  break;
898  } break;
899  case 5:
900  switch (n)
901  {
902  case 1:
903  R = 10.0*r2*r2*r-12.0*r2*r+3.0*r;;
904  break;
905  case 3:
906  R = 5.0*r2*r2*r-4.0*r2*r;
907  break;
908  case 5:
909  R = r2*r2*r;
910  break;
911  }break;
912  }
913 
914  // Spherical harmonic
915  double Y=0.0;
916 
917  switch (l2)
918  {
919  case 0:
920  Y = (1.0/2.0)*sqrt(1.0/PI);
921  break;
922  case 1:
923  switch (m)
924  {
925  case -1:
926  Y = sqrt(3.0/(4.0*PI))*yr;
927  break;
928  case 0:
929  Y = sqrt(3.0/(4.0*PI))*zr;
930  break;
931  case 1:
932  Y = sqrt(3.0/(4.0*PI))*xr;
933  break;
934  } break;
935  case 2:
936  switch (m)
937  {
938  case -2:
939  Y = sqrt(15.0/(4.0*PI))*xr*yr;
940  break;
941  case -1:
942  Y = sqrt(15.0/(4.0*PI))*zr*yr;
943  break;
944  case 0:
945  Y = sqrt(5.0/(16.0*PI))*(-xr2-yr2+2.0*zr2);
946  break;
947  case 1:
948  Y = sqrt(15.0/(4.0*PI))*xr*zr;
949  break;
950  case 2:
951  Y = sqrt(15.0/(16.0*PI))*(xr2-yr2);
952  break;
953  } break;
954  case 3:
955  switch (m)
956  {
957  case -3:
958  Y = sqrt(35.0/(16.0*2.0*PI))*yr*(3.0*xr2-yr2);
959  break;
960  case -2:
961  Y = sqrt(105.0/(4.0*PI))*zr*yr*xr;
962  break;
963  case -1:
964  Y = sqrt(21.0/(16.0*2.0*PI))*yr*(4.0*zr2-xr2-yr2);
965  break;
966  case 0:
967  Y = sqrt(7.0/(16.0*PI))*zr*(2.0*zr2-3.0*xr2-3.0*yr2);
968  break;
969  case 1:
970  Y = sqrt(21.0/(16.0*2.0*PI))*xr*(4.0*zr2-xr2-yr2);
971  break;
972  case 2:
973  Y = sqrt(105.0/(16.0*PI))*zr*(xr2-yr2);
974  break;
975  case 3:
976  Y = sqrt(35.0/(16.0*2.0*PI))*xr*(xr2-3.0*yr2);
977  break;
978  } break;
979  case 4:
980  switch (m)
981  {
982  case -4:
983  Y = sqrt((35.0*9.0)/(16.0*PI))*yr*xr*(xr2-yr2);
984  break;
985  case -3:
986  Y = sqrt((9.0*35.0)/(16.0*2.0*PI))*yr*zr*(3.0*xr2-yr2);
987  break;
988  case -2:
989  Y = sqrt((9.0*5.0)/(16.0*PI))*yr*xr*(7.0*zr2-(xr2+yr2+zr2));
990  break;
991  case -1:
992  Y = sqrt((9.0*5.0)/(16.0*2.0*PI))*yr*zr*(7.0*zr2-3.0*(xr2+yr2+zr2));
993  break;
994  case 0:
995  Y = sqrt(9.0/(16.0*16.0*PI))*(35.0*zr2*zr2-30.0*zr2+3.0);
996  break;
997  case 1:
998  Y = sqrt((9.0*5.0)/(16.0*2.0*PI))*xr*zr*(7.0*zr2-3.0*(xr2+yr2+zr2));
999  break;
1000  case 2:
1001  Y = sqrt((9.0*5.0)/(8.0*8.0*PI))*(xr2-yr2)*(7.0*zr2-(xr2+yr2+zr2));
1002  break;
1003  case 3:
1004  Y = sqrt((9.0*35.0)/(16.0*2.0*PI))*xr*zr*(xr2-3.0*yr2);
1005  break;
1006  case 4:
1007  Y = sqrt((9.0*35.0)/(16.0*16.0*PI))*(xr2*(xr2-3.0*yr2)-yr2*(3.0*xr2-yr2));
1008  break;
1009  } break;
1010  case 5:
1011  switch (m)
1012  {
1013  case -5:
1014  Y = (3.0/16.0)*sqrt(77.0/(2.0*PI))*sint2*sint2*sinth*sin(5.0*phi);
1015  break;
1016  case -4:
1017  Y = (3.0/8.0)*sqrt(385.0/(2.0*PI))*sint2*sint2*sin(4.0*phi);
1018  break;
1019  case -3:
1020  Y = (1.0/16.0)*sqrt(385.0/(2.0*PI))*sint2*sinth*(9.0*cost2-1.0)*sin(3.0*phi);
1021  break;
1022  case -2:
1023  Y = (1.0/4.0)*sqrt(1155.0/(4.0*PI))*sint2*(3.0*cost2*costh-costh)*sin(2.0*phi);
1024  break;
1025  case -1:
1026  Y = (1.0/8.0)*sqrt(165.0/(4.0*PI))*sinth*(21.0*cost2*cost2-14.0*cost2+1)*sin(phi);
1027  break;
1028  case 0:
1029  Y = (1.0/16.0)*sqrt(11.0/PI)*(63.0*cost2*cost2*costh-70.0*cost2*costh+15.0*costh);
1030  break;
1031  case 1:
1032  Y = (1.0/8.0)*sqrt(165.0/(4.0*PI))*sinth*(21.0*cost2*cost2-14.0*cost2+1)*cos(phi);
1033  break;
1034  case 2:
1035  Y = (1.0/4.0)*sqrt(1155.0/(4.0*PI))*sint2*(3.0*cost2*costh-costh)*cos(2.0*phi);
1036  break;
1037  case 3:
1038  Y = (1.0/16.0)*sqrt(385.0/(2.0*PI))*sint2*sinth*(9.0*cost2-1.0)*cos(3.0*phi);
1039  break;
1040  case 4:
1041  Y = (3.0/8.0)*sqrt(385.0/(2.0*PI))*sint2*sint2*cos(4.0*phi);
1042  break;
1043  case 5:
1044  Y = (3.0/16.0)*sqrt(77.0/(2.0*PI))*sint2*sint2*sinth*cos(5.0*phi);
1045  break;
1046  }break;
1047  }
1048 
1049  return R*Y;
1050 }
1051 #endif
1052 
1053 template<int L1, int L2>
1054 double ZernikeSphericalHarmonics(int n, int m, double xr, double yr, double zr, double r)
1055 {
1056  // General variables
1057  double r2=r*r;
1058  double xr2=xr*xr;
1059  double yr2=yr*yr;
1060  double zr2=zr*zr;
1061 
1062  //Variables needed for l>=5
1063  double tht=0.0;
1064  double phi=0.0;
1065  double costh=0.0;
1066  double sinth=0.0;
1067  double costh2=0.0;
1068  double sinth2=0.0;
1069  double cosph=0.0;
1070  double cosph2=0.0;
1071  if (L2>=5)
1072  {
1073  tht = atan2(yr,xr);
1074  phi = atan2(zr,sqrt(xr2 + yr2));
1075  sinth = sin(abs(m)*phi); costh = cos(tht); cosph = cos(abs(m)*phi);
1076  sinth2 = sinth*sinth; costh2 = costh*costh; cosph2 = cosph * cosph;
1077  }
1078 
1079  // Zernike polynomial
1080  double R=0.0;
1081 
1082  switch (L1)
1083  {
1084  case 0:
1085  R = std::sqrt(3.0);
1086  break;
1087  case 1:
1088  R = std::sqrt(5.0)*r;
1089  break;
1090  case 2:
1091  switch (n)
1092  {
1093  case 0:
1094  R = -0.5*std::sqrt(7.0)*(2.5*(1-2*r2)+0.5);
1095  break;
1096  case 2:
1097  R = std::sqrt(7.0)*r2;
1098  break;
1099  } break;
1100  case 3:
1101  switch (n)
1102  {
1103  case 1:
1104  R = -1.5*r*(3.5*(1-2*r2)+1.5);
1105  break;
1106  case 3:
1107  R = 3.0*r2*r;
1108  } break;
1109  case 4:
1110  switch (n)
1111  {
1112  case 0:
1113  R = std::sqrt(11.0)*((63.0*r2*r2/8.0)-(35.0*r2/4.0)+(15.0/8.0));
1114  break;
1115  case 2:
1116  R = -0.5*std::sqrt(11.0)*r2*(4.5*(1-2*r2)+2.5);
1117  break;
1118  case 4:
1119  R = std::sqrt(11.0)*r2*r2;
1120  break;
1121  } break;
1122  case 5:
1123  switch (n)
1124  {
1125  case 1:
1126  R = std::sqrt(13.0)*r*((99.0*r2*r2/8.0)-(63.0*r2/4.0)+(35.0/8.0));
1127  break;
1128  case 3:
1129  R = -0.5*std::sqrt(13.0)*r2*r*(5.5*(1-2*r2)+3.5);
1130  break;
1131  case 5:
1132  R = std::sqrt(13.0)*r2*r2*r;
1133  break;
1134  } break;
1135  case 6:
1136  switch (n)
1137  {
1138  case 0:
1139  R = 103.8 * r2 * r2 * r2 - 167.7 * r2 * r2 + 76.25 * r2 - 8.472;
1140  break;
1141  case 2:
1142  R = 69.23 * r2 * r2 * r2 - 95.86 * r2 * r2 + 30.5 * r2;
1143  break;
1144  case 4:
1145  R = 25.17 * r2 * r2 * r2 - 21.3 * r2 * r2;
1146  break;
1147  case 6:
1148  R = 3.873 * r2 * r2 * r2;
1149  break;
1150  } break;
1151  case 7:
1152  switch (n)
1153  {
1154  case 1:
1155  R = 184.3 * r2 * r2 * r2 * r - 331.7 * r2 * r2 * r + 178.6 * r2 * r - 27.06 * r;
1156  break;
1157  case 3:
1158  R = 100.5 * r2 * r2 * r2 * r - 147.4 * r2 * r2 * r + 51.02 * r2 * r;
1159  break;
1160  case 5:
1161  R = 30.92 * r2 * r2 * r2 * r - 26.8 * r2 * r2 * r;
1162  break;
1163  case 7:
1164  R = 4.123 * r2 * r2 * r2 * r;
1165  break;
1166  } break;
1167  case 8:
1168  switch (n)
1169  {
1170  case 0:
1171  R = 413.9*r2*r2*r2*r2 - 876.5*r2*r2*r2 + 613.6*r2*r2 - 157.3*r2 + 10.73;
1172  break;
1173  case 2:
1174  R = 301.0*r2*r2*r2*r2 - 584.4*r2*r2*r2 + 350.6*r2*r2 - 62.93*r2;
1175  break;
1176  case 4:
1177  R = 138.9*r2*r2*r2*r2 - 212.5*r2*r2*r2 + 77.92*r2*r2;
1178  break;
1179  case 6:
1180  R = 37.05*r2*r2*r2*r2 - 32.69*r2*r2*r2;
1181  break;
1182  case 8:
1183  R = 4.359*r2*r2*r2*r2;
1184  break;
1185  } break;
1186  case 9:
1187  switch (n)
1188  {
1189  case 1:
1190  R = 751.6*r2*r2*r2*r2*r - 1741.0*r2*r2*r2*r + 1382.0*r2*r2*r - 430.0*r2*r + 41.35*r;
1191  break;
1192  case 3:
1193  R = 462.6*r2*r2*r2*r2*r - 949.5*r2*r2*r2*r + 614.4*r2*r2*r - 122.9*r2*r;
1194  break;
1195  case 5:
1196  R = 185.0*r2*r2*r2*r2*r - 292.1*r2*r2*r2*r + 111.7*r2*r2*r;
1197  break;
1198  case 7:
1199  R = 43.53*r2*r2*r2*r2*r - 38.95*r2*r2*r2*r;
1200  break;
1201  case 9:
1202  R = 4.583*r2*r2*r2*r2*r;
1203  break;
1204  } break;
1205  case 10:
1206  switch (n)
1207  {
1208  case 0:
1209  R = 1652.0*r2*r2*r2*r2*r2 - 4326.0*r2*r2*r2*r2 + 4099.0*r2*r2*r2 - 1688.0*r2*r2 + 281.3*r2 - 12.98;
1210  break;
1211  case 2:
1212  R = 1271.0*r2*r2*r2*r2*r2 - 3147.0*r2*r2*r2*r2 + 2732.0*r2*r2*r2 - 964.4*r2*r2 + 112.5*r2;
1213  break;
1214  case 4:
1215  R = 677.7*r2*r2*r2*r2*r2 - 1452.0*r2*r2*r2*r2 + 993.6*r2*r2*r2 - 214.3*r2*r2;
1216  break;
1217  case 6:
1218  R = 239.2*r2*r2*r2*r2*r2 - 387.3*r2*r2*r2*r2 + 152.9*r2*r2*r2;
1219  break;
1220  case 8:
1221  R = 50.36*r2*r2*r2*r2*r2 - 45.56*r2*r2*r2*r2;
1222  break;
1223  case 10:
1224  R = 4.796*r2*r2*r2*r2*r2;
1225  break;
1226  } break;
1227  case 11:
1228  switch (n)
1229  {
1230  case 1:
1231  R = r*-5.865234375E+1+(r*r*r)*8.7978515625E+2-(r*r*r*r*r)*4.2732421875E+3+(r*r*r*r*r*r*r)*9.0212890625E+3-(r*r*r*r*r*r*r*r*r)*8.61123046875E+3+pow(r,1.1E+1)*3.04705078125E+3;
1232  break;
1233  case 3:
1234  R = (r*r*r)*2.513671875E+2-(r*r*r*r*r)*1.89921875E+3+(r*r*r*r*r*r*r)*4.920703125E+3-(r*r*r*r*r*r*r*r*r)*5.29921875E+3+pow(r,1.1E+1)*2.0313671875E+3;
1235  break;
1236  case 5:
1237  R = (r*r*r*r*r)*-3.453125E+2+(r*r*r*r*r*r*r)*1.5140625E+3-(r*r*r*r*r*r*r*r*r)*2.1196875E+3+pow(r,1.1E+1)*9.559375E+2;
1238  break;
1239  case 7:
1240  R = (r*r*r*r*r*r*r)*2.01875E+2-(r*r*r*r*r*r*r*r*r)*4.9875E+2+pow(r,1.1E+1)*3.01875E+2;
1241  break;
1242  case 9:
1243  R = (r*r*r*r*r*r*r*r*r)*-5.25E+1+pow(r,1.1E+1)*5.75E+1;
1244  break;
1245  case 11:
1246  R = pow(r,1.1E+1)*5.0;
1247  break;
1248  } break;
1249  case 12:
1250  switch (n)
1251  {
1252  case 0:
1253  R = (r*r)*-4.57149777110666E+2+(r*r*r*r)*3.885773105442524E+3-(r*r*r*r*r*r)*1.40627979054153E+4+(r*r*r*r*r*r*r*r)*2.460989633446932E+4-pow(r,1.0E+1)*2.05828223888278E+4+pow(r,1.2E+1)*6.597058457955718E+3+1.523832590368693E+1;
1254  break;
1255  case 2:
1256  R = (r*r)*-1.828599108443595E+2+(r*r*r*r)*2.220441774539649E+3-(r*r*r*r*r*r)*9.375198603600264E+3+(r*r*r*r*r*r*r*r)*1.789810642504692E+4-pow(r,1.0E+1)*1.583294029909372E+4+pow(r,1.2E+1)*5.277646766364574E+3;
1257  break;
1258  case 4:
1259  R = (r*r*r*r)*4.934315054528415E+2-(r*r*r*r*r*r)*3.409163128584623E+3+(r*r*r*r*r*r*r*r)*8.260664503872395E+3-pow(r,1.0E+1)*8.444234826177359E+3+pow(r,1.2E+1)*3.104498097866774E+3;
1260  break;
1261  case 6:
1262  R = (r*r*r*r*r*r)*-5.244866351671517E+2+(r*r*r*r*r*r*r*r)*2.202843867704272E+3-pow(r,1.0E+1)*2.98031817394495E+3+pow(r,1.2E+1)*1.307157093837857E+3;
1263  break;
1264  case 8:
1265  R = (r*r*r*r*r*r*r*r)*2.591581020820886E+2-pow(r,1.0E+1)*6.274354050420225E+2+pow(r,1.2E+1)*3.734734553815797E+2;
1266  break;
1267  case 10:
1268  R = pow(r,1.0E+1)*-5.975575286115054E+1+pow(r,1.2E+1)*6.49519052838441E+1;
1269  break;
1270  case 12:
1271  R = pow(r,1.2E+1)*5.19615242270811;
1272  break;
1273  } break;
1274  case 13:
1275  switch (n)
1276  {
1277  case 1:
1278  R = r*7.896313435467891E+1-(r*r*r)*1.610847940832376E+3+(r*r*r*r*r)*1.093075388422608E+4-(r*r*r*r*r*r*r)*3.400678986203671E+4+(r*r*r*r*r*r*r*r*r)*5.332882955634594E+4-pow(r,1.1E+1)*4.102217658185959E+4+pow(r,1.3E+1)*1.230665297454596E+4;
1279  break;
1280  case 3:
1281  R = (r*r*r)*-4.602422688100487E+2+(r*r*r*r*r)*4.858112837433815E+3-(r*r*r*r*r*r*r)*1.854915810656548E+4+(r*r*r*r*r*r*r*r*r)*3.281774126553535E+4-pow(r,1.1E+1)*2.734811772125959E+4+pow(r,1.3E+1)*8.687049158513546E+3;
1282  break;
1283  case 5:
1284  R = (r*r*r*r*r)*8.832932431697845E+2-(r*r*r*r*r*r*r)*5.707433263555169E+3+(r*r*r*r*r*r*r*r*r)*1.312709650617838E+4-pow(r,1.1E+1)*1.286970245704055E+4+pow(r,1.3E+1)*4.572131136059761E+3;
1285  break;
1286  case 7:
1287  R = (r*r*r*r*r*r*r)*-7.60991101808846E+2+(r*r*r*r*r*r*r*r*r)*3.088728589691222E+3-pow(r,1.1E+1)*4.064116565383971E+3+pow(r,1.3E+1)*1.741764242306352E+3;
1288  break;
1289  case 9:
1290  R = (r*r*r*r*r*r*r*r*r)*3.251293252306059E+2-pow(r,1.1E+1)*7.741174410264939E+2+pow(r,1.3E+1)*4.54373280601576E+2;
1291  break;
1292  case 11:
1293  R = pow(r,1.1E+1)*-6.731456008902751E+1+pow(r,1.3E+1)*7.269972489634529E+1;
1294  break;
1295  case 13:
1296  R = pow(r,1.3E+1)*5.385164807128604;
1297  break;
1298  } break;
1299  case 14:
1300  switch (n)
1301  {
1302  case 0:
1303  R = (r*r)*6.939451623205096E+2-(r*r*r*r)*7.910974850460887E+3+(r*r*r*r*r*r)*3.955487425231934E+4-(r*r*r*r*r*r*r*r)*1.010846786448956E+5+pow(r,1.0E+1)*1.378427436065674E+5-pow(r,1.2E+1)*9.542959172773361E+4+pow(r,1.4E+1)*2.63567443819046E+4-1.749441585683962E+1;
1304  break;
1305  case 2:
1306  R = (r*r)*2.775780649287626E+2-(r*r*r*r)*4.520557057410479E+3+(r*r*r*r*r*r)*2.636991616821289E+4-(r*r*r*r*r*r*r*r)*7.351612992358208E+4+pow(r,1.0E+1)*1.060328796973228E+5-pow(r,1.2E+1)*7.634367338204384E+4+pow(r,1.4E+1)*2.170555419689417E+4;
1307  break;
1308  case 4:
1309  R = (r*r*r*r)*-1.004568234980106E+3+(r*r*r*r*r*r)*9.589060424804688E+3-(r*r*r*r*r*r*r*r)*3.393052150321007E+4+pow(r,1.0E+1)*5.655086917197704E+4-pow(r,1.2E+1)*4.490804316592216E+4+pow(r,1.4E+1)*1.370877107170224E+4;
1310  break;
1311  case 6:
1312  R = (r*r*r*r*r*r)*1.475240065354854E+3-(r*r*r*r*r*r*r*r)*9.04813906750083E+3+pow(r,1.0E+1)*1.99591302959919E+4-pow(r,1.2E+1)*1.8908649754107E+4+pow(r,1.4E+1)*6.527986224621534E+3;
1313  break;
1314  case 8:
1315  R = (r*r*r*r*r*r*r*r)*-1.064486949119717E+3+pow(r,1.0E+1)*4.201922167569399E+3-pow(r,1.2E+1)*5.402471358314157E+3+pow(r,1.4E+1)*2.270603904217482E+3;
1316  break;
1317  case 10:
1318  R = pow(r,1.0E+1)*4.001830635787919E+2-pow(r,1.2E+1)*9.395602362267673E+2+pow(r,1.4E+1)*5.44944937011227E+2;
1319  break;
1320  case 12:
1321  R = pow(r,1.2E+1)*-7.516481889830902E+1+pow(r,1.4E+1)*8.073258326109499E+1;
1322  break;
1323  case 14:
1324  R = pow(r,1.4E+1)*5.567764362829621;
1325  break;
1326  } break;
1327  case 15:
1328  switch (n)
1329  {
1330  case 1:
1331  R = r*-1.022829477079213E+2+(r*r*r)*2.720726409032941E+3-(r*r*r*r*r)*2.448653768128157E+4+(r*r*r*r*r*r*r)*1.042945123462677E+5-(r*r*r*r*r*r*r*r*r)*2.370329826049805E+5+pow(r,1.1E+1)*2.953795629386902E+5-pow(r,1.3E+1)*1.903557183384895E+5+pow(r,1.5E+1)*4.958846444106102E+4;
1332  break;
1333  case 3:
1334  R = (r*r*r)*7.773504025805742E+2-(r*r*r*r*r)*1.088290563613176E+4+(r*r*r*r*r*r*r)*5.688791582524776E+4-(r*r*r*r*r*r*r*r*r)*1.458664508337975E+5+pow(r,1.1E+1)*1.969197086257935E+5-pow(r,1.3E+1)*1.343687423563004E+5+pow(r,1.5E+1)*3.653886853551865E+4;
1335  break;
1336  case 5:
1337  R = (r*r*r*r*r)*-1.978710115659982E+3+(r*r*r*r*r*r*r)*1.750397410005331E+4-(r*r*r*r*r*r*r*r*r)*5.834658033359051E+4+pow(r,1.1E+1)*9.266809817695618E+4-pow(r,1.3E+1)*7.072039071393013E+4+pow(r,1.5E+1)*2.08793534488678E+4;
1338  break;
1339  case 7:
1340  R = (r*r*r*r*r*r*r)*2.333863213345408E+3-(r*r*r*r*r*r*r*r*r)*1.372860713732243E+4+pow(r,1.1E+1)*2.926360995060205E+4-pow(r,1.3E+1)*2.694110122436285E+4+pow(r,1.5E+1)*9.077979760378599E+3;
1341  break;
1342  case 9:
1343  R = (r*r*r*r*r*r*r*r*r)*-1.445116540770978E+3+pow(r,1.1E+1)*5.57402094297111E+3-pow(r,1.3E+1)*7.028113362878561E+3+pow(r,1.5E+1)*2.90495352332294E+3;
1344  break;
1345  case 11:
1346  R = pow(r,1.1E+1)*4.846974733015522E+2-pow(r,1.3E+1)*1.124498138058931E+3+pow(r,1.5E+1)*6.455452274046838E+2;
1347  break;
1348  case 13:
1349  R = pow(r,1.3E+1)*-8.329615837475285E+1+pow(r,1.5E+1)*8.904072102135979E+1;
1350  break;
1351  case 15:
1352  R = pow(r,1.5E+1)*5.744562646534177;
1353  break;
1354  } break;
1355  default: break;
1356  }
1357 
1358  // Spherical harmonic
1359  double Y=0.0;
1360 
1361  switch (L2)
1362  {
1363  case 0:
1364  Y = (1.0/2.0)*sqrt(1.0/PI);
1365  break;
1366  case 1:
1367  switch (m)
1368  {
1369  case -1:
1370  Y = sqrt(3.0/(4.0*PI))*yr;
1371  break;
1372  case 0:
1373  Y = sqrt(3.0/(4.0*PI))*zr;
1374  break;
1375  case 1:
1376  Y = sqrt(3.0/(4.0*PI))*xr;
1377  break;
1378  } break;
1379  case 2:
1380  switch (m)
1381  {
1382  case -2:
1383  Y = sqrt(15.0/(4.0*PI))*xr*yr;
1384  break;
1385  case -1:
1386  Y = sqrt(15.0/(4.0*PI))*zr*yr;
1387  break;
1388  case 0:
1389  Y = sqrt(5.0/(16.0*PI))*(-xr2-yr2+2.0*zr2);
1390  break;
1391  case 1:
1392  Y = sqrt(15.0/(4.0*PI))*xr*zr;
1393  break;
1394  case 2:
1395  Y = sqrt(15.0/(16.0*PI))*(xr2-yr2);
1396  break;
1397  } break;
1398  case 3:
1399  switch (m)
1400  {
1401  case -3:
1402  Y = sqrt(35.0/(16.0*2.0*PI))*yr*(3.0*xr2-yr2);
1403  break;
1404  case -2:
1405  Y = sqrt(105.0/(4.0*PI))*zr*yr*xr;
1406  break;
1407  case -1:
1408  Y = sqrt(21.0/(16.0*2.0*PI))*yr*(4.0*zr2-xr2-yr2);
1409  break;
1410  case 0:
1411  Y = sqrt(7.0/(16.0*PI))*zr*(2.0*zr2-3.0*xr2-3.0*yr2);
1412  break;
1413  case 1:
1414  Y = sqrt(21.0/(16.0*2.0*PI))*xr*(4.0*zr2-xr2-yr2);
1415  break;
1416  case 2:
1417  Y = sqrt(105.0/(16.0*PI))*zr*(xr2-yr2);
1418  break;
1419  case 3:
1420  Y = sqrt(35.0/(16.0*2.0*PI))*xr*(xr2-3.0*yr2);
1421  break;
1422  } break;
1423  case 4:
1424  switch (m)
1425  {
1426  case -4:
1427  Y = sqrt((35.0*9.0)/(16.0*PI))*yr*xr*(xr2-yr2);
1428  break;
1429  case -3:
1430  Y = sqrt((9.0*35.0)/(16.0*2.0*PI))*yr*zr*(3.0*xr2-yr2);
1431  break;
1432  case -2:
1433  Y = sqrt((9.0*5.0)/(16.0*PI))*yr*xr*(7.0*zr2-(xr2+yr2+zr2));
1434  break;
1435  case -1:
1436  Y = sqrt((9.0*5.0)/(16.0*2.0*PI))*yr*zr*(7.0*zr2-3.0*(xr2+yr2+zr2));
1437  break;
1438  case 0:
1439  Y = sqrt(9.0/(16.0*16.0*PI))*(35.0*zr2*zr2-30.0*zr2+3.0);
1440  break;
1441  case 1:
1442  Y = sqrt((9.0*5.0)/(16.0*2.0*PI))*xr*zr*(7.0*zr2-3.0*(xr2+yr2+zr2));
1443  break;
1444  case 2:
1445  Y = sqrt((9.0*5.0)/(8.0*8.0*PI))*(xr2-yr2)*(7.0*zr2-(xr2+yr2+zr2));
1446  break;
1447  case 3:
1448  Y = sqrt((9.0*35.0)/(16.0*2.0*PI))*xr*zr*(xr2-3.0*yr2);
1449  break;
1450  case 4:
1451  Y = sqrt((9.0*35.0)/(16.0*16.0*PI))*(xr2*(xr2-3.0*yr2)-yr2*(3.0*xr2-yr2));
1452  break;
1453  } break;
1454  case 5:
1455  switch (m)
1456  {
1457  case -5:
1458  Y = (3.0/16.0)*sqrt(77.0/(2.0*PI))*sinth2*sinth2*sinth*sin(5.0*phi);
1459  break;
1460  case -4:
1461  Y = (3.0/8.0)*sqrt(385.0/(2.0*PI))*sinth2*sinth2*sin(4.0*phi);
1462  break;
1463  case -3:
1464  Y = (1.0/16.0)*sqrt(385.0/(2.0*PI))*sinth2*sinth*(9.0*costh2-1.0)*sin(3.0*phi);
1465  break;
1466  case -2:
1467  Y = (1.0/4.0)*sqrt(1155.0/(4.0*PI))*sinth2*(3.0*costh2*costh-costh)*sin(2.0*phi);
1468  break;
1469  case -1:
1470  Y = (1.0/8.0)*sqrt(165.0/(4.0*PI))*sinth*(21.0*costh2*costh2-14.0*costh2+1)*sin(phi);
1471  break;
1472  case 0:
1473  Y = (1.0/16.0)*sqrt(11.0/PI)*(63.0*costh2*costh2*costh-70.0*costh2*costh+15.0*costh);
1474  break;
1475  case 1:
1476  Y = (1.0/8.0)*sqrt(165.0/(4.0*PI))*sinth*(21.0*costh2*costh2-14.0*costh2+1)*cos(phi);
1477  break;
1478  case 2:
1479  Y = (1.0/4.0)*sqrt(1155.0/(4.0*PI))*sinth2*(3.0*costh2*costh-costh)*cos(2.0*phi);
1480  break;
1481  case 3:
1482  Y = (1.0/16.0)*sqrt(385.0/(2.0*PI))*sinth2*sinth*(9.0*costh2-1.0)*cos(3.0*phi);
1483  break;
1484  case 4:
1485  Y = (3.0/8.0)*sqrt(385.0/(2.0*PI))*sinth2*sinth2*cos(4.0*phi);
1486  break;
1487  case 5:
1488  Y = (3.0/16.0)*sqrt(77.0/(2.0*PI))*sinth2*sinth2*sinth*cos(5.0*phi);
1489  break;
1490  }break;
1491  case 6:
1492  switch (m)
1493  {
1494  case -6:
1495  Y = -0.6832*sinth*pow(costh2 - 1.0, 3);
1496  break;
1497  case -5:
1498  Y = 2.367*costh*sinth*pow(1.0 - 1.0*costh2, 2.5);
1499  break;
1500  case -4:
1501  Y = 0.001068*sinth*(5198.0*costh2 - 472.5)*pow(costh2 - 1.0, 2);
1502  break;
1503  case -3:
1504  Y = -0.005849*sinth*pow(1.0 - 1.0*costh2, 1.5)*(- 1732.0*costh2*costh + 472.5*costh);
1505  break;
1506  case -2:
1507  Y = -0.03509*sinth*(costh2 - 1.0)*(433.1*costh2*costh2 - 236.2*costh2 + 13.12);
1508  break;
1509  case -1:
1510  Y = 0.222*sinth*pow(1.0 - 1.0*costh2, 0.5)*(86.62*costh2*costh2*costh - 78.75*costh2*costh + 13.12*costh);
1511  break;
1512  case 0:
1513  Y = 14.68*costh2*costh2*costh2 - 20.02*costh2*costh2 + 6.675*costh2 - 0.3178;
1514  break;
1515  case 1:
1516  Y = 0.222*cosph*pow(1.0 - 1.0*costh2, 0.5)*(86.62*costh2*costh2*costh - 78.75*costh2*costh + 13.12*costh);
1517  break;
1518  case 2:
1519  Y = -0.03509*cosph*(costh2 - 1.0)*(433.1*costh2*costh2 - 236.2*costh2 + 13.12);
1520  break;
1521  case 3:
1522  Y = -0.005849*cosph*pow(1.0 - 1.0*costh2, 1.5)*(-1732.0*costh2*costh + 472.5*costh);
1523  break;
1524  case 4:
1525  Y = 0.001068*cosph*(5198.0*costh2 - 472.5)*pow(costh2 - 1.0, 2);
1526  break;
1527  case 5:
1528  Y = 2.367*costh*cosph*pow(1.0 - 1.0*costh2, 2.5);
1529  break;
1530  case 6:
1531  Y = -0.6832*cosph*pow(costh2 - 1.0, 3);
1532  break;
1533  }break;
1534  case 7:
1535  switch (m)
1536  {
1537  case -7:
1538  Y = 0.7072*sinth*pow(1.0 - 1.0*costh2, 3.5);
1539  break;
1540  case -6:
1541  Y = -2.646*costh*sinth*pow(costh2 - 1.0, 3);
1542  break;
1543  case -5:
1544  Y = 9.984e-5*sinth*pow(1.0 - 1.0*costh2, 2.5)*(67570.0*costh2 - 5198.0);
1545  break;
1546  case -4:
1547  Y = -0.000599*sinth*pow(costh2 - 1.0, 2)*(-22520.0*costh2*costh + 5198.0*costh);
1548  break;
1549  case -3:
1550  Y = 0.003974*sinth*pow(1.0 - 1.0*costh2, 1.5)*(5631.0*costh2*costh2 - 2599.0*costh2 + 118.1);
1551  break;
1552  case -2:
1553  Y = -0.0281*sinth*(costh2 - 1.0)*(1126.0*costh2*costh2*costh - 866.2*costh2*costh + 118.1*costh);
1554  break;
1555  case -1:
1556  Y = 0.2065*sinth*pow(1.0 - 1.0*costh2, 0.5)*(187.7*costh2*costh2*costh2 - 216.6*costh2*costh2 + 59.06*costh2 - 2.188);
1557  break;
1558  case 0:
1559  Y = 29.29*costh2*costh2*costh2*costh - 47.32*costh2*costh2*costh + 21.51*costh2*costh - 2.39*costh;
1560  break;
1561  case 1:
1562  Y = 0.2065*cosph*pow(1.0 - 1.0*costh2, 0.5)*(187.7*costh2*costh2*costh2 - 216.6*costh2*costh2 + 59.06*costh2 - 2.188);
1563  break;
1564  case 2:
1565  Y = -0.0281*cosph*(costh2 - 1.0)*(1126.0*costh2*costh2*costh - 866.2*costh2*costh + 118.1*costh);
1566  break;
1567  case 3:
1568  Y = 0.003974*cosph*pow(1.0 - 1.0*costh2, 1.5)*(5631.0*costh2*costh2 - 2599.0*costh2 + 118.1);
1569  break;
1570  case 4:
1571  Y = -0.000599*cosph*pow(costh2 - 1.0, 2)*(- 22520.0*costh2*costh + 5198.0*costh);
1572  break;
1573  case 5:
1574  Y = 9.984e-5*cosph*pow(1.0 - 1.0*costh2, 2.5)*(67570.0*costh2 - 5198.0);
1575  break;
1576  case 6:
1577  Y = -2.646*cosph*costh*pow(costh2 - 1.0, 3);
1578  break;
1579  case 7:
1580  Y = 0.7072*cosph*pow(1.0 - 1.0*costh2, 3.5);
1581  break;
1582  }break;
1583  case 8:
1584  switch (m)
1585  {
1586  case -8:
1587  Y = sinth*pow(costh2-1.0,4.0)*7.289266601746931E-1;
1588  break;
1589  case -7:
1590  Y = costh*sinth*pow(costh2*-1.0+1.0,7.0/2.0)*2.915706640698772;
1591  break;
1592  case -6:
1593  Y = sinth*(costh2*1.0135125E+6-6.75675E+4)*pow(costh2-1.0,3.0)*-7.878532816224526E-6;
1594  break;
1595  case -5:
1596  Y = sinth*pow(costh2*-1.0+1.0,5.0/2.0)*(costh*6.75675E+4-(costh2*costh)*3.378375E+5)*-5.105872826582925E-56;
1597  break;
1598  case -4:
1599  Y = sinth*pow(costh2-1.0,2.0)*(costh2*-3.378375E+4+(costh2*costh2)*8.4459375E+4+1.299375E+3)*3.681897256448963E-4;
1600  break;
1601  case -3:
1602  Y = sinth*pow(costh2*-1.0+1.0,3.0/2.0)*(costh*1.299375E+3-(costh*costh2)*1.126125E+4+(costh*costh2*costh2)*1.6891875E+4)*2.851985351334463E-3;
1603  break;
1604  case -2:
1605  Y = sinth*(costh2-1.0)*(costh2*6.496875E+2-(costh2*costh2)*2.8153125E+3+(costh2*costh2*costh2)*2.8153125E+3-1.96875E+1)*-2.316963852365461E-2;
1606  break;
1607  case -1:
1608  Y = sinth*sqrt(costh2*-1.0+1.0)*(costh*1.96875E+1-(costh*costh2)*2.165625E+2+(costh*costh2*costh2)*5.630625E+2-(costh*costh2*costh2*costh2)*4.021875E+2)*-1.938511038201796E-1;;
1609  break;
1610  case 0:
1611  Y = costh2*-1.144933081936324E+1+(costh2*costh2)*6.297131950652692E+1-(costh2*costh2*costh2)*1.091502871445846E+2+(costh2*costh2*costh2*costh2)*5.847336811327841E+1+3.180369672045344E-1;
1612  break;
1613  case 1:
1614  Y = cosph*sqrt(costh2*-1.0+1.0)*(costh*1.96875E+1-(costh*costh2)*2.165625E+2+(costh*costh2*costh2)*5.630625E+2-(costh*costh2*costh2*costh2)*4.021875E+2)*-1.938511038201796E-1;;
1615  break;
1616  case 2:
1617  Y = cosph*(costh2-1.0)*(costh2*6.496875E+2-(costh2*costh2)*2.8153125E+3+(costh2*costh2*costh2)*2.8153125E+3-1.96875E+1)*-2.316963852365461E-2;
1618  break;
1619  case 3:
1620  Y = cosph*pow(costh2*-1.0+1.0,3.0/2.0)*(costh*1.299375E+3-(costh*costh2)*1.126125E+4+(costh*costh2*costh2)*1.6891875E+4)*2.851985351334463E-3;
1621  break;
1622  case 4:
1623  Y = cosph*pow(costh2-1.0,2.0)*(costh2*-3.378375E+4+(costh2*costh2)*8.4459375E+4+1.299375E+3)*3.681897256448963E-4;
1624  break;
1625  case 5:
1626  Y = cosph*pow(costh2*-1.0+1.0,5.0/2.0)*(costh*6.75675E+4-(costh*costh2)*3.378375E+5)*-5.105872826582925E-5;
1627  break;
1628  case 6:
1629  Y = cosph*(costh2*1.0135125E+6-6.75675E+4)*pow(costh2-1.0,3.0)*-7.878532816224526E-6;
1630  break;
1631  case 7:
1632  Y = cosph*costh*pow(costh2*-1.0+1.0,7.0/2.0)*2.915706640698772;
1633  break;
1634  case 8:
1635  Y = cosph*pow(costh2-1.0,4.0)*7.289266601746931E-1;
1636  break;
1637  }break;
1638  case 9:
1639  switch (m)
1640  {
1641  case -9:
1642  Y = sinth*pow(costh2*-1.0+1.0,9.0/2.0)*7.489009518540115E-1;
1643  break;
1644  case -8:
1645  Y = costh*sinth*pow(costh2-1.0,4.0)*3.17731764895143;
1646  break;
1647  case -7:
1648  Y = sinth*pow(costh2*-1.0+1.0,7.0/2.0)*(costh2*1.72297125E+7-1.0135125E+6)*5.376406125665728E-7;
1649  break;
1650  case -6:
1651  Y = sinth*(costh*1.0135125E+6-(costh*costh2)*5.7432375E+6)*pow(costh2-1.0,3.0)*3.724883428715686E-6;
1652  break;
1653  case -5:
1654  Y = sinth*pow(costh2*-1.0+1.0,5.0/2.0)*(costh2*-5.0675625E+5+(costh2*costh2)*1.435809375E+6+1.6891875E+4)*2.885282297193648E-5;
1655  break;
1656  case -4:
1657  Y = sinth*pow(costh2-1.0,2.0)*(costh*1.6891875E+4-(costh*costh2)*1.6891875E+5+(costh*costh2*costh2)*2.87161875E+5)*2.414000363328839E-4;
1658  break;
1659  case -3:
1660  Y = sinth*pow(costh2*-1.0+1.0,3.0/2.0)*((costh2)*8.4459375E+3-(costh2*costh2)*4.22296875E+4+(costh2*costh2*costh2)*4.78603125E+4-2.165625E+2)*2.131987394015766E-3;
1661  break;
1662  case -2:
1663  Y = sinth*(costh2-1.0)*(costh*2.165625E+2-(costh*costh2)*2.8153125E+3+(costh*costh2*costh2)*8.4459375E+3-(costh*costh2*costh2*costh2)*6.8371875E+3)*1.953998722751749E-2;
1664  break;
1665  case -1:
1666  Y = sinth*sqrt(costh2*-1.0+1.0)*(costh2*-1.0828125E+2+(costh2*costh2)*7.03828125E+2-(costh2*costh2*costh2)*1.40765625E+3+(costh2*costh2*costh2*costh2)*8.546484375E+2+2.4609375)*1.833013280775049E-1;
1667  break;
1668  case 0:
1669  Y = costh*3.026024588281871-(costh*costh2)*4.438169396144804E+1+(costh*costh2*costh2)*1.730886064497754E+2-(costh*costh2*costh2*costh2)*2.472694377852604E+2+(costh*costh2*costh2*costh2*costh2)*1.167661233986728E+2;
1670  break;
1671  case 1:
1672  Y = cosph*sqrt(costh2*-1.0+1.0)*(costh2*-1.0828125E+2+(costh2*costh2)*7.03828125E+2-(costh2*costh2*costh2)*1.40765625E+3+(costh2*costh2*costh2*costh2)*8.546484375E+2+2.4609375)*1.833013280775049E-1;
1673  break;
1674  case 2:
1675  Y = cosph*(costh2-1.0)*(costh*2.165625E+2-(costh*costh2)*2.8153125E+3+(costh*costh2*costh2)*8.4459375E+3-(costh*costh2*costh2*costh2)*6.8371875E+3)*1.953998722751749E-2;
1676  break;
1677  case 3:
1678  Y = cosph*pow(costh2*-1.0+1.0,3.0/2.0)*(costh2*8.4459375E+3-(costh2*costh2)*4.22296875E+4+(costh2*costh2*costh2)*4.78603125E+4-2.165625E+2)*2.131987394015766E-3;
1679  break;
1680  case 4:
1681  Y = cosph*pow(costh2-1.0,2.0)*(costh*1.6891875E+4-(costh*costh2)*1.6891875E+5+(costh*costh2*costh2)*2.87161875E+5)*2.414000363328839E-4;
1682  break;
1683  case 5:
1684  Y = cosph*pow(costh2*-1.0+1.0,5.0/2.0)*(costh2*-5.0675625E+5+(costh2*costh2)*1.435809375E+6+1.6891875E+4)*2.885282297193648E-5;
1685  break;
1686  case 6:
1687  Y = cosph*(costh*1.0135125E+6-(costh*costh2)*5.7432375E+6)*pow(costh2-1.0,3.0)*3.724883428715686E-6;
1688  break;
1689  case 7:
1690  Y = cosph*pow(costh2*-1.0+1.0,7.0/2.0)*(costh2*1.72297125E+7-1.0135125E+6)*5.376406125665728E-7;
1691  break;
1692  case 8:
1693  Y = cosph*costh*pow(costh2-1.0,4.0)*3.17731764895143;
1694  break;
1695  case 9:
1696  Y = cosph*pow(costh2*-1.0+1.0,9.0/2.0)*7.489009518540115E-1;
1697  break;
1698  }break;
1699  default: break;
1700  }
1701 
1702  return R*Y;
1703 }
1704 
1705 double ZernikeSphericalHarmonics(int l1, int n, int l2, int m, double xr, double yr, double zr, double r) {
1706  switch (l1) {
1707  case 0 : {
1708  switch (l2) {
1709  case 0: return ZernikeSphericalHarmonics<0, 0>(n, m, xr, yr, zr, r);
1710  case 1: return ZernikeSphericalHarmonics<0, 1>(n, m, xr, yr, zr, r);
1711  case 2: return ZernikeSphericalHarmonics<0, 2>(n, m, xr, yr, zr, r);
1712  case 3: return ZernikeSphericalHarmonics<0, 3>(n, m, xr, yr, zr, r);
1713  case 4: return ZernikeSphericalHarmonics<0, 4>(n, m, xr, yr, zr, r);
1714  case 5: return ZernikeSphericalHarmonics<0, 5>(n, m, xr, yr, zr, r);
1715  case 6: return ZernikeSphericalHarmonics<0, 6>(n, m, xr, yr, zr, r);
1716  case 7: return ZernikeSphericalHarmonics<0, 7>(n, m, xr, yr, zr, r);
1717  case 8: return ZernikeSphericalHarmonics<0, 8>(n, m, xr, yr, zr, r);
1718  case 9: return ZernikeSphericalHarmonics<0, 9>(n, m, xr, yr, zr, r);
1719  default: break;
1720  }
1721  } break;
1722  case 1 : {
1723  switch (l2) {
1724  case 0: return ZernikeSphericalHarmonics<1, 0>(n, m, xr, yr, zr, r);
1725  case 1: return ZernikeSphericalHarmonics<1, 1>(n, m, xr, yr, zr, r);
1726  case 2: return ZernikeSphericalHarmonics<1, 2>(n, m, xr, yr, zr, r);
1727  case 3: return ZernikeSphericalHarmonics<1, 3>(n, m, xr, yr, zr, r);
1728  case 4: return ZernikeSphericalHarmonics<1, 4>(n, m, xr, yr, zr, r);
1729  case 5: return ZernikeSphericalHarmonics<1, 5>(n, m, xr, yr, zr, r);
1730  case 6: return ZernikeSphericalHarmonics<1, 6>(n, m, xr, yr, zr, r);
1731  case 7: return ZernikeSphericalHarmonics<1, 7>(n, m, xr, yr, zr, r);
1732  case 8: return ZernikeSphericalHarmonics<1, 8>(n, m, xr, yr, zr, r);
1733  case 9: return ZernikeSphericalHarmonics<1, 9>(n, m, xr, yr, zr, r);
1734  default: break;
1735  }
1736  } break;
1737  case 2 : {
1738  switch (l2) {
1739  case 0: return ZernikeSphericalHarmonics<2, 0>(n, m, xr, yr, zr, r);
1740  case 1: return ZernikeSphericalHarmonics<2, 1>(n, m, xr, yr, zr, r);
1741  case 2: return ZernikeSphericalHarmonics<2, 2>(n, m, xr, yr, zr, r);
1742  case 3: return ZernikeSphericalHarmonics<2, 3>(n, m, xr, yr, zr, r);
1743  case 4: return ZernikeSphericalHarmonics<2, 4>(n, m, xr, yr, zr, r);
1744  case 5: return ZernikeSphericalHarmonics<2, 5>(n, m, xr, yr, zr, r);
1745  case 6: return ZernikeSphericalHarmonics<2, 6>(n, m, xr, yr, zr, r);
1746  case 7: return ZernikeSphericalHarmonics<2, 7>(n, m, xr, yr, zr, r);
1747  case 8: return ZernikeSphericalHarmonics<2, 8>(n, m, xr, yr, zr, r);
1748  case 9: return ZernikeSphericalHarmonics<2, 9>(n, m, xr, yr, zr, r);
1749  default: break;
1750  }
1751  } break;
1752  case 3 : {
1753  switch (l2) {
1754  case 0: return ZernikeSphericalHarmonics<3, 0>(n, m, xr, yr, zr, r);
1755  case 1: return ZernikeSphericalHarmonics<3, 1>(n, m, xr, yr, zr, r);
1756  case 2: return ZernikeSphericalHarmonics<3, 2>(n, m, xr, yr, zr, r);
1757  case 3: return ZernikeSphericalHarmonics<3, 3>(n, m, xr, yr, zr, r);
1758  case 4: return ZernikeSphericalHarmonics<3, 4>(n, m, xr, yr, zr, r);
1759  case 5: return ZernikeSphericalHarmonics<3, 5>(n, m, xr, yr, zr, r);
1760  case 6: return ZernikeSphericalHarmonics<3, 6>(n, m, xr, yr, zr, r);
1761  case 7: return ZernikeSphericalHarmonics<3, 7>(n, m, xr, yr, zr, r);
1762  case 8: return ZernikeSphericalHarmonics<3, 8>(n, m, xr, yr, zr, r);
1763  case 9: return ZernikeSphericalHarmonics<3, 9>(n, m, xr, yr, zr, r);
1764  default: break;
1765  }
1766  } break;
1767  case 4 : {
1768  switch (l2) {
1769  case 0: return ZernikeSphericalHarmonics<4, 0>(n, m, xr, yr, zr, r);
1770  case 1: return ZernikeSphericalHarmonics<4, 1>(n, m, xr, yr, zr, r);
1771  case 2: return ZernikeSphericalHarmonics<4, 2>(n, m, xr, yr, zr, r);
1772  case 3: return ZernikeSphericalHarmonics<4, 3>(n, m, xr, yr, zr, r);
1773  case 4: return ZernikeSphericalHarmonics<4, 4>(n, m, xr, yr, zr, r);
1774  case 5: return ZernikeSphericalHarmonics<4, 5>(n, m, xr, yr, zr, r);
1775  case 6: return ZernikeSphericalHarmonics<4, 6>(n, m, xr, yr, zr, r);
1776  case 7: return ZernikeSphericalHarmonics<4, 7>(n, m, xr, yr, zr, r);
1777  case 8: return ZernikeSphericalHarmonics<4, 8>(n, m, xr, yr, zr, r);
1778  case 9: return ZernikeSphericalHarmonics<4, 9>(n, m, xr, yr, zr, r);
1779  default: break;
1780  }
1781  } break;
1782  case 5 : {
1783  switch (l2) {
1784  case 0: return ZernikeSphericalHarmonics<5, 0>(n, m, xr, yr, zr, r);
1785  case 1: return ZernikeSphericalHarmonics<5, 1>(n, m, xr, yr, zr, r);
1786  case 2: return ZernikeSphericalHarmonics<5, 2>(n, m, xr, yr, zr, r);
1787  case 3: return ZernikeSphericalHarmonics<5, 3>(n, m, xr, yr, zr, r);
1788  case 4: return ZernikeSphericalHarmonics<5, 4>(n, m, xr, yr, zr, r);
1789  case 5: return ZernikeSphericalHarmonics<5, 5>(n, m, xr, yr, zr, r);
1790  case 6: return ZernikeSphericalHarmonics<5, 6>(n, m, xr, yr, zr, r);
1791  case 7: return ZernikeSphericalHarmonics<5, 7>(n, m, xr, yr, zr, r);
1792  case 8: return ZernikeSphericalHarmonics<5, 8>(n, m, xr, yr, zr, r);
1793  case 9: return ZernikeSphericalHarmonics<5, 9>(n, m, xr, yr, zr, r);
1794  default: break;
1795  }
1796  } break;
1797  case 6 : {
1798  switch (l2) {
1799  case 0: return ZernikeSphericalHarmonics<6, 0>(n, m, xr, yr, zr, r);
1800  case 1: return ZernikeSphericalHarmonics<6, 1>(n, m, xr, yr, zr, r);
1801  case 2: return ZernikeSphericalHarmonics<6, 2>(n, m, xr, yr, zr, r);
1802  case 3: return ZernikeSphericalHarmonics<6, 3>(n, m, xr, yr, zr, r);
1803  case 4: return ZernikeSphericalHarmonics<6, 4>(n, m, xr, yr, zr, r);
1804  case 5: return ZernikeSphericalHarmonics<6, 5>(n, m, xr, yr, zr, r);
1805  case 6: return ZernikeSphericalHarmonics<6, 6>(n, m, xr, yr, zr, r);
1806  case 7: return ZernikeSphericalHarmonics<6, 7>(n, m, xr, yr, zr, r);
1807  case 8: return ZernikeSphericalHarmonics<6, 8>(n, m, xr, yr, zr, r);
1808  case 9: return ZernikeSphericalHarmonics<6, 9>(n, m, xr, yr, zr, r);
1809  default: break;
1810  }
1811  } break;
1812  case 7 : {
1813  switch (l2) {
1814  case 0: return ZernikeSphericalHarmonics<7, 0>(n, m, xr, yr, zr, r);
1815  case 1: return ZernikeSphericalHarmonics<7, 1>(n, m, xr, yr, zr, r);
1816  case 2: return ZernikeSphericalHarmonics<7, 2>(n, m, xr, yr, zr, r);
1817  case 3: return ZernikeSphericalHarmonics<7, 3>(n, m, xr, yr, zr, r);
1818  case 4: return ZernikeSphericalHarmonics<7, 4>(n, m, xr, yr, zr, r);
1819  case 5: return ZernikeSphericalHarmonics<7, 5>(n, m, xr, yr, zr, r);
1820  case 6: return ZernikeSphericalHarmonics<7, 6>(n, m, xr, yr, zr, r);
1821  case 7: return ZernikeSphericalHarmonics<7, 7>(n, m, xr, yr, zr, r);
1822  case 8: return ZernikeSphericalHarmonics<7, 8>(n, m, xr, yr, zr, r);
1823  case 9: return ZernikeSphericalHarmonics<7, 9>(n, m, xr, yr, zr, r);
1824  default: break;
1825  }
1826  } break;
1827  case 8 : {
1828  switch (l2) {
1829  case 0: return ZernikeSphericalHarmonics<8, 0>(n, m, xr, yr, zr, r);
1830  case 1: return ZernikeSphericalHarmonics<8, 1>(n, m, xr, yr, zr, r);
1831  case 2: return ZernikeSphericalHarmonics<8, 2>(n, m, xr, yr, zr, r);
1832  case 3: return ZernikeSphericalHarmonics<8, 3>(n, m, xr, yr, zr, r);
1833  case 4: return ZernikeSphericalHarmonics<8, 4>(n, m, xr, yr, zr, r);
1834  case 5: return ZernikeSphericalHarmonics<8, 5>(n, m, xr, yr, zr, r);
1835  case 6: return ZernikeSphericalHarmonics<8, 6>(n, m, xr, yr, zr, r);
1836  case 7: return ZernikeSphericalHarmonics<8, 7>(n, m, xr, yr, zr, r);
1837  case 8: return ZernikeSphericalHarmonics<8, 8>(n, m, xr, yr, zr, r);
1838  case 9: return ZernikeSphericalHarmonics<8, 9>(n, m, xr, yr, zr, r);
1839  default: break;
1840  }
1841  } break;
1842  case 9 : {
1843  switch (l2) {
1844  case 0: return ZernikeSphericalHarmonics<9, 0>(n, m, xr, yr, zr, r);
1845  case 1: return ZernikeSphericalHarmonics<9, 1>(n, m, xr, yr, zr, r);
1846  case 2: return ZernikeSphericalHarmonics<9, 2>(n, m, xr, yr, zr, r);
1847  case 3: return ZernikeSphericalHarmonics<9, 3>(n, m, xr, yr, zr, r);
1848  case 4: return ZernikeSphericalHarmonics<9, 4>(n, m, xr, yr, zr, r);
1849  case 5: return ZernikeSphericalHarmonics<9, 5>(n, m, xr, yr, zr, r);
1850  case 6: return ZernikeSphericalHarmonics<9, 6>(n, m, xr, yr, zr, r);
1851  case 7: return ZernikeSphericalHarmonics<9, 7>(n, m, xr, yr, zr, r);
1852  case 8: return ZernikeSphericalHarmonics<9, 8>(n, m, xr, yr, zr, r);
1853  case 9: return ZernikeSphericalHarmonics<9, 9>(n, m, xr, yr, zr, r);
1854  default: break;
1855  }
1856  } break;
1857  case 10 : {
1858  switch (l2) {
1859  case 0: return ZernikeSphericalHarmonics<10, 0>(n, m, xr, yr, zr, r);
1860  case 1: return ZernikeSphericalHarmonics<10, 1>(n, m, xr, yr, zr, r);
1861  case 2: return ZernikeSphericalHarmonics<10, 2>(n, m, xr, yr, zr, r);
1862  case 3: return ZernikeSphericalHarmonics<10, 3>(n, m, xr, yr, zr, r);
1863  case 4: return ZernikeSphericalHarmonics<10, 4>(n, m, xr, yr, zr, r);
1864  case 5: return ZernikeSphericalHarmonics<10, 5>(n, m, xr, yr, zr, r);
1865  case 6: return ZernikeSphericalHarmonics<10, 6>(n, m, xr, yr, zr, r);
1866  case 7: return ZernikeSphericalHarmonics<10, 7>(n, m, xr, yr, zr, r);
1867  case 8: return ZernikeSphericalHarmonics<10, 8>(n, m, xr, yr, zr, r);
1868  case 9: return ZernikeSphericalHarmonics<10, 9>(n, m, xr, yr, zr, r);
1869  default: break;
1870  }
1871  } break;
1872  case 11 : {
1873  switch (l2) {
1874  case 0: return ZernikeSphericalHarmonics<11, 0>(n, m, xr, yr, zr, r);
1875  case 1: return ZernikeSphericalHarmonics<11, 1>(n, m, xr, yr, zr, r);
1876  case 2: return ZernikeSphericalHarmonics<11, 2>(n, m, xr, yr, zr, r);
1877  case 3: return ZernikeSphericalHarmonics<11, 3>(n, m, xr, yr, zr, r);
1878  case 4: return ZernikeSphericalHarmonics<11, 4>(n, m, xr, yr, zr, r);
1879  case 5: return ZernikeSphericalHarmonics<11, 5>(n, m, xr, yr, zr, r);
1880  case 6: return ZernikeSphericalHarmonics<11, 6>(n, m, xr, yr, zr, r);
1881  case 7: return ZernikeSphericalHarmonics<11, 7>(n, m, xr, yr, zr, r);
1882  case 8: return ZernikeSphericalHarmonics<11, 8>(n, m, xr, yr, zr, r);
1883  case 9: return ZernikeSphericalHarmonics<11, 9>(n, m, xr, yr, zr, r);
1884  default: break;
1885  }
1886  } break;
1887  case 12 : {
1888  switch (l2) {
1889  case 0: return ZernikeSphericalHarmonics<12, 0>(n, m, xr, yr, zr, r);
1890  case 1: return ZernikeSphericalHarmonics<12, 1>(n, m, xr, yr, zr, r);
1891  case 2: return ZernikeSphericalHarmonics<12, 2>(n, m, xr, yr, zr, r);
1892  case 3: return ZernikeSphericalHarmonics<12, 3>(n, m, xr, yr, zr, r);
1893  case 4: return ZernikeSphericalHarmonics<12, 4>(n, m, xr, yr, zr, r);
1894  case 5: return ZernikeSphericalHarmonics<12, 5>(n, m, xr, yr, zr, r);
1895  case 6: return ZernikeSphericalHarmonics<12, 6>(n, m, xr, yr, zr, r);
1896  case 7: return ZernikeSphericalHarmonics<12, 7>(n, m, xr, yr, zr, r);
1897  case 8: return ZernikeSphericalHarmonics<12, 8>(n, m, xr, yr, zr, r);
1898  case 9: return ZernikeSphericalHarmonics<12, 9>(n, m, xr, yr, zr, r);
1899  default: break;
1900  }
1901  } break;
1902  case 13 : {
1903  switch (l2) {
1904  case 0: return ZernikeSphericalHarmonics<13, 0>(n, m, xr, yr, zr, r);
1905  case 1: return ZernikeSphericalHarmonics<13, 1>(n, m, xr, yr, zr, r);
1906  case 2: return ZernikeSphericalHarmonics<13, 2>(n, m, xr, yr, zr, r);
1907  case 3: return ZernikeSphericalHarmonics<13, 3>(n, m, xr, yr, zr, r);
1908  case 4: return ZernikeSphericalHarmonics<13, 4>(n, m, xr, yr, zr, r);
1909  case 5: return ZernikeSphericalHarmonics<13, 5>(n, m, xr, yr, zr, r);
1910  case 6: return ZernikeSphericalHarmonics<13, 6>(n, m, xr, yr, zr, r);
1911  case 7: return ZernikeSphericalHarmonics<13, 7>(n, m, xr, yr, zr, r);
1912  case 8: return ZernikeSphericalHarmonics<13, 8>(n, m, xr, yr, zr, r);
1913  case 9: return ZernikeSphericalHarmonics<13, 9>(n, m, xr, yr, zr, r);
1914  default: break;
1915  }
1916  } break;
1917  case 14 : {
1918  switch (l2) {
1919  case 0: return ZernikeSphericalHarmonics<14, 0>(n, m, xr, yr, zr, r);
1920  case 1: return ZernikeSphericalHarmonics<14, 1>(n, m, xr, yr, zr, r);
1921  case 2: return ZernikeSphericalHarmonics<14, 2>(n, m, xr, yr, zr, r);
1922  case 3: return ZernikeSphericalHarmonics<14, 3>(n, m, xr, yr, zr, r);
1923  case 4: return ZernikeSphericalHarmonics<14, 4>(n, m, xr, yr, zr, r);
1924  case 5: return ZernikeSphericalHarmonics<14, 5>(n, m, xr, yr, zr, r);
1925  case 6: return ZernikeSphericalHarmonics<14, 6>(n, m, xr, yr, zr, r);
1926  case 7: return ZernikeSphericalHarmonics<14, 7>(n, m, xr, yr, zr, r);
1927  case 8: return ZernikeSphericalHarmonics<14, 8>(n, m, xr, yr, zr, r);
1928  case 9: return ZernikeSphericalHarmonics<14, 9>(n, m, xr, yr, zr, r);
1929  default: break;
1930  }
1931  } break;
1932  case 15 : {
1933  switch (l2) {
1934  case 0: return ZernikeSphericalHarmonics<15, 0>(n, m, xr, yr, zr, r);
1935  case 1: return ZernikeSphericalHarmonics<15, 1>(n, m, xr, yr, zr, r);
1936  case 2: return ZernikeSphericalHarmonics<15, 2>(n, m, xr, yr, zr, r);
1937  case 3: return ZernikeSphericalHarmonics<15, 3>(n, m, xr, yr, zr, r);
1938  case 4: return ZernikeSphericalHarmonics<15, 4>(n, m, xr, yr, zr, r);
1939  case 5: return ZernikeSphericalHarmonics<15, 5>(n, m, xr, yr, zr, r);
1940  case 6: return ZernikeSphericalHarmonics<15, 6>(n, m, xr, yr, zr, r);
1941  case 7: return ZernikeSphericalHarmonics<15, 7>(n, m, xr, yr, zr, r);
1942  case 8: return ZernikeSphericalHarmonics<15, 8>(n, m, xr, yr, zr, r);
1943  case 9: return ZernikeSphericalHarmonics<15, 9>(n, m, xr, yr, zr, r);
1944  default: break;
1945  }
1946  } break;
1947  default: break;
1948  }
1949  REPORT_ERROR(ERR_ARG_INCORRECT, "ZernikeSphericalHarmonics not supported for l1 = " + std::to_string(l1) + " and l2 = " + std::to_string(l2));
1950 }
1951 
1952 #ifdef NEVERDEFINED
1953 double ALegendreSphericalHarmonics(int l, int m, double xr, double yr, double zr, double r)
1954 {
1955  // General variables
1956  double rp=2.0*r-1.0;
1957  double rp2=rp*rp,xr2=xr*xr,yr2=yr*yr,zr2=zr*zr;
1958  double pol=sqrt(1.0-rp2);
1959  double pol2=pol*pol;
1960 
1961  // Associated Legendre polynomial
1962  double R=0.0;
1963 
1964  switch (l)
1965  {
1966  case 0:
1967  R=1.0/2.0;
1968  break;
1969  case 1:
1970  switch (m)
1971  {
1972  case -1:
1973  R=(1.0/4.0)*pol;
1974  break;
1975  case 0:
1976  R=(1.0/2.0)*pol;
1977  break;
1978  case 1:
1979  R=-(1.0/2.0)*pol;
1980  break;
1981  }break;
1982  case 2:
1983  switch (m)
1984  {
1985  case -2:
1986  R=(1.0/(2.0*24.0))*3.0*pol2;
1987  break;
1988  case -1:
1989  R=(1.0/12.0)*rp*pol;
1990  break;
1991  case 0:
1992  R=(1.0/4.0)*(3*rp2-1.0);
1993  break;
1994  case 1:
1995  R=-3.0*rp*pol;
1996  break;
1997  case 2:
1998  R=3.0*pol2;
1999  }break;
2000  case 3:
2001  switch (m)
2002  {
2003  case -3:
2004  R=(1.0/(2.0*720.0))*15.0*pol2*pol;
2005  break;
2006  case -2:
2007  R=(1.0/240.0)*15.0*rp*pol2;
2008  break;
2009  case -1:
2010  R=(1.0/24.0)*(3.0/2.0)*(5.0*rp2-1.0)*pol;
2011  break;
2012  case 0:
2013  R=(1.0/4.0)*(5.0*rp2*rp-3.0*rp);
2014  break;
2015  case 1:
2016  R=-(3.0/4.0)*(5.0*rp2-1.0)*pol;
2017  break;
2018  case 2:
2019  R=(15.0/2.0)*rp*pol2;
2020  break;
2021  case 3:
2022  R=-15.0*pol2*pol;
2023  break;
2024  }break;
2025  case 4:
2026  switch (m)
2027  {
2028  case -4:
2029  R=(1.0/(2.0*40320.0))*105.0*pol2*pol2;
2030  break;
2031  case -3:
2032  R=(1.0/(5040.0*2.0))*105.0*rp*pol2*pol;
2033  break;
2034  case -2:
2035  R=(1.0/(2.0*360.0))*(15.0/2.0)*(7.0*rp2-1.0)*pol2;
2036  break;
2037  case -1:
2038  R=(1.0/40.0)*(5.0/2.0)*(7.0*rp2*rp-3.0*rp)*pol;
2039  break;
2040  case 0:
2041  R=(1.0/16.0)*(35.0*rp2*rp2-30.0*rp2+3.0);
2042  break;
2043  case 1:
2044  R=-(5.0/4.0)*(7.0*rp2*rp-3*rp)*pol;
2045  break;
2046  case 2:
2047  R=(15.0/4.0)*(7.0*rp2-1)*pol2;
2048  break;
2049  case 3:
2050  R=-(105.0/2.0)*rp*pol2*pol;
2051  break;
2052  case 4:
2053  R=(105.0/2.0)*pol2*pol2;
2054  }break;
2055  }
2056 
2057  // Spherical harmonic
2058  double Y=0.0;
2059 
2060  switch (l)
2061  {
2062  case 0:
2063  Y = (1.0/2.0)*sqrt(1.0/PI);
2064  break;
2065  case 1:
2066  switch (m)
2067  {
2068  case -1:
2069  Y = sqrt(3.0/(4.0*PI))*yr;
2070  break;
2071  case 0:
2072  Y = sqrt(3.0/(4.0*PI))*zr;
2073  break;
2074  case 1:
2075  Y = sqrt(3.0/(4.0*PI))*xr;
2076  break;
2077  } break;
2078  case 2:
2079  switch (m)
2080  {
2081  case -2:
2082  Y = sqrt(15.0/(4.0*PI))*xr*yr;
2083  break;
2084  case -1:
2085  Y = sqrt(15.0/(4.0*PI))*zr*yr;
2086  break;
2087  case 0:
2088  Y = sqrt(5.0/(16.0*PI))*(-xr2-yr2+2.0*zr2);
2089  break;
2090  case 1:
2091  Y = sqrt(15.0/(4.0*PI))*xr*zr;
2092  break;
2093  case 2:
2094  Y = sqrt(15.0/(16.0*PI))*(xr2-yr2);
2095  break;
2096  } break;
2097  case 3:
2098  switch (m)
2099  {
2100  case -3:
2101  Y = sqrt(35.0/(16.0*2.0*PI))*yr*(3.0*xr2-yr2);
2102  break;
2103  case -2:
2104  Y = sqrt(105.0/(4.0*PI))*zr*yr*xr;
2105  break;
2106  case -1:
2107  Y = sqrt(21.0/(16.0*2.0*PI))*yr*(4.0*zr2-xr2-yr2);
2108  break;
2109  case 0:
2110  Y = sqrt(7.0/(16.0*PI))*zr*(2.0*zr2-3.0*xr2-3.0*yr2);
2111  break;
2112  case 1:
2113  Y = sqrt(21.0/(16.0*2.0*PI))*xr*(4.0*zr2-xr2-yr2);
2114  break;
2115  case 2:
2116  Y = sqrt(105.0/(16.0*PI))*zr*(xr2-yr2);
2117  break;
2118  case 3:
2119  Y = sqrt(35.0/(16.0*2.0*PI))*xr*(xr2-3.0*yr2);
2120  break;
2121  } break;
2122  case 4:
2123  switch (m)
2124  {
2125  case -4:
2126  Y = sqrt((35.0*9.0)/(16.0*PI))*yr*xr*(xr2-yr2);
2127  break;
2128  case -3:
2129  Y = sqrt((9.0*35.0)/(16.0*2.0*PI))*yr*zr*(3.0*xr2-yr2);
2130  break;
2131  case -2:
2132  Y = sqrt((9.0*5.0)/(16.0*PI))*yr*xr*(7.0*zr2-(xr2+yr2+zr2));
2133  break;
2134  case -1:
2135  Y = sqrt((9.0*5.0)/(16.0*2.0*PI))*yr*zr*(7.0*zr2-3.0*(xr2+yr2+zr2));
2136  break;
2137  case 0:
2138  Y = sqrt(9.0/(16.0*16.0*PI))*(35.0*zr2*zr2-30.0*zr2+3.0);
2139  break;
2140  case 1:
2141  Y = sqrt((9.0*5.0)/(16.0*2.0*PI))*xr*zr*(7.0*zr2-3.0*(xr2+yr2+zr2));
2142  break;
2143  case 2:
2144  Y = sqrt((9.0*5.0)/(8.0*8.0*PI))*(xr2-yr2)*(7.0*zr2-(xr2+yr2+zr2));
2145  break;
2146  case 3:
2147  Y = sqrt((9.0*35.0)/(16.0*2.0*PI))*xr*zr*(xr2-3.0*yr2);
2148  break;
2149  case 4:
2150  Y = sqrt((9.0*35.0)/(16.0*16.0*PI))*(xr2*(xr2-3.0*yr2)-yr2*(3.0*xr2-yr2));
2151  break;
2152  } break;
2153  }
2154 
2155  return R*Y;
2156 }
2157 #endif
2158 
2159 void spherical_index2lnm(int idx, int &l1, int &n, int &l2, int &m, int max_l1)
2160 {
2161  auto numR = static_cast<int>(std::floor((4+4*max_l1+std::pow(max_l1,2))/4));
2162  double aux_id = std::floor(idx-(idx/numR)*numR);
2163  l1 = static_cast<int>(std::floor((1.0+std::sqrt(1.0+4.0*aux_id))/2.0) + std::floor((2.0+2.0*std::sqrt(aux_id))/2.0) - 2.0);
2164  n = static_cast<int>(std::ceil((4.0*aux_id - l1*(l1+2.0))/2.0));
2165  l2 = static_cast<int>(std::floor(std::sqrt(std::floor(idx/numR))));
2166  m = static_cast<int>(std::floor(idx/numR)-l2*(l2+1));
2167 }
2168 
2169 int spherical_lnm2index(int l1, int n, int l2, int m, int max_l1)
2170 {
2171  auto numR = static_cast<int>(std::floor((4+4*max_l1+std::pow(max_l1,2))/4));
2172  int id_SH = l2*(l2+1)+m;
2173  auto id_R = static_cast<int>(std::floor((2*n + l1*(l1 + 2))/4));
2174  int id_Z = id_SH*numR+id_R;
2175  return id_Z;
2176 }
2177 
2178 template<typename T>
2179 void solveBySVD(const Matrix2D< T >& A, const Matrix1D< T >& b,
2180  Matrix1D< double >& result, double tolerance)
2181 {
2182  if (A.Xdim() == 0)
2183  REPORT_ERROR(ERR_MATRIX_EMPTY, "Solve: Matrix is empty");
2184 
2185  if (A.Xdim() != A.Ydim())
2186  REPORT_ERROR(ERR_MATRIX_SIZE, "Solve: Matrix is not squared");
2187 
2188  if (A.Xdim() != b.size())
2189  REPORT_ERROR(ERR_MATRIX_SIZE, "Solve: Different sizes of Matrix and Vector");
2190 
2191  if (b.isRow())
2192  REPORT_ERROR(ERR_MATRIX_DIM, "Solve: Not correct vector shape");
2193 
2194  // First perform de single value decomposition
2195  // Xmipp interface that calls to svdcmp of numerical recipes
2196  Matrix2D< double > u;
2197  Matrix2D< double > v;
2198  Matrix1D< double > w;
2199  svdcmp(A, u, w, v);
2200 
2201  // Here is checked if eigenvalues of the svd decomposition are acceptable
2202  // If a value is lower than tolerance, the it's zeroed, as this increases
2203  // the precision of the routine.
2204  for (int i = 0; i < w.size(); i++)
2205  if (w(i) < tolerance)
2206  w(i) = 0;
2207 
2208  // Set size of matrices
2209  result.resize(b.size());
2210 
2211  // Xmipp interface that calls to svdksb of numerical recipes
2212  Matrix1D< double > bd;
2213  typeCast(b, bd);
2214  svbksb(u, w, v, bd, result);
2215 }
2216 
2217 template<typename T>
2218 void solve(const Matrix2D<T>& A, const Matrix2D<T>& b, Matrix2D<T>& result)
2219 {
2220  if (A.Xdim() == 0)
2221  REPORT_ERROR(ERR_MATRIX_EMPTY, "Solve: Matrix is empty");
2222 
2223  if (A.Xdim() != A.Ydim())
2224  REPORT_ERROR(ERR_MATRIX_SIZE, "Solve: Matrix is not squared");
2225 
2226  if (A.Ydim() != b.Ydim())
2227  REPORT_ERROR(ERR_MATRIX_SIZE, "Solve: Different sizes of A and b");
2228 
2229  // Solve
2230  result = b;
2231  Matrix2D<T> Aux = A;
2232  gaussj(Aux.adaptForNumericalRecipes2(), Aux.Ydim(),
2233  result.adaptForNumericalRecipes2(), b.Xdim());
2234 }
2235 
2236 // explicit instantiation
2237 template void solveBySVD<double>(Matrix2D<double> const&, Matrix1D<double> const&, Matrix1D<double>&, double);
2238 template void solve<double>(Matrix2D<double> const&, Matrix2D<double> const&, Matrix2D<double>&);
FOR_ALL_ELEMENTS_IN_MATRIX2D
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
Definition: matrix2d.h:104
cholsl
void cholsl(double *a, int n, double *p, double *b, double *x)
Definition: numerical_recipes.cpp:9370
xi
double xi
Definition: numerical_recipes.cpp:4635
Matrix2D::Xdim
size_t Xdim() const
Definition: matrix2d.h:575
DESolver::probability
double probability
Definition: numerical_tools.h:388
numerical_tools.h
quadraticProgramming_grcn32
void quadraticProgramming_grcn32(int nparam, int j, double *gradgj, void *cd)
Definition: numerical_tools.cpp:206
min
void min(Image< double > &op1, const Image< double > &op2)
Definition: image_operate.cpp:118
MultidimArray::resize
void resize(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim, bool copy=true)
Definition: multidim_array.cpp:983
xmipp_memory.h
DESolver::population
double * population
Definition: numerical_tools.h:396
numerical_recipes.h
nineq
void nineq
Definition: numerical_recipes.cpp:4411
spherical_lnm2index
int spherical_lnm2index(int l1, int n, int l2, int m, int max_l1)
Definition: numerical_tools.cpp:2169
RowVector
#define RowVector(a, b)
Definition: numerical_tools.cpp:384
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
DESolver::bestSolution
double * bestSolution
Definition: numerical_tools.h:394
stBest1Exp
constexpr int stBest1Exp
Definition: numerical_tools.h:262
REPORT_ERROR
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
solve< double >
template void solve< double >(Matrix2D< double > const &, Matrix2D< double > const &, Matrix2D< double > &)
grobfd
void grobfd()
DESolver::~DESolver
virtual ~DESolver(void)
Destructor.
Definition: numerical_tools.cpp:402
Matrix1D::size
size_t size() const
Definition: matrix1d.h:508
Matrix2D::fromVector
void fromVector(const Matrix1D< T > &op1)
Definition: matrix2d.cpp:803
ERR_MATRIX_SIZE
Problem with matrix size.
Definition: xmipp_error.h:152
MultidimArray::initRandom
void initRandom(double op1, double op2, RandomMode mode=RND_UNIFORM)
Definition: multidim_array.cpp:821
L2
#define L2
Definition: cuda_volume_deform_sph_defines.h:17
c
doublereal * c
Definition: numerical_recipes.cpp:2230
g
doublereal * g
Definition: numerical_recipes.cpp:2473
bl
double * bl
Definition: numerical_recipes.cpp:4414
cd
void * cd
Definition: numerical_recipes.cpp:4416
quadraticProgramming_grob32
void quadraticProgramming_grob32(int nparam, double *x, double *gradfj, void *cd)
Definition: numerical_tools.cpp:192
xmipp_funcs.h
grad
doublereal * grad
Definition: numerical_recipes.cpp:2473
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
buffer
HBITMAP buffer
Definition: svm-toy.cpp:37
stBest2Exp
constexpr int stBest2Exp
Definition: numerical_tools.h:265
CDAB::A
Matrix2D< double > A
Definition: numerical_tools.cpp:160
quadraticProgramming_cntr32
void quadraticProgramming_cntr32(int nparam, int j, double *x, double *gj, void *cd)
Definition: numerical_tools.cpp:182
ERR_MEM_NOTENOUGH
There is not enough memory for allocation.
Definition: xmipp_error.h:166
choldc
void choldc(double *a, int n, double *p)
Definition: numerical_recipes.cpp:9347
nparam
static void nparam
Definition: numerical_recipes.cpp:5153
ERR_MATRIX_EMPTY
The matrix is empty.
Definition: xmipp_error.h:151
Matrix2D::toVector
void toVector(Matrix1D< T > &op1) const
Definition: matrix2d.cpp:832
Matrix2D::inv
void inv(Matrix2D< T > &result) const
Definition: matrix2d.cpp:663
w
doublereal * w
Definition: numerical_recipes.cpp:2477
abs
void abs(Image< double > &op)
Definition: image_operate.cpp:219
stBest1Bin
constexpr int stBest1Bin
Definition: numerical_tools.h:267
udelta
double udelta
Definition: numerical_recipes.cpp:4413
powellOptimizer
void powellOptimizer(Matrix1D< double > &p, int i0, int n, double(*f)(double *x, void *), void *prm, double ftol, double &fret, int &iter, const Matrix1D< double > &steps, bool show)
Definition: numerical_tools.cpp:44
nfsr
void nfsr
Definition: numerical_recipes.cpp:4411
stRand2Bin
constexpr int stRand2Bin
Definition: numerical_tools.h:271
Matrix2D::transpose
Matrix2D< T > transpose() const
Definition: matrix2d.cpp:1314
DESolver::Solve
virtual bool Solve(int maxGenerations)
Definition: numerical_tools.cpp:488
mesh_pts
void * mesh_pts
Definition: numerical_recipes.cpp:4411
CDAB::B
Matrix2D< double > B
Definition: numerical_tools.cpp:161
iter
glob_prnt iter
Definition: numerical_recipes.cpp:4698
Matrix2D::adaptForNumericalRecipes2
T * adaptForNumericalRecipes2() const
Definition: matrix2d.h:875
eps
cmache_1 eps
Definition: numerical_recipes.cpp:2273
DESolver::calcTrialSolution
StrategyFunction calcTrialSolution
Definition: numerical_tools.h:386
iprint
integer * iprint
Definition: numerical_recipes.cpp:2231
nineqn
void nineqn
Definition: numerical_recipes.cpp:4411
x
doublereal * x
Definition: numerical_recipes.cpp:2230
i
#define i
Definition: numerical_recipes.cpp:2493
k
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
gaussj
void gaussj(T *a, int n, T *b, int m)
Definition: numerical_recipes.h:266
d
doublereal * d
Definition: numerical_recipes.cpp:2230
randomPermutation
void randomPermutation(int N, MultidimArray< int > &result)
Definition: numerical_tools.cpp:33
DESolver::EnergyFunction
virtual double EnergyFunction(double testSolution[], bool &bAtSolution)=0
DESolver::trialEnergy
double trialEnergy
Definition: numerical_tools.h:390
rnd_unif
double rnd_unif()
Definition: xmipp_funcs.cpp:440
bu
double * bu
Definition: numerical_recipes.cpp:4414
solve
void solve(const Matrix2D< T > &A, const Matrix2D< T > &b, Matrix2D< T > &result)
Definition: numerical_tools.cpp:2218
CDAB::C
Matrix2D< double > C
Definition: numerical_tools.cpp:158
FOR_ALL_ELEMENTS_IN_MATRIX1D
#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
Definition: matrix1d.h:72
evaluateQuadratic
void evaluateQuadratic(const Matrix1D< double > &x, const Matrix1D< double > &c, const Matrix2D< double > &H, double &val, Matrix1D< double > &grad)
Definition: numerical_tools.cpp:126
miter
void miter
Definition: numerical_recipes.cpp:4411
b
doublereal * b
Definition: numerical_recipes.cpp:2230
nf
void nf
Definition: numerical_recipes.cpp:4411
solveBySVD< double >
template void solveBySVD< double >(Matrix2D< double > const &, Matrix1D< double > const &, Matrix1D< double > &, double)
quadraticProgramming
void quadraticProgramming(const Matrix2D< double > &C, const Matrix1D< double > &d, const Matrix2D< double > &A, const Matrix1D< double > &b, const Matrix2D< double > &Aeq, const Matrix1D< double > &beq, Matrix1D< double > &bl, Matrix1D< double > &bu, Matrix1D< double > &x)
Definition: numerical_tools.cpp:222
bigbnd
double bigbnd
Definition: numerical_recipes.cpp:4413
cfsqp
void cfsqp(nparam, nf, nfsr, nineqn, nineq, neqn, neq, ncsrl, ncsrn, mesh_pts, mode, iprint, miter, inform, bigbnd, eps, epseqn, udelta, bl, bu, x, f, g, lambda, obj, constr, gradob, gradcn, cd) int nparam
lambda
double * lambda
Definition: numerical_recipes.cpp:4414
Matrix1D::isRow
int isRow() const
Definition: matrix1d.h:520
ERR_MATRIX_DIM
Problem with matrix dimensions.
Definition: xmipp_error.h:150
in
int in
Definition: image_find_center.cpp:56
f
double * f
Definition: numerical_recipes.cpp:4414
ERR_ARG_INCORRECT
Incorrect argument received.
Definition: xmipp_error.h:113
DESolver::DESolver
DESolver(int dim, int popSize)
Empty constructor.
Definition: numerical_tools.cpp:387
Matrix1D::resize
void resize(size_t Xdim, bool copy=true)
Definition: matrix1d.h:410
ncsrn
void ncsrn
Definition: numerical_recipes.cpp:4411
neqn
void neqn
Definition: numerical_recipes.cpp:4411
stBest2Bin
constexpr int stBest2Bin
Definition: numerical_tools.h:270
cif::normalize
quaternion_type< T > normalize(quaternion_type< T > q)
Definition: point.cpp:278
max
void max(Image< double > &op1, const Image< double > &op2)
Definition: image_operate.cpp:128
ABS
#define ABS(x)
Definition: xmipp_macros.h:142
stRand2Exp
constexpr int stRand2Exp
Definition: numerical_tools.h:266
powell
void powell(double *p, double *xi, int n, double ftol, int &iter, double &fret, double(*func)(double *, void *), void *prm, bool show)
Definition: numerical_recipes.cpp:877
z
double z
Definition: numerical_recipes.cpp:8490
Matrix2D::Ydim
size_t Ydim() const
Definition: matrix2d.h:584
DESolver::Setup
void Setup(double min[], double max[], int deStrategy, double diffScale, double crossoverProb)
Setup() must be called before solve to set min, max, strategy etc.
Definition: numerical_tools.cpp:417
CDAB::D
Matrix2D< double > D
Definition: numerical_tools.cpp:159
matrixOperation_Ax
void matrixOperation_Ax(const Matrix2D< double > &A, const Matrix1D< double > &x, Matrix1D< double > &y)
Definition: matrix2d.cpp:424
ncsrl
void ncsrl
Definition: numerical_recipes.cpp:4411
mode
void mode
Definition: numerical_recipes.cpp:4411
quadraticProgramming_obj32
void quadraticProgramming_obj32(int nparam, int j, double *x, double *fj, void *cd)
Definition: numerical_tools.cpp:169
Element
#define Element(a, b, c)
Definition: numerical_tools.cpp:383
DESolver::trialSolution
double * trialSolution
Definition: numerical_tools.h:393
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
DESolver::bestEnergy
double bestEnergy
Definition: numerical_tools.h:391
svdcmp
void svdcmp(const Matrix2D< T > &a, Matrix2D< double > &u, Matrix1D< double > &w, Matrix2D< double > &v)
Definition: matrix2d.cpp:125
j
#define j
Definition: numerical_recipes.cpp:2493
steps
double steps
Definition: numerical_recipes.cpp:5166
checkRandomness
double checkRandomness(const std::string &sequence)
Definition: numerical_tools.cpp:806
m
int m
Definition: image_find_center.cpp:57
stRand1Bin
constexpr int stRand1Bin
Definition: numerical_tools.h:268
FOR_ALL_ELEMENTS_IN_ARRAY1D
#define FOR_ALL_ELEMENTS_IN_ARRAY1D(v)
Definition: multidim_array_base.h:554
regularizedLeastSquare
void regularizedLeastSquare(const Matrix2D< double > &A, const Matrix1D< double > &d, double lambda, const Matrix2D< double > &G, Matrix1D< double > &x)
Definition: numerical_tools.cpp:342
stRandToBest1Exp
constexpr int stRandToBest1Exp
Definition: numerical_tools.h:264
MATRIX1D_ARRAY
#define MATRIX1D_ARRAY(v)
Definition: matrix1d.h:58
DESolver::scale
double scale
Definition: numerical_tools.h:387
nnls
int nnls(double *a, int m, int n, double *b, double *x, double *rnorm, double *wp, double *zzp, int *indexp)
Definition: numerical_recipes.cpp:1165
typeCast
void typeCast(const Matrix1D< T1 > &v1, Matrix1D< T2 > &v2)
Definition: matrix1d.h:1227
neq
void neq
Definition: numerical_recipes.cpp:4411
r3
float r3
Definition: image_find_center.cpp:54
ZernikeSphericalHarmonics
double ZernikeSphericalHarmonics(int n, int m, double xr, double yr, double zr, double r)
Definition: numerical_tools.cpp:1054
stRandToBest1Bin
constexpr int stRandToBest1Bin
Definition: numerical_tools.h:269
GaussianInterpolator::initialize
void initialize(double xmax, int N, bool normalize=true)
Definition: numerical_tools.cpp:65
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626
inform
void * inform
Definition: numerical_recipes.cpp:4411
DESolver::popEnergy
double * popEnergy
Definition: numerical_tools.h:395
u
doublereal * u
Definition: numerical_recipes.cpp:2230
DESolver::SelectSamples
void SelectSamples(int candidate, int *r1, int *r2=nullptr, int *r3=nullptr, int *r4=nullptr, int *r5=nullptr)
Definition: numerical_tools.cpp:749
r2
float r2
Definition: image_find_center.cpp:54
prm
ProgClassifyCL2D * prm
Definition: mpi_classify_CL2D.cpp:35
solveNonNegative
double solveNonNegative(const Matrix2D< double > &C, const Matrix1D< double > &d, Matrix1D< double > &result)
Definition: numerical_tools.cpp:82
cif::to_string
std::string to_string(bond_type bondType)
Definition: compound.cpp:43
stRand1Exp
constexpr int stRand1Exp
Definition: numerical_tools.h:263
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
randomize_random_generator
unsigned int randomize_random_generator()
Definition: xmipp_funcs.cpp:407
solveViaCholesky
void solveViaCholesky(const Matrix2D< double > &A, const Matrix1D< double > &b, Matrix1D< double > &result)
Definition: numerical_tools.cpp:111
PI
#define PI
Definition: tools.h:43
Matrix1D::adaptForNumericalRecipes
T * adaptForNumericalRecipes() const
Definition: matrix1d.h:844
spherical_index2lnm
void spherical_index2lnm(int idx, int &l1, int &n, int &l2, int &m, int max_l1)
Definition: numerical_tools.cpp:2159
MATRIX2D_ARRAY
#define MATRIX2D_ARRAY(m)
Definition: matrix2d.h:89
Matrix1D::initConstant
void initConstant(T val)
Definition: matrix1d.cpp:83
n
int * n
Definition: numerical_recipes.cpp:2229
solveBySVD
void solveBySVD(const Matrix2D< T > &A, const Matrix1D< T > &b, Matrix1D< double > &result, double tolerance)
Definition: numerical_tools.cpp:2179
MultidimArray::indexSort
void indexSort(MultidimArray< int > &indx) const
Definition: multidim_array.cpp:882
L1
#define L1
Definition: cuda_volume_deform_sph_defines.h:14
grcnfd
void grcnfd()
svbksb
void svbksb(Matrix2D< double > &u, Matrix1D< double > &w, Matrix2D< double > &v, Matrix1D< double > &b, Matrix1D< double > &x)
Definition: matrix2d.cpp:174
leastSquare
void leastSquare(const Matrix2D< double > &C, const Matrix1D< double > &d, const Matrix2D< double > &A, const Matrix1D< double > &b, const Matrix2D< double > &Aeq, const Matrix1D< double > &beq, Matrix1D< double > &bl, Matrix1D< double > &bu, Matrix1D< double > &x)
Definition: numerical_tools.cpp:322
CopyVector
#define CopyVector(a, b)
Definition: numerical_tools.cpp:385
r1
float r1
Definition: image_find_center.cpp:54
Generated by   doxygen 1.8.13