Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:47:11

0001 #include <cassert>
0002 
0003 #include "CalibCalorimetry/HcalAlgos/interface/ConstantStepOdeSolver.h"
0004 
0005 inline static double interpolateLinear(const double x, const double f0, const double f1) {
0006   return f0 * (1.0 - x) + f1 * x;
0007 }
0008 
0009 double ConstantStepOdeSolver::getPeakTime(const unsigned which) const {
0010   if (which >= dim_)
0011     throw cms::Exception("In ConstantStepOdeSolver::getPeakTime: index out of range");
0012   if (runLen_ < 3)
0013     throw cms::Exception("In ConstantStepOdeSolver::getPeakTime: not enough data");
0014 
0015   const double* hbuf = &historyBuffer_[which];
0016   double maxval = hbuf[0];
0017   unsigned maxind = 0;
0018   for (unsigned i = 1; i < runLen_; ++i)
0019     if (hbuf[dim_ * i] > maxval) {
0020       maxval = hbuf[dim_ * i];
0021       maxind = i;
0022     }
0023   if (maxind == 0U)
0024     return 0.0;
0025   if (maxind == runLen_ - 1U)
0026     return dt_ * maxind;
0027   const double l = hbuf[dim_ * (maxind - 1U)];
0028   const double r = hbuf[dim_ * (maxind + 1U)];
0029   if (l < maxval || r < maxval)
0030     return dt_ * (maxind + (l - r) / 2.0 / (l + r - 2.0 * maxval));
0031   else
0032     return dt_ * maxind;
0033 }
0034 
0035 double ConstantStepOdeSolver::getIntegrated(const unsigned which, const unsigned idx) const {
0036   if (which >= dim_ || idx >= runLen_)
0037     throw cms::Exception("In ConstantStepOdeSolver::getIntegrated: index out of range");
0038   if (lastIntegrated_ != which)
0039     (const_cast<ConstantStepOdeSolver*>(this))->integrateCoordinate(which);
0040   return chargeBuffer_[idx];
0041 }
0042 
0043 void ConstantStepOdeSolver::integrateCoordinate(const unsigned which) {
0044   if (runLen_ < 4)
0045     throw cms::Exception("In ConstantStepOdeSolver::integrateCoordinate: not enough data");
0046   if (chargeBuffer_.size() < runLen_)
0047     chargeBuffer_.resize(runLen_);
0048   double* integ = &chargeBuffer_[0];
0049   const double* coord = &historyBuffer_[which];
0050 
0051   integ[0] = 0.0;
0052   integ[1] = coord[dim_ * 0] * (3.0 / 8.0) + coord[dim_ * 1] * (19.0 / 24.0) + coord[dim_ * 2] * (-5.0 / 24.0) +
0053              coord[dim_ * 3] * (1.0 / 24.0);
0054   long double sum = integ[1];
0055   const unsigned rlenm1 = runLen_ - 1U;
0056   for (unsigned i = 2; i < rlenm1; ++i) {
0057     sum += (coord[dim_ * (i - 2U)] * (-1.0 / 24.0) + coord[dim_ * (i - 1U)] * (13.0 / 24.0) +
0058             coord[dim_ * i] * (13.0 / 24.0) + coord[dim_ * (i + 1U)] * (-1.0 / 24.0));
0059     integ[i] = sum;
0060   }
0061   sum += (coord[dim_ * rlenm1] * (3.0 / 8.0) + coord[dim_ * (rlenm1 - 1U)] * (19.0 / 24.0) +
0062           coord[dim_ * (rlenm1 - 2U)] * (-5.0 / 24.0) + coord[dim_ * (rlenm1 - 3U)] * (1.0 / 24.0));
0063   integ[rlenm1] = sum;
0064   if (dt_ != 1.0)
0065     for (unsigned i = 1; i < runLen_; ++i)
0066       integ[i] *= dt_;
0067   lastIntegrated_ = which;
0068 }
0069 
0070 void ConstantStepOdeSolver::writeHistory(std::ostream& os, const double dt, const bool cubic) const {
0071   os << "# " << this->methodName() << '\n';
0072   if (dim_ && runLen_) {
0073     if (dt == dt_) {
0074       for (unsigned ipt = 0; ipt < runLen_; ++ipt) {
0075         os << getTime(ipt);
0076         for (unsigned which = 0; which < dim_; ++which)
0077           os << ' ' << getCoordinate(which, ipt);
0078         os << '\n';
0079       }
0080     } else {
0081       const double tmax = lastMaxT();
0082       for (unsigned i = 0;; ++i) {
0083         const double t = i * dt;
0084         if (t > tmax)
0085           break;
0086         os << t;
0087         for (unsigned which = 0; which < dim_; ++which)
0088           os << ' ' << interpolateCoordinate(which, t, cubic);
0089         os << '\n';
0090       }
0091     }
0092   }
0093 }
0094 
0095 void ConstantStepOdeSolver::writeIntegrated(std::ostream& os,
0096                                             const unsigned which,
0097                                             const double dt,
0098                                             const bool cubic) const {
0099   os << "# " << this->methodName() << '\n';
0100   if (dim_ && runLen_) {
0101     if (dt == dt_) {
0102       for (unsigned ipt = 0; ipt < runLen_; ++ipt)
0103         os << getTime(ipt) << ' ' << getIntegrated(which, ipt) << '\n';
0104     } else {
0105       const double tmax = lastMaxT();
0106       for (unsigned i = 0;; ++i) {
0107         const double t = i * dt;
0108         if (t > tmax)
0109           break;
0110         os << t << ' ' << interpolateIntegrated(which, t, cubic) << '\n';
0111       }
0112     }
0113   }
0114 }
0115 
0116 ConstantStepOdeSolver::ConstantStepOdeSolver(const ConstantStepOdeSolver& r)
0117     : rhs_(nullptr),
0118       dt_(r.dt_),
0119       dim_(r.dim_),
0120       runLen_(r.runLen_),
0121       lastIntegrated_(r.lastIntegrated_),
0122       historyBuffer_(r.historyBuffer_),
0123       chargeBuffer_(r.chargeBuffer_) {
0124   if (r.rhs_)
0125     rhs_ = r.rhs_->clone();
0126 }
0127 
0128 ConstantStepOdeSolver& ConstantStepOdeSolver::operator=(const ConstantStepOdeSolver& r) {
0129   if (this != &r) {
0130     delete rhs_;
0131     rhs_ = nullptr;
0132     dt_ = r.dt_;
0133     dim_ = r.dim_;
0134     runLen_ = r.runLen_;
0135     lastIntegrated_ = r.lastIntegrated_;
0136     historyBuffer_ = r.historyBuffer_;
0137     chargeBuffer_ = r.chargeBuffer_;
0138     if (r.rhs_)
0139       rhs_ = r.rhs_->clone();
0140   }
0141   return *this;
0142 }
0143 
0144 void ConstantStepOdeSolver::truncateCoordinate(const unsigned which, const double minValue, const double maxValue) {
0145   if (which >= dim_)
0146     throw cms::Exception("In ConstantStepOdeSolver::truncateCoordinate: index out of range");
0147   if (minValue > maxValue)
0148     throw cms::Exception("In ConstantStepOdeSolver::truncateCoordinate: invalid truncation range");
0149 
0150   double* buf = &historyBuffer_[which];
0151   for (unsigned i = 0; i < runLen_; ++i) {
0152     if (buf[dim_ * i] < minValue)
0153       buf[dim_ * i] = minValue;
0154     else if (buf[dim_ * i] > maxValue)
0155       buf[dim_ * i] = maxValue;
0156   }
0157 }
0158 
0159 double ConstantStepOdeSolver::interpolateCoordinate(const unsigned which, const double t, const bool cubic) const {
0160   if (which >= dim_)
0161     throw cms::Exception("In ConstantStepOdeSolver::interpolateCoordinate: index out of range");
0162   if (runLen_ < 2U || (cubic && runLen_ < 4U))
0163     throw cms::Exception("In ConstantStepOdeSolver::interpolateCoordinate: not enough data");
0164   const double maxt = runLen_ ? dt_ * (runLen_ - 1U) : 0.0;
0165   if (t < 0.0 || t > maxt)
0166     throw cms::Exception("In ConstantStepOdeSolver::interpolateCoordinate: time out of range");
0167 
0168   const double* arr = &historyBuffer_[0];
0169   if (t == 0.0)
0170     return arr[which];
0171   else if (t == maxt)
0172     return arr[which + dim_ * (runLen_ - 1U)];
0173 
0174   // Translate time into timestep units
0175   const double tSteps = t / dt_;
0176   unsigned nLow = tSteps;
0177   if (nLow >= runLen_ - 1)
0178     nLow = runLen_ - 2;
0179   double x = tSteps - nLow;
0180 
0181   if (cubic) {
0182     unsigned i0 = 0;
0183     if (nLow == runLen_ - 2) {
0184       i0 = nLow - 2U;
0185       x += 2.0;
0186     } else if (nLow) {
0187       i0 = nLow - 1U;
0188       x += 1.0;
0189     }
0190     const double* base = arr + (which + dim_ * i0);
0191     return interpolateLinear(x * (3.0 - x) / 2.0,
0192                              interpolateLinear(x / 3.0, base[0], base[dim_ * 3]),
0193                              interpolateLinear(x - 1.0, base[dim_], base[dim_ * 2]));
0194   } else
0195     return interpolateLinear(x, arr[which + dim_ * nLow], arr[which + dim_ * (nLow + 1U)]);
0196 }
0197 
0198 double ConstantStepOdeSolver::interpolateIntegrated(const unsigned which, const double t, const bool cubic) const {
0199   if (which >= dim_)
0200     throw cms::Exception("In ConstantStepOdeSolver::interpolateIntegrated: index out of range");
0201   if (runLen_ < 2U || (cubic && runLen_ < 4U))
0202     throw cms::Exception("In ConstantStepOdeSolver::interpolateIntegrated: not enough data");
0203   const double maxt = runLen_ ? dt_ * (runLen_ - 1U) : 0.0;
0204   if (t < 0.0 || t > maxt)
0205     throw cms::Exception("In ConstantStepOdeSolver::interpolateIntegrated: time out of range");
0206   if (lastIntegrated_ != which)
0207     (const_cast<ConstantStepOdeSolver*>(this))->integrateCoordinate(which);
0208 
0209   const double* buf = &chargeBuffer_[0];
0210   if (t == 0.0)
0211     return buf[0];
0212   else if (t == maxt)
0213     return buf[runLen_ - 1U];
0214 
0215   // Translate time into timestep units
0216   const double tSteps = t / dt_;
0217   unsigned nLow = tSteps;
0218   if (nLow >= runLen_ - 1)
0219     nLow = runLen_ - 2;
0220   double x = tSteps - nLow;
0221 
0222   if (cubic) {
0223     unsigned i0 = 0;
0224     if (nLow == runLen_ - 2) {
0225       i0 = nLow - 2U;
0226       x += 2.0;
0227     } else if (nLow) {
0228       i0 = nLow - 1U;
0229       x += 1.0;
0230     }
0231     const double* base = buf + i0;
0232     return interpolateLinear(x * (3.0 - x) / 2.0,
0233                              interpolateLinear(x / 3.0, base[0], base[3]),
0234                              interpolateLinear(x - 1.0, base[1], base[2]));
0235   } else
0236     return interpolateLinear(x, buf[nLow], buf[nLow + 1U]);
0237 }
0238 
0239 void ConstantStepOdeSolver::setHistory(const double dt, const double* data, const unsigned dim, const unsigned runLen) {
0240   const unsigned len = dim * runLen;
0241   if (!len)
0242     return;
0243   if (dt <= 0.0)
0244     throw cms::Exception("In ConstantStepOdeSolver::setHistory: can not run backwards in time");
0245   assert(data);
0246   const unsigned sz = dim * (runLen + 1U);
0247   if (historyBuffer_.size() < sz)
0248     historyBuffer_.resize(sz);
0249   dt_ = dt;
0250   dim_ = dim;
0251   runLen_ = runLen;
0252   lastIntegrated_ = dim_;
0253   double* arr = &historyBuffer_[0];
0254   for (unsigned i = 0; i < len; ++i)
0255     *arr++ = *data++;
0256   for (unsigned i = 0; i < dim; ++i)
0257     *arr++ = 0.0;
0258 }
0259 
0260 void ConstantStepOdeSolver::run(const double* initialConditions,
0261                                 const unsigned lenInitialConditions,
0262                                 const double dt,
0263                                 const unsigned nSteps) {
0264   if (!nSteps || !lenInitialConditions)
0265     return;
0266   if (!rhs_)
0267     throw cms::Exception("In ConstantStepOdeSolver::run: ODE right hand side has not been set");
0268   if (dt <= 0.0)
0269     throw cms::Exception("In ConstantStepOdeSolver::run: can not run backwards in time");
0270   assert(initialConditions);
0271   dt_ = dt;
0272   dim_ = lenInitialConditions;
0273   lastIntegrated_ = dim_;
0274   runLen_ = nSteps + 1U;
0275   const unsigned sz = (runLen_ + 1U) * dim_;
0276   if (historyBuffer_.size() < sz)
0277     historyBuffer_.resize(sz);
0278   double* arr = &historyBuffer_[0];
0279   for (unsigned i = 0; i < lenInitialConditions; ++i)
0280     arr[i] = initialConditions[i];
0281   double* stepBuffer = arr + runLen_ * dim_;
0282 
0283   for (unsigned i = 0; i < nSteps; ++i, arr += lenInitialConditions) {
0284     const double t = i * dt;
0285     this->step(t, dt, arr, lenInitialConditions, stepBuffer);
0286     double* next = arr + lenInitialConditions;
0287     for (unsigned i = 0; i < lenInitialConditions; ++i)
0288       next[i] = arr[i] + stepBuffer[i];
0289   }
0290 }
0291 
0292 void EulerOdeSolver::step(
0293     const double t, const double dt, const double* x, const unsigned lenX, double* coordIncrement) const {
0294   rhs_->calc(t, x, lenX, coordIncrement);
0295   for (unsigned i = 0; i < lenX; ++i)
0296     coordIncrement[i] *= dt;
0297 }
0298 
0299 void RK2::step(const double t, const double dt, const double* x, const unsigned lenX, double* coordIncrement) const {
0300   const double halfstep = dt / 2.0;
0301   if (buf_.size() < lenX)
0302     buf_.resize(lenX);
0303   double* midpoint = &buf_[0];
0304   rhs_->calc(t, x, lenX, midpoint);
0305   for (unsigned i = 0; i < lenX; ++i) {
0306     midpoint[i] *= halfstep;
0307     midpoint[i] += x[i];
0308   }
0309   rhs_->calc(t + halfstep, midpoint, lenX, coordIncrement);
0310   for (unsigned i = 0; i < lenX; ++i)
0311     coordIncrement[i] *= dt;
0312 }
0313 
0314 void RK4::step(const double t, const double dt, const double* x, const unsigned lenX, double* coordIncrement) const {
0315   const double halfstep = dt / 2.0;
0316   if (buf_.size() < 4 * lenX)
0317     buf_.resize(4 * lenX);
0318   double* k1x = &buf_[0];
0319   double* k2x = k1x + lenX;
0320   double* k3x = k2x + lenX;
0321   double* k4x = k3x + lenX;
0322   rhs_->calc(t, x, lenX, k1x);
0323   for (unsigned i = 0; i < lenX; ++i)
0324     coordIncrement[i] = x[i] + halfstep * k1x[i];
0325   rhs_->calc(t + halfstep, coordIncrement, lenX, k2x);
0326   for (unsigned i = 0; i < lenX; ++i)
0327     coordIncrement[i] = x[i] + halfstep * k2x[i];
0328   rhs_->calc(t + halfstep, coordIncrement, lenX, k3x);
0329   for (unsigned i = 0; i < lenX; ++i)
0330     coordIncrement[i] = x[i] + dt * k3x[i];
0331   rhs_->calc(t + dt, coordIncrement, lenX, k4x);
0332   for (unsigned i = 0; i < lenX; ++i)
0333     coordIncrement[i] = dt / 6.0 * (k1x[i] + 2.0 * (k2x[i] + k3x[i]) + k4x[i]);
0334 }