29#ifndef quantlib_cubic_interpolation_hpp
30#define quantlib_cubic_interpolation_hpp
157 template <
class I1,
class I2>
164 Real leftConditionValue,
166 Real rightConditionValue) {
167 impl_ = ext::shared_ptr<Interpolation::Impl>(
new
174 rightConditionValue));
198 template <
class I1,
class I2>
211 template <
class I1,
class I2>
224 template <
class I1,
class I2>
237 template <
class I1,
class I2>
250 template <
class I1,
class I2>
263 template <
class I1,
class I2>
276 template <
class I1,
class I2>
289 template <
class I1,
class I2>
302 template <
class I1,
class I2>
315 template <
class I1,
class I2>
331 bool monotonic =
false,
334 Real leftConditionValue = 0.0,
337 Real rightConditionValue = 0.0)
341 template <
class I1,
class I2>
344 const I2& yBegin)
const {
362 template <
class I1,
class I2>
372 Real leftConditionValue,
374 Real rightConditionValue)
377 Cubic::requiredPoints),
387 "Lagrange boundary condition requires at least "
388 "4 points (" << (xEnd-xBegin) <<
" are given)");
394 for (
Size i=0; i<
n_-1; ++i) {
401 for (
Size i=1; i<
n_-1; ++i) {
424 QL_FAIL(
"this end condition is not implemented yet");
435 QL_FAIL(
"unknown end condition");
456 QL_FAIL(
"this end condition is not implemented yet");
467 QL_FAIL(
"unknown end condition");
474 for (
Size i=0; i<
n_-2; ++i) {
477 T_[i][i+2]=
dx_[i+1]/6.0;
480 for (
Size i=0; i<
n_-2; ++i) {
482 S_[i][i+1]=-(1.0/
dx_[i+1]+1.0/
dx_[i]);
483 S_[i][i+2]=1.0/
dx_[i+1];
489 for (
Size i=0; i<
n_-2; ++i)
500 for (
Size i=1; i<
n_-1; ++i) {
501 Q_[i][i-1]=7.0/8*1.0/(
n_-1)*
dx_[i-1]*
dx_[i-1]*
dx_[i-1];
512 for (
Size i=0; i<
n_-1; ++i)
513 tmp_[i]=(Y_[i+1]-Y_[i])/
dx_[i]-(2.0*D_[i]+D_[i+1])*
dx_[i]/6.0;
518 for (
Size i=0; i<
n_-2; ++i) {
521 T_[i][i+2]=
dx_[i+1]/6.0;
524 for (
Size i=0; i<
n_-2; ++i) {
526 S_[i][i+1]=-(1.0/
dx_[i+1]+1.0/
dx_[i]);
527 S_[i][i+2]=1.0/
dx_[i+1];
533 for (
Size i=0; i<
n_-2; ++i)
542 Q_[0][0]=1.0/(
n_-1)*
dx_[0];
543 Q_[0][1]=1.0/2*1.0/(
n_-1)*
dx_[0];
544 for (
Size i=1; i<
n_-1; ++i) {
545 Q_[i][i-1]=1.0/2*1.0/(
n_-1)*
dx_[i-1];
546 Q_[i][i]=1.0/(
n_-1)*
dx_[i]+1.0/(
n_-1)*
dx_[i-1];
547 Q_[i][i+1]=1.0/2*1.0/(
n_-1)*
dx_[i];
556 for (
Size i=0; i<
n_-1; ++i)
557 tmp_[i]=(Y_[i+1]-Y_[i])/
dx_[i]-(2.0*D_[i]+D_[i+1])*
dx_[i]/6.0;
565 QL_FAIL(
"FourthOrder not implemented yet");
569 for (
Size i=1; i<
n_-1; ++i)
577 for (
Size i=1; i<
n_-1; ++i) {
578 Real Smin = std::min(
S_[i-1],
S_[i]);
579 Real Smax = std::max(
S_[i-1],
S_[i]);
580 if(Smax+2.0*Smin == 0){
583 else if (Smin*Smax == 0)
589 tmp_[i] = 3.0*Smin*Smax/(Smax+2.0*Smin);
598 for (
Size i=2; i<
n_-2; ++i) {
599 if ((
S_[i-2]==
S_[i-1]) && (
S_[i]!=
S_[i+1]))
601 else if ((
S_[i-2]!=
S_[i-1]) && (
S_[i]==
S_[i+1]))
603 else if (
S_[i]==
S_[i-1])
605 else if ((
S_[i-2]==
S_[i-1]) && (
S_[i-1]!=
S_[i]) && (
S_[i]==
S_[i+1]))
608 tmp_[i] = (std::abs(
S_[i+1]-
S_[i])*
S_[i-1]+std::abs(
S_[i-1]-
S_[i-2])*
S_[i])/(std::abs(
S_[i+1]-
S_[i])+std::abs(
S_[i-1]-
S_[i-2]));
615 for (
Size i=1; i<
n_-1; ++i) {
616 if (
S_[i-1]*
S_[i]<0.0)
623 tmp_[i] = 2.0/(1.0/
S_[i-1]+1.0/
S_[i]);
631 for (
Size i=1; i<
n_-1; ++i) {
634 if (
S_[i-1]*
S_[i]<=0.0)
640 tmp_[i] = (w1+w2)/(w1/
S_[i-1]+w2/
S_[i]);
647 else if (
S_[0]*
S_[1]<0) {
648 if (std::fabs(
tmp_[0])>std::fabs(3*
S_[0])) {
658 if (std::fabs(
tmp_[
n_-1])>std::fabs(3*
S_[
n_-2])) {
675 for (
Size i=0; i<
n_; ++i) {
678 correction =
tmp_[i]/std::fabs(
tmp_[i]) *
679 std::min<Real>(std::fabs(
tmp_[i]),
680 std::fabs(3.0*
S_[0]));
684 if (correction!=
tmp_[i]) {
685 tmp_[i] = correction;
688 }
else if (i==
n_-1) {
690 correction =
tmp_[i]/std::fabs(
tmp_[i]) *
691 std::min<Real>(std::fabs(
tmp_[i]),
692 std::fabs(3.0*
S_[
n_-2]));
696 if (correction!=
tmp_[i]) {
697 tmp_[i] = correction;
703 M = 3.0 * std::min(std::min(std::fabs(
S_[i-1]),
707 if ((
S_[i-1]-
S_[i-2])*(
S_[i]-
S_[i-1])>0.0) {
711 if (pm*pd>0.0 && pm*(
S_[i-1]-
S_[i-2])>0.0) {
712 M = std::max<Real>(M, 1.5*std::min(
713 std::fabs(pm),std::fabs(pd)));
718 if ((
S_[i]-
S_[i-1])*(
S_[i+1]-
S_[i])>0.0) {
721 if (pm*pu>0.0 && -pm*(
S_[i]-
S_[i-1])>0.0) {
722 M = std::max<Real>(M, 1.5*std::min(
723 std::fabs(pm),std::fabs(pu)));
727 if (
tmp_[i]*pm>0.0) {
728 correction =
tmp_[i]/std::fabs(
tmp_[i]) *
729 std::min(std::fabs(
tmp_[i]), M);
733 if (correction!=
tmp_[i]) {
734 tmp_[i] = correction;
743 for (
Size i=0; i<
n_-1; ++i) {
750 for (
Size i=1; i<
n_-1; ++i) {
754 (
a_[i-1]/2.0 +
dx_[i-1] *
755 (
b_[i-1]/3.0 +
dx_[i-1] *
c_[i-1]/4.0)));
778 return 2.0*
b_[j] + 6.0*
c_[j]*
dx_;
793 return (-((((a-c)*(
b-c)*(c-x)*z-(a-
d)*(
b-
d)*(
d-x)*w)*(a-x+
b-x)
794 +((a-c)*(
b-c)*z-(a-
d)*(
b-
d)*w)*(a-x)*(
b-x))*(a-
b)+
795 ((a-c)*(a-
d)*
v-(
b-c)*(
b-
d)*u)*(c-
d)*(c-x)*(
d-x)
796 +((a-c)*(a-
d)*(a-x)*
v-(
b-c)*(
b-
d)*(
b-x)*u)
798 ((a-
b)*(a-c)*(a-
d)*(
b-c)*(
b-
d)*(c-
d));
AkimaCubicInterpolation(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
1-D array used in linear algebra.
Cubic interpolation factory and traits
CubicInterpolation::DerivativeApprox da_
CubicInterpolation::BoundaryCondition rightType_
CubicInterpolation::BoundaryCondition leftType_
static const Size requiredPoints
Cubic(CubicInterpolation::DerivativeApprox da=CubicInterpolation::Kruger, bool monotonic=false, CubicInterpolation::BoundaryCondition leftCondition=CubicInterpolation::SecondDerivative, Real leftConditionValue=0.0, CubicInterpolation::BoundaryCondition rightCondition=CubicInterpolation::SecondDerivative, Real rightConditionValue=0.0)
Interpolation interpolate(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin) const
Cubic interpolation between discrete points.
const std::vector< Real > & cCoefficients() const
const std::vector< bool > & monotonicityAdjustments() const
const std::vector< Real > & bCoefficients() const
CubicInterpolation(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, CubicInterpolation::DerivativeApprox da, bool monotonic, CubicInterpolation::BoundaryCondition leftCond, Real leftConditionValue, CubicInterpolation::BoundaryCondition rightCond, Real rightConditionValue)
const std::vector< Real > & aCoefficients() const
const std::vector< Real > & primitiveConstants() const
@ Kruger
Kruger approximation (local, monotonic, non-linear)
@ SplineOM2
Overshooting minimization 2nd derivative.
@ Parabolic
Parabolic approximation (local, non-monotonic, linear)
@ Harmonic
Weighted harmonic mean approximation (local, monotonic, non-linear)
@ FourthOrder
Fourth-order approximation (local, non-monotonic, linear)
@ SplineOM1
Overshooting minimization 1st derivative.
@ Akima
Akima approximation (local, non-monotonic, non-linear)
@ FritschButland
Fritsch-Butland approximation (local, monotonic, non-linear)
@ SecondDerivative
Match value of second derivative at end.
@ Periodic
Match first and second derivative at either end.
@ NotAKnot
Make second(-last) point an inactive knot.
@ FirstDerivative
Match value of end-slope.
const detail::CoefficientHolder & coeffs() const
CubicNaturalSpline(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
CubicSplineOvershootingMinimization1(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
CubicSplineOvershootingMinimization2(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
FritschButlandCubic(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
HarmonicCubic(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
basic template implementation
Size locate(Real x) const
templateImpl(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, const int requiredPoints=2)
base class for 1-D interpolations.
ext::shared_ptr< Impl > impl_
KrugerCubic(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
Matrix used in linear algebra.
MonotonicCubicNaturalSpline(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
MonotonicParabolic(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
Parabolic(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
Base implementation for tridiagonal operator.
Array solveFor(const Array &rhs) const
solve linear system for a given right-hand side
void setLastRow(Real, Real)
void setFirstRow(Real, Real)
void setMidRow(Size, Real, Real, Real)
std::vector< bool > monotonicityAdjustments_
virtual ~CoefficientHolder()=default
std::vector< Real > primitiveConst_
CoefficientHolder(Size n)
CubicInterpolation::DerivativeApprox da_
Real cubicInterpolatingPolynomialDerivative(Real a, Real b, Real c, Real d, Real u, Real v, Real w, Real z, Real x) const
Real value(Real x) const override
CubicInterpolation::BoundaryCondition rightType_
Real derivative(Real x) const override
CubicInterpolation::BoundaryCondition leftType_
CubicInterpolationImpl(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, CubicInterpolation::DerivativeApprox da, bool monotonic, CubicInterpolation::BoundaryCondition leftCondition, Real leftConditionValue, CubicInterpolation::BoundaryCondition rightCondition, Real rightConditionValue)
Real primitive(Real x) const override
Real secondDerivative(Real x) const override
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
#define QL_FAIL(message)
throw an error (possibly with file and line information)
ext::function< Real(Real)> b
std::size_t Size
size of a container
base class for 1-D interpolations
matrix used in linear algebra.
Matrix transpose(const Matrix &m)
Matrix inverse(const Matrix &m)
ext::shared_ptr< BlackVolTermStructure > v