27 #ifndef OPM_SPLINE_HPP
28 #define OPM_SPLINE_HPP
89 template<
class Scalar>
93 typedef std::vector<Scalar> Vector;
127 Scalar y0, Scalar y1,
128 Scalar m0, Scalar m1)
129 {
set(x0, x1, y0, y1, m0, m1); }
139 template <
class ScalarArrayX,
class ScalarArrayY>
141 const ScalarArrayX& x,
142 const ScalarArrayY& y,
144 bool sortInputs =
true)
145 { this->
setXYArrays(nSamples, x, y, splineType, sortInputs); }
154 template <
class Po
intArray>
156 const PointArray& points,
158 bool sortInputs =
true)
168 template <
class ScalarContainer>
170 const ScalarContainer& y,
172 bool sortInputs =
true)
181 template <
class Po
intContainer>
184 bool sortInputs =
true)
197 template <
class ScalarArray>
199 const ScalarArray& x,
200 const ScalarArray& y,
203 bool sortInputs =
true)
204 { this->
setXYArrays(nSamples, x, y, m0, m1, sortInputs); }
215 template <
class Po
intArray>
217 const PointArray& points,
220 bool sortInputs =
true)
232 template <
class ScalarContainerX,
class ScalarContainerY>
234 const ScalarContainerY& y,
237 bool sortInputs =
true)
248 template <
class Po
intContainer>
252 bool sortInputs =
true)
266 void set(Scalar x0, Scalar x1,
267 Scalar y0, Scalar y1,
268 Scalar m0, Scalar m1)
321 template <
class ScalarArrayX,
class ScalarArrayY>
323 const ScalarArrayX& x,
324 const ScalarArrayY& y,
325 Scalar m0, Scalar m1,
326 bool sortInputs =
true)
328 assert(nSamples > 1);
331 for (
size_t i = 0; i < nSamples; ++i) {
355 template <
class ScalarContainerX,
class ScalarContainerY>
357 const ScalarContainerY& y,
358 Scalar m0, Scalar m1,
359 bool sortInputs =
true)
361 assert(x.size() == y.size());
362 assert(x.size() > 1);
366 std::copy(x.begin(), x.end(), xPos_.begin());
367 std::copy(y.begin(), y.end(), yPos_.begin());
389 template <
class Po
intArray>
391 const PointArray& points,
394 bool sortInputs =
true)
398 assert(nSamples > 1);
401 for (
size_t i = 0; i < nSamples; ++i) {
402 xPos_[i] = points[i][0];
403 yPos_[i] = points[i][1];
426 template <
class XYContainer>
430 bool sortInputs =
true)
434 assert(points.size() > 1);
437 typename XYContainer::const_iterator it = points.begin();
438 typename XYContainer::const_iterator endIt = points.end();
439 for (
size_t i = 0; it != endIt; ++i, ++it) {
468 template <
class XYContainer>
472 bool sortInputs =
true)
476 typename XYContainer::const_iterator it = points.begin();
477 typename XYContainer::const_iterator endIt = points.end();
478 for (
unsigned i = 0; it != endIt; ++i, ++it) {
479 xPos_[i] = std::get<0>(*it);
480 yPos_[i] = std::get<1>(*it);
509 template <
class ScalarArrayX,
class ScalarArrayY>
511 const ScalarArrayX& x,
512 const ScalarArrayY& y,
514 bool sortInputs =
true)
516 assert(nSamples > 1);
519 for (
size_t i = 0; i < nSamples; ++i) {
529 if (splineType == Periodic)
531 else if (splineType == Natural)
533 else if (splineType == Monotonic)
536 throw std::runtime_error(
"Spline type "+std::to_string(
int(splineType))+
" not supported at this place");
550 template <
class ScalarContainerX,
class ScalarContainerY>
552 const ScalarContainerY& y,
554 bool sortInputs =
true)
556 assert(x.size() == y.size());
557 assert(x.size() > 1);
560 std::copy(x.begin(), x.end(), xPos_.begin());
561 std::copy(y.begin(), y.end(), yPos_.begin());
568 if (splineType == Periodic)
570 else if (splineType == Natural)
572 else if (splineType == Monotonic)
575 throw std::runtime_error(
"Spline type "+std::to_string(
int(splineType))+
" not supported at this place");
590 template <
class Po
intArray>
592 const PointArray& points,
594 bool sortInputs =
true)
598 assert(nSamples > 1);
601 for (
size_t i = 0; i < nSamples; ++i) {
602 xPos_[i] = points[i][0];
603 yPos_[i] = points[i][1];
611 if (splineType == Periodic)
613 else if (splineType == Natural)
615 else if (splineType == Monotonic)
618 throw std::runtime_error(
"Spline type "+std::to_string(
int(splineType))+
" not supported at this place");
632 template <
class XYContainer>
635 bool sortInputs =
true)
639 assert(points.size() > 1);
642 typename XYContainer::const_iterator it = points.begin();
643 typename XYContainer::const_iterator endIt = points.end();
644 for (
size_t i = 0; it != endIt; ++ i, ++it) {
654 if (splineType == Periodic)
656 else if (splineType == Natural)
658 else if (splineType == Monotonic)
661 throw std::runtime_error(
"Spline type "+std::to_string(
int(splineType))+
" not supported at this place");
678 template <
class XYContainer>
681 bool sortInputs =
true)
685 typename XYContainer::const_iterator it = points.begin();
686 typename XYContainer::const_iterator endIt = points.end();
687 for (
unsigned i = 0; it != endIt; ++i, ++it) {
688 xPos_[i] = std::get<0>(*it);
689 yPos_[i] = std::get<1>(*it);
697 if (splineType == Periodic)
699 else if (splineType == Natural)
701 else if (splineType == Monotonic)
704 throw std::runtime_error(
"Spline type "+std::to_string(
int(splineType))+
" not supported at this place");
710 template <
class Evaluation>
718 {
return xPos_.size(); }
723 Scalar
xAt(
size_t sampleIdx)
const
724 {
return x_(sampleIdx); }
730 {
return y_(sampleIdx); }
749 void printCSV(Scalar xi0, Scalar xi1,
size_t k, std::ostream& os = std::cout)
const
751 Scalar x0 = std::min(xi0, xi1);
752 Scalar x1 = std::max(xi0, xi1);
754 for (
size_t i = 0; i <= k; ++i) {
755 double x = i*(x1 - x0)/k + x0;
756 double x_p1 = x + (x1 - x0)/k;
763 y = (x -
x_(0))*dy_dx +
y_(0);
764 mono = (dy_dx>0)?1:-1;
766 else if (x >
x_(n)) {
768 y = (x -
x_(n))*dy_dx +
y_(n);
769 mono = (dy_dx>0)?1:-1;
772 throw std::runtime_error(
"The sampling points given to a spline must be sorted by their x value!");
781 os << x <<
" " << y <<
" " << dy_dx <<
" " << mono <<
"\n";
796 template <
class Evaluation>
797 Evaluation
eval(
const Evaluation& x,
bool extrapolate =
false)
const
799 if (!extrapolate && !
applies(x))
800 throw NumericalIssue(
"Tried to evaluate a spline outside of its range");
805 Scalar m = evalDerivative_(
xAt(0), 0);
807 return y0 + m*(x -
xAt(0));
809 else if (x >
xAt(
static_cast<size_t>(
static_cast<long int>(
numSamples()) - 1))) {
810 Scalar m = evalDerivative_(
xAt(
static_cast<size_t>(
numSamples() - 1)),
813 return y0 + m*(x -
xAt(
static_cast<size_t>(
numSamples() - 1)));
817 return eval_(x, segmentIdx_(scalarValue(x)));
832 template <
class Evaluation>
835 if (!extrapolate && !
applies(x))
836 throw NumericalIssue(
"Tried to evaluate the derivative of a spline outside of its range");
841 return evalDerivative_(
xAt(0), 0);
846 return evalDerivative_(x, segmentIdx_(scalarValue(x)));
861 template <
class Evaluation>
864 if (!extrapolate && !
applies(x))
865 throw NumericalIssue(
"Tried to evaluate the second derivative of a spline outside of its range");
866 else if (extrapolate)
869 return evalDerivative2_(x, segmentIdx_(scalarValue(x)));
884 template <
class Evaluation>
887 if (!extrapolate && !
applies(x))
888 throw NumericalIssue(
"Tried to evaluate the third derivative of a spline outside of its range");
889 else if (extrapolate)
892 return evalDerivative3_(x, segmentIdx_(scalarValue(x)));
901 template <
class Evaluation>
905 const Evaluation& d)
const
916 template <
class Evaluation>
921 const Evaluation& d)
const
925 Evaluation tmpSol[3], sol = 0;
927 size_t iFirst = segmentIdx_(x0);
928 size_t iLast = segmentIdx_(x1);
929 for (
size_t i = iFirst; i <= iLast; ++i)
937 throw std::runtime_error(
"Spline has more than one intersection");
942 throw std::runtime_error(
"Spline has no intersection");
956 [[maybe_unused]]
bool extrapolate =
false)
const
958 assert(std::abs(x0 - x1) > 1e-30);
969 Scalar m = evalDerivative_(
xAt(0), 0);
970 if (std::abs(m) < 1e-20)
975 size_t i = segmentIdx_(x0);
976 if (
x_(i + 1) >= x1) {
979 monotonic_(i, x0, x1, r);
985 monotonic_(i, x0,
x_(i+1), r);
990 size_t iEnd = segmentIdx_(x1);
991 for (; i < iEnd - 1; ++i) {
992 monotonic_(i,
x_(i),
x_(i + 1), r);
1001 assert(extrapolate);
1005 return (r < 0 || r==3)?-1:0;
1007 return (r > 0 || r==3)?1:0;
1013 monotonic_(iEnd,
x_(iEnd), x1, r);
1035 bool operator ()(
unsigned idxA,
unsigned idxB)
const
1036 {
return x_.at(idxA) < x_.at(idxB); }
1038 const std::vector<Scalar>& x_;
1049 std::vector<unsigned> idxVector(n);
1050 for (
unsigned i = 0; i < n; ++i)
1056 std::sort(idxVector.begin(), idxVector.end(), cmp);
1059 std::vector<Scalar> tmpX(n), tmpY(n);
1060 for (
size_t i = 0; i < idxVector.size(); ++ i) {
1061 tmpX[i] = xPos_[idxVector[i]];
1062 tmpY[i] = yPos_[idxVector[i]];
1076 for (
unsigned i = 0; i <= (n - 1)/2; ++i) {
1077 std::swap(xPos_[i], xPos_[n - i - 1]);
1078 std::swap(yPos_[i], yPos_[n - i - 1]);
1087 xPos_.resize(nSamples);
1088 yPos_.resize(nSamples);
1089 slopeVec_.resize(nSamples);
1107 M.
solve(moments, d);
1128 M.
solve(moments, d);
1152 M.
solve(moments, d);
1155 for (
int i =
static_cast<int>(
numSamples()) - 2; i >= 0; --i) {
1156 unsigned ui =
static_cast<unsigned>(i);
1157 moments[ui+1] = moments[ui];
1171 template <
class DestVector,
class SourceVector>
1174 const SourceVector& srcX,
1175 const SourceVector& srcY,
1178 assert(nSamples >= 2);
1182 for (
unsigned i = 0; i < nSamples; ++i) {
1184 if (srcX[0] > srcX[nSamples - 1])
1185 idx = nSamples - i - 1;
1186 destX[i] = srcX[idx];
1187 destY[i] = srcY[idx];
1191 template <
class DestVector,
class ListIterator>
1192 void assignFromArrayList_(DestVector& destX,
1194 const ListIterator& srcBegin,
1195 const ListIterator& srcEnd,
1198 assert(nSamples >= 2);
1201 ListIterator it = srcBegin;
1203 bool reverse =
false;
1204 if ((*srcBegin)[0] > (*it)[0])
1209 for (
unsigned i = 0; it != srcEnd; ++i, ++it) {
1212 idx = nSamples - i - 1;
1213 destX[idx] = (*it)[0];
1214 destY[idx] = (*it)[1];
1225 template <
class DestVector,
class ListIterator>
1228 ListIterator srcBegin,
1229 ListIterator srcEnd,
1232 assert(nSamples >= 2);
1238 ListIterator it = srcBegin;
1240 bool reverse =
false;
1241 if (std::get<0>(*srcBegin) > std::get<0>(*it))
1246 for (
unsigned i = 0; it != srcEnd; ++i, ++it) {
1249 idx = nSamples - i - 1;
1250 destX[idx] = std::get<0>(*it);
1251 destY[idx] = std::get<1>(*it);
1259 template <
class Vector,
class Matrix>
1267 d[0] = 6/
h_(1) * ( (
y_(1) -
y_(0))/
h_(1) - m0);
1276 (m1 - (
y_(n) -
y_(n - 1))/
h_(n));
1283 template <
class Vector,
class Matrix>
1293 for (
size_t i = 1; i < n; ++i) {
1294 Scalar lambda_i =
h_(i + 1) / (
h_(i) +
h_(i + 1));
1295 Scalar mu_i = 1 - lambda_i;
1297 6 / (
h_(i) +
h_(i + 1))
1299 ( (
y_(i + 1) -
y_(i))/
h_(i + 1) - (
y_(i) -
y_(i - 1))/
h_(i));
1303 M[i][i + 1] = lambda_i;
1308 Scalar lambda_0 = 0;
1329 template <
class Matrix,
class Vector>
1338 assert(M.
rows() == n);
1341 for (
size_t i = 2; i < n; ++i) {
1342 Scalar lambda_i =
h_(i + 1) / (
h_(i) +
h_(i + 1));
1343 Scalar mu_i = 1 - lambda_i;
1345 6 / (
h_(i) +
h_(i + 1))
1347 ( (
y_(i + 1) -
y_(i))/
h_(i + 1) - (
y_(i) -
y_(i - 1))/
h_(i));
1351 M[i-1][i] = lambda_i;
1355 Scalar lambda_n =
h_(1) / (
h_(n) +
h_(1));
1356 Scalar lambda_1 =
h_(2) / (
h_(1) +
h_(2));;
1357 Scalar mu_1 = 1 - lambda_1;
1358 Scalar mu_n = 1 - lambda_n;
1377 M[n-1][0] = lambda_n;
1391 template <
class Vector>
1397 std::vector<Scalar> delta(n);
1398 for (
size_t k = 0; k < n - 1; ++k)
1399 delta[k] = (
y_(k + 1) -
y_(k))/(
x_(k + 1) -
x_(k));
1402 for (
size_t k = 1; k < n - 1; ++k)
1403 slopes[k] = (delta[k - 1] + delta[k])/2;
1404 slopes[0] = delta[0];
1405 slopes[n - 1] = delta[n - 2];
1408 for (
size_t k = 0; k < n - 1; ++k) {
1409 if (std::abs(delta[k]) < 1e-50) {
1417 Scalar alpha = slopes[k] / delta[k];
1418 Scalar beta = slopes[k + 1] / delta[k];
1420 if (alpha < 0 || (k > 0 && slopes[k] / delta[k - 1] < 0)) {
1424 else if (alpha*alpha + beta*beta > 3*3) {
1425 Scalar tau = 3.0/std::sqrt(alpha*alpha + beta*beta);
1426 slopes[k] = tau*alpha*delta[k];
1427 slopes[k + 1] = tau*beta*delta[k];
1440 template <
class MomentsVector,
class SlopeVector>
1451 Scalar h = this->
h_(n - 1);
1456 (
y_(n - 1) -
y_(n - 2))/h
1458 h/6*(moments[n-1] - moments[n - 2]);
1463 moments[n - 1] * x*x / (2 * h)
1469 for (
size_t i = 0; i < n - 1; ++ i) {
1472 Scalar h_i = this->
h_(i + 1);
1477 (
y_(i+1) -
y_(i))/h_i
1479 h_i/6*(moments[i+1] - moments[i]);
1482 - moments[i] * x_i1*x_i1 / (2 * h_i)
1489 slopes[n - 1] = mRight;
1495 template <
class Evaluation>
1496 Evaluation eval_(
const Evaluation& x,
size_t i)
const
1499 Scalar delta =
h_(i + 1);
1500 Evaluation t = (x -
x_(i))/delta;
1504 + h10_(t) *
slope_(i)*delta
1505 + h01_(t) *
y_(i + 1)
1506 + h11_(t) *
slope_(i + 1)*delta;
1511 template <
class Evaluation>
1512 Evaluation evalDerivative_(
const Evaluation& x,
size_t i)
const
1515 Scalar delta =
h_(i + 1);
1516 Evaluation t = (x -
x_(i))/delta;
1517 Evaluation alpha = 1 / delta;
1521 (h00_prime_(t) *
y_(i)
1522 + h10_prime_(t) *
slope_(i)*delta
1523 + h01_prime_(t) *
y_(i + 1)
1524 + h11_prime_(t) *
slope_(i + 1)*delta);
1529 template <
class Evaluation>
1530 Evaluation evalDerivative2_(
const Evaluation& x,
size_t i)
const
1533 Scalar delta =
h_(i + 1);
1534 Evaluation t = (x -
x_(i))/delta;
1535 Evaluation alpha = 1 / delta;
1539 *(h00_prime2_(t) *
y_(i)
1540 + h10_prime2_(t) *
slope_(i)*delta
1541 + h01_prime2_(t) *
y_(i + 1)
1542 + h11_prime2_(t) *
slope_(i + 1)*delta);
1547 template <
class Evaluation>
1548 Evaluation evalDerivative3_(
const Evaluation& x,
size_t i)
const
1551 Scalar delta =
h_(i + 1);
1552 Evaluation t = (x -
x_(i))/delta;
1553 Evaluation alpha = 1 / delta;
1557 *(h00_prime3_(t)*
y_(i)
1558 + h10_prime3_(t)*
slope_(i)*delta
1559 + h01_prime3_(t)*
y_(i + 1)
1560 + h11_prime3_(t)*
slope_(i + 1)*delta);
1564 template <
class Evaluation>
1565 Evaluation h00_(
const Evaluation& t)
const
1566 {
return (2*t - 3)*t*t + 1; }
1568 template <
class Evaluation>
1569 Evaluation h10_(
const Evaluation& t)
const
1570 {
return ((t - 2)*t + 1)*t; }
1572 template <
class Evaluation>
1573 Evaluation h01_(
const Evaluation& t)
const
1574 {
return (-2*t + 3)*t*t; }
1576 template <
class Evaluation>
1577 Evaluation h11_(
const Evaluation& t)
const
1578 {
return (t - 1)*t*t; }
1581 template <
class Evaluation>
1582 Evaluation h00_prime_(
const Evaluation& t)
const
1583 {
return (3*2*t - 2*3)*t; }
1585 template <
class Evaluation>
1586 Evaluation h10_prime_(
const Evaluation& t)
const
1587 {
return (3*t - 2*2)*t + 1; }
1589 template <
class Evaluation>
1590 Evaluation h01_prime_(
const Evaluation& t)
const
1591 {
return (-3*2*t + 2*3)*t; }
1593 template <
class Evaluation>
1594 Evaluation h11_prime_(
const Evaluation& t)
const
1595 {
return (3*t - 2)*t; }
1598 template <
class Evaluation>
1599 Evaluation h00_prime2_(
const Evaluation& t)
const
1600 {
return 2*3*2*t - 2*3; }
1602 template <
class Evaluation>
1603 Evaluation h10_prime2_(
const Evaluation& t)
const
1604 {
return 2*3*t - 2*2; }
1606 template <
class Evaluation>
1607 Evaluation h01_prime2_(
const Evaluation& t)
const
1608 {
return -2*3*2*t + 2*3; }
1610 template <
class Evaluation>
1611 Evaluation h11_prime2_(
const Evaluation& t)
const
1612 {
return 2*3*t - 2; }
1615 template <
class Evaluation>
1616 Scalar h00_prime3_(
const Evaluation&)
const
1619 template <
class Evaluation>
1620 Scalar h10_prime3_(
const Evaluation&)
const
1623 template <
class Evaluation>
1624 Scalar h01_prime3_(
const Evaluation&)
const
1627 template <
class Evaluation>
1628 Scalar h11_prime3_(
const Evaluation&)
const
1639 int monotonic_(
size_t i, Scalar x0, Scalar x1,
int& r)
const
1646 if (std::abs(a) < 1e-20 && std::abs(b) < 1e-20 && std::abs(c) < 1e-20)
1649 Scalar disc = b*b - 4*a*c;
1653 if (x0*(x0*a + b) + c > 0) {
1654 r = (r==3 || r == 1)?1:0;
1658 r = (r==3 || r == -1)?-1:0;
1662 disc = std::sqrt(disc);
1663 Scalar xE1 = (-b + disc)/(2*a);
1664 Scalar xE2 = (-b - disc)/(2*a);
1666 if (std::abs(disc) < 1e-30) {
1668 if (std::abs(xE1 - x0) < 1e-30)
1673 if (x0*(x0*a + b) + c > 0) {
1674 r = (r==3 || r == 1)?1:0;
1678 r = (r==3 || r == -1)?-1:0;
1682 if ((x0 < xE1 && xE1 < x1) ||
1683 (x0 < xE2 && xE2 < x1))
1692 if (x0*(x0*a + b) + c > 0) {
1693 r = (r==3 || r == 1)?1:0;
1697 r = (r==3 || r == -1)?-1:0;
1706 template <
class Evaluation>
1709 const Evaluation& a,
1710 const Evaluation& b,
1711 const Evaluation& c,
1712 const Evaluation& d,
1713 Scalar x0 = -1e30, Scalar x1 = 1e30)
const
1721 x0 = std::max(
x_(segIdx), x0);
1722 x1 = std::min(
x_(segIdx+1), x1);
1726 for (
unsigned j = 0; j < n; ++j) {
1727 if (x0 <= sol[j] && sol[j] <= x1) {
1736 size_t segmentIdx_(Scalar x)
const
1742 while (iLow + 1 < iHigh) {
1743 size_t i = (iLow + iHigh) / 2;
1755 Scalar
h_(
size_t i)
const
1757 assert(
x_(i) >
x_(i-1));
1759 return x_(i) -
x_(i - 1);
1765 Scalar
x_(
size_t i)
const
1766 {
return xPos_[i]; }
1771 Scalar
y_(
size_t i)
const
1772 {
return yPos_[i]; }
1779 {
return slopeVec_[i]; }
1783 Scalar a_(
size_t i)
const
1784 {
return evalDerivative3_(Scalar(0.0), i)/6.0; }
1788 Scalar b_(
size_t i)
const
1789 {
return evalDerivative2_(Scalar(0.0), i)/2.0; }
1793 Scalar c_(
size_t i)
const
1794 {
return evalDerivative_(Scalar(0.0), i); }
1798 Scalar d_(
size_t i)
const
1799 {
return eval_(Scalar(0.0), i); }
Provides the opm-material specific exception classes.
Provides free functions to invert polynomials of degree 1, 2 and 3.
unsigned invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:148
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left.
Definition: Exceptions.hpp:46
Class implementing cubic splines.
Definition: Spline.hpp:91
size_t numSamples() const
Return the number of (x, y) values.
Definition: Spline.hpp:717
void makeFullSystem_(Matrix &M, Vector &d, Scalar m0, Scalar m1)
Make the linear system of equations Mx = d which results in the moments of the full spline.
Definition: Spline.hpp:1260
void setArrayOfPoints(size_t nSamples, const PointArray &points, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points of a natural spline using a C-style array.
Definition: Spline.hpp:591
void setXYArrays(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using C-style arrays.
Definition: Spline.hpp:322
Scalar valueAt(size_t sampleIdx) const
Return the x value of a given sampling point.
Definition: Spline.hpp:729
void makeNaturalSpline_()
Create a natural spline from the already set sampling points.
Definition: Spline.hpp:1118
void sortInput_()
Sort the sample points in ascending order of their x value.
Definition: Spline.hpp:1044
Spline(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:140
Spline(const PointContainer &points, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:249
void setSlopesFromMoments_(SlopeVector &slopes, const MomentsVector &moments)
Convert the moments at the sample points to slopes.
Definition: Spline.hpp:1441
void setContainerOfTuples(const XYContainer &points, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using a STL-compatible container of ...
Definition: Spline.hpp:469
void makePeriodicSystem_(Matrix &M, Vector &d)
Make the linear system of equations Mx = d which results in the moments of the periodic spline.
Definition: Spline.hpp:1330
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition: Spline.hpp:711
Spline(const ScalarContainer &x, const ScalarContainer &y, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:169
Spline(Scalar x0, Scalar x1, Scalar y0, Scalar y1, Scalar m0, Scalar m1)
Convenience constructor for a full spline with just two sampling points.
Definition: Spline.hpp:126
Evaluation intersect(const Evaluation &a, const Evaluation &b, const Evaluation &c, const Evaluation &d) const
Find the intersections of the spline with a cubic polynomial in the whole interval,...
Definition: Spline.hpp:902
Scalar y_(size_t i) const
Returns the y coordinate of the i-th sampling point.
Definition: Spline.hpp:1771
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's derivative at a given position.
Definition: Spline.hpp:833
SplineType
The type of the spline to be created.
Definition: Spline.hpp:101
Spline(const PointContainer &points, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:182
int monotonic(Scalar x0, Scalar x1, [[maybe_unused]] bool extrapolate=false) const
Returns 1 if the spline is monotonically increasing, -1 if the spline is mononously decreasing and 0 ...
Definition: Spline.hpp:955
void makeNaturalSystem_(Matrix &M, Vector &d)
Make the linear system of equations Mx = d which results in the moments of the natural spline.
Definition: Spline.hpp:1284
Scalar h_(size_t i) const
Returns x[i] - x[i - 1].
Definition: Spline.hpp:1755
void set(Scalar x0, Scalar x1, Scalar y0, Scalar y1, Scalar m0, Scalar m1)
Set the sampling points and the boundary slopes of the spline with two sampling points.
Definition: Spline.hpp:266
Spline(const ScalarContainerX &x, const ScalarContainerY &y, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:233
void setContainerOfPoints(const XYContainer &points, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points of a natural spline using a STL-compatible container of array-like objects.
Definition: Spline.hpp:633
Evaluation evalThirdDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's third derivative at a given position.
Definition: Spline.hpp:885
void printCSV(Scalar xi0, Scalar xi1, size_t k, std::ostream &os=std::cout) const
Prints k tuples of the format (x, y, dx/dy, isMonotonic) to stdout.
Definition: Spline.hpp:749
Spline(size_t nSamples, const ScalarArray &x, const ScalarArray &y, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:198
void assignSamplingPoints_(DestVector &destX, DestVector &destY, const SourceVector &srcX, const SourceVector &srcY, unsigned nSamples)
Set the sampling point vectors.
Definition: Spline.hpp:1172
Scalar xAt(size_t sampleIdx) const
Return the x value of a given sampling point.
Definition: Spline.hpp:723
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the spline as interval.
Definition: Spline.hpp:1022
void assignFromTupleList_(DestVector &destX, DestVector &destY, ListIterator srcBegin, ListIterator srcEnd, unsigned nSamples)
Set the sampling points.
Definition: Spline.hpp:1226
Spline(size_t nSamples, const PointArray &points, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:216
Scalar slope_(size_t i) const
Returns the slope (i.e.
Definition: Spline.hpp:1778
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Spline.hpp:797
size_t intersectSegment_(Evaluation *sol, size_t segIdx, const Evaluation &a, const Evaluation &b, const Evaluation &c, const Evaluation &d, Scalar x0=-1e30, Scalar x1=1e30) const
Find all the intersections of a segment of the spline with a cubic polynomial within a specified inte...
Definition: Spline.hpp:1707
Scalar x_(size_t i) const
Returns the y coordinate of the i-th sampling point.
Definition: Spline.hpp:1765
void setArrayOfPoints(size_t nSamples, const PointArray &points, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using a C-style array.
Definition: Spline.hpp:390
void makeFullSpline_(Scalar m0, Scalar m1)
Create a natural spline from the already set sampling points.
Definition: Spline.hpp:1097
Evaluation evalSecondDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's second derivative at a given position.
Definition: Spline.hpp:862
void makeMonotonicSpline_(Vector &slopes)
Create a monotonic spline from the already set sampling points.
Definition: Spline.hpp:1392
void setNumSamples_(size_t nSamples)
Resizes the internal vectors to store the sample points.
Definition: Spline.hpp:1085
void setXYArrays(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points natural spline using C-style arrays.
Definition: Spline.hpp:510
Spline()
Default constructor for a spline.
Definition: Spline.hpp:113
void reverseSamplingPoints_()
Reverse order of the elements in the arrays which contain the sampling points.
Definition: Spline.hpp:1072
void setContainerOfTuples(const XYContainer &points, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points of a natural spline using a STL-compatible container of tuple-like objects.
Definition: Spline.hpp:679
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points of a natural spline using STL-compatible containers.
Definition: Spline.hpp:551
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using STL-compatible containers.
Definition: Spline.hpp:356
void makePeriodicSpline_()
Create a periodic spline from the already set sampling points.
Definition: Spline.hpp:1139
void setContainerOfPoints(const XYContainer &points, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using a STL-compatible container of ...
Definition: Spline.hpp:427
Evaluation intersectInterval(Scalar x0, Scalar x1, const Evaluation &a, const Evaluation &b, const Evaluation &c, const Evaluation &d) const
Find the intersections of the spline with a cubic polynomial in a sub-interval of the spline,...
Definition: Spline.hpp:917
Spline(size_t nSamples, const PointArray &points, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:155
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left.
Definition: TridiagonalMatrix.hpp:51
size_t rows() const
Return the number of rows of the matrix.
Definition: TridiagonalMatrix.hpp:140
void solve(XVector &x, const BVector &b) const
Calculate the solution for a linear system of equations.
Definition: TridiagonalMatrix.hpp:714
Helper class needed to sort the input sampling points.
Definition: Spline.hpp:1030