Line data    Source code

       1             : #include <cmath>
       2             : #include <string>
       3             : #include <vector>
       4             : 
       5             : #include "data.hpp"
       6             : #include "exceptions.hpp"
       7             : #include "matrix.hpp"
       8             : #include "qrfactorization.hpp"
       9             : #include "vectorf.hpp"
      10             : 
      11          13 : void constructSixTermSystem(const vectorf<double> &coeff, const int nUse,
      12             :                             matrix<double> &W, vectorf<double> &b)
      13             : {
      14             :   // TODO Check coeff.get_size, nUse, and W.get_rows are compatible?
      15             : 
      16          65 :   for (int i{1}, k{(int)coeff.get_size() - nUse}; i <= nUse; i++, k++)
      17             :   {
      18          52 :     W(i, 1) = 2.0 * coeff(k);
      19          52 :     W(i, 2) = 2.0 * (k - 1) * coeff(k);
      20          52 :     W(i, 3) = -2.0 * coeff(k - 1);
      21          52 :     W(i, 4) = -(k - 2) * coeff(k - 1);
      22          52 :     b(i) = k * coeff(k + 1);
      23             :   }
      24          13 : }
      25             : 
      26          14 : void testBeta4(double beta4)
      27             : {
      28          14 :   if (beta4 < 0)
      29             :   {
      30             :     std::string message =
      31          10 :         "Unconstrained optimization lead to Sqrt of negative number: " +
      32          15 :         std::to_string(beta4) + " \n";
      33           5 :     throw sayMessage(message);
      34             :   }
      35           9 : }
      36             : 
      37          10 : void testRCSix(double rc)
      38             : {
      39          10 :   if (std::isnan(rc))
      40             :   {
      41             :     std::string message =
      42           2 :         "The radius of convergence is infinity, which is highly unlikely\n";
      43           1 :     throw sayMessage(message);
      44             :   }
      45           9 : }
      46             : 
      47          12 : void testCosTheta(double cosTheta)
      48             : {
      49          12 :   if (std::isnan(cosTheta))
      50             :   {
      51             :     std::string message =
      52             :         "Unconstrained optimization lead to infinite CosTheta which is not in "
      53           2 :         "[-1, 1]\n";
      54           1 :     throw sayMessage(message);
      55             :   }
      56          11 :   if ((cosTheta < -1.0) || (cosTheta > 1.0))
      57             :   {
      58          10 :     std::string message = "Unconstrained optimization lead to CosTheta [" +
      59          15 :                           std::to_string(cosTheta) + "] not in [-1, 1]\n";
      60           5 :     throw sayMessage(message);
      61             :   }
      62           6 : }
      63             : 
      64          10 : double testSingularityOrder(double singularityOrder1, double singularityOrder2)
      65             : {
      66          10 :   double order{0.0};
      67          10 :   if (std::isnan(singularityOrder1) && std::isnan(singularityOrder2))
      68             :   {
      69             :     std::string message =
      70           2 :         "Unconstrained optimization lead to NaN for Order of Singularity\n";
      71           1 :     throw sayMessage(message);
      72             :   }
      73             : 
      74           9 :   if (std::isnan(singularityOrder1) && !std::isnan(singularityOrder2))
      75           3 :     order = singularityOrder2;
      76             : 
      77           9 :   if (!std::isnan(singularityOrder1) && std::isnan(singularityOrder2))
      78           1 :     order = singularityOrder1;
      79             : 
      80           9 :   if (!std::isnan(singularityOrder1) && !std::isnan(singularityOrder2))
      81           5 :     order = (singularityOrder1 + singularityOrder2) / 2.0;
      82             : 
      83           9 :   return order;
      84             : }
      85             : 
      86             : // The six-term-test of Chang and Corliss
      87             : //
      88             : // The method returns the absolute value of the penultimate equation to nUse
      89             : // parameter, which must be at least 4, evaluated at the solution to the least
      90             : // squares optimal solution.
      91             : //
      92             : // The computation is successful whenever the value return is acceptably small,
      93             : // otherwise the algorithm is said to detect that the coefficients do not
      94             : // resemble pair of complex conjugate poles.
      95             : //
      96             : // The calling subroutine should maintain the invarient that coeff.size()
      97             : // must be sufficiently large, say 10.
      98          13 : double sixterm(const std::vector<double> &coeff, const double &scale,
      99             :                double &rc, double &order)
     100             : {
     101          13 :   int nUse = SIXTERM_NUSE;
     102          13 :   int m{nUse}, n{4};
     103          26 :   matrix<double> W(m, n);
     104          26 :   vectorf<double> beta(n);
     105          26 :   vectorf<double> b(m);
     106          20 :   vectorf<double> tc(coeff);
     107             : 
     108          13 :   constructSixTermSystem(tc, nUse, W, b);
     109             : 
     110             :   // Solve W beta = b for beta
     111          13 :   qr(m, n, W, b, beta);
     112             : 
     113             :   // Interpret the variables found from Least Squares Optimization Problem
     114          13 :   double hOverRc{0.0}, cosTheta{0.0}, singularityOrder1{0.0},
     115          13 :       singularityOrder2{0.0};
     116             : 
     117             :   // Evaluate h/Rc
     118          13 :   testBeta4(beta(4));
     119           9 :   hOverRc = sqrt(beta(4));
     120             : 
     121             :   // Evaluate Rc
     122           9 :   rc = scale / hOverRc;
     123           9 :   testRCSix(rc);
     124             : 
     125             :   // Evaluate Cos(Theta)
     126           9 :   cosTheta = beta(2) / hOverRc;
     127           9 :   testCosTheta(cosTheta);
     128             : 
     129             :   // Evaluate Order of the Singularity
     130           6 :   singularityOrder1 = beta(1) / beta(2);
     131           6 :   singularityOrder2 = beta(3) / beta(4);
     132           6 :   order = testSingularityOrder(singularityOrder1, singularityOrder2);
     133             : 
     134             :   // Compare order
     135             :   // printf( "Abs of difference between two computations for order:
     136             :   // [%22.16f]\n",
     137             :   //     fabs(singularityOrder1-singularityOrder2));
     138             :   // TODO Add this error to checked error on return?
     139             : 
     140             :   // Evaluate previous W equation at the solution of the least squares
     141             :   // optimization problem. Return the backward error, the absolute value of this
     142             :   // evaluation.
     143             :   //
     144             :   // TODO Use all equations?
     145           6 :   nUse++;
     146           6 :   int k = coeff.size() - nUse;
     147             :   double check =
     148           6 :       k * tc(k + 1) -
     149           6 :       ((2.0 * tc(k)) * beta(1) + (2.0 * (k - 1) * tc(k)) * beta(2) +
     150           6 :        (-2.0 * tc(k - 1)) * beta(3) + (-(k - 2) * tc(k - 1)) * beta(4));
     151             : 
     152          12 :   return fabs(check);
     153             : }