Program Listing for File modified_fire.hpp

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

#ifndef _PELE_MODIFIED_FIRE_H__
#define _PELE_MODIFIED_FIRE_H__

#include <fstream>
#include <iomanip>
#include <iostream>
#include <ostream>
#include <stdexcept>

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

namespace pele {

class MODIFIED_FIRE : public ODEBasedOptimizer {
 private:
  double _dtstart, _dt, _dtmax, _maxstep, _Nmin, _finc, _fdec, _fa, _astart, _a,
      _fold, _ifnorm, _vnorm;
  pele::Array<double> _v, _dx, _xold, _gold;
  size_t _fire_iter_number, _N;
  bool _stepback;
  inline void _ForwardEuler_integration();
  inline void _VelocityVerlet_integration();
#if PRINT_TO_FILE == 1
  std::ofstream trajectory_file;
  std::ofstream time_file;
  std::ofstream gradient_file;
#endif

 public:
  MODIFIED_FIRE(
      std::shared_ptr<pele::BasePotential> potential, pele::Array<double> &x0,
      double dtstart, double dtmax, double maxstep, size_t Nmin = 5,
      double finc = 1.1, double fdec = 0.5, double fa = 0.99,
      double astart = 0.1, double tol = 1e-4, bool stepback = true,
      bool save_trajectory = false, int iterations_before_save = 1,
      StopCriterionType stop_criterion_type = StopCriterionType::GRADIENT);
  virtual ~MODIFIED_FIRE() {
#if PRINT_TO_FILE == 1
    trajectory_file.close();
#endif
  }


  void one_iteration();


  void initialize_func_gradient();

  void set_func_gradient(double f, Array<double> grad);

  inline void reset(Array<double> &x0);
};

inline void MODIFIED_FIRE::reset(Array<double> &x0) {
  // arrays are already wrapped by the integrator, must not wrap them again,
  // just update their values, dont's use array.copy()
  //  or will wrap a copy to the array
  if (!func_initialized_) {
    initialize_func_gradient();
  }
  iter_number_ = 0;
  x_.assign(x0);
  f_ = potential_->get_energy_gradient(x_, g_);
  nfev_ = 1;
  // fire specific
  _fire_iter_number = 0;
  _dt = _dtstart;
  _a = _astart;
  _fold = f_;
  _xold.assign(x_);
  _gold.assign(g_);
  for (size_t k = 0; k < x_.size(); ++k) {
    _v[k] = -g_[k] * _dt;
  }
  _ifnorm = 1. / norm(g_);
  _vnorm = norm(_v);
  gradient_norm_ = 1. / (_ifnorm * sqrt(_N));
}

inline void MODIFIED_FIRE::_VelocityVerlet_integration() {
  /* the minuses in the following expressions are due to the fact that
   * the gradients rather than the forces appear in the expression
   */
  if (_dt == 0.0) {
    throw std::logic_error("time step cannot be zero");
  }
  for (size_t i = 0; i < _N; ++i) {
    _v[i] -=
        0.5 * _dt * (_gold[i] + g_[i]);  // update velocity assumes all masses 1
    _dx[i] =
        _dt *
        (_v[i] -
         0.5 * _dt * g_[i]);  // build displacement vector, assumes all masses 1
  }
  _gold.assign(g_);  // save gradient as old g
  double normdx = compute_pot_norm(_dx);

  if (normdx > _maxstep) {
    _dx *= (_maxstep /
            normdx);  // resize displacement vector is greater than _maxstep
  }

  x_ += _dx;

  f_ = potential_->get_energy_gradient(x_, g_);  // update gradient
}

inline void MODIFIED_FIRE::_ForwardEuler_integration() {
  /* the minuses in the following expressions are due to the fact that
   * the gradients rather than the forces appear in the expression
   */
  for (size_t i = 0; i < _N;
       ++i) {  // this was after get_energy_gradient, moved for testing
    _v[i] -= _dt * g_[i];  // update velocity, assumes all masses are 1
    _dx[i] = _dt * _v[i];  // build displacement vector
  }

  _gold.assign(g_);  // save gradient as old g
  double normdx = compute_pot_norm(_dx);

  if (normdx > _maxstep) {
    _dx *= (_maxstep /
            normdx);  // resize displacement vector is greater than _maxstep
  }

  x_ += _dx;

  f_ = potential_->get_energy_gradient(x_, g_);  // update gradient
}

inline void MODIFIED_FIRE::one_iteration() {
#if OPTIMIZER_DEBUG_LEVEL > 2
  std::cout << "func initialized" << func_initialized_ << std::endl;
#endif
  time_ += _dt;
  if (!func_initialized_) {
    initialize_func_gradient();
  }
  nfev_ += 1;
  iter_number_ += 1;
  _fire_iter_number +=
      1;  // this is different from iter_number_ which does not get reset

  // save old configuration in case next step has P < 0
  _fold = f_;        // set f_ old before integration step
  _xold.assign(x_);  // save x as xold, gold saved in integrator (because
                     // velocity verlet needs it, if vv is used)

  /*equation written in this conditional statement _v = (1- _a)*_v + _a * funit
   * * vnorm*/

  for (size_t i = 0; i < _N; ++i) {
    _v[i] = (1. - _a) * _v[i] - _a * g_[i] * _ifnorm * _vnorm;
  }

  /*run MD*/
  this->_ForwardEuler_integration();
  // this->_VelocityVerlet_integration();

  double P = -1 * dot(_v, g_);

  if (P > 0) {
    if (_fire_iter_number > _Nmin) {
      _dt = std::min(_dt * _finc, _dtmax);
      _a *= _fa;
    }

    _ifnorm = 1. / norm(g_);
    _vnorm = norm(_v);
    gradient_norm_ = 1. / (_ifnorm * sqrt(_N));  // update rms
  } else {
    _dt *= _fdec;
    _a = _astart;
    _fire_iter_number = 0;
    _v.assign(0);

    // reset position and gradient to the one before the step (core of modified
    // fire) reset velocity to initial (0)
    if (_stepback == true) {
      f_ = _fold;
      // reset position and gradient to the one before the step (core of
      // modified fire) reset velocity to initial (0)
      x_.assign(_xold);
      g_.assign(_gold);
    }
    if (verbosity_ > 1) {
      std::cout << "warning: step direction was uphill.  inverting\n";
    }
  }
  step_norm_ = norm(x_ - _xold);

  // print some status information
  if ((iprint_ > 0) && (iter_number_ % iprint_ == 0)) {
    std::cout << "fire: " << iter_number_ << " fire_iter_number "
              << _fire_iter_number << " dt " << _dt << " a " << _a << " P " << P
              << " vnorm " << _vnorm << " E " << f_ << " rms " << gradient_norm_
              << " nfev " << nfev_ << "\n";
  }
#if PRINT_TO_FILE == 1
  trajectory_file << std::setprecision(17) << x_;
  time_file << std::setprecision(17) << time_ << std::endl;
  gradient_file << std::setprecision(17) << g_;
#endif
}

}  // namespace pele

#endif  // #ifndef _PELE_MODIFIED_FIRE_H__