File indexing completed on 2024-04-06 11:58:05
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
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
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 }
