Program Listing for File newton.hpp

↰ Return to documentation for file (source/pele/newton.hpp)

// Newton method that ignores singular while calculating the inverse
// of a matrix. Modified to remove rattlers.

#ifndef _PELE_NEWTON_H_
#define _PELE_NEWTON_H_

#include <Eigen/Dense>
#include <Eigen/src/Core/Matrix.h>
#include <array>
#include <cstddef>
#include <cstdint>
#include <pele/vecn.hpp>

#include "array.hpp"
#include "base_potential.hpp"
#include "eigen_interface.hpp"
#include "pele/backtracking.hpp"
#include "pele/optimizer.hpp"

namespace pele {

class Newton : public GradientOptimizer {
 private:
  Eigen::MatrixXd _hessian;
  Eigen::VectorXd _gradient;
  Eigen::VectorXd _x;
  Eigen::VectorXd _x_old;
  Eigen::VectorXd _gradient_old;
  Eigen::VectorXd _step;
  BacktrackingLineSearch _line_search;
  double _energy;
  double _tolerance;
  double _threshold;
  int _nhev;
  bool _use_rattler_mask;
  Array<uint8_t> _not_rattlers;
  bool _jammed;
  bool _rattlers_found;

 public:
  Newton(std::shared_ptr<BasePotential> potential,
         const pele::Array<double> &x0, double tol = 1e-6,
         double threshold = std::numeric_limits<double>::epsilon(),
         bool use_rattler_mask = false);

  virtual ~Newton() {}

  void one_iteration();
  double compute_step();

  // Trust region stuff
  double compute_trust_region_step();
  Array<double> solve_subproblem();
  double ev(Eigen::MatrixXd, Array<double>);
  double get_rho_k(Array<double> p_k);
  double delta_k;
  double eta = 0.1;
  Array<double> p_k;
  Array<double> x_k;
  double delta_max = 0.1;
  double prev_e;

  // TODO: implement
  // void reset(pele::Array<double> x0);
  inline int get_nhev() const { return _nhev; }
  inline Eigen::VectorXd get_step() const { return _step; }
  void reset(Array<double> x);
  void get_rattler_details(Array<uint8_t> &not_rattlers, bool &jammed) {
    if (!_rattlers_found) {
      throw std::runtime_error(" Rattler based method not called");
    }
    not_rattlers = not_rattlers;
    jammed = _jammed;
  }
  // Warning: does not pass info about whether we're near a minimum or not.
  bool is_jammed() {
    if (_rattlers_found) {
      return _jammed;
    } else {
      return false;
    }
  }

  void compute_func_gradient(Array<double> x, double &func,
                             Array<double> gradient);
};

}  // namespace pele

#endif  // !_PELE_NEWTON_H_