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_(0), dt_(0.0), dim_(0), runLen_(0), lastIntegrated_(0) {}

    inline ConstantStepOdeSolver(const AbsODERHS& rhs) :
        rhs_(0), 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 {return "Euler";}

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


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

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

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

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

    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 {return "4th order Runge-Kutta";}

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

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

#endif // CalibCalorimetry_HcalAlgos_ConstantStepOdeSolver_h_