[ceph.git] / ceph / src / boost / boost / math / interpolators / detail / cardinal_quintic_b_spline_detail.hpp

// Copyright Nick Thompson, 2019
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP
#define BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP
#include <cmath>
#include <vector>
#include <utility>
#include <boost/math/special_functions/cardinal_b_spline.hpp>

namespace boost{ namespace math{ namespace interpolators{ namespace detail{


template <class Real>
class cardinal_quintic_b_spline_detail
{
public:
    // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
    // y[0] = y(a), y[n -1] = y(b), step_size = (b - a)/(n -1).

    cardinal_quintic_b_spline_detail(const Real* const y,
                                     size_t n,
                                     Real t0 /* initial time, left endpoint */,
                                     Real h  /*spacing, stepsize*/,
                                     std::pair<Real, Real> left_endpoint_derivatives,
                                     std::pair<Real, Real> right_endpoint_derivatives)
    {
        static_assert(!std::is_integral<Real>::value, "The quintic B-spline interpolator only works with floating point types.");
        if (h <= 0) {
            throw std::logic_error("Spacing must be > 0.");
        }
        m_inv_h = 1/h;
        m_t0 = t0;

        if (n < 8) {
            throw std::logic_error("The quintic B-spline interpolator requires at least 8 points.");
        }

        using std::isnan;
        // This interpolator has error of order h^6, so the derivatives should be estimated with the same error.
        // See: https://en.wikipedia.org/wiki/Finite_difference_coefficient
        if (isnan(left_endpoint_derivatives.first)) {
            Real tmp = -49*y[0]/20 + 6*y[1] - 15*y[2]/2 + 20*y[3]/3 - 15*y[4]/4 + 6*y[5]/5 - y[6]/6;
            left_endpoint_derivatives.first = tmp/h;
        }
        if (isnan(right_endpoint_derivatives.first)) {
            Real tmp = 49*y[n-1]/20 - 6*y[n-2] + 15*y[n-3]/2 - 20*y[n-4]/3 + 15*y[n-5]/4 - 6*y[n-6]/5 + y[n-7]/6;
            right_endpoint_derivatives.first = tmp/h;
        }
        if(isnan(left_endpoint_derivatives.second)) {
            Real tmp = 469*y[0]/90 - 223*y[1]/10 + 879*y[2]/20 - 949*y[3]/18 + 41*y[4] - 201*y[5]/10 + 1019*y[6]/180 - 7*y[7]/10;
            left_endpoint_derivatives.second = tmp/(h*h);
        }
        if (isnan(right_endpoint_derivatives.second)) {
            Real tmp = 469*y[n-1]/90 - 223*y[n-2]/10 + 879*y[n-3]/20 - 949*y[n-4]/18 + 41*y[n-5] - 201*y[n-6]/10 + 1019*y[n-7]/180 - 7*y[n-8]/10;
            right_endpoint_derivatives.second = tmp/(h*h);
        }

        // This is really challenging my mental limits on by-hand row reduction.
        // I debated bringing in a dependency on a sparse linear solver, but given that that would cause much agony for users I decided against it.

        std::vector<Real> rhs(n+4);
        rhs[0] = 20*y[0] - 12*h*left_endpoint_derivatives.first +  2*h*h*left_endpoint_derivatives.second;
        rhs[1] = 60*y[0] - 12*h*left_endpoint_derivatives.first;
        for (size_t i = 2; i < n + 2; ++i) {
            rhs[i] = 120*y[i-2];
        }
        rhs[n+2] = 60*y[n-1] + 12*h*right_endpoint_derivatives.first;
        rhs[n+3] = 20*y[n-1] + 12*h*right_endpoint_derivatives.first +  2*h*h*right_endpoint_derivatives.second;

        std::vector<Real> diagonal(n+4, 66);
        diagonal[0] = 1;
        diagonal[1] = 18;
        diagonal[n+2] = 18;
        diagonal[n+3] = 1;

        std::vector<Real> first_superdiagonal(n+4, 26);
        first_superdiagonal[0] = 10;
        first_superdiagonal[1] = 33;
        first_superdiagonal[n+2] = 1;
        // There is one less superdiagonal than diagonal; make sure that if we read it, it shows up as a bug:
        first_superdiagonal[n+3] = std::numeric_limits<Real>::quiet_NaN();

        std::vector<Real> second_superdiagonal(n+4, 1);
        second_superdiagonal[0] = 9;
        second_superdiagonal[1] = 8;
        second_superdiagonal[n+2] = std::numeric_limits<Real>::quiet_NaN();
        second_superdiagonal[n+3] = std::numeric_limits<Real>::quiet_NaN();

        std::vector<Real> first_subdiagonal(n+4, 26);
        first_subdiagonal[0] = std::numeric_limits<Real>::quiet_NaN();
        first_subdiagonal[1] = 1;
        first_subdiagonal[n+2] = 33;
        first_subdiagonal[n+3] = 10;

        std::vector<Real> second_subdiagonal(n+4, 1);
        second_subdiagonal[0] = std::numeric_limits<Real>::quiet_NaN();
        second_subdiagonal[1] = std::numeric_limits<Real>::quiet_NaN();
        second_subdiagonal[n+2] = 8;
        second_subdiagonal[n+3] = 9;


        for (size_t i = 0; i < n+2; ++i) {
            Real di = diagonal[i];
            diagonal[i] = 1;
            first_superdiagonal[i] /= di;
            second_superdiagonal[i] /= di;
            rhs[i] /= di;

            // Eliminate first subdiagonal:
            Real nfsub = -first_subdiagonal[i+1];
            // Superfluous:
            first_subdiagonal[i+1] /= nfsub;
            // Not superfluous:
            diagonal[i+1] /= nfsub;
            first_superdiagonal[i+1] /= nfsub;
            second_superdiagonal[i+1] /= nfsub;
            rhs[i+1] /= nfsub;

            diagonal[i+1] += first_superdiagonal[i];
            first_superdiagonal[i+1] += second_superdiagonal[i];
            rhs[i+1] += rhs[i];
            // Superfluous, but clarifying:
            first_subdiagonal[i+1] = 0;

            // Eliminate second subdiagonal:
            Real nssub = -second_subdiagonal[i+2];
            first_subdiagonal[i+2] /= nssub;
            diagonal[i+2] /= nssub;
            first_superdiagonal[i+2] /= nssub;
            second_superdiagonal[i+2] /= nssub;
            rhs[i+2] /= nssub;

            first_subdiagonal[i+2] += first_superdiagonal[i];
            diagonal[i+2] += second_superdiagonal[i];
            rhs[i+2] += rhs[i];
            // Superfluous, but clarifying:
            second_subdiagonal[i+2] = 0;
        }

        // Eliminate last subdiagonal:
        Real dnp2 = diagonal[n+2];
        diagonal[n+2] = 1;
        first_superdiagonal[n+2] /= dnp2;
        rhs[n+2] /= dnp2;
        Real nfsubnp3 = -first_subdiagonal[n+3];
        diagonal[n+3] /= nfsubnp3;
        rhs[n+3] /= nfsubnp3;

        diagonal[n+3] += first_superdiagonal[n+2];
        rhs[n+3] += rhs[n+2];

        m_alpha.resize(n + 4, std::numeric_limits<Real>::quiet_NaN());

        m_alpha[n+3] = rhs[n+3]/diagonal[n+3];
        m_alpha[n+2] = rhs[n+2] - first_superdiagonal[n+2]*m_alpha[n+3];
        for (int64_t i = int64_t(n+1); i >= 0; --i) {
            m_alpha[i] = rhs[i] - first_superdiagonal[i]*m_alpha[i+1] - second_superdiagonal[i]*m_alpha[i+2];
        }

    }

