in lib/maths/time_series/CAdaptiveBucketing.cc [552:687]
bool CAdaptiveBucketing::knots(core_t::TTime time,
common::CSplineTypes::EBoundaryCondition boundary,
TDoubleVec& knots,
TDoubleVec& values,
TDoubleVec& variances) const {
knots.clear();
values.clear();
variances.clear();
std::size_t n{m_Centres.size()};
for (std::size_t i = 0; i < n; ++i) {
if (this->bucketCount(i) > 0.0) {
double wide{3.0 * (m_Endpoints[n] - m_Endpoints[0]) / static_cast<double>(n)};
LOG_TRACE(<< "period " << m_Endpoints[n] - m_Endpoints[0]
<< ", # buckets = " << n << ", wide = " << wide);
// We get two points for each wide bucket but at most
// one third of the buckets can be wide. In this case
// we have 2 * n/3 + 2*n/3 knot points.
knots.reserve(4 * n / 3);
values.reserve(4 * n / 3);
variances.reserve(4 * n / 3);
double c{m_Centres[i]};
double c0{c - m_Endpoints[0]};
knots.push_back(m_Endpoints[0]);
values.push_back(this->predict(i, time, c));
variances.push_back(this->variance(i));
for (/**/; i < n; ++i) {
if (this->bucketCount(i) > 0.0) {
double a{m_Endpoints[i]};
double b{m_Endpoints[i + 1]};
double min{(9.0 * a + b) / 10.0};
double max{(a + 9.0 * b) / 10.0};
c = m_Centres[i];
double m{this->predict(i, time, c)};
double v{this->variance(i)};
if (b - a > wide) {
knots.push_back(std::max(c - (b - a) / 4.0, min));
values.push_back(m);
variances.push_back(v);
knots.push_back(std::min(c + (b - a) / 4.0, max));
values.push_back(m);
variances.push_back(v);
} else {
knots.push_back(std::min(std::max(c, min), max));
values.push_back(m);
variances.push_back(v);
}
}
}
switch (boundary) {
case common::CSplineTypes::E_Natural:
case common::CSplineTypes::E_ParabolicRunout:
knots.push_back(m_Endpoints[n]);
values.push_back(values.back());
variances.push_back(variances.back());
break;
case common::CSplineTypes::E_Periodic:
// We search for the value in the last and next period which
// are adjacent. Note that values need not be the same at the
// start and end of the period because the gradient can vary,
// but we expect them to be continuous.
for (std::size_t j = n - 1; j > 0; --j) {
if (this->bucketCount(j) > 0.0) {
double alpha{m_Endpoints[n] - m_Centres[j]};
double beta{c0};
double Z{alpha + beta};
if (Z == 0.0) {
alpha = beta = 0.5;
} else {
alpha /= Z;
beta /= Z;
}
double lastPeriodValue{
this->predict(j, time, m_Centres[j] - m_Endpoints[n])};
double lastPeriodVariance{this->variance(j)};
knots[0] = m_Endpoints[0];
values[0] = alpha * values[0] + beta * lastPeriodValue;
variances[0] = alpha * variances[0] + beta * lastPeriodVariance;
break;
}
}
for (std::size_t j = 0; j < n; ++j) {
if (this->bucketCount(j) > 0.0) {
double alpha{m_Centres[j]};
double beta{m_Endpoints[n] - knots.back()};
double Z{alpha + beta};
if (Z == 0.0) {
alpha = beta = 0.5;
} else {
alpha /= Z;
beta /= Z;
}
double nextPeriodValue{
this->predict(j, time, m_Endpoints[n] + m_Centres[j])};
double nextPeriodVariance{this->variance(j)};
values.push_back(alpha * values.back() + beta * nextPeriodValue);
variances.push_back(alpha * variances.back() + beta * nextPeriodVariance);
knots.push_back(m_Endpoints[n]);
break;
}
}
break;
}
}
}
if (knots.size() > 2) {
// If the distance between knot points becomes too small the
// spline can become poorly conditioned. We can safely discard
// knot points which are very close to one another.
TSizeVec indices(knots.size());
std::iota(indices.begin(), indices.end(), 0);
indices.erase(std::remove_if(indices.begin() + 1, indices.end() - 1,
[&knots](std::size_t i) {
return knots[i] - knots[i - 1] < 1.0 ||
knots[i + 1] - knots[i] < 1.0;
}),
indices.end() - 1);
if (indices.size() < knots.size()) {
for (std::size_t i = 0; i < indices.size(); ++i) {
knots[i] = knots[indices[i]];
values[i] = values[indices[i]];
variances[i] = variances[indices[i]];
}
knots.resize(indices.size());
values.resize(indices.size());
variances.resize(indices.size());
}
}
return knots.size() >= 2;
}