CubicSplineClass/example.cpp at master · joshhooker/CubicSplineClass · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#include <array>
#include <chrono>
#include <cmath>
#include <iostream>
#include <string>
#include <vector>

#include "CubicSpline.h"

// Timing helper
// returns wall time in nanoseconds for a callable
template<typename Fn>
double time_ns(Fn&& fn) {
    const auto t0 = std::chrono::high_resolution_clock::now();
    fn();
    const auto t1 = std::chrono::high_resolution_clock::now();
    return static_cast<double>(
        std::chrono::duration_cast<std::chrono::nanoseconds>(t1 - t0).count());
}

// Evaluation Loop
// returns sum of |true - spline|
static std::vector<double> make_query_points(double lo, double hi, double step) {
    std::vector<double> pts;
    for (double v = lo; v <= hi; v += step)
        pts.push_back(v);
    return pts;
}

template<typename Spline>
static double eval_error(const Spline& f, const std::vector<double>& pts) {
    double err = 0.0;
    for (double x : pts)
        err += std::fabs(std::sin(2.0 * std::cos(3.0 * x)) - f(x));
    return err / static_cast<double>(pts.size());
}

// Benchmark
// setup + eval, averaged over `iter` iterations.
struct Result {
    double setup_us;
    double eval_us;
    double mean_err;
};

template<typename MakeSpline>
Result benchmark(MakeSpline make_spline,
                 const std::vector<double>& query_pts,
                 const int iter,
                 const int warmup = 5) {
    double total_setup = 0.0;
    double total_eval  = 0.0;
    double last_err    = 0.0;

    for (int i = -warmup; i < iter; i++) {
        double setup_ns = time_ns([&]{ make_spline(); });

        auto f = make_spline();  // rebuild for eval (same cost, already timed above)
        double eval_ns  = time_ns([&]{ last_err = eval_error(f, query_pts); });

        if (i >= 0) {
            total_setup += setup_ns;
            total_eval  += eval_ns;
        }
    }

    const double scale = 1.0 / (iter * 1000.0);  // ns -> us
    return { total_setup * scale, total_eval * scale, last_err };
}

static void print_result(const std::string& label, const Result& r,
                         std::size_t n_queries) {
    std::cout << "  " << label << "\n"
              << "    Setup:  " << r.setup_us  << " us/call\n"
              << "    Eval:   " << r.eval_us   << " us / "
              << n_queries << " queries\n"
              << "    Error:  " << r.mean_err  << " (mean |true - spline|)\n\n";
}

int main() {

    std::vector<double> xv, yv;
    std::array<double, 110> xArr{}, yArr{};
    double xC[110], yC[110];

    int j = 0;
    for (double v = -2.5; v <= 8.5; v += 0.1) {
        double fv = std::sin(2.0 * std::cos(3.0 * v));
        xv.push_back(v);
        yv.push_back(fv);
        if (j < 110) {
            xArr[j] = v; yArr[j] = fv;
            xC[j]   = v; yC[j]   = fv;
            j++;
        }
    }

    // Points to evaluate splines
    const auto query_pts = make_query_points(-1.0, 5.0, 0.0034);

    constexpr int iter   = 1000;
    constexpr int warmup = 10;

    std::cout << "Initial points in spline: " << xv.size()
              << "   Queries per eval: " << query_pts.size()
              << "   Iterations: " << iter << "\n\n";

    // std::vector
    {
        auto r = benchmark([&]{ return CubicSpline(xv, yv); },
                           query_pts, iter, warmup);
        print_result("std::vector", r, query_pts.size());
    }

    // std::array
    {
        auto r = benchmark([&]{ return CubicSpline(xArr, yArr); },
                           query_pts, iter, warmup);
        print_result("std::array", r, query_pts.size());
    }

    // c array
    {
        auto r = benchmark([&]{ return CubicSpline(xC, yC); },
                           query_pts, iter, warmup);
        print_result("C array", r, query_pts.size());
    }

    return 0;
}