Program Listing for File lbfgs.hpp

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

#ifndef _PELE_LBFGS_H__
#define _PELE_LBFGS_H__

#include <fstream>
#include <iostream>
#include <memory>
#include <stdexcept>
#include <vector>

#include "array.hpp"
#include "base_potential.hpp"
#include "optimizer.hpp"

namespace pele {

class LBFGS : public GradientOptimizer {
 private:
  int M_;
  double max_f_rise_;
  bool use_relative_f_;

  // places to store the lbfgs memory
  Array<double> s_;
  Array<double> y_;
  Array<double> rho_;
  double H0_;
  int k_;

  Array<double> alpha;

  Array<double> xold;
  Array<double> gold;
  Array<double> step;
  double
      inv_sqrt_size;
#if PRINT_TO_FILE == 1
  std::ofstream trajectory_file;
#endif

 public:
  LBFGS(std::shared_ptr<pele::BasePotential> potential,
        const pele::Array<double> x0, double tol = 1e-4, int M = 4,
        bool save_trajectory = false, int iterations_before_save = 1,
        StopCriterionType stop_criterion = StopCriterionType::GRADIENT);
  virtual ~LBFGS() {
#if PRINT_TO_FILE == 1
    trajectory_file.close();
#endif
  }

  void one_iteration();

  // functions for setting the parameters
  inline void set_H0(double H0) {
    if (iter_number_ > 0) {
      std::cout << "warning: setting H0 after the first iteration.\n";
    }
    H0_ = H0;
  }
  inline void set_max_f_rise(double max_f_rise) { max_f_rise_ = max_f_rise; }

  inline void set_use_relative_f(int use_relative_f) {
    use_relative_f_ = (bool)use_relative_f;
  }

  // functions for accessing the results
  inline double get_H0() const { return H0_; }

  inline Array<double> get_s() const { return s_.copy(); }
  inline Array<double> get_y() const { return y_.copy(); }
  inline Array<double> get_rho() const { return rho_.copy(); }
  inline int get_k() const { return k_; }

  void set_state(Array<double> const &s, Array<double> const &y,
                 Array<double> const &rho, double H0, int k) {
    if (s.size() != s_.size() || y.size() != y_.size() ||
        rho.size() != rho_.size()) {
      throw std::invalid_argument(
          "LBFGS::set_state: memory arrays have the wrong size");
    }
    s_.assign(s);
    y_.assign(y);
    rho_.assign(rho);
    H0_ = H0;
    k_ = k;
  }

  virtual void reset(pele::Array<double> &x0);

 private:
  void update_memory(Array<double> xold, Array<double> gold, Array<double> xnew,
                     Array<double> gnew);

  void compute_lbfgs_step(Array<double> step);

  double backtracking_linesearch(Array<double> step);
};
}  // namespace pele

#endif