Line data    Source code

       1             : /**
       2             :  * \file SimplifiedVectorizedNewtonRaphson.h
       3             :  */
       4             : 
       5             : #ifndef ATK_UTILITY_SIMPLIFIEDVECTORIZEDNEWTONRAPHSON_H
       6             : #define ATK_UTILITY_SIMPLIFIEDVECTORIZEDNEWTONRAPHSON_H
       7             : 
       8             : #include <ATK/config.h>
       9             : 
      10             : #include <Eigen/Dense>
      11             : 
      12             : #include <cmath>
      13             : #include <limits>
      14             : 
      15             : namespace ATK
      16             : {
      17             :   /// Simplified Vectorized Newton Raphson optimizer
      18             :   /*!
      19             :    * A NR optimizer, 10 iterations max, for stabilized computation
      20             :    */
      21             :   template<typename Function, int size, int max_iterations=10>
      22             :   class SimplifiedVectorizedNewtonRaphson
      23             :   {
      24             :     using DataType = typename Function::DataType;
      25             :     using Vector = Eigen::Matrix<DataType, size, 1>;
      26             : 
      27             :     Function function;
      28             :     
      29             :     DataType precision;
      30             :     DataType maxstep{10};
      31             :     Vector y0{Vector::Zero()};
      32             :     
      33             :   public:
      34             :     /*!
      35             :      * @brief Constructs the optimizer
      36             :      * @param function is the function that we will try to optimize.
      37             :      *   It is a functor taking x[n-1], x[n], y[n-1] and an estimate y[n], returning the value of the cost function and its derivative according to y[n]
      38             :      * @param precision is the precision that the optimizer will try to achieve. By default uses $$\\sqrt{\\epsilon_{Datatype}}$$
      39             :      */
      40             :     SimplifiedVectorizedNewtonRaphson(Function&& function, DataType precision = 0)
      41             :     :function(std::move(function)), precision(precision)
      42             :     {
      43             :       if(precision == 0)
      44             :       {
      45             :         this->precision = std::sqrt(std::numeric_limits<DataType>::epsilon());
      46             :       }
      47             :     }
      48             : 
      49          12 :     SimplifiedVectorizedNewtonRaphson(Function&& function, Vector&& y0, DataType precision = 0)
      50          12 :       :function(std::move(function)), precision(precision), y0(y0)
      51             :     {
      52          12 :       if (precision == 0)
      53             :       {
      54          12 :         this->precision = std::sqrt(std::numeric_limits<DataType>::epsilon());
      55             :       }
      56          12 :     }
      57             : 
      58             : 
      59             :     SimplifiedVectorizedNewtonRaphson(const SimplifiedVectorizedNewtonRaphson&) = delete;
      60             :     SimplifiedVectorizedNewtonRaphson& operator=(const SimplifiedVectorizedNewtonRaphson&) = delete;
      61             : 
      62             :     /// Optimize the function and sets its internal state
      63          12 :     Vector optimize()
      64             :     {
      65          12 :       Vector y1 = y0;
      66             :       gsl::index j;
      67         171 :       for(j = 0; j < max_iterations; ++j)
      68             :       {
      69         169 :         Vector cx = function(y1);
      70         169 :         auto maximum = cx.maxCoeff();
      71         169 :         auto minimum = cx.minCoeff();
      72         169 :         auto r = std::max(maximum, -minimum);
      73         169 :         if(r > maxstep)
      74             :         {
      75          75 :           cx *= maxstep/r;
      76             :         }
      77             : 
      78         169 :         auto yk = y1 - cx;
      79         169 :         if((cx.array().abs() < precision).all())
      80             :         {
      81          10 :           return yk;
      82             :         }
      83         159 :         y1 = yk;
      84             :       }
      85           2 :       return y1;
      86             :     }
      87             :   };
      88             : }
      89             : 
      90             : #endif