simulatedAnnealing.hpp

Go to the documentation of this file.

#ifndef JULIAN_SIMULATEDANNEALING_HPP
#define JULIAN_SIMULATEDANNEALING_HPP

#include <cmath>
#include <vector>
#include <utils/smartPointer.hpp>
#include <mathematics/RNGs/uniformRNG.hpp>
#include <mathematics/distributions/randomVariable.hpp>

namespace julian {
  class SimulatedAnnealing {
  public:
    SimulatedAnnealing(SmartPointer<UniformRNG> rng, SmartPointer<RandomVariable> rand):
      step_(rand), rng_(rng) {}


    void setCoolingSchedule(const std::vector<double>& cooling_schedule);
    void setLinearCooling(double Tstart, double Tend, double param);
    void setExponentialCooling(double Tstart, double Tend, double param);

    template<typename F, typename PA>
    double calculate(F f,PA prob_accept,double x_initial, int iters_per_t);    
    template<typename F, typename PA>
    std::vector<double> calculate(F f,PA prob_accept,std::vector<double> x_initial, int iters_per_t);

    double takeStep(double x);
    std::vector<double> takeStep(std::vector<double> x);
    ~SimulatedAnnealing(){}
  private:
    SmartPointer<RandomVariable> step_; 
    SmartPointer<UniformRNG> rng_; 
    std::vector<double> cooling_schedule_; 
  };
  
  template<typename F, typename PA>
  double SimulatedAnnealing::calculate(F f,PA prob_accept, double x_initial, int iters_per_t) { 
    double x = x_initial;
    double x_best = x;
    double x_new;
    
    double E = f(x_initial); 
    double E_best = E;
    double E_new;
    
    for (auto T : cooling_schedule_) {
      for (int iter = 0; iter < iters_per_t; iter++) {
    x_new = takeStep(x);
    E_new = f(x_new);
    if (E_new <= E_best) {
      x_best = x_new;
      E_best = E_new;
    }
    
    if (E_new < E) {
      x = x_new;
      E = E_new;
    } else if ( rng_->getRandom() < prob_accept(E, E_new, T)) { 
      x = x_new;
      E = E_new;
    }
      }
    }
    return x_best;
  }

  template<typename F, typename PA>
  std::vector<double> SimulatedAnnealing::calculate(F f, PA prob_accept, std::vector<double> x_initial, int iters_per_t) {
    std::vector<double> x = x_initial;
    std::vector<double> x_best = x;
    std::vector<double> x_new;
    
    double E = f(x_initial); 
    double E_best = E;
    double E_new;
        
    for (auto T : cooling_schedule_) {
      for (int iter = 0; iter < iters_per_t; iter++) {
    x_new = takeStep(x);
    E_new = f(x_new);
    
    if (E_new <= E_best) {
      x_best = x_new;
      E_best = E_new;
    }
    
    if (E_new < E) {
      x = x_new;
      E = E_new;
    } else if ( rng_->getRandom() < prob_accept(E, E_new, T)) { 
      x = x_new;
      E = E_new;
    }
      }
    }
    
    return x_best;
  }
}  // namespace julian

#endif /* SIMULATEDANNEALING_HPP */