pele
Python energy landscape explorer
/home/js850/projects/pele/source/modified_fire.cpp
00001 #include "pele/modified_fire.h"
00002 
00003 using std::vector;
00004 using std::cout;
00005 
00006 namespace pele{
00042 MODIFIED_FIRE::MODIFIED_FIRE(std::shared_ptr<pele::BasePotential> potential,
00043         pele::Array<double>& x0, double dtstart, double dtmax, double maxstep,
00044         size_t Nmin, double finc, double fdec, double fa, double astart, double
00045         tol, bool stepback)
00046     : GradientOptimizer(potential,x0,tol), //call GradientOptimizer constructor
00047       _dtstart(dtstart), _dt(dtstart),
00048       _dtmax(dtmax), _maxstep(maxstep), _Nmin(Nmin),
00049       _finc(finc), _fdec(fdec), _fa(fa),
00050       _astart(astart), _a(astart), _fold(f_),
00051       _ifnorm(0),_vnorm(0),
00052       _v(x0.size(),0), _dx(x0.size()), _xold(x0.copy()),_gold(g_.copy()),
00053       _fire_iter_number(0), _N(x_.size()),
00054       _stepback(stepback)
00055 {}
00056 
00061 void MODIFIED_FIRE::initialize_func_gradient()
00062 {
00063     nfev_ += 1;                     //this accounts for the energy evaluation done by the integrator
00064     f_ = potential_->get_energy_gradient(x_, g_);
00065     _fold = f_;
00066     for (size_t k=0; k<x_.size();++k) { //set initial velocities (using forward Euler)
00067         _v[k] = -g_[k]*_dt;
00068     }
00069     _ifnorm = 1./norm(g_);
00070     _vnorm = norm(_v);
00071     rms_ = 1. / (_ifnorm*sqrt(_N));
00072     func_initialized_ = true;
00073 }
00074 
00075 void MODIFIED_FIRE::set_func_gradient(double f, Array<double> grad)
00076 {
00077     if (grad.size() != g_.size()) {
00078         throw std::invalid_argument("the gradient has the wrong size");
00079     }
00080     if (iter_number_ > 0){
00081         cout << "warning: setting f and grad after the first iteration.  this is dangerous.\n";
00082     }
00083 
00084     // copy the function and gradient
00085     f_ = f;
00086     _fold = f_;
00087     g_.assign(grad);
00088     for(size_t k=0; k<x_.size();++k) { //set initial velocities (using forward Euler)
00089         _v[k] = -g_[k]*_dt;
00090     }
00091     _ifnorm = 1./norm(g_);
00092     _vnorm = norm(_v);
00093     rms_ = 1. / (_ifnorm*sqrt(_N));
00094     func_initialized_ = true;
00095 }
00096 }
 All Classes Namespaces Functions Variables Typedefs