Greenstadt.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco, Mattia Piazza and Enrico Bertolazzi.                       *
 *                                                                                               *
 * The Optimist project is distributed under the BSD 2-Clause License.                           *
 *                                                                                               *
 * Davide Stocco                          Mattia Piazza                        Enrico Bertolazzi *
 * University of Trento               University of Trento                  University of Trento *
 * davide.stocco@unitn.it            mattia.piazza@unitn.it           enrico.bertolazzi@unitn.it *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#pragma once
 
#ifndef OPTIMIST_ROOTFINDER_GREENSTADT_HH
#define OPTIMIST_ROOTFINDER_GREENSTADT_HH
 
#include "Optimist/RootFinder/QuasiNewton.hh"
 
namespace Optimist
{
  namespace RootFinder
  {
 
    /*\
     |    ____                         _            _ _
     |   / ___|_ __ ___  ___ _ __  ___| |_ __ _  __| | |_
     |  | |  _| '__/ _ \/ _ \ '_ \/ __| __/ _` |/ _` | __|
     |  | |_| | | |  __/  __/ | | \__ \ || (_| | (_| | |_
     |   \____|_|  \___|\___|_| |_|___/\__\__,_|\__,_|\__|
     |
    \*/

    template <typename Real, Integer N>
    class Greenstadt : public QuasiNewton<Real, N, Greenstadt<Real, N>>
    {
    public:
      static constexpr bool requires_function{true};
      static constexpr bool requires_first_derivative{true};
      static constexpr bool requires_second_derivative{false};
 
      OPTIMIST_BASIC_CONSTANTS(Real)
 
      using typename QuasiNewton<Real, N, Greenstadt<Real, N>>::Vector;
      using typename QuasiNewton<Real, N, Greenstadt<Real, N>>::Matrix;

      using Method = enum class Method : Integer {ONE = 1, TWO = 2};
 
    private:
      Method m_method{Method::ONE}; 
 
    public:
      Greenstadt() {}

      std::string name_impl() const
      {
        std::ostringstream os;
        os << "Greenstadt";
        if (this->m_method == Method::ONE) {
          os << "1";
        } else if (this->m_method == Method::TWO) {
          os << "2";
        }
        return os.str();
      }

      Method method() const {return this->m_method;}

      void method(Method t_method) {this->m_method = t_method;}

      void enable_one_method() {this->m_method = Method::ONE;}

      void enable_two_method() {this->m_method = Method::TWO;}

      void set_method(Method t_method) {this->m_method = t_method;}

      void update_impl(
        Vector const & /*delta_x_old*/, Vector const & /*delta_function_old*/, Matrix const & jacobian_old,
        Vector const & delta_x_new,     Vector const & delta_function_new,     Vector const & function_new,
        Matrix       & jacobian_new
      ) {
        if (this->m_method == Method::ONE) {
          // Greenstadt's 1st method
          // J1 = J0 - (J0*DF1-DX1)/(C'*DF1)*C', where C = F1;
          jacobian_new = jacobian_old - (jacobian_old*delta_function_new-delta_x_new)/(function_new.transpose()*delta_function_new)*function_new.transpose();
        } else if (this->m_method == Method::TWO) {
          // Greenstadt's 2nd method
          // J1 = J0 - (J0*DF1-DX1)/(C'*DF1)*C', where C  = J0'*J0*DF1;
          Vector C(jacobian_old.transpose()*jacobian_old*delta_function_new);
          jacobian_new = jacobian_old - (jacobian_old*delta_function_new-delta_x_new)/(C.transpose()*delta_function_new)*C.transpose();
        }
      }
 
    }; // class Greenstadt
 
  } // namespace RootFinder
 
} // namespace Optimist
 
#endif // OPTIMIST_ROOTFINDER_GREENSTADT_HH