generic-library/boost
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
boost/libs/math/example/daubechies_wavelets/bench.cpp
Edouard DUPIN d0115b733d [DEV] add v1.76.0
2021-10-05 21:37:46 +02:00

542 lines
16 KiB
C++
Raw Permalink Blame History

/*
* Copyright Nick Thompson, 2020
* 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 <cmath>
#include <random>
#include <benchmark/benchmark.h>
#include <boost/random/uniform_real_distribution.hpp>
#include <boost/math/special_functions/daubechies_scaling.hpp>
#include <boost/math/quadrature/wavelet_transforms.hpp>
#include <boost/math/interpolators/cubic_hermite.hpp>
#include <boost/math/interpolators/detail/quintic_hermite_detail.hpp>
#include <boost/math/interpolators/detail/septic_hermite_detail.hpp>
double exponential(benchmark::IterationCount j)
{
return std::pow(2, j);
}
template<typename Real, int p>
void DyadicGrid(benchmark::State & state)
{
int j = state.range(0);
size_t s = 0;
for (auto _ : state)
{
auto v = boost::math::daubechies_scaling_dyadic_grid<Real, 4, 0>(j);
benchmark::DoNotOptimize(v[0]);
s = v.size();
}
state.counters["RAM"] = s*sizeof(Real);
state.SetComplexityN(state.range(0));
}
BENCHMARK_TEMPLATE(DyadicGrid, double, 4)->DenseRange(3, 22, 1)->Unit(benchmark::kMillisecond)->Complexity(exponential);
//BENCHMARK_TEMPLATE(DyadicGrid, double, 8)->DenseRange(3, 22, 1)->Unit(benchmark::kMillisecond)->Complexity(exponential);
//BENCHMARK_TEMPLATE(DyadicGrid, double, 11)->DenseRange(3,22,1)->Unit(benchmark::kMillisecond)->Complexity(exponential);
uint64_t s[2] = { 0x41, 0x29837592 };
static inline uint64_t rotl(const uint64_t x, int k) {
return (x << k) | (x >> (64 - k));
}
uint64_t next(void) {
const uint64_t s0 = s[0];
uint64_t s1 = s[1];
const uint64_t result = s0 + s1;
s1 ^= s0;
s[0] = rotl(s0, 55) ^ s1 ^ (s1 << 14); // a, b
s[1] = rotl(s1, 36); // c
return result;
}
double uniform() {
return next()*(1.0/18446744073709551616.0);
}
template<typename Real, int p>
void ScalingEvaluation(benchmark::State & state)
{
auto phi = boost::math::daubechies_scaling<Real, p>();
Real x = 0;
Real step = std::numeric_limits<Real>::epsilon();
for (auto _ : state)
{
benchmark::DoNotOptimize(phi(x));
x += step;
}
}
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 2);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 3);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 4);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 5);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 6);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 7);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 8);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 9);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 10);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 11);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 12);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 13);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 14);
BENCHMARK_TEMPLATE(ScalingEvaluation, double, 15);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 2);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 3);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 4);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 5);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 6);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 7);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 8);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 9);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 10);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 11);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 12);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 13);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 14);
BENCHMARK_TEMPLATE(ScalingEvaluation, float, 15);
template<typename Real, int p>
void ScalingConstructor(benchmark::State & state)
{
for (auto _ : state)
{
benchmark::DoNotOptimize(boost::math::daubechies_scaling<Real, p>());
}
}
BENCHMARK_TEMPLATE(ScalingConstructor, float, 2)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, double, 2)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, long double, 2)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, float, 3)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, double, 3)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, long double, 3)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, float, 4)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, double, 4)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, long double, 4)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, float, 5)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, double, 5)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, long double, 5)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, float, 11)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, double, 11)->Unit(benchmark::kMillisecond);
BENCHMARK_TEMPLATE(ScalingConstructor, long double, 11)->Unit(benchmark::kMillisecond);
template<typename Real>
void CubicHermite(benchmark::State & state)
{
using boost::math::interpolators::cubic_hermite;
auto n = state.range(0);
std::vector<Real> x(n);
std::vector<Real> y(n);
std::vector<Real> dydx(n);
std::random_device rd;
boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
x[0] = dis(rd);
y[0] = dis(rd);
dydx[0] = dis(rd);
for (size_t i = 1; i < y.size(); ++i)
{
x[i] = x[i-1] + dis(rd);
y[i] = dis(rd);
dydx[i] = dis(rd);
}
Real x0 = x.front();
Real xf = x.back();
auto qh = cubic_hermite(std::move(x), std::move(y), std::move(dydx));
Real t = x0;
Real step = uniform()*(xf-x0)/2048;
for (auto _ : state)
{
benchmark::DoNotOptimize(qh(t));
t += step;
if (t >= xf)
{
t = x0;
step = uniform()*(xf-x0)/2048;
}
}
state.SetComplexityN(state.range(0));
}
BENCHMARK_TEMPLATE(CubicHermite, double)->RangeMultiplier(2)->Range(1<<8, 1<<20)->Complexity(benchmark::oLogN);
template<typename Real>
void CardinalCubicHermite(benchmark::State & state)
{
using boost::math::interpolators::detail::cardinal_cubic_hermite_detail;
auto n = state.range(0);
std::vector<Real> y(n);
std::vector<Real> dydx(n);
std::random_device rd;
boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
for (size_t i = 0; i < y.size(); ++i)
{
y[i] = uniform();
dydx[i] = uniform();
}
Real dx = Real(1)/Real(8);
Real x0 = 0;
Real xf = x0 + (y.size()-1)*dx;
auto qh = cardinal_cubic_hermite_detail(std::move(y), std::move(dydx), x0, dx);
Real x = x0;
Real step = uniform()*(xf-x0)/2048;
for (auto _ : state)
{
benchmark::DoNotOptimize(qh.unchecked_evaluation(x));
x += step;
if (x >= xf)
{
x = x0;
step = uniform()*(xf-x0)/2048;
}
}
state.SetComplexityN(state.range(0));
}
template<typename Real>
void CardinalCubicHermiteAOS(benchmark::State & state)
{
auto n = state.range(0);
std::vector<std::array<Real, 2>> dat(n);
std::random_device rd;
boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
for (size_t i = 0; i < dat.size(); ++i)
{
dat[i][0] = uniform();
dat[i][1] = uniform();
}
using boost::math::interpolators::detail::cardinal_cubic_hermite_detail_aos;
Real dx = Real(1)/Real(8);
Real x0 = 0;
Real xf = x0 + (dat.size()-1)*dx;
auto qh = cardinal_cubic_hermite_detail_aos(std::move(dat), x0, dx);
Real x = 0;
Real step = uniform()*(xf-x0)/2048;
for (auto _ : state)
{
benchmark::DoNotOptimize(qh.unchecked_evaluation(x));
x += step;
if (x >= xf)
{
x = x0;
step = uniform()*(xf-x0)/2048;
}
}
state.SetComplexityN(state.range(0));
}
BENCHMARK_TEMPLATE(CardinalCubicHermiteAOS, double)->RangeMultiplier(2)->Range(1<<8, 1<<21)->Complexity(benchmark::o1);
BENCHMARK_TEMPLATE(CardinalCubicHermite, double)->RangeMultiplier(2)->Range(1<<8, 1<<21)->Complexity(benchmark::o1);
template<class Real>
void SineEvaluation(benchmark::State& state)
{
std::default_random_engine gen;
std::uniform_real_distribution<Real> x_dis(0, 3.14159);
Real x = x_dis(gen);
for (auto _ : state)
{
benchmark::DoNotOptimize(std::sin(x));
x += std::numeric_limits<Real>::epsilon();
}
}
BENCHMARK_TEMPLATE(SineEvaluation, float);
BENCHMARK_TEMPLATE(SineEvaluation, double);
BENCHMARK_TEMPLATE(SineEvaluation, long double);
template<class Real>
void ExpEvaluation(benchmark::State& state)
{
std::default_random_engine gen;
std::uniform_real_distribution<Real> x_dis(0, 3.14159);
Real x = x_dis(gen);
for (auto _ : state)
{
benchmark::DoNotOptimize(std::exp(x));
x += std::numeric_limits<Real>::epsilon();
}
}
BENCHMARK_TEMPLATE(ExpEvaluation, float);
BENCHMARK_TEMPLATE(ExpEvaluation, double);
BENCHMARK_TEMPLATE(ExpEvaluation, long double);
template<class Real>
void PowEvaluation(benchmark::State& state)
{
std::default_random_engine gen;
std::uniform_real_distribution<Real> x_dis(0, 3.14159);
Real x = x_dis(gen);
for (auto _ : state)
{
benchmark::DoNotOptimize(std::pow(x, x+1));
x += std::numeric_limits<Real>::epsilon();
}
}
BENCHMARK_TEMPLATE(PowEvaluation, float);
BENCHMARK_TEMPLATE(PowEvaluation, double);
BENCHMARK_TEMPLATE(PowEvaluation, long double);
template<typename Real>
void CardinalQuinticHermite(benchmark::State & state)
{
using boost::math::interpolators::detail::cardinal_quintic_hermite_detail;
auto n = state.range(0);
std::vector<Real> y(n);
std::vector<Real> dydx(n);
std::vector<Real> d2ydx2(n);
std::random_device rd;
boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
for (size_t i = 0; i < y.size(); ++i)
{
y[i] = uniform();
dydx[i] = uniform();
d2ydx2[i] = uniform();
}
Real dx = Real(1)/Real(8);
Real x0 = 0;
Real xf = x0 + (y.size()-1)*dx;
auto qh = cardinal_quintic_hermite_detail(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx);
Real x = 0;
Real step = uniform()*(xf-x0)/2048;
for (auto _ : state)
{
benchmark::DoNotOptimize(qh.unchecked_evaluation(x));
x += step;
if (x >= xf)
{
x = x0;
step = uniform()*(xf-x0)/2048;
}
}
state.SetComplexityN(state.range(0));
}
template<typename Real>
void CardinalQuinticHermiteAOS(benchmark::State & state)
{
auto n = state.range(0);
std::vector<std::array<Real, 3>> dat(n);
std::random_device rd;
boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
for (size_t i = 0; i < dat.size(); ++i)
{
dat[i][0] = uniform();
dat[i][1] = uniform();
dat[i][2] = uniform();
}
using boost::math::interpolators::detail::cardinal_quintic_hermite_detail_aos;
Real dx = Real(1)/Real(8);
Real x0 = 0;
Real xf = x0 + (dat.size()-1)*dx;
auto qh = cardinal_quintic_hermite_detail_aos(std::move(dat), x0, dx);
Real x = x0;
Real step = uniform()*(xf-x0)/2048;
for (auto _ : state) {
benchmark::DoNotOptimize(qh.unchecked_evaluation(x));
x += step;
if (x >= xf)
{
x = x0;
step = uniform()*(xf-x0)/2048;
}
}
state.SetComplexityN(state.range(0));
}
BENCHMARK_TEMPLATE(CardinalQuinticHermiteAOS, double)->RangeMultiplier(2)->Range(1<<8, 1<<22)->Complexity(benchmark::o1);
BENCHMARK_TEMPLATE(CardinalQuinticHermite, double)->RangeMultiplier(2)->Range(1<<8, 1<<22)->Complexity(benchmark::o1);
template<typename Real>
void SepticHermite(benchmark::State & state)
{
using boost::math::interpolators::detail::septic_hermite_detail;
auto n = state.range(0);
std::vector<Real> x(n);
std::vector<Real> y(n);
std::vector<Real> dydx(n);
std::vector<Real> d2ydx2(n);
std::vector<Real> d3ydx3(n);
std::random_device rd;
boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
Real x0 = dis(rd);
x[0] = x0;
for (decltype(n) i = 1; i < n; ++i)
{
x[i] = x[i-1] + dis(rd);
}
for (size_t i = 0; i < y.size(); ++i)
{
y[i] = dis(rd);
dydx[i] = dis(rd);
d2ydx2[i] = dis(rd);
d3ydx3[i] = dis(rd);
}
Real xf = x.back();
auto sh = septic_hermite_detail(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
Real t = x0;
for (auto _ : state)
{
benchmark::DoNotOptimize(sh(t));
t += xf/128;
if (t >= xf)
{
t = x0;
}
}
state.SetComplexityN(state.range(0));
}
BENCHMARK_TEMPLATE(SepticHermite, double)->RangeMultiplier(2)->Range(1<<8, 1<<20)->Complexity();
template<typename Real>
void CardinalSepticHermite(benchmark::State & state)
{
using boost::math::interpolators::detail::cardinal_septic_hermite_detail;
auto n = state.range(0);
std::vector<Real> y(n);
std::vector<Real> dydx(n);
std::vector<Real> d2ydx2(n);
std::vector<Real> d3ydx3(n);
std::random_device rd;
boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
for (size_t i = 0; i < y.size(); ++i)
{
y[i] = dis(rd);
dydx[i] = dis(rd);
d2ydx2[i] = dis(rd);
d3ydx3[i] = dis(rd);
}
Real dx = Real(1)/Real(8);
Real x0 = 0;
Real xf = x0 + (y.size()-1)*dx;
auto sh = cardinal_septic_hermite_detail(std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), x0, dx);
Real x = 0;
for (auto _ : state)
{
benchmark::DoNotOptimize(sh.unchecked_evaluation(x));
x += xf/128;
if (x >= xf)
{
x = x0;
}
}
state.SetComplexityN(state.range(0));
}
BENCHMARK_TEMPLATE(CardinalSepticHermite, double)->RangeMultiplier(2)->Range(1<<8, 1<<20)->Complexity();
template<typename Real>
void CardinalSepticHermiteAOS(benchmark::State & state)
{
using boost::math::interpolators::detail::cardinal_septic_hermite_detail_aos;
auto n = state.range(0);
std::vector<std::array<Real, 4>> data(n);
std::random_device rd;
boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
for (size_t i = 0; i < data.size(); ++i)
{
for (size_t j = 0; j < 4; ++j)
{
data[i][j] = dis(rd);
}
}
Real dx = Real(1)/Real(8);
Real x0 = 0;
Real xf = x0 + (data.size()-1)*dx;
auto sh = cardinal_septic_hermite_detail_aos(std::move(data), x0, dx);
Real x = 0;
for (auto _ : state)
{
benchmark::DoNotOptimize(sh.unchecked_evaluation(x));
x += xf/128;
if (x >= xf)
{
x = x0;
}
}
state.SetComplexityN(state.range(0));
}
BENCHMARK_TEMPLATE(CardinalSepticHermiteAOS, double)->RangeMultiplier(2)->Range(1<<8, 1<<20)->Complexity();
template<typename Real, int p>
void WaveletTransform(benchmark::State & state)
{
auto psi = boost::math::daubechies_wavelet<Real, p>();
auto f = [](Real t) {
return std::exp(-t*t);
};
auto Wf = boost::math::quadrature::daubechies_wavelet_transform(f, psi);
for (auto _ : state)
{
Real s = 1 + uniform();
Real t = uniform();
benchmark::DoNotOptimize(Wf(s, t));
}
}
BENCHMARK_TEMPLATE(WaveletTransform, float, 3);
BENCHMARK_TEMPLATE(WaveletTransform, float, 4);
BENCHMARK_TEMPLATE(WaveletTransform, float, 5);
BENCHMARK_TEMPLATE(WaveletTransform, float, 6);
BENCHMARK_TEMPLATE(WaveletTransform, float, 7);
BENCHMARK_TEMPLATE(WaveletTransform, float, 8);
BENCHMARK_TEMPLATE(WaveletTransform, float, 9);
BENCHMARK_TEMPLATE(WaveletTransform, float, 10);
BENCHMARK_TEMPLATE(WaveletTransform, float, 11);
BENCHMARK_TEMPLATE(WaveletTransform, float, 12);
BENCHMARK_TEMPLATE(WaveletTransform, float, 13);
BENCHMARK_TEMPLATE(WaveletTransform, float, 14);
BENCHMARK_TEMPLATE(WaveletTransform, float, 15);
BENCHMARK_TEMPLATE(WaveletTransform, float, 16);
BENCHMARK_TEMPLATE(WaveletTransform, float, 17);
BENCHMARK_TEMPLATE(WaveletTransform, float, 18);
BENCHMARK_TEMPLATE(WaveletTransform, float, 19);
BENCHMARK_TEMPLATE(WaveletTransform, double, 4);
BENCHMARK_TEMPLATE(WaveletTransform, double, 8);
BENCHMARK_TEMPLATE(WaveletTransform, double, 12);
BENCHMARK_TEMPLATE(WaveletTransform, double, 15);
BENCHMARK_TEMPLATE(WaveletTransform, double, 19);
BENCHMARK_MAIN();
Reference in New Issue View Git Blame Copy Permalink