Program Listing for File cosine_gradient_descent.hpp

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

#ifndef _PELE_COS_GD_H__
#define _PELE_COS_GD_H__

#include <memory>
#include <pele/vecn.hpp>
#include <vector>

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

namespace pele {

class CosineGradientDescent : public ODEBasedOptimizer {
 private:
  Array<double> xold;
  Array<double> gold;
  Array<double> step;
  double cos_sim_tol_;
  double dt_;
  double dt_initial_;
  int back_track_n_;
                      //  if step keeps satisfying the condition for
                      //  back_track_n_ times,
                      // i.e we're not making progress, multiply the step by
                      // back_track_n_
  const int max_back_track = 100;
  int n_success_without_backtrack = 0;

 public:
  CosineGradientDescent(std::shared_ptr<pele::BasePotential> potential,
                        const pele::Array<double> x0, double tol = 1e-4,
                        double cos_sim_tol = 1e-2, double dt_initial = 1e-5,
                        int back_track_n = 5, bool save_trajectory = true,
                        int iterations_before_save = 1)
      : ODEBasedOptimizer(potential, x0, tol, save_trajectory,
                          iterations_before_save),
        xold(x_.size(), 0),
        gold(x_.size(), 0),
        step(x_.size(), 0),
        cos_sim_tol_(cos_sim_tol),
        dt_(dt_initial),
        dt_initial_(dt_initial),
        back_track_n_(back_track_n) {
    std::cout << "cosine gradient descent initialized" << std::endl;
    std::cout << "with parameters: " << std::endl;
    std::cout << "tol: " << tol << std::endl;
    std::cout << "cos_sim_tol: " << cos_sim_tol << std::endl;
    std::cout << "dt_initial: " << dt_initial << std::endl;
    std::cout << "back_track_n: " << back_track_n << std::endl;
    std::cout << "save_trajectory: " << save_trajectory << std::endl;
    std::cout << "iterations_before_save: " << iterations_before_save
              << std::endl;
    if (back_track_n_ < 1) {
      throw std::invalid_argument(
          "back_track_n must be greater than or equal to 1");
    }
  }

  virtual ~CosineGradientDescent() {}

  void one_iteration() {
    if (!func_initialized_) {
      initialize_func_gradient();
    }
    // make a copy of the position and gradient
    xold.assign(x_);
    //  really need to refactor and clean this up
    gold.assign(g_);
    // Gradient defines step direction
#if OPTIMIZER_DEBUG_LEVEL > 2
    std::cout << "cosine gradient descent: " << iter_number_ << " E " << f_
              << " rms " << gradient_norm_ << " nfev " << nfev_ << std::endl;
#endif
#if OPTIMIZER_DEBUG_LEVEL > 3
    std::cout << "x" << x_ << std::endl;
    std::cout << g_ << std::endl;
#endif
    // if the cosine similarity is less than the tolerance, reset the step size

    for (int i = 0; i < max_back_track; i++) {
      x_.assign(xold - dt_ * g_);
      f_ = potential_->get_energy_gradient(x_, g_);
      nfev_ += 1;
      double cos_sim_eps = abs(1 - dot(gold, g_) / (norm(gold) * norm(g_)));
#if OPTIMIZER_DEBUG_LEVEL > 2
      std::cout << "cosine similarity: " << cos_sim_eps << std::endl;
#endif
      if (cos_sim_eps < cos_sim_tol_) {
        if (i == 0) {
          n_success_without_backtrack += 1;
        }
        break;
      }
      // backtracking failed
      n_success_without_backtrack = 0;
      dt_ = dt_ / back_track_n_;
    }
    time_ += dt_;
    step.assign(x_ - xold);
    if (n_success_without_backtrack == back_track_n_) {
      dt_ = dt_ * back_track_n_;
      n_success_without_backtrack = 0;
    }
    step_norm_ = norm(step);
    gradient_norm_ = norm(g_);
    // print some status information
    if ((iprint_ > 0) && (iter_number_ % iprint_ == 0)) {
      std::cout << "steepest descent: " << iter_number_ << " E " << f_
                << " rms " << gradient_norm_ << " nfev " << nfev_
                << " step norm " << step_norm_ << std::endl;
    }
    iter_number_ += 1;
  }
  virtual void reset(pele::Array<double> &x0) {
    if (x0.size() != x_.size()) {
      throw std::invalid_argument(
          "The number of degrees of freedom (x0.size()) cannot change when "
          "calling reset()");
    }
    iter_number_ = 0;
    nfev_ = 0;
    n_success_without_backtrack = 0;
    x_.assign(x0);
    dt_ = dt_initial_;
    initialize_func_gradient();
    time_ = 0;
#if OPTIMIZER_DEBUG_LEVEL > 2
    std::cout << "resetting optimizer" << std::endl;
#endif
  }
};
}  // namespace pele

#endif