in include/maths/common/CSpline.h [400:560]
bool interpolate(const TDoubleVec& knots, const TDoubleVec& values, EBoundaryCondition boundary) {
if (knots.size() < 2) {
LOG_ERROR(<< "Insufficient knot points supplied");
return false;
}
if (knots.size() != values.size()) {
LOG_ERROR(<< "Number knots not equal to number of values: "
<< " knots = " << knots << " values = " << values);
return false;
}
TNonConstKnots oldKnots;
TNonConstValues oldValues;
oldKnots.swap(this->knotsRef());
oldValues.swap(this->valuesRef());
this->knotsRef().assign(knots.begin(), knots.end());
this->valuesRef().assign(values.begin(), values.end());
// If two knots are equal to working precision then we
// de-duplicate and use the average of the values at the
// duplicate knots. The exit condition from this loop is
// rather subtle and ensures that "last" always points
// to the last element in the reduced knot set.
std::size_t last = std::numeric_limits<std::size_t>::max();
std::size_t n = this->knots().size();
for (std::size_t i = 1; i <= n; ++i) {
std::size_t i_ = i - 1;
double knot = this->knots()[i_];
for (/**/; i < n && this->knots()[i] == knot; ++i) {
}
if (i - i_ > 1) {
TMeanAccumulator value;
for (std::size_t j = i_; j < i; ++j) {
value.add(this->values()[j]);
}
this->valuesRef()[i_] = CBasicStatistics::mean(value);
}
if (++last != i_) {
this->knotsRef()[last] = this->knots()[i_];
this->valuesRef()[last] = this->values()[i_];
}
}
this->knotsRef().erase(this->knotsRef().begin() + last + 1,
this->knotsRef().end());
this->valuesRef().erase(this->valuesRef().begin() + last + 1,
this->valuesRef().end());
n = this->knots().size();
LOG_TRACE(<< "knots = " << this->knots());
LOG_TRACE(<< "values = " << this->values());
if (this->knots().size() < 2) {
LOG_ERROR(<< "Insufficient distinct knot points supplied");
this->knotsRef().swap(oldKnots);
this->valuesRef().swap(oldValues);
return false;
}
switch (m_Type) {
case E_Linear:
// Curvatures are all zero and we don't bother to store them.
break;
case E_Cubic: {
this->curvaturesRef().clear();
this->curvaturesRef().reserve(n);
// Construct the diagonals: a is the subdiagonal, b is the
// main diagonal and c is the superdiagonal.
TDoubleVec a;
TDoubleVec b;
TDoubleVec c;
a.reserve(n - 1);
b.reserve(n);
c.reserve(n - 1);
double h = this->knots()[1] - this->knots()[0];
double h_ = this->knots()[n - 1] - this->knots()[n - 2];
switch (boundary) {
case E_Natural:
b.push_back(1.0);
c.push_back(0.0);
this->curvaturesRef().push_back(0.0);
break;
case E_ParabolicRunout:
b.push_back(1.0);
c.push_back(-1.0);
this->curvaturesRef().push_back(0.0);
break;
case E_Periodic:
b.push_back(2.0 * (h + h_));
c.push_back(h - 1.0);
this->curvaturesRef().push_back(
6.0 * ((this->values()[1] - this->values()[0]) / h -
(this->values()[0] - this->values()[n - 2]) / h_));
break;
}
for (std::size_t i = 1; i + 1 < n; ++i) {
h_ = h;
h = this->knots()[i + 1] - this->knots()[i];
a.push_back(h_);
b.push_back(2.0 * (h + h_));
c.push_back(h);
this->curvaturesRef().push_back(
6.0 * ((this->values()[i + 1] - this->values()[i]) / h -
(this->values()[i] - this->values()[i - 1]) / h_));
}
h_ = h;
h = this->knots()[1] - this->knots()[0];
switch (boundary) {
case E_Natural:
a.push_back(0.0);
b.push_back(1.0);
this->curvaturesRef().push_back(0.0);
if (!spline_detail::solveTridiagonal(a, b, c, this->curvaturesRef())) {
LOG_ERROR(<< "Failed to calculate curvatures");
return false;
}
break;
case E_ParabolicRunout:
a.push_back(-1.0);
b.push_back(1.0);
this->curvaturesRef().push_back(0.0);
if (!spline_detail::solveTridiagonal(a, b, c, this->curvaturesRef())) {
LOG_ERROR(<< "Failed to calculate curvatures");
return false;
}
break;
case E_Periodic: {
a.push_back(h_ * (1.0 - h));
b.push_back(2.0 * (h + h_));
TDoubleVec u(n, 0.0);
u[0] = 1.0;
u[n - 1] = h;
TDoubleVec v(n, 0.0);
v[1] = 1.0;
v[n - 2] = h_;
this->curvaturesRef().push_back(
6.0 * ((this->values()[1] - this->values()[n - 1]) / h -
(this->values()[n - 1] - this->values()[n - 2]) / h_));
if (!spline_detail::solvePeturbedTridiagonal(
a, b, c, u, v, this->curvaturesRef())) {
LOG_ERROR(<< "Failed to calculate curvatures");
return false;
}
} break;
}
break;
}
}
return true;
}