    Real operator()(Real t) const {
        using std::ceil;
        using std::floor;
        using boost::math::cardinal_b_spline;
        // tf = t0 + (n-1)*h
        // alpha.size() = n+4
        if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) {
            const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work.";
            throw std::domain_error(err_msg);
        }
        Real x = (t-m_t0)*m_inv_h;
        // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5.
        // TODO: Zero pad m_alpha so that only the domain check is necessary.
        int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1)));
        int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) );
        Real s = 0;
        for (int64_t j = j_min; j <= j_max; ++j) {
            // TODO: Use Cox 1972 to generate all integer translates of B5 simultaneously.
            s += m_alpha[j]*cardinal_b_spline<5, Real>(x - j + 2);
        }
        return s;
    }

    Real prime(Real t) const {
        using std::ceil;
        using std::floor;
        using boost::math::cardinal_b_spline_prime;
        if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) {
            const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work.";
            throw std::domain_error(err_msg);
        }
        Real x = (t-m_t0)*m_inv_h;
        // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5
        int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1)));
        int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) );
        Real s = 0;
        for (int64_t j = j_min; j <= j_max; ++j) {
            s += m_alpha[j]*cardinal_b_spline_prime<5, Real>(x - j + 2);
        }
        return s*m_inv_h;

    }

    Real double_prime(Real t) const {
        using std::ceil;
        using std::floor;
        using boost::math::cardinal_b_spline_double_prime;
        if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) {
            const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work.";
            throw std::domain_error(err_msg);
        }
        Real x = (t-m_t0)*m_inv_h;
        // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5
        int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1)));
        int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) );
        Real s = 0;
        for (int64_t j = j_min; j <= j_max; ++j) {
            s += m_alpha[j]*cardinal_b_spline_double_prime<5, Real>(x - j + 2);
        }
        return s*m_inv_h*m_inv_h;
    }


    Real t_max() const {
        return m_t0 + (m_alpha.size()-5)/m_inv_h;
    }

