libs/math/test/lanczos_smoothing_test.cpp

   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 #define BOOST_TEST_MODULE lanczos_smoothing_test
   8
   9 #include <random>
  10 #include <array>
  11 #include <boost/range.hpp>
  12 #include <boost/numeric/ublas/vector.hpp>
  13 #include <boost/math/constants/constants.hpp>
  14 #include <boost/test/included/unit_test.hpp>
  15 #include <boost/test/tools/floating_point_comparison.hpp>
  16 #include <boost/math/differentiation/lanczos_smoothing.hpp>
  17 #include <boost/multiprecision/cpp_bin_float.hpp>
  18 #include <boost/math/special_functions/next.hpp> // for float_distance
  19 #include <boost/math/tools/condition_numbers.hpp>
  20
  21 using std::abs;
  22 using std::pow;
  23 using std::sqrt;
  24 using std::sin;
  25 using boost::math::constants::two_pi;
  26 using boost::multiprecision::cpp_bin_float_50;
  27 using boost::multiprecision::cpp_bin_float_100;
  28 using boost::math::differentiation::discrete_lanczos_derivative;
  29 using boost::math::differentiation::detail::discrete_legendre;
  30 using boost::math::differentiation::detail::interior_velocity_filter;
  31 using boost::math::differentiation::detail::boundary_velocity_filter;
  32 using boost::math::tools::summation_condition_number;
  33
  34 template<class Real>
  35 void test_dlp_norms()
  36 {
  37     std::cout << "Testing Discrete Legendre Polynomial norms on type " << typeid(Real).name() << "\n";
  38     Real tol = std::numeric_limits<Real>::epsilon();
  39     auto dlp = discrete_legendre<Real>(1, Real(0));
  40     BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(0), 3, tol);
  41     BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(1), 2, tol);
  42     dlp = discrete_legendre<Real>(2, Real(0));
  43     BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(0), Real(5)/Real(2), tol);
  44     BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(1), Real(5)/Real(4), tol);
  45     BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(2), Real(3*3*7)/Real(pow(2,6)), 2*tol);
  46     dlp = discrete_legendre<Real>(200, Real(0));
  47     for(size_t r = 0; r < 10; ++r)
  48     {
  49         Real calc = dlp.norm_sq(r);
  50         Real expected = Real(2)/Real(2*r+1);
  51         // As long as r << n, ||q_r||^2 -> 2/(2r+1) as n->infty
  52         BOOST_CHECK_CLOSE_FRACTION(calc, expected, 0.05);
  53     }
  54
  55 }
  56
  57 template<class Real>
  58 void test_dlp_evaluation()
  59 {
  60     std::cout << "Testing evaluation of Discrete Legendre polynomials on type " << typeid(Real).name() << "\n";
  61     Real tol = std::numeric_limits<Real>::epsilon();
  62     size_t n = 25;
  63     Real x = 0.72;
  64     auto dlp = discrete_legendre<Real>(n, x);
  65     Real q0 = dlp(x, 0);
  66     BOOST_TEST(q0 == 1);
  67     Real q1 = dlp(x, 1);
  68     BOOST_TEST(q1 == x);
  69     Real q2 = dlp(x, 2);
  70     int N = 2*n+1;
  71     Real expected = 0.5*(3*x*x - Real(N*N - 1)/Real(4*n*n));
  72     BOOST_CHECK_CLOSE_FRACTION(q2, expected, tol);
  73     Real q3 = dlp(x, 3);
  74     expected = (x/3)*(5*expected - (Real(N*N - 4))/(2*n*n));
  75     BOOST_CHECK_CLOSE_FRACTION(q3, expected, 2*tol);
  76
  77     // q_r(x) is even for even r, and odd for odd r:
  78     for (size_t n = 8; n < 22; ++n)
  79     {
  80         dlp = discrete_legendre<Real>(n, x);
  81         for(size_t r = 2; r <= n; ++r)
  82         {
  83             if (r & 1)
  84             {
  85                 Real q1 = dlp(x, r);
  86                 Real q2 = -dlp(-x, r);
  87                 BOOST_CHECK_CLOSE_FRACTION(q1, q2, tol);
  88             }
  89             else
  90             {
  91                 Real q1 = dlp(x, r);
  92                 Real q2 = dlp(-x, r);
  93                 BOOST_CHECK_CLOSE_FRACTION(q1, q2, tol);
  94             }
  95
  96             Real l2_sq = 0;
  97             for (int j = -(int)n; j <= (int) n; ++j)
  98             {
  99                 Real y = Real(j)/Real(n);
 100                 Real term = dlp(y, r);
 101                 l2_sq += term*term;
 102             }
 103             l2_sq /= n;
 104             Real l2_sq_expected = dlp.norm_sq(r);
 105             BOOST_CHECK_CLOSE_FRACTION(l2_sq, l2_sq_expected, 20*tol);
 106         }
 107     }
 108 }
 109
 110 template<class Real>
 111 void test_dlp_next()
 112 {
 113     std::cout << "Testing Discrete Legendre polynomial 'next' function on type " << typeid(Real).name() << "\n";
 114     Real tol = std::numeric_limits<Real>::epsilon();
 115
 116     for(size_t n = 2; n < 20; ++n)
 117     {
 118         for(Real x = -1; x <= 1; x += 0.1)
 119         {
 120             auto dlp = discrete_legendre<Real>(n, x);
 121             for (size_t k = 2; k < n; ++k)
 122             {
 123                 BOOST_CHECK_CLOSE(dlp.next(), dlp(x, k), tol);
 124             }
 125
 126             dlp = discrete_legendre<Real>(n, x);
 127             for (size_t k = 2; k < n; ++k)
 128             {
 129                 BOOST_CHECK_CLOSE(dlp.next_prime(), dlp.prime(x, k), tol);
 130             }
 131         }
 132     }
 133 }
 134
 135
 136 template<class Real>
 137 void test_dlp_derivatives()
 138 {
 139     std::cout << "Testing Discrete Legendre polynomial derivatives on type " << typeid(Real).name() << "\n";
 140     Real tol = 10*std::numeric_limits<Real>::epsilon();
 141     int n = 25;
 142     Real x = 0.72;
 143     auto dlp = discrete_legendre<Real>(n, x);
 144     Real q0p = dlp.prime(x, 0);
 145     BOOST_TEST(q0p == 0);
 146     Real q1p = dlp.prime(x, 1);
 147     BOOST_TEST(q1p == 1);
 148     Real q2p = dlp.prime(x, 2);
 149     Real expected = 3*x;
 150     BOOST_CHECK_CLOSE_FRACTION(q2p, expected, tol);
 151 }
 152
 153 template<class Real>
 154 void test_dlp_second_derivative()
 155 {
 156     std::cout << "Testing Discrete Legendre polynomial derivatives on type " << typeid(Real).name() << "\n";
 157     int n = 25;
 158     Real x = Real(1)/Real(3);
 159     auto dlp = discrete_legendre<Real>(n, x);
 160     Real q2pp = dlp.next_dbl_prime();
 161     BOOST_TEST(q2pp == 3);
 162 }
 163
 164
 165 template<class Real>
 166 void test_interior_velocity_filter()
 167 {
 168     using boost::math::constants::half;
 169     std::cout << "Testing interior filter on type " << typeid(Real).name() << "\n";
 170     Real tol = std::numeric_limits<Real>::epsilon();
 171     for(int n = 1; n < 10; ++n)
 172     {
 173         for (int p = 1; p < n; p += 2)
 174         {
 175             auto f = interior_velocity_filter<Real>(n,p);
 176             // Since we only store half the filter coefficients,
 177             // we need to reindex the moment sums:
 178             auto cond = summation_condition_number<Real>(0);
 179             for (size_t j = 0; j < f.size(); ++j)
 180             {
 181                 cond += j*f[j];
 182             }
 183             BOOST_CHECK_CLOSE_FRACTION(cond.sum(), half<Real>(), 2*cond()*tol);
 184
 185             for (int l = 3; l <= p; l += 2)
 186             {
 187                 cond = summation_condition_number<Real>(0);
 188                 for (size_t j = 0; j < f.size() - 1; ++j)
 189                 {
 190                     cond += pow(Real(j), l)*f[j];
 191                 }
 192                 Real expected = -pow(Real(f.size() - 1), l)*f[f.size()-1];
 193                 BOOST_CHECK_CLOSE_FRACTION(expected, cond.sum(), 7*cond()*tol);
 194             }
 195             //std::cout << "(n,p) = (" << n  << "," << p << ") = {";
 196             //for (auto & x : f)
 197             //{
 198             //    std::cout << x << ", ";
 199             //}
 200             //std::cout << "}\n";
 201         }
 202     }
 203 }
 204
 205 template<class Real>
 206 void test_interior_lanczos()
 207 {
 208     std::cout << "Testing interior Lanczos on type " << typeid(Real).name() << "\n";
 209     Real tol = std::numeric_limits<Real>::epsilon();
 210     std::vector<Real> v(500);
 211     std::fill(v.begin(), v.end(), 7);
 212
 213     for (size_t n = 1; n < 10; ++n)
 214     {
 215         for (size_t p = 2; p < 2*n; p += 2)
 216         {
 217             auto dld = discrete_lanczos_derivative(Real(0.1), n, p);
 218             for (size_t m = n; m < v.size() - n; ++m)
 219             {
 220                 Real dvdt = dld(v, m);
 221                 BOOST_CHECK_SMALL(dvdt, tol);
 222             }
 223             auto dvdt = dld(v);
 224             for (size_t m = n; m < v.size() - n; ++m)
 225             {
 226                 BOOST_CHECK_SMALL(dvdt[m], tol);
 227             }
 228         }
 229     }
 230
 231
 232     for(size_t i = 0; i < v.size(); ++i)
 233     {
 234         v[i] = 7*i+8;
 235     }
 236
 237     for (size_t n = 1; n < 10; ++n)
 238     {
 239         for (size_t p = 2; p < 2*n; p += 2)
 240         {
 241             auto dld = discrete_lanczos_derivative(Real(1), n, p);
 242             for (size_t m = n; m < v.size() - n; ++m)
 243             {
 244                 Real dvdt = dld(v, m);
 245                 BOOST_CHECK_CLOSE_FRACTION(dvdt, 7, 2000*tol);
 246             }
 247             auto dvdt = dld(v);
 248             for (size_t m = n; m < v.size() - n; ++m)
 249             {
 250                 BOOST_CHECK_CLOSE_FRACTION(dvdt[m], 7, 2000*tol);
 251             }
 252         }
 253     }
 254
 255     //std::random_device rd{};
 256     //auto seed = rd();
 257     //std::cout << "Seed = " << seed << "\n";
 258     std::mt19937 gen(4172378669);
 259     std::normal_distribution<> dis{0, 0.01};
 260     for (size_t i = 0; i < v.size(); ++i)
 261     {
 262         v[i] = 7*i+8 + dis(gen);
 263     }
 264
 265     for (size_t n = 1; n < 10; ++n)
 266     {
 267         for (size_t p = 2; p < 2*n; p += 2)
 268         {
 269             auto dld = discrete_lanczos_derivative(Real(1), n, p);
 270             for (size_t m = n; m < v.size() - n; ++m)
 271             {
 272                 BOOST_CHECK_CLOSE_FRACTION(dld(v, m), Real(7), Real(0.0042));
 273             }
 274         }
 275     }
 276
 277
 278     for (size_t i = 0; i < v.size(); ++i)
 279     {
 280         v[i] = 15*i*i + 7*i+8 + dis(gen);
 281     }
 282
 283     for (size_t n = 1; n < 10; ++n)
 284     {
 285         for (size_t p = 2; p < 2*n; p += 2)
 286         {
 287             auto dld = discrete_lanczos_derivative(Real(1), n, p);
 288             for (size_t m = n; m < v.size() - n; ++m)
 289             {
 290                 BOOST_CHECK_CLOSE_FRACTION(dld(v,m), Real(30*m + 7), Real(0.00008));
 291             }
 292         }
 293     }
 294
 295     std::normal_distribution<> dis1{0, 0.0001};
 296     Real omega = Real(1)/Real(16);
 297     for (size_t i = 0; i < v.size(); ++i)
 298     {
 299         v[i] = sin(i*omega) + dis1(gen);
 300     }
 301
 302     for (size_t n = 10; n < 20; ++n)
 303     {
 304         for (size_t p = 3; p < 100 && p < n/2; p += 2)
 305         {
 306             auto dld = discrete_lanczos_derivative(Real(1), n, p);
 307
 308             for (size_t m = n; m < v.size() - n && m < n + 10; ++m)
 309             {
 310                 BOOST_CHECK_CLOSE_FRACTION(dld(v,m), omega*cos(omega*m), Real(0.03));
 311             }
 312         }
 313     }
 314 }
 315
 316 template<class Real>
 317 void test_boundary_velocity_filters()
 318 {
 319     std::cout << "Testing boundary filters on type " << typeid(Real).name() << "\n";
 320     Real tol = std::numeric_limits<Real>::epsilon();
 321     for(int n = 1; n < 5; ++n)
 322     {
 323         for (int p = 1; p < 2*n+1; ++p)
 324         {
 325             for (int s = -n; s <= n; ++s)
 326             {
 327                 auto f = boundary_velocity_filter<Real>(n, p, s);
 328                 // Sum is zero:
 329                 auto cond = summation_condition_number<Real>(0);
 330                 for (size_t i = 0; i < f.size() - 1; ++i)
 331                 {
 332                     cond += f[i];
 333                 }
 334
 335                 BOOST_CHECK_CLOSE_FRACTION(cond.sum(), -f[f.size()-1], 6*cond()*tol);
 336
 337                 cond = summation_condition_number<Real>(0);
 338                 for (size_t k = 0; k < f.size(); ++k)
 339                 {
 340                     Real j = Real(k) - Real(n);
 341                     // note the shifted index here:
 342                     cond += (j-s)*f[k];
 343                 }
 344                 BOOST_CHECK_CLOSE_FRACTION(cond.sum(), 1, 6*cond()*tol);
 345
 346
 347                 for (int l = 2; l <= p; ++l)
 348                 {
 349                     cond = summation_condition_number<Real>(0);
 350                     for (size_t k = 0; k < f.size() - 1; ++k)
 351                     {
 352                         Real j = Real(k) - Real(n);
 353                         // The condition number of this sum is infinite!
 354                         // No need to get to worked up about the tolerance.
 355                         cond += pow(j-s, l)*f[k];
 356                     }
 357
 358                     Real expected = -pow(Real(f.size()-1) - Real(n) - Real(s), l)*f[f.size()-1];
 359                     if (expected == 0)
 360                     {
 361                         BOOST_CHECK_SMALL(cond.sum(), cond()*tol);
 362                     }
 363                     else
 364                     {
 365                         BOOST_CHECK_CLOSE_FRACTION(expected, cond.sum(), 200*cond()*tol);
 366                     }
 367                 }
 368
 369                 //std::cout << "(n,p,s) = ("<< n << ", " << p << "," << s << ") = {";
 370                 //for (auto & x : f)
 371                 //{
 372                 //    std::cout << x << ", ";
 373                 //}
 374                 //std::cout << "}\n";*/
 375             }
 376         }
 377     }
 378 }
 379
 380 template<class Real>
 381 void test_boundary_lanczos()
 382 {
 383     std::cout << "Testing Lanczos boundary on type " << typeid(Real).name() << "\n";
 384     Real tol = std::numeric_limits<Real>::epsilon();
 385     std::vector<Real> v(500, 7);
 386
 387     for (size_t n = 1; n < 10; ++n)
 388     {
 389         for (size_t p = 2; p < 2*n; ++p)
 390         {
 391             auto lsd = discrete_lanczos_derivative(Real(0.0125), n, p);
 392             for (size_t m = 0; m < n; ++m)
 393             {
 394                 Real dvdt = lsd(v,m);
 395                 BOOST_CHECK_SMALL(dvdt, 4*sqrt(tol));
 396             }
 397             for (size_t m = v.size() - n; m < v.size(); ++m)
 398             {
 399                 Real dvdt = lsd(v,m);
 400                 BOOST_CHECK_SMALL(dvdt, 4*sqrt(tol));
 401             }
 402         }
 403     }
 404
 405     for(size_t i = 0; i < v.size(); ++i)
 406     {
 407         v[i] = 7*i+8;
 408     }
 409
 410     for (size_t n = 3; n < 10; ++n)
 411     {
 412         for (size_t p = 2; p < 2*n; ++p)
 413         {
 414             auto lsd = discrete_lanczos_derivative(Real(1), n, p);
 415             for (size_t m = 0; m < n; ++m)
 416             {
 417                 Real dvdt = lsd(v,m);
 418                 BOOST_CHECK_CLOSE_FRACTION(dvdt, 7, sqrt(tol));
 419             }
 420
 421             for (size_t m = v.size() - n; m < v.size(); ++m)
 422             {
 423                 Real dvdt = lsd(v,m);
 424                 BOOST_CHECK_CLOSE_FRACTION(dvdt, 7, 4*sqrt(tol));
 425             }
 426         }
 427     }
 428
 429     for (size_t i = 0; i < v.size(); ++i)
 430     {
 431         v[i] = 15*i*i + 7*i+8;
 432     }
 433
 434     for (size_t n = 1; n < 10; ++n)
 435     {
 436         for (size_t p = 2; p < 2*n; ++p)
 437         {
 438             auto lsd = discrete_lanczos_derivative(Real(1), n, p);
 439             for (size_t m = 0; m < v.size(); ++m)
 440             {
 441                 BOOST_CHECK_CLOSE_FRACTION(lsd(v,m), 30*m+7, 30*sqrt(tol));
 442             }
 443         }
 444     }
 445
 446     // Demonstrate that the boundary filters are also denoising:
 447     //std::random_device rd{};
 448     //auto seed = rd();
 449     //std::cout << "seed = " << seed << "\n";
 450     std::mt19937 gen(311354333);
 451     std::normal_distribution<> dis{0, 0.01};
 452     for (size_t i = 0; i < v.size(); ++i)
 453     {
 454         v[i] += dis(gen);
 455     }
 456
 457     for (size_t n = 1; n < 10; ++n)
 458     {
 459         for (size_t p = 2; p < n; ++p)
 460         {
 461             auto lsd = discrete_lanczos_derivative(Real(1), n, p);
 462             for (size_t m = 0; m < v.size(); ++m)
 463             {
 464                 BOOST_CHECK_CLOSE_FRACTION(lsd(v,m), 30*m+7, 0.005);
 465             }
 466             auto dvdt = lsd(v);
 467             for (size_t m = 0; m < v.size(); ++m)
 468             {
 469                 BOOST_CHECK_CLOSE_FRACTION(dvdt[m], 30*m+7, 0.005);
 470             }
 471         }
 472     }
 473 }
 474
 475 template<class Real>
 476 void test_acceleration_filters()
 477 {
 478     Real eps = std::numeric_limits<Real>::epsilon();
 479     for (size_t n = 1; n < 5; ++n)
 480     {
 481         for(size_t p = 3; p <= 2*n; ++p)
 482         {
 483             for(int64_t s = -int64_t(n); s <= 0; ++s)
 484             {
 485                 auto g = boost::math::differentiation::detail::acceleration_filter<long double>(n,p,s);
 486
 487                 std::vector<Real> f(g.size());
 488                 for (size_t i = 0; i < g.size(); ++i)
 489                 {
 490                     f[i] = static_cast<Real>(g[i]);
 491                 }
 492
 493                 auto cond = summation_condition_number<Real>(0);
 494
 495                 for (size_t i = 0; i < f.size() - 1; ++i)
 496                 {
 497                     cond += f[i];
 498                 }
 499                 BOOST_CHECK_CLOSE_FRACTION(cond.sum(), -f[f.size()-1], 10*cond()*eps);
 500
 501
 502                 cond = summation_condition_number<Real>(0);
 503                 for (size_t k = 0; k < f.size() -1; ++k)
 504                 {
 505                     Real j = Real(k) - Real(n);
 506                     cond += (j-s)*f[k];
 507                 }
 508                 Real expected = -(Real(f.size()-1)- Real(n) - s)*f[f.size()-1];
 509                 BOOST_CHECK_CLOSE_FRACTION(cond.sum(), expected, 10*cond()*eps);
 510
 511                 cond = summation_condition_number<Real>(0);
 512                 for (size_t k = 0; k < f.size(); ++k)
 513                 {
 514                     Real j = Real(k) - Real(n);
 515                     cond += (j-s)*(j-s)*f[k];
 516                 }
 517                 BOOST_CHECK_CLOSE_FRACTION(cond.sum(), 2, 100*cond()*eps);
 518                 // See unlabelled equation in McDevitt, 2012, just after equation 26:
 519                 // It appears that there is an off-by-one error in that equation, since p + 1 moments don't vanish, only p.
 520                 // This test is itself suspect; the condition number of the moment sum is infinite.
 521                 // So the *slightest* error in the filter gets amplified by the test; in terms of the
 522                 // behavior of the actual filter, it's not a big deal.
 523                 for (size_t l = 3; l <= p; ++l)
 524                 {
 525                     cond = summation_condition_number<Real>(0);
 526                     for (size_t k = 0; k < f.size() - 1; ++k)
 527                     {
 528                         Real j = Real(k) - Real(n);
 529                         cond += pow((j-s), l)*f[k];
 530                     }
 531                     Real expected = -pow(Real(f.size()- 1 - n -s), l)*f[f.size()-1];
 532                     BOOST_CHECK_CLOSE_FRACTION(cond.sum(), expected, 1000*cond()*eps);
 533                 }
 534             }
 535         }
 536     }
 537 }
 538
 539 template<class Real>
 540 void test_lanczos_acceleration()
 541 {
 542     Real eps = std::numeric_limits<Real>::epsilon();
 543     std::vector<Real> v(100, 7);
 544     auto lanczos = discrete_lanczos_derivative<Real, 2>(Real(1), 4, 3);
 545     for (size_t i = 0; i < v.size(); ++i)
 546     {
 547         BOOST_CHECK_SMALL(lanczos(v, i), eps);
 548     }
 549
 550     for(size_t i = 0; i < v.size(); ++i)
 551     {
 552         v[i] = 7*i + 6;
 553     }
 554     for (size_t i = 0; i < v.size(); ++i)
 555     {
 556         BOOST_CHECK_SMALL(lanczos(v,i), 200*eps);
 557     }
 558
 559     for(size_t i = 0; i < v.size(); ++i)
 560     {
 561         v[i] = 7*i*i + 9*i + 6;
 562     }
 563     for (size_t i = 0; i < v.size(); ++i)
 564     {
 565         BOOST_CHECK_CLOSE_FRACTION(lanczos(v, i), 14, 1500*eps);
 566     }
 567
 568     // Now add noise, and kick up the smoothing of the Lanzcos derivative (increase n):
 569     //std::random_device rd{};
 570     //auto seed = rd();
 571     //std::cout << "seed = " << seed << "\n";
 572     size_t seed = 2507134629;
 573     std::mt19937 gen(seed);
 574     Real std_dev = 0.1;
 575     std::normal_distribution<Real> dis{0, std_dev};
 576     for (size_t i = 0; i < v.size(); ++i)
 577     {
 578         v[i] += dis(gen);
 579     }
 580     lanczos = discrete_lanczos_derivative<Real, 2>(Real(1), 18, 3);
 581     auto w = lanczos(v);
 582     for (size_t i = 0; i < v.size(); ++i)
 583     {
 584         BOOST_CHECK_CLOSE_FRACTION(w[i], 14, std_dev/200);
 585     }
 586 }
 587
 588 template<class Real>
 589 void test_rescaling()
 590 {
 591     std::cout << "Test rescaling on type " << typeid(Real).name() << "\n";
 592     Real tol = std::numeric_limits<Real>::epsilon();
 593     std::vector<Real> v(500);
 594     for(size_t i = 0; i < v.size(); ++i)
 595     {
 596         v[i] = 7*i*i + 9*i + 6;
 597     }
 598     std::vector<Real> dvdt1(500);
 599     std::vector<Real> dvdt2(500);
 600     auto lanczos1 = discrete_lanczos_derivative(Real(1));
 601     auto lanczos2 = discrete_lanczos_derivative(Real(2));
 602
 603     lanczos1(v, dvdt1);
 604     lanczos2(v, dvdt2);
 605
 606     for(size_t i = 0; i < v.size(); ++i)
 607     {
 608         BOOST_CHECK_CLOSE_FRACTION(dvdt1[i], 2*dvdt2[i], tol);
 609     }
 610
 611     auto lanczos3 = discrete_lanczos_derivative<Real, 2>(Real(1));
 612     auto lanczos4 = discrete_lanczos_derivative<Real, 2>(Real(2));
 613
 614
 615     std::vector<Real> dv2dt21(500);
 616     std::vector<Real> dv2dt22(500);
 617
 618     for(size_t i = 0; i < v.size(); ++i)
 619     {
 620         BOOST_CHECK_CLOSE_FRACTION(dv2dt21[i], 4*dv2dt22[i], tol);
 621     }
 622 }
 623
 624 template<class Real>
 625 void test_data_representations()
 626 {
 627     std::cout << "Test rescaling on type " << typeid(Real).name() << "\n";
 628     Real tol = 150*std::numeric_limits<Real>::epsilon();
 629     std::array<Real, 500> v;
 630     for(size_t i = 0; i < v.size(); ++i)
 631     {
 632         v[i] = 9*i + 6;
 633     }
 634     std::array<Real, 500> dvdt;
 635     auto lanczos = discrete_lanczos_derivative(Real(1));
 636
 637     lanczos(v, dvdt);
 638
 639     for(size_t i = 0; i < v.size(); ++i)
 640     {
 641         BOOST_CHECK_CLOSE_FRACTION(dvdt[i], 9, tol);
 642     }
 643
 644     boost::numeric::ublas::vector<Real> w(500);
 645     boost::numeric::ublas::vector<Real> dwdt(500);
 646     for(size_t i = 0; i < w.size(); ++i)
 647     {
 648         w[i] = 9*i + 6;
 649     }
 650
 651     lanczos(w, dwdt);
 652
 653     for(size_t i = 0; i < v.size(); ++i)
 654     {
 655         BOOST_CHECK_CLOSE_FRACTION(dwdt[i], 9, tol);
 656     }
 657
 658     auto v1 = boost::make_iterator_range(v.begin(), v.end());
 659     auto v2 = boost::make_iterator_range(dvdt.begin(), dvdt.end());
 660     lanczos(v1, v2);
 661
 662     for(size_t i = 0; i < v2.size(); ++i)
 663     {
 664         BOOST_CHECK_CLOSE_FRACTION(v2[i], 9, tol);
 665     }
 666
 667     auto lanczos2 = discrete_lanczos_derivative<Real, 2>(Real(1));
 668
 669     lanczos2(v1, v2);
 670
 671     for(size_t i = 0; i < v2.size(); ++i)
 672     {
 673         BOOST_CHECK_SMALL(v2[i], 10*tol);
 674     }
 675
 676 }
 677
 678 BOOST_AUTO_TEST_CASE(lanczos_smoothing_test)
 679 {
 680     test_dlp_second_derivative<double>();
 681     test_dlp_norms<double>();
 682     test_dlp_evaluation<double>();
 683     test_dlp_derivatives<double>();
 684     test_dlp_next<double>();
 685
 686     // Takes too long!
 687     //test_dlp_norms<cpp_bin_float_50>();
 688     test_boundary_velocity_filters<double>();
 689     test_boundary_velocity_filters<long double>();
 690     // Takes too long!
 691     //test_boundary_velocity_filters<cpp_bin_float_50>();
 692     test_boundary_lanczos<double>();
 693     test_boundary_lanczos<long double>();
 694     // Takes too long!
 695     //test_boundary_lanczos<cpp_bin_float_50>();
 696
 697     test_interior_velocity_filter<double>();
 698     test_interior_velocity_filter<long double>();
 699     test_interior_lanczos<double>();
 700
 701     test_acceleration_filters<double>();
 702
 703     test_lanczos_acceleration<float>();
 704     test_lanczos_acceleration<double>();
 705
 706     test_rescaling<double>();
 707     test_data_representations<double>();
 708 }