ConstantStepOdeSolver.h

Go to the documentation of this file.

#ifndef CalibCalorimetry_HcalAlgos_ConstantStepOdeSolver_h_
#define CalibCalorimetry_HcalAlgos_ConstantStepOdeSolver_h_

#include <vector>
#include <iostream>
#include "FWCore/Utilities/interface/Exception.h"

#include "CalibCalorimetry/HcalAlgos/interface/AbsODERHS.h"

//
// ODE solver with a constant time step. The derived classes are supposed
// to implement concrete ODE solving algorithms (Runge-Kutta, etc).
//
class ConstantStepOdeSolver {
public:
  inline ConstantStepOdeSolver() : rhs_(nullptr), dt_(0.0), dim_(0), runLen_(0), lastIntegrated_(0) {}

  inline ConstantStepOdeSolver(const AbsODERHS& rhs)
      : rhs_(nullptr), dt_(0.0), dim_(0), runLen_(0), lastIntegrated_(0) {
    rhs_ = rhs.clone();
  }

  // The copy constructor and the assignment operator are explicitly provided
  ConstantStepOdeSolver(const ConstantStepOdeSolver& r);
  ConstantStepOdeSolver& operator=(const ConstantStepOdeSolver& r);

  inline virtual ~ConstantStepOdeSolver() { delete rhs_; }

  // Access the equation right hand side
  inline void setRHS(const AbsODERHS& rhs) {
    delete rhs_;
    rhs_ = rhs.clone();
  }
  inline const AbsODERHS* getRHS() const { return rhs_; }
  inline AbsODERHS* getRHS() { return rhs_; }

  // Inspectors (will be valid after at least one "run" call)
  inline unsigned lastDim() const { return dim_; }
  inline unsigned lastRunLength() const { return runLen_; }
  inline double lastDeltaT() const { return dt_; }
  inline double lastMaxT() const { return runLen_ ? dt_ * (runLen_ - 1U) : 0.0; }

  inline double getTime(const unsigned idx) const {
    if (idx >= runLen_)
      throw cms::Exception("In ConstantStepOdeSolver::getTime: index out of range");
    return idx * dt_;
  }

  inline double getCoordinate(const unsigned which, const unsigned idx) const {
    if (which >= dim_ || idx >= runLen_)
      throw cms::Exception("In ConstantStepOdeSolver::getCoordinate: index out of range");
    return historyBuffer_[dim_ * idx + which];
  }

  // Integrate the node with the given number and get
  // the value of the integral at the given history point
  double getIntegrated(unsigned which, unsigned idx) const;

  // Truncate some coordinate
  void truncateCoordinate(unsigned which, double minValue, double maxValue);

  // Linear interpolation methods will be used in case the "cubic"
  // argument is false, and cubic in case it is true
  double interpolateCoordinate(unsigned which, double t, bool cubic = false) const;

  // Interpolate the integrated node
  double interpolateIntegrated(unsigned which, double t, bool cubic = false) const;

  // Get the time of the peak position
  double getPeakTime(unsigned which) const;

  // Solve the ODE and remember the history
  void run(const double* initialConditions, unsigned lenConditions, double dt, unsigned nSteps);

  // Set the history from some external source. The size
  // of the data array should be at least dim*runLen.
  void setHistory(double dt, const double* data, unsigned dim, unsigned runLen);

  // Write out the history
  void writeHistory(std::ostream& os, double dt, bool cubic = false) const;

  // Write out the integrated node
  void writeIntegrated(std::ostream& os, unsigned which, double dt, bool cubic = false) const;

  // The following method must be overriden by derived classes
  virtual const char* methodName() const = 0;

protected:
  AbsODERHS* rhs_;

private:
  // The following method must be overriden by derived classes
  virtual void step(double t, double dt, const double* x, unsigned lenX, double* coordIncrement) const = 0;

  // The following integration corresponds to the cubic
  // interpolation of the coordinate
  void integrateCoordinate(const unsigned which);

  double dt_;
  unsigned dim_;
  unsigned runLen_;
  unsigned lastIntegrated_;

  std::vector<double> historyBuffer_;
  std::vector<double> chargeBuffer_;
};

// Dump the coordinate history as it was collected
inline std::ostream& operator<<(std::ostream& os, const ConstantStepOdeSolver& s) {
  s.writeHistory(os, s.lastDeltaT());
  return os;
}

// A few concrete ODE solvers
class EulerOdeSolver : public ConstantStepOdeSolver {
public:
  inline EulerOdeSolver() : ConstantStepOdeSolver() {}

  inline explicit EulerOdeSolver(const AbsODERHS& rhs) : ConstantStepOdeSolver(rhs) {}

  inline const char* methodName() const override { return "Euler"; }

private:
  void step(double t, double dt, const double* x, unsigned lenX, double* coordIncrement) const override;
};

class RK2 : public ConstantStepOdeSolver {
public:
  inline RK2() : ConstantStepOdeSolver() {}

  inline explicit RK2(const AbsODERHS& rhs) : ConstantStepOdeSolver(rhs) {}

  inline const char* methodName() const override { return "2nd order Runge-Kutta"; }

private:
  void step(double t, double dt, const double* x, unsigned lenX, double* coordIncrement) const override;

  mutable std::vector<double> buf_;
};

class RK4 : public ConstantStepOdeSolver {
public:
  inline RK4() : ConstantStepOdeSolver() {}

  inline explicit RK4(const AbsODERHS& rhs) : ConstantStepOdeSolver(rhs) {}

  inline const char* methodName() const override { return "4th order Runge-Kutta"; }

private:
  void step(double t, double dt, const double* x, unsigned lenX, double* coordIncrement) const override;

  mutable std::vector<double> buf_;
};

#endif  // CalibCalorimetry_HcalAlgos_ConstantStepOdeSolver_h_