Line data    Source code

       1             : #include "Opt/QuadraticSolver.hpp"
       2             : 
       3             : #include <algorithm>
       4             : #include <cmath>
       5             : #include <stdexcept>
       6             : 
       7             : using namespace std;
       8             : using namespace Opt;
       9             : 
      10             : typedef UnivariateSolver::Problem Problem;
      11             : typedef UnivariateSolver::Var     Var;
      12             : 
      13           0 : void QuadraticSolver::addInitialSolution(Var v) {
      14           0 :     initialSols.push_back(v);
      15           0 : }
      16             : 
      17           0 : void QuadraticSolver::clearInitialSolutions() {
      18           0 :     initialSols.clear();
      19           0 : }
      20             : 
      21           0 : void QuadraticSolver::setStopCriteria(Var e) {
      22           0 :     epsilon = e;
      23           0 : }
      24             : 
      25           0 : Var QuadraticSolver::solve(Problem f) {
      26           0 :     if(initialSols.size() < 3) {
      27           0 :         throw logic_error("QuadraticSolver requires at least 3 initial solutions");
      28             :     }
      29             : 
      30           0 :     vector<pair<Var, Var>> sols;
      31           0 :     sols.reserve(initialSols.size());
      32             : 
      33           0 :     for(const Var &x: initialSols)
      34           0 :         sols.emplace_back(f(x), x);
      35             : 
      36           0 :     while(sols.size() > 3)
      37           0 :         sols.pop_back();
      38           0 :     sort(sols.begin(), sols.end());
      39             : 
      40           0 :     sols.reserve(4);
      41             : 
      42           0 :     Var xPrev = sols.begin()->second;
      43           0 :     Var x     = sols.rbegin()->second;
      44             : 
      45           0 :     while(fabs(xPrev - x) > epsilon) {
      46           0 :         xPrev = x;
      47             : 
      48           0 :         const auto &[z1, x1] = sols[0];
      49           0 :         const auto &[z2, x2] = sols[1];
      50           0 :         const auto &[z3, x3] = sols[2];
      51           0 :         const Var
      52           0 :             xx1 = x1 * x1,
      53           0 :             xx2 = x2 * x2,
      54           0 :             xx3 = x3 * x3;
      55             : 
      56             :         // clang-format off
      57           0 :         x = 0.5 * (
      58           0 :             (xx2-xx3)*z1 + (xx3-xx1)*z2 + (xx1-xx2)*z3
      59           0 :         )/(
      60           0 :             ( x2- x3)*z1 + ( x3- x1)*z2 + ( x1- x2)*z3
      61             :         );
      62             :         // clang-format on
      63             : 
      64           0 :         sols.emplace_back(f(x), x);
      65           0 :         sort(sols.begin(), sols.end());
      66           0 :         while(sols.size() > 3)
      67           0 :             sols.pop_back();
      68             :     }
      69             : 
      70           0 :     return x;
      71             : }