[ceph.git] / ceph / src / boost / libs / math / test / cardinal_b_spline_test.cpp

/*
 * 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)
 */

#include "math_unit_test.hpp"
#include <numeric>
#include <utility>
#include <random>
#include <cmath>
#include <boost/math/special_functions/cardinal_b_spline.hpp>
#include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>
#ifdef BOOST_HAS_FLOAT128
#include <boost/multiprecision/float128.hpp>
using boost::multiprecision::float128;
#endif

using std::abs;
using boost::math::cardinal_b_spline;
using boost::math::cardinal_b_spline_prime;
using boost::math::forward_cardinal_b_spline;
using boost::math::cardinal_b_spline_double_prime;


template<class Real>
void test_box()
{
    Real t = cardinal_b_spline<0>(Real(1.1));
    Real expected = 0;
    CHECK_ULP_CLOSE(expected, t, 0);
    CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<0>(Real(1.1)), 0);

    t = cardinal_b_spline<0>(Real(-1.1));
    expected = 0;
    CHECK_ULP_CLOSE(expected, t, 0);

    Real h = Real(1)/Real(256);
    for (t = -Real(1)/Real(2)+h; t < Real(1)/Real(2); t += h)
    {
        expected = 1;
        CHECK_ULP_CLOSE(expected, cardinal_b_spline<0>(t), 0);
        expected = 0;
        CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<0>(Real(1.1)), 0);
    }

    for (t = h; t < 1; t += h)
    {
        expected = 1;
        CHECK_ULP_CLOSE(expected, forward_cardinal_b_spline<0>(t), 0);
    }
}

template<class Real>
void test_hat()
{
    Real t = cardinal_b_spline<1>(Real(2.1));
    Real expected = 0;
    CHECK_ULP_CLOSE(expected, t, 0);

    t = cardinal_b_spline<1>(Real(-2.1));
    expected = 0;
    CHECK_ULP_CLOSE(expected, t, 0);

    Real h = Real(1)/Real(256);
    for (t = -1; t <= 1; t += h)
    {
        expected = 1-abs(t);
        if(!CHECK_ULP_CLOSE(expected, cardinal_b_spline<1>(t), 0) )
        {
            std::cerr << "  Problem at t = " << t << "\n";
        }
        if (t == -Real(1)) {
            if (!CHECK_ULP_CLOSE(Real(1)/Real(2), cardinal_b_spline_prime<1>(t), 0)) {
                std::cout << "  Problem at t = " << t << "\n";
            }
        }
        else if (t == Real(1)) {
            CHECK_ULP_CLOSE(-Real(1)/Real(2), cardinal_b_spline_prime<1>(t), 0);
        }
        else if (t < 0) {
            CHECK_ULP_CLOSE(Real(1), cardinal_b_spline_prime<1>(t), 0);
        }
        else if (t == 0) {
            CHECK_ULP_CLOSE(Real(0), cardinal_b_spline_prime<1>(t), 0);
        }
        else if (t > 0) {
            CHECK_ULP_CLOSE(Real(-1), cardinal_b_spline_prime<1>(t), 0);
        }
    }

    for (t = 0; t < 2; t += h)
    {
        expected = 1 - abs(t-1);
        CHECK_ULP_CLOSE(expected, forward_cardinal_b_spline<1>(t), 0);
    }
}

template<class Real>
void test_quadratic()
{
    using std::abs;
    auto b2 = [](Real x) {
        Real absx = abs(x);
        if (absx >= 3/Real(2)) {
            return Real(0);
        }
        if (absx >= 1/Real(2)) {
            Real t = absx - 3/Real(2);
            return t*t/2;
        }
        Real t1 = absx - 1/Real(2);
        Real t2 = absx + 1/Real(2);
        return (2-t1*t1 -t2*t2)/2;
    };

    auto b2_prime = [&](Real x)->Real {
        Real absx = abs(x);
        Real signx  = 1;
        if (x < 0) {
            signx = -1;
        }
        if (absx >= 3/Real(2)) {
            return Real(0);
        }
        if (absx >= 1/Real(2)) {
            return (absx - 3/Real(2))*signx;
        }
        return -2*absx*signx;
    };


    Real h = 1/Real(256);
    for (Real t = -5; t <= 5; t += h) {
        Real expected = b2(t);
        CHECK_ULP_CLOSE(expected, cardinal_b_spline<2>(t), 0);
        expected = b2_prime(t);

        if (!CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<2>(t), 0))
        {
            std::cerr << "  Problem at t = " << t << "\n";
        }

    }
}

