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          19 : void constructThreeTermSystem(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         152 :   for (int i{1}, k{(int)coeff.get_size() - nUse}; i <= nUse; i++, k++)
      17             :   {
      18         133 :     W(i, 1) = (k - 1) * coeff(k);
      19         133 :     W(i, 2) = coeff(k);
      20         133 :     b(i) = k * coeff(k + 1);
      21             :   }
      22          19 : }
      23             : 
      24          20 : void testRCThree(double rc)
      25             : {
      26          20 :   if (std::isnan(rc))
      27             :   {
      28             :     std::string message =
      29           2 :         "Unconstrained optimization lead to NaN for Radius of Convergence\n";
      30           1 :     throw sayMessage(message);
      31             :   }
      32          19 : }
      33             : 
      34          20 : void testOrder(double order)
      35             : {
      36          20 :   if (std::isnan(order))
      37             :   {
      38             :     std::string message =
      39           2 :         "Unconstrained optimization lead to NaN for Order of Singularity\n";
      40           1 :     throw sayMessage(message);
      41             :   }
      42          19 : }
      43             : 
      44             : // The three-term-test of Chang and Corliss
      45          19 : double threeterm(const std::vector<double> &coeff, const double &scale,
      46             :                  double &rc, double &order)
      47             : {
      48          19 :   int nUse = THREETERM_NUSE;
      49          19 :   int m{nUse}, n{2};
      50          38 :   matrix<double> W(m, n);
      51          38 :   vectorf<double> beta(n);
      52          38 :   vectorf<double> b(m);
      53          19 :   vectorf<double> tc(coeff);
      54             : 
      55          19 :   constructThreeTermSystem(tc, nUse, W, b);
      56             : 
      57             :   // Solve W beta = b for beta
      58          19 :   qr(m, n, W, b, beta);
      59             : 
      60             :   // Interpret the variables found from Least Squares Optimization Problem
      61          19 :   double hOverRc = beta(1);
      62          19 :   rc = scale / hOverRc;
      63          19 :   testRCThree(rc);
      64             : 
      65          19 :   order = beta(2) / beta(1);
      66          19 :   testOrder(order);
      67             : 
      68             :   // Evaluate previous W equation at the solution of the least squares
      69             :   // optimization problem. Return the backward error, the absolute value of this
      70             :   // evaluation.
      71             :   //
      72             :   // TODO Use all equations?
      73          19 :   nUse++;
      74          19 :   int k = coeff.size() - nUse;
      75          19 :   double check = k * tc(k + 1) - ((k - 1) * tc(k) * beta(1) + tc(k) * beta(2));
      76             : 
      77          38 :   return fabs(check);
      78             : }