private:
    std::vector<Real> m_alpha;
    Real m_inv_h;
    Real m_t0;
};

}}}}
#endif
Commit	Line	Data
92f5a8d4 TL	1	// Copyright Nick Thompson, 2019
	2	// Use, modification and distribution are subject to the
	3	// Boost Software License, Version 1.0.
	4	// (See accompanying file LICENSE_1_0.txt
	5	// or copy at http://www.boost.org/LICENSE_1_0.txt)
	6
	7	#ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP
	8	#define BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP
	9	#include <cmath>
	10	#include <vector>
	11	#include <utility>
	12	#include <boost/math/special_functions/cardinal_b_spline.hpp>
	13
	14	namespace boost{ namespace math{ namespace interpolators{ namespace detail{
	15
	16
	17	template <class Real>
	18	class cardinal_quintic_b_spline_detail
	19	{
	20	public:
	21	// If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
	22	// y[0] = y(a), y[n -1] = y(b), step_size = (b - a)/(n -1).
	23
	24	cardinal_quintic_b_spline_detail(const Real* const y,
	25	size_t n,
	26	Real t0 /* initial time, left endpoint */,
	27	Real h /spacing, stepsize/,
	28	std::pair<Real, Real> left_endpoint_derivatives,
	29	std::pair<Real, Real> right_endpoint_derivatives)
	30	{
	31	static_assert(!std::is_integral<Real>::value, "The quintic B-spline interpolator only works with floating point types.");
	32	if (h <= 0) {
	33	throw std::logic_error("Spacing must be > 0.");
	34	}
	35	m_inv_h = 1/h;
	36	m_t0 = t0;
	37
	38	if (n < 8) {
f67539c2	39	throw std::logic_error("The quintic B-spline interpolator requires at least 8 points.");
92f5a8d4 TL	40	}
	41
	42	using std::isnan;
	43	// This interpolator has error of order h^6, so the derivatives should be estimated with the same error.
	44	// See: https://en.wikipedia.org/wiki/Finite_difference_coefficient
	45	if (isnan(left_endpoint_derivatives.first)) {
	46	Real tmp = -49y[0]/20 + 6y[1] - 15y[2]/2 + 20y[3]/3 - 15y[4]/4 + 6y[5]/5 - y[6]/6;
	47	left_endpoint_derivatives.first = tmp/h;
	48	}
	49	if (isnan(right_endpoint_derivatives.first)) {
	50	Real tmp = 49y[n-1]/20 - 6y[n-2] + 15y[n-3]/2 - 20y[n-4]/3 + 15y[n-5]/4 - 6y[n-6]/5 + y[n-7]/6;
	51	right_endpoint_derivatives.first = tmp/h;
	52	}
	53	if(isnan(left_endpoint_derivatives.second)) {
	54	Real tmp = 469y[0]/90 - 223y[1]/10 + 879y[2]/20 - 949y[3]/18 + 41y[4] - 201y[5]/10 + 1019y[6]/180 - 7y[7]/10;
	55	left_endpoint_derivatives.second = tmp/(h*h);
	56	}
	57	if (isnan(right_endpoint_derivatives.second)) {
	58	Real tmp = 469y[n-1]/90 - 223y[n-2]/10 + 879y[n-3]/20 - 949y[n-4]/18 + 41y[n-5] - 201y[n-6]/10 + 1019y[n-7]/180 - 7y[n-8]/10;
	59	right_endpoint_derivatives.second = tmp/(h*h);
	60	}
	61
	62	// This is really challenging my mental limits on by-hand row reduction.
	63	// I debated bringing in a dependency on a sparse linear solver, but given that that would cause much agony for users I decided against it.
	64
	65	std::vector<Real> rhs(n+4);
	66	rhs[0] = 20y[0] - 12hleft_endpoint_derivatives.first + 2hhleft_endpoint_derivatives.second;
	67	rhs[1] = 60y[0] - 12h*left_endpoint_derivatives.first;
	68	for (size_t i = 2; i < n + 2; ++i) {
	69	rhs[i] = 120*y[i-2];
	70	}
	71	rhs[n+2] = 60y[n-1] + 12h*right_endpoint_derivatives.first;
	72	rhs[n+3] = 20y[n-1] + 12hright_endpoint_derivatives.first + 2hhright_endpoint_derivatives.second;
	73
	74	std::vector<Real> diagonal(n+4, 66);
	75	diagonal[0] = 1;
	76	diagonal[1] = 18;
	77	diagonal[n+2] = 18;
	78	diagonal[n+3] = 1;
	79
	80	std::vector<Real> first_superdiagonal(n+4, 26);
	81	first_superdiagonal[0] = 10;
	82	first_superdiagonal[1] = 33;
	83	first_superdiagonal[n+2] = 1;
	84	// There is one less superdiagonal than diagonal; make sure that if we read it, it shows up as a bug:
	85	first_superdiagonal[n+3] = std::numeric_limits<Real>::quiet_NaN();
	86
	87	std::vector<Real> second_superdiagonal(n+4, 1);
	88	second_superdiagonal[0] = 9;
	89	second_superdiagonal[1] = 8;
	90	second_superdiagonal[n+2] = std::numeric_limits<Real>::quiet_NaN();
	91	second_superdiagonal[n+3] = std::numeric_limits<Real>::quiet_NaN();
	92
	93	std::vector<Real> first_subdiagonal(n+4, 26);
	94	first_subdiagonal[0] = std::numeric_limits<Real>::quiet_NaN();
	95	first_subdiagonal[1] = 1;
	96	first_subdiagonal[n+2] = 33;
	97	first_subdiagonal[n+3] = 10;
	98
	99	std::vector<Real> second_subdiagonal(n+4, 1);
	100	second_subdiagonal[0] = std::numeric_limits<Real>::quiet_NaN();
	101	second_subdiagonal[1] = std::numeric_limits<Real>::quiet_NaN();
	102	second_subdiagonal[n+2] = 8;
	103	second_subdiagonal[n+3] = 9;
104
105
106	for (size_t i = 0; i < n+2; ++i) {
107	Real di = diagonal[i];
108	diagonal[i] = 1;
109	first_superdiagonal[i] /= di;
110	second_superdiagonal[i] /= di;
111	rhs[i] /= di;
112
113	// Eliminate first subdiagonal:
114	Real nfsub = -first_subdiagonal[i+1];
115	// Superfluous:
116	first_subdiagonal[i+1] /= nfsub;
117	// Not superfluous:
118	diagonal[i+1] /= nfsub;
119	first_superdiagonal[i+1] /= nfsub;
120	second_superdiagonal[i+1] /= nfsub;
121	rhs[i+1] /= nfsub;
122
123	diagonal[i+1] += first_superdiagonal[i];
124	first_superdiagonal[i+1] += second_superdiagonal[i];
125	rhs[i+1] += rhs[i];
126	// Superfluous, but clarifying:
127	first_subdiagonal[i+1] = 0;
128
129	// Eliminate second subdiagonal:
130	Real nssub = -second_subdiagonal[i+2];
131	first_subdiagonal[i+2] /= nssub;
132	diagonal[i+2] /= nssub;
133	first_superdiagonal[i+2] /= nssub;
134	second_superdiagonal[i+2] /= nssub;
135	rhs[i+2] /= nssub;
136
137	first_subdiagonal[i+2] += first_superdiagonal[i];
138	diagonal[i+2] += second_superdiagonal[i];
139	rhs[i+2] += rhs[i];
140	// Superfluous, but clarifying:
141	second_subdiagonal[i+2] = 0;
142	}
143
144	// Eliminate last subdiagonal:
145	Real dnp2 = diagonal[n+2];
146	diagonal[n+2] = 1;
147	first_superdiagonal[n+2] /= dnp2;
148	rhs[n+2] /= dnp2;
149	Real nfsubnp3 = -first_subdiagonal[n+3];
150	diagonal[n+3] /= nfsubnp3;
151	rhs[n+3] /= nfsubnp3;
152
153	diagonal[n+3] += first_superdiagonal[n+2];
154	rhs[n+3] += rhs[n+2];
155
156	m_alpha.resize(n + 4, std::numeric_limits<Real>::quiet_NaN());
157
158	m_alpha[n+3] = rhs[n+3]/diagonal[n+3];
159	m_alpha[n+2] = rhs[n+2] - first_superdiagonal[n+2]*m_alpha[n+3];
160	for (int64_t i = int64_t(n+1); i >= 0; --i) {
161	m_alpha[i] = rhs[i] - first_superdiagonal[i]m_alpha[i+1] - second_superdiagonal[i]m_alpha[i+2];
162	}
163
164	}
165
166	Real operator()(Real t) const {
167	using std::ceil;
168	using std::floor;
169	using boost::math::cardinal_b_spline;
170	// tf = t0 + (n-1)*h
171	// alpha.size() = n+4
172	if (t < m_t0 \|\| t > m_t0 + (m_alpha.size()-5)/m_inv_h) {
173	const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work.";
174	throw std::domain_error(err_msg);
175	}
176	Real x = (t-m_t0)*m_inv_h;
177	// Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5.
178	// TODO: Zero pad m_alpha so that only the domain check is necessary.
179	int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1)));
180	int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) );
181	Real s = 0;
182	for (int64_t j = j_min; j <= j_max; ++j) {
183	// TODO: Use Cox 1972 to generate all integer translates of B5 simultaneously.
184	s += m_alpha[j]*cardinal_b_spline<5, Real>(x - j + 2);
185	}
186	return s;
187	}
188
189	Real prime(Real t) const {
190	using std::ceil;
191	using std::floor;
192	using boost::math::cardinal_b_spline_prime;
193	if (t < m_t0 \|\| t > m_t0 + (m_alpha.size()-5)/m_inv_h) {
194	const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work.";
195	throw std::domain_error(err_msg);
196	}
197	Real x = (t-m_t0)*m_inv_h;
198	// Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5
199	int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1)));
200	int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) );
201	Real s = 0;
202	for (int64_t j = j_min; j <= j_max; ++j) {
203	s += m_alpha[j]*cardinal_b_spline_prime<5, Real>(x - j + 2);
204	}
205	return s*m_inv_h;
206
207	}
208
209	Real double_prime(Real t) const {
210	using std::ceil;
211	using std::floor;
212	using boost::math::cardinal_b_spline_double_prime;
213	if (t < m_t0 \|\| t > m_t0 + (m_alpha.size()-5)/m_inv_h) {
214	const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work.";
215	throw std::domain_error(err_msg);
216	}
217	Real x = (t-m_t0)*m_inv_h;
218	// Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5
219	int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1)));
220	int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) );
221	Real s = 0;
222	for (int64_t j = j_min; j <= j_max; ++j) {
223	s += m_alpha[j]*cardinal_b_spline_double_prime<5, Real>(x - j + 2);
224	}
225	return sm_inv_hm_inv_h;
226	}
227
228
229	Real t_max() const {
230	return m_t0 + (m_alpha.size()-5)/m_inv_h;
231	}
232
233	private:
234	std::vector<Real> m_alpha;
235	Real m_inv_h;
236	Real m_t0;
237	};
238
239	}}}}
240	#endif