template<class Real>
void test_cubic()
{
    Real expected = Real(2)/Real(3);
    Real computed = cardinal_b_spline<3, Real>(0);
    CHECK_ULP_CLOSE(expected, computed, 0);

    expected = Real(1)/Real(6);
    computed = cardinal_b_spline<3, Real>(1);
    CHECK_ULP_CLOSE(expected, computed, 0);

    expected = Real(0);
    computed = cardinal_b_spline<3, Real>(2);
    CHECK_ULP_CLOSE(expected, computed, 0);

    Real h = 1/Real(256);
    for (Real t = -4; t <= 4; t += h) {
        expected = boost::math::detail::b3_spline_prime<Real>(t);
        computed = cardinal_b_spline_prime<3>(t);
        CHECK_ULP_CLOSE(expected, computed, 0);
        expected = boost::math::detail::b3_spline_double_prime<Real>(t);
        computed = cardinal_b_spline_double_prime<3>(t);
        if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
            std::cerr << "  Problem at t = " << t << "\n";
        }
    }
}

template<class Real>
void test_quintic()
{
  Real expected = Real(11)/Real(20);
  Real computed = cardinal_b_spline<5, Real>(0);
  CHECK_ULP_CLOSE(expected, computed, 0);

  expected = Real(13)/Real(60);
  computed = cardinal_b_spline<5, Real>(1);
  CHECK_ULP_CLOSE(expected, computed, 1);

  expected = Real(1)/Real(120);
  computed = cardinal_b_spline<5, Real>(2);
  CHECK_ULP_CLOSE(expected, computed, 0);

  expected = Real(0);
  computed = cardinal_b_spline<5, Real>(3);
  CHECK_ULP_CLOSE(expected, computed, 0);

}

template<unsigned n, typename Real>
void test_b_spline_derivatives()
{
    Real h = 1/Real(256);
    Real supp = (n+Real(1))/Real(2);
    for (Real t = -supp - 1; t <= supp+1; t+= h)
    {
        Real expected = cardinal_b_spline<n-1>(t+Real(1)/Real(2)) - cardinal_b_spline<n-1>(t - Real(1)/Real(2));
        Real computed = cardinal_b_spline_prime<n>(t);
        CHECK_MOLLIFIED_CLOSE(expected, computed, std::numeric_limits<Real>::epsilon());

        expected = cardinal_b_spline<n-2>(t+1) - 2*cardinal_b_spline<n-2>(t) + cardinal_b_spline<n-2>(t-1);
        computed = cardinal_b_spline_double_prime<n>(t);
        CHECK_MOLLIFIED_CLOSE(expected, computed, 2*std::numeric_limits<Real>::epsilon());
    }
}

template<unsigned n, typename Real>
void test_partition_of_unity()
{
  std::mt19937 gen(323723);
  Real supp = (n+1.0)/2.0;
  std::uniform_real_distribution<Real> dis(-supp, -supp+1);

  for(size_t i = 0; i < 500; ++i) {
    Real x = dis(gen);
    Real one = 0;
    while (x < supp) {
        one += cardinal_b_spline<n>(x);
        x += 1;
    }
    if(!CHECK_ULP_CLOSE(Real(1), one, n)) {
      std::cerr << "  Partition of unity failure at n = " << n << "\n";
    }
  }
}


int main()
{
    test_box<float>();
    test_box<double>();
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    test_box<long double>();
#endif

    test_hat<float>();
    test_hat<double>();
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    test_hat<long double>();
#endif

    test_quadratic<float>();
    test_quadratic<double>();
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    test_quadratic<long double>();
#endif

    test_cubic<float>();
    test_cubic<double>();
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    test_cubic<long double>();
#endif

    test_quintic<float>();
    test_quintic<double>();
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    test_quintic<long double>();
#endif

    test_partition_of_unity<1, double>();
    test_partition_of_unity<2, double>();
    test_partition_of_unity<3, double>();
    test_partition_of_unity<4, double>();
    test_partition_of_unity<5, double>();
    test_partition_of_unity<6, double>();

    test_b_spline_derivatives<3, double>();
    test_b_spline_derivatives<4, double>();
    test_b_spline_derivatives<5, double>();
    test_b_spline_derivatives<6, double>();
    test_b_spline_derivatives<7, double>();
    test_b_spline_derivatives<8, double>();
    test_b_spline_derivatives<9, double>();

#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    test_b_spline_derivatives<3, long double>();
    test_b_spline_derivatives<4, long double>();
    test_b_spline_derivatives<5, long double>();
    test_b_spline_derivatives<6, long double>();
    test_b_spline_derivatives<7, long double>();
    test_b_spline_derivatives<8, long double>();
    test_b_spline_derivatives<9, long double>();
#endif

#ifdef BOOST_HAS_FLOAT128
    test_box<float128>();
    test_hat<float128>();
    test_quadratic<float128>();
    test_cubic<float128>();
    test_quintic<float128>();
#endif

    return boost::math::test::report_errors();
}
Commit	Line	Data
92f5a8d4 TL	1	/*
	2	* Copyright Nick Thompson, 2019
	3	* Use, modification and distribution are subject to the
	4	* Boost Software License, Version 1.0. (See accompanying file
	5	* LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
	6	*/
	7
	8	#include "math_unit_test.hpp"
	9	#include <numeric>
	10	#include <utility>
	11	#include <random>
	12	#include <cmath>
	13	#include <boost/math/special_functions/cardinal_b_spline.hpp>
	14	#include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>
	15	#ifdef BOOST_HAS_FLOAT128
	16	#include <boost/multiprecision/float128.hpp>
	17	using boost::multiprecision::float128;
	18	#endif
	19
	20	using std::abs;
	21	using boost::math::cardinal_b_spline;
	22	using boost::math::cardinal_b_spline_prime;
	23	using boost::math::forward_cardinal_b_spline;
	24	using boost::math::cardinal_b_spline_double_prime;
	25
	26
	27	template<class Real>
	28	void test_box()
	29	{
	30	Real t = cardinal_b_spline<0>(Real(1.1));
	31	Real expected = 0;
	32	CHECK_ULP_CLOSE(expected, t, 0);
	33	CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<0>(Real(1.1)), 0);
	34
	35	t = cardinal_b_spline<0>(Real(-1.1));
	36	expected = 0;
	37	CHECK_ULP_CLOSE(expected, t, 0);
	38
	39	Real h = Real(1)/Real(256);
	40	for (t = -Real(1)/Real(2)+h; t < Real(1)/Real(2); t += h)
	41	{
	42	expected = 1;
	43	CHECK_ULP_CLOSE(expected, cardinal_b_spline<0>(t), 0);
	44	expected = 0;
	45	CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<0>(Real(1.1)), 0);
	46	}
	47
	48	for (t = h; t < 1; t += h)
	49	{
	50	expected = 1;
	51	CHECK_ULP_CLOSE(expected, forward_cardinal_b_spline<0>(t), 0);
	52	}
	53	}
	54
	55	template<class Real>
	56	void test_hat()
	57	{
	58	Real t = cardinal_b_spline<1>(Real(2.1));
	59	Real expected = 0;
	60	CHECK_ULP_CLOSE(expected, t, 0);
	61
	62	t = cardinal_b_spline<1>(Real(-2.1));
	63	expected = 0;
	64	CHECK_ULP_CLOSE(expected, t, 0);
65
66	Real h = Real(1)/Real(256);
67	for (t = -1; t <= 1; t += h)
68	{
69	expected = 1-abs(t);
70	if(!CHECK_ULP_CLOSE(expected, cardinal_b_spline<1>(t), 0) )
71	{
72	std::cerr << " Problem at t = " << t << "\n";
73	}
74	if (t == -Real(1)) {
75	if (!CHECK_ULP_CLOSE(Real(1)/Real(2), cardinal_b_spline_prime<1>(t), 0)) {
76	std::cout << " Problem at t = " << t << "\n";
77	}
78	}
79	else if (t == Real(1)) {
80	CHECK_ULP_CLOSE(-Real(1)/Real(2), cardinal_b_spline_prime<1>(t), 0);
81	}
82	else if (t < 0) {
83	CHECK_ULP_CLOSE(Real(1), cardinal_b_spline_prime<1>(t), 0);
84	}
85	else if (t == 0) {
86	CHECK_ULP_CLOSE(Real(0), cardinal_b_spline_prime<1>(t), 0);
87	}
88	else if (t > 0) {
89	CHECK_ULP_CLOSE(Real(-1), cardinal_b_spline_prime<1>(t), 0);
90	}
91	}
92
93	for (t = 0; t < 2; t += h)
94	{
95	expected = 1 - abs(t-1);
96	CHECK_ULP_CLOSE(expected, forward_cardinal_b_spline<1>(t), 0);
97	}
98	}
99
100	template<class Real>
101	void test_quadratic()
102	{
103	using std::abs;
104	auto b2 = [](Real x) {
105	Real absx = abs(x);
106	if (absx >= 3/Real(2)) {
107	return Real(0);
108	}
109	if (absx >= 1/Real(2)) {
110	Real t = absx - 3/Real(2);
111	return t*t/2;
112	}
113	Real t1 = absx - 1/Real(2);
114	Real t2 = absx + 1/Real(2);
115	return (2-t1t1 -t2t2)/2;
116	};
117
118	auto b2_prime = [&](Real x)->Real {
119	Real absx = abs(x);
120	Real signx = 1;
121	if (x < 0) {
122	signx = -1;
123	}
124	if (absx >= 3/Real(2)) {
125	return Real(0);
126	}
127	if (absx >= 1/Real(2)) {
128	return (absx - 3/Real(2))*signx;
129	}
130	return -2absxsignx;
131	};
132
133
134	Real h = 1/Real(256);
135	for (Real t = -5; t <= 5; t += h) {
136	Real expected = b2(t);
137	CHECK_ULP_CLOSE(expected, cardinal_b_spline<2>(t), 0);
138	expected = b2_prime(t);
139
140	if (!CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<2>(t), 0))
141	{
142	std::cerr << " Problem at t = " << t << "\n";
143	}
144
145	}
146	}
147
148	template<class Real>
149	void test_cubic()
150	{
151	Real expected = Real(2)/Real(3);
152	Real computed = cardinal_b_spline<3, Real>(0);
153	CHECK_ULP_CLOSE(expected, computed, 0);
154
155	expected = Real(1)/Real(6);
156	computed = cardinal_b_spline<3, Real>(1);
157	CHECK_ULP_CLOSE(expected, computed, 0);
158
159	expected = Real(0);
160	computed = cardinal_b_spline<3, Real>(2);
161	CHECK_ULP_CLOSE(expected, computed, 0);
162
163	Real h = 1/Real(256);
164	for (Real t = -4; t <= 4; t += h) {
165	expected = boost::math::detail::b3_spline_prime<Real>(t);
166	computed = cardinal_b_spline_prime<3>(t);
167	CHECK_ULP_CLOSE(expected, computed, 0);
168	expected = boost::math::detail::b3_spline_double_prime<Real>(t);
169	computed = cardinal_b_spline_double_prime<3>(t);
170	if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
171	std::cerr << " Problem at t = " << t << "\n";
172	}
173	}
174	}
175
176	template<class Real>
177	void test_quintic()
178	{
179	Real expected = Real(11)/Real(20);
180	Real computed = cardinal_b_spline<5, Real>(0);
181	CHECK_ULP_CLOSE(expected, computed, 0);
182
183	expected = Real(13)/Real(60);
184	computed = cardinal_b_spline<5, Real>(1);
185	CHECK_ULP_CLOSE(expected, computed, 1);
186
187	expected = Real(1)/Real(120);
188	computed = cardinal_b_spline<5, Real>(2);
189	CHECK_ULP_CLOSE(expected, computed, 0);
190
191	expected = Real(0);
192	computed = cardinal_b_spline<5, Real>(3);
193	CHECK_ULP_CLOSE(expected, computed, 0);
194
195	}
196
197	template<unsigned n, typename Real>
198	void test_b_spline_derivatives()
199	{
200	Real h = 1/Real(256);
201	Real supp = (n+Real(1))/Real(2);
202	for (Real t = -supp - 1; t <= supp+1; t+= h)
203	{
204	Real expected = cardinal_b_spline<n-1>(t+Real(1)/Real(2)) - cardinal_b_spline<n-1>(t - Real(1)/Real(2));
205	Real computed = cardinal_b_spline_prime<n>(t);
206	CHECK_MOLLIFIED_CLOSE(expected, computed, std::numeric_limits<Real>::epsilon());
207
208	expected = cardinal_b_spline<n-2>(t+1) - 2*cardinal_b_spline<n-2>(t) + cardinal_b_spline<n-2>(t-1);
209	computed = cardinal_b_spline_double_prime<n>(t);
210	CHECK_MOLLIFIED_CLOSE(expected, computed, 2*std::numeric_limits<Real>::epsilon());
211	}
212	}
213
214	template<unsigned n, typename Real>
215	void test_partition_of_unity()
216	{
217	std::mt19937 gen(323723);
218	Real supp = (n+1.0)/2.0;
219	std::uniform_real_distribution<Real> dis(-supp, -supp+1);
220
221	for(size_t i = 0; i < 500; ++i) {
222	Real x = dis(gen);
223	Real one = 0;
224	while (x < supp) {
225	one += cardinal_b_spline<n>(x);
226	x += 1;
227	}
228	if(!CHECK_ULP_CLOSE(Real(1), one, n)) {
229	std::cerr << " Partition of unity failure at n = " << n << "\n";
230	}
231	}
232	}
233
234
235	int main()
236	{
237	test_box<float>();
238	test_box<double>();
1e59de90	239	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
92f5a8d4	240	test_box<long double>();
1e59de90	241	#endif
92f5a8d4 TL	242
	243	test_hat<float>();
	244	test_hat<double>();
1e59de90	245	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
92f5a8d4	246	test_hat<long double>();
1e59de90	247	#endif
92f5a8d4 TL	248
	249	test_quadratic<float>();
	250	test_quadratic<double>();
1e59de90	251	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
92f5a8d4	252	test_quadratic<long double>();
1e59de90	253	#endif
92f5a8d4 TL	254
	255	test_cubic<float>();
	256	test_cubic<double>();
1e59de90	257	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
92f5a8d4	258	test_cubic<long double>();
1e59de90	259	#endif
92f5a8d4 TL	260
	261	test_quintic<float>();
	262	test_quintic<double>();
1e59de90	263	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
92f5a8d4	264	test_quintic<long double>();
1e59de90	265	#endif
92f5a8d4 TL	266
	267	test_partition_of_unity<1, double>();
	268	test_partition_of_unity<2, double>();
	269	test_partition_of_unity<3, double>();
	270	test_partition_of_unity<4, double>();
	271	test_partition_of_unity<5, double>();
	272	test_partition_of_unity<6, double>();
	273
	274	test_b_spline_derivatives<3, double>();
	275	test_b_spline_derivatives<4, double>();
	276	test_b_spline_derivatives<5, double>();
	277	test_b_spline_derivatives<6, double>();
	278	test_b_spline_derivatives<7, double>();
	279	test_b_spline_derivatives<8, double>();
	280	test_b_spline_derivatives<9, double>();
	281
1e59de90	282	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
92f5a8d4 TL	283	test_b_spline_derivatives<3, long double>();
	284	test_b_spline_derivatives<4, long double>();
	285	test_b_spline_derivatives<5, long double>();
	286	test_b_spline_derivatives<6, long double>();
	287	test_b_spline_derivatives<7, long double>();
	288	test_b_spline_derivatives<8, long double>();
	289	test_b_spline_derivatives<9, long double>();
1e59de90	290	#endif
92f5a8d4 TL	291
	292	#ifdef BOOST_HAS_FLOAT128
	293	test_box<float128>();
	294	test_hat<float128>();
	295	test_quadratic<float128>();
	296	test_cubic<float128>();
	297	test_quintic<float128>();
	298	#endif
	299
	300	return boost::math::test::report_errors();
	301	}