HcalAlgos/interface/ConstantStepOdeSolver.h

0001 #ifndef CalibCalorimetry_HcalAlgos_ConstantStepOdeSolver_h_
0002 #define CalibCalorimetry_HcalAlgos_ConstantStepOdeSolver_h_
0003
0004 #include <vector>
0005 #include <iostream>
0006 #include "FWCore/Utilities/interface/Exception.h"
0007
0008 #include "CalibCalorimetry/HcalAlgos/interface/AbsODERHS.h"
0009
0010 //
0011 // ODE solver with a constant time step. The derived classes are supposed
0012 // to implement concrete ODE solving algorithms (Runge-Kutta, etc).
0013 //
0014 class ConstantStepOdeSolver {
0015 public:
0016   inline ConstantStepOdeSolver() : rhs_(nullptr), dt_(0.0), dim_(0), runLen_(0), lastIntegrated_(0) {}
0017
0018   inline ConstantStepOdeSolver(const AbsODERHS& rhs)
0019       : rhs_(nullptr), dt_(0.0), dim_(0), runLen_(0), lastIntegrated_(0) {
0020     rhs_ = rhs.clone();
0021   }
0022
0023   // The copy constructor and the assignment operator are explicitly provided
0024   ConstantStepOdeSolver(const ConstantStepOdeSolver& r);
0025   ConstantStepOdeSolver& operator=(const ConstantStepOdeSolver& r);
0026
0027   inline virtual ~ConstantStepOdeSolver() { delete rhs_; }
0028
0029   // Access the equation right hand side
0030   inline void setRHS(const AbsODERHS& rhs) {
0031     delete rhs_;
0032     rhs_ = rhs.clone();
0033   }
0034   inline const AbsODERHS* getRHS() const { return rhs_; }
0035   inline AbsODERHS* getRHS() { return rhs_; }
0036
0037   // Inspectors (will be valid after at least one "run" call)
0038   inline unsigned lastDim() const { return dim_; }
0039   inline unsigned lastRunLength() const { return runLen_; }
0040   inline double lastDeltaT() const { return dt_; }
0041   inline double lastMaxT() const { return runLen_ ? dt_ * (runLen_ - 1U) : 0.0; }
0042
0043   inline double getTime(const unsigned idx) const {
0044     if (idx >= runLen_)
0045       throw cms::Exception("In ConstantStepOdeSolver::getTime: index out of range");
0046     return idx * dt_;
0047   }
0048
0049   inline double getCoordinate(const unsigned which, const unsigned idx) const {
0050     if (which >= dim_ || idx >= runLen_)
0051       throw cms::Exception("In ConstantStepOdeSolver::getCoordinate: index out of range");
0052     return historyBuffer_[dim_ * idx + which];
0053   }
0054
0055   // Integrate the node with the given number and get
0056   // the value of the integral at the given history point
0057   double getIntegrated(unsigned which, unsigned idx) const;
0058
0059   // Truncate some coordinate
0060   void truncateCoordinate(unsigned which, double minValue, double maxValue);
0061
0062   // Linear interpolation methods will be used in case the "cubic"
0063   // argument is false, and cubic in case it is true
0064   double interpolateCoordinate(unsigned which, double t, bool cubic = false) const;
0065
0066   // Interpolate the integrated node
0067   double interpolateIntegrated(unsigned which, double t, bool cubic = false) const;
0068
0069   // Get the time of the peak position
0070   double getPeakTime(unsigned which) const;
0071
0072   // Solve the ODE and remember the history
0073   void run(const double* initialConditions, unsigned lenConditions, double dt, unsigned nSteps);
0074
0075   // Set the history from some external source. The size
0076   // of the data array should be at least dim*runLen.
0077   void setHistory(double dt, const double* data, unsigned dim, unsigned runLen);
0078
0079   // Write out the history
0080   void writeHistory(std::ostream& os, double dt, bool cubic = false) const;
0081
0082   // Write out the integrated node
0083   void writeIntegrated(std::ostream& os, unsigned which, double dt, bool cubic = false) const;
0084
0085   // The following method must be overriden by derived classes
0086   virtual const char* methodName() const = 0;
0087
0088 protected:
0089   AbsODERHS* rhs_;
0090
0091 private:
0092   // The following method must be overriden by derived classes
0093   virtual void step(double t, double dt, const double* x, unsigned lenX, double* coordIncrement) const = 0;
0094
0095   // The following integration corresponds to the cubic
0096   // interpolation of the coordinate
0097   void integrateCoordinate(const unsigned which);
0098
0099   double dt_;
0100   unsigned dim_;
0101   unsigned runLen_;
0102   unsigned lastIntegrated_;
0103
0104   std::vector<double> historyBuffer_;
0105   std::vector<double> chargeBuffer_;
0106 };
0107
0108 // Dump the coordinate history as it was collected
0109 inline std::ostream& operator<<(std::ostream& os, const ConstantStepOdeSolver& s) {
0110   s.writeHistory(os, s.lastDeltaT());
0111   return os;
0112 }
0113
0114 // A few concrete ODE solvers
0115 class EulerOdeSolver : public ConstantStepOdeSolver {
0116 public:
0117   inline EulerOdeSolver() : ConstantStepOdeSolver() {}
0118
0119   inline explicit EulerOdeSolver(const AbsODERHS& rhs) : ConstantStepOdeSolver(rhs) {}
0120
0121   inline const char* methodName() const override { return "Euler"; }
0122
0123 private:
0124   void step(double t, double dt, const double* x, unsigned lenX, double* coordIncrement) const override;
0125 };
0126
0127 class RK2 : public ConstantStepOdeSolver {
0128 public:
0129   inline RK2() : ConstantStepOdeSolver() {}
0130
0131   inline explicit RK2(const AbsODERHS& rhs) : ConstantStepOdeSolver(rhs) {}
0132
0133   inline const char* methodName() const override { return "2nd order Runge-Kutta"; }
0134
0135 private:
0136   void step(double t, double dt, const double* x, unsigned lenX, double* coordIncrement) const override;
0137
0138   mutable std::vector<double> buf_;
0139 };
0140
0141 class RK4 : public ConstantStepOdeSolver {
0142 public:
0143   inline RK4() : ConstantStepOdeSolver() {}
0144
0145   inline explicit RK4(const AbsODERHS& rhs) : ConstantStepOdeSolver(rhs) {}
0146
0147   inline const char* methodName() const override { return "4th order Runge-Kutta"; }
0148
0149 private:
0150   void step(double t, double dt, const double* x, unsigned lenX, double* coordIncrement) const override;
0151
0152   mutable std::vector<double> buf_;
0153 };
0154
0155 #endif  // CalibCalorimetry_HcalAlgos_ConstantStepOdeSolver_h_