libs/math/example/cardinal_cubic_b_spline_example.cpp

   1 // Copyright Nicholas Thompson 2017.
   2 // Copyright Paul A. Bristow 2017.
   3 // Copyright John Maddock 2017.
   4
   5 // Distributed under the Boost Software License, Version 1.0.
   6 // (See accompanying file LICENSE_1_0.txt or
   7 // copy at http://www.boost.org/LICENSE_1_0.txt).
   8
   9 #include <iostream>
  10 #include <limits>
  11 #include <vector>
  12 #include <algorithm>
  13 #include <iomanip>
  14 #include <iterator>
  15 #include <cmath>
  16 #include <random>
  17 #include <cstdint>
  18 #include <boost/random/uniform_real_distribution.hpp>
  19 #include <boost/math/tools/roots.hpp>
  20
  21 //[cubic_b_spline_example
  22
  23 /*`This example demonstrates how to use the cubic b spline interpolator for regularly spaced data.
  24 */
  25 #include <boost/math/interpolators/cardinal_cubic_b_spline.hpp>
  26
  27 int main()
  28 {
  29     // We begin with an array of samples:
  30     std::vector<double> v(500);
  31     // And decide on a stepsize:
  32     double step = 0.01;
  33
  34     // Initialize the vector with a function we'd like to interpolate:
  35     for (size_t i = 0; i < v.size(); ++i)
  36     {
  37         v[i] = sin(i*step);
  38     }
  39     // We could define an arbitrary start time, but for now we'll just use 0:
  40     boost::math::interpolators::cardinal_cubic_b_spline<double> spline(v.data(), v.size(), 0 /* start time */, step);
  41
  42     // Now we can evaluate the spline wherever we please.
  43     std::mt19937 gen;
  44     boost::random::uniform_real_distribution<double> absissa(0, v.size()*step);
  45     for (size_t i = 0; i < 10; ++i)
  46     {
  47         double x = absissa(gen);
  48         std::cout << "sin(" << x << ") = " << sin(x) << ", spline interpolation gives " << spline(x) << std::endl;
  49         std::cout << "cos(" << x << ") = " << cos(x) << ", spline derivative interpolation gives " << spline.prime(x) << std::endl;
  50     }
  51
  52     // The next example is less trivial:
  53     // We will try to figure out when the population of the United States crossed 100 million.
  54     // Since the census is taken every 10 years, the data is equally spaced, so we can use the cubic b spline.
  55     // Data taken from https://en.wikipedia.org/wiki/United_States_Census
  56     // We'll start at the year 1860:
  57     double t0 = 1860;
  58     double time_step = 10;
  59     std::vector<double> population{31443321,  /* 1860 */
  60                                    39818449,  /* 1870 */
  61                                    50189209,  /* 1880 */
  62                                    62947714,  /* 1890 */
  63                                    76212168,  /* 1900 */
  64                                    92228496,  /* 1910 */
  65                                    106021537, /* 1920 */
  66                                    122775046, /* 1930 */
  67                                    132164569, /* 1940 */
  68                                    150697361, /* 1950 */
  69                                    179323175};/* 1960 */
  70
  71     // An eyeball estimate indicates that the population crossed 100 million around 1915.
  72     // Let's see what interpolation says:
  73     boost::math::interpolators::cardinal_cubic_b_spline<double> p(population.data(), population.size(), t0, time_step);
  74
  75     // Now create a function which has a zero at p = 100,000,000:
  76     auto f = [=](double t){ return p(t) - 100000000; };
  77
  78     // Boost includes a bisection algorithm, which is robust, though not as fast as some others
  79     // we provide, but let's try that first.  We need a termination condition for it, which
  80     // takes the two endpoints of the range and returns either true (stop) or false (keep going),
  81     // we could use a predefined one such as boost::math::tools::eps_tolerance<double>, but that
  82     // won't stop until we have full double precision which is overkill, since we just need the
  83     // endpoint to yield the same month.  While we're at it, we'll keep track of the number of
  84     // iterations required too, though this is strictly optional:
  85
  86     auto termination = [](double left, double right)
  87     {
  88        double left_month = std::round((left - std::floor(left)) * 12 + 1);
  89        double right_month = std::round((right - std::floor(right)) * 12 + 1);
  90        return (left_month == right_month) && (std::floor(left) == std::floor(right));
  91     };
  92     std::uintmax_t iterations = 1000;
  93     auto result =  boost::math::tools::bisect(f, 1910.0, 1920.0, termination, iterations);
  94     auto time = result.first;  // termination condition ensures that both endpoints yield the same result
  95     auto month = std::round((time - std::floor(time))*12  + 1);
  96     auto year = std::floor(time);
  97     std::cout << "The population of the United States surpassed 100 million on the ";
  98     std::cout << month << "th month of " << year << std::endl;
  99     std::cout << "Found in " << iterations << " iterations" << std::endl;
 100
 101     // Since the cubic B spline offers the first derivative, we could equally have used Newton iterations,
 102     // this takes "number of bits correct" as a termination condition - 20 should be plenty for what we need,
 103     // and once again, we track how many iterations are taken:
 104
 105     auto f_n = [=](double t) { return std::make_pair(p(t) - 100000000, p.prime(t)); };
 106     iterations = 1000;
 107     time = boost::math::tools::newton_raphson_iterate(f_n, 1910.0, 1900.0, 2000.0, 20, iterations);
 108     month = std::round((time - std::floor(time))*12  + 1);
 109     year = std::floor(time);
 110     std::cout << "The population of the United States surpassed 100 million on the ";
 111     std::cout << month << "th month of " << year << std::endl;
 112     std::cout << "Found in " << iterations << " iterations" << std::endl;
 113
 114 }
 115
 116 //] [/cubic_b_spline_example]
 117
 118 //[cubic_b_spline_example_out
 119 /*` Program output is:
 120 [pre
 121 sin(4.07362) = -0.802829, spline interpolation gives - 0.802829
 122 cos(4.07362) = -0.596209, spline derivative interpolation gives - 0.596209
 123 sin(0.677385) = 0.626758, spline interpolation gives 0.626758
 124 cos(0.677385) = 0.779214, spline derivative interpolation gives 0.779214
 125 sin(4.52896) = -0.983224, spline interpolation gives - 0.983224
 126 cos(4.52896) = -0.182402, spline derivative interpolation gives - 0.182402
 127 sin(4.17504) = -0.85907, spline interpolation gives - 0.85907
 128 cos(4.17504) = -0.511858, spline derivative interpolation gives - 0.511858
 129 sin(0.634934) = 0.593124, spline interpolation gives 0.593124
 130 cos(0.634934) = 0.805111, spline derivative interpolation gives 0.805111
 131 sin(4.84434) = -0.991307, spline interpolation gives - 0.991307
 132 cos(4.84434) = 0.131567, spline derivative interpolation gives 0.131567
 133 sin(4.56688) = -0.989432, spline interpolation gives - 0.989432
 134 cos(4.56688) = -0.144997, spline derivative interpolation gives - 0.144997
 135 sin(1.10517) = 0.893541, spline interpolation gives 0.893541
 136 cos(1.10517) = 0.448982, spline derivative interpolation gives 0.448982
 137 sin(3.1618) = -0.0202022, spline interpolation gives - 0.0202022
 138 cos(3.1618) = -0.999796, spline derivative interpolation gives - 0.999796
 139 sin(1.54084) = 0.999551, spline interpolation gives 0.999551
 140 cos(1.54084) = 0.0299566, spline derivative interpolation gives 0.0299566
 141 The population of the United States surpassed 100 million on the 11th month of 1915
 142 Found in 12 iterations
 143 The population of the United States surpassed 100 million on the 11th month of 1915
 144 Found in 3 iterations
 145 ]
 146 */
 147 //]