ceph/src/boost/libs/math/example/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/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::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     // An eye
  57     // We'll start at the year 1860:
  58     double t0 = 1860;
  59     double time_step = 10;
  60     std::vector<double> population{31443321,  /* 1860 */
  61                                    39818449,  /* 1870 */
  62                                    50189209,  /* 1880 */
  63                                    62947714,  /* 1890 */
  64                                    76212168,  /* 1900 */
  65                                    92228496,  /* 1910 */
  66                                    106021537, /* 1920 */
  67                                    122775046, /* 1930 */
  68                                    132164569, /* 1940 */
  69                                    150697361, /* 1950 */
  70                                    179323175};/* 1960 */
  71
  72     // An eyeball estimate indicates that the population crossed 100 million around 1915.
  73     // Let's see what interpolation says:
  74     boost::math::cubic_b_spline<double> p(population.data(), population.size(), t0, 10);
  75
  76     // Now create a function which has a zero at p = 100,000,000:
  77     auto f = [=](double t){ return p(t) - 100000000; };
  78
  79     // Boost includes a bisection algorithm, which is robust, though not as fast as some others
  80     // we provide, but lets try that first.  We need a termination condition for it, which
  81     // takes the two endpoints of the range and returns either true (stop) or false (keep going),
  82     // we could use a predefined one such as boost::math::tools::eps_tolerance<double>, but that
  83     // won't stop until we have full double precision which is overkill, since we just need the
  84     // endpoint to yield the same month.  While we're at it, we'll keep track of the number of
  85     // iterations required too, though this is strictly optional:
  86
  87     auto termination = [](double left, double right)
  88     {
  89        double left_month = std::round((left - std::floor(left)) * 12 + 1);
  90        double right_month = std::round((right - std::floor(right)) * 12 + 1);
  91        return (left_month == right_month) && (std::floor(left) == std::floor(right));
  92     };
  93     std::uintmax_t iterations = 1000;
  94     auto result =  boost::math::tools::bisect(f, 1910.0, 1920.0, termination, iterations);
  95     auto time = result.first;  // termination condition ensures that both endpoints yield the same result
  96     auto month = std::round((time - std::floor(time))*12  + 1);
  97     auto year = std::floor(time);
  98     std::cout << "The population of the United States surpassed 100 million on the ";
  99     std::cout << month << "th month of " << year << std::endl;
 100     std::cout << "Found in " << iterations << " iterations" << std::endl;
 101
 102     // Since the cubic B spline offers the first derivative, we could equally have used Newton iterations,
 103     // this takes "number of bits correct" as a termination condition - 20 should be plenty for what we need,
 104     // and once again, we track how many iterations are taken:
 105
 106     auto f_n = [=](double t) { return std::make_pair(p(t) - 100000000, p.prime(t)); };
 107     iterations = 1000;
 108     time = boost::math::tools::newton_raphson_iterate(f_n, 1910.0, 1900.0, 2000.0, 20, iterations);
 109     month = std::round((time - std::floor(time))*12  + 1);
 110     year = std::floor(time);
 111     std::cout << "The population of the United States surpassed 100 million on the ";
 112     std::cout << month << "th month of " << year << std::endl;
 113     std::cout << "Found in " << iterations << " iterations" << std::endl;
 114
 115 }
 116
 117 //] [/cubic_b_spline_example]
 118
 119 //[cubic_b_spline_example_out
 120 /*` Program output is:
 121 [pre
 122 sin(4.07362) = -0.802829, spline interpolation gives - 0.802829
 123 cos(4.07362) = -0.596209, spline derivative interpolation gives - 0.596209
 124 sin(0.677385) = 0.626758, spline interpolation gives 0.626758
 125 cos(0.677385) = 0.779214, spline derivative interpolation gives 0.779214
 126 sin(4.52896) = -0.983224, spline interpolation gives - 0.983224
 127 cos(4.52896) = -0.182402, spline derivative interpolation gives - 0.182402
 128 sin(4.17504) = -0.85907, spline interpolation gives - 0.85907
 129 cos(4.17504) = -0.511858, spline derivative interpolation gives - 0.511858
 130 sin(0.634934) = 0.593124, spline interpolation gives 0.593124
 131 cos(0.634934) = 0.805111, spline derivative interpolation gives 0.805111
 132 sin(4.84434) = -0.991307, spline interpolation gives - 0.991307
 133 cos(4.84434) = 0.131567, spline derivative interpolation gives 0.131567
 134 sin(4.56688) = -0.989432, spline interpolation gives - 0.989432
 135 cos(4.56688) = -0.144997, spline derivative interpolation gives - 0.144997
 136 sin(1.10517) = 0.893541, spline interpolation gives 0.893541
 137 cos(1.10517) = 0.448982, spline derivative interpolation gives 0.448982
 138 sin(3.1618) = -0.0202022, spline interpolation gives - 0.0202022
 139 cos(3.1618) = -0.999796, spline derivative interpolation gives - 0.999796
 140 sin(1.54084) = 0.999551, spline interpolation gives 0.999551
 141 cos(1.54084) = 0.0299566, spline derivative interpolation gives 0.0299566
 142 The population of the United States surpassed 100 million on the 11th month of 1915
 143 Found in 12 iterations
 144 The population of the United States surpassed 100 million on the 11th month of 1915
 145 Found in 3 iterations
 146 ]
 147 */
 148 //]