diff --git a/src/CubicInterpolation/Axis.h b/src/CubicInterpolation/Axis.h index 2e76da0..a58ef34 100644 --- a/src/CubicInterpolation/Axis.h +++ b/src/CubicInterpolation/Axis.h @@ -81,7 +81,7 @@ template std::ostream &operator<<(std::ostream &out, const Axis /** * @brief Exponential axis to interpolate over many different scales. As basis * the euler number choosen by default if no further specified. Nodes will - * distributed in form of \f$ \text{low} \cdot (\exp(n \cdot \text{stepsize} + \log(2)) - 1) \f$ + * distributed in form of \f$ \exp(n \cdot \text{stepsize}) \f$ */ template class ExpAxis : public Axis { @@ -109,6 +109,32 @@ template class ExpAxis : public Axis { T back_derive(T t) const final; }; +/** + * @brief Exponential axis to interpolate over many different scales. As basis + * the euler number choosen by default if no further specified. Nodes will + * distributed in form of \f$ \text{low} \cdot (\exp(n \cdot \text{stepsize} + \log(2)) - + * 1) \f$ + */ +template class ExpM1Axis : public Axis { + + void print(std::ostream &os) const { + os << "ExpM1Axis("; + Axis::print(os); + os << ")"; + } + +public: + ExpM1Axis(T _low, T _high, T _stepsize = 1); + + ExpM1Axis(T _low, T _high, size_t _nodes); + + T transform(T x) const final; + T back_transform(T t) const final; + + T derive(T x) const final; + T back_derive(T t) const final; +}; + /** * @brief Linear axis to describe data which varify not to much in orders of * magnitudes. It's the fastest axis evaluation. diff --git a/src/detail/Axis.cxx b/src/detail/Axis.cxx index 4ff09bd..23279d0 100644 --- a/src/detail/Axis.cxx +++ b/src/detail/Axis.cxx @@ -11,29 +11,64 @@ template unsigned int Axis::required_nodes() const { return std::round(transform(high) + 1u); } +// +// Exp Axis +// + template ExpAxis::ExpAxis(T _low, T _high, T _stepsize) : Axis(_low, _high, _stepsize) {} template ExpAxis::ExpAxis(T _low, T _high, size_t _nodes) : Axis(_low, _high, 0.f) { - this->stepsize = (std::log1p(_high / _low) - M_LN2) / static_cast(_nodes - 1); + this->stepsize = std::log(_high / _low) / static_cast(_nodes - 1); } template T ExpAxis::transform(T x) const { - return (std::log1p(x / this->low) - M_LN2) / this->stepsize; + return std::log(x / this->low) / this->stepsize; } template T ExpAxis::back_transform(T t) const { - return this->low * std::expm1(t * this->stepsize + M_LN2); + return this->low * std::exp(t * this->stepsize); } template T ExpAxis::derive(T x) const { - return 1. / (this->stepsize * (x + this->low)); + return 1. / (x * this->stepsize); } template T ExpAxis::back_derive(T t) const { + return this->low * this->stepsize * std::exp(t * this->stepsize); +} + +// +// ExpM1 Axis +// + +template +ExpM1Axis::ExpM1Axis(T _low, T _high, T _stepsize) : Axis(_low, _high, _stepsize) {} + +template +ExpM1Axis::ExpM1Axis(T _low, T _high, size_t _nodes) : Axis(_low, _high, 0.f) { + this->stepsize = (std::log1p(_high / _low) - M_LN2) / static_cast(_nodes - 1); +} + +template T ExpM1Axis::transform(T x) const { + return (std::log1p(x / this->low) - M_LN2) / this->stepsize; +} + +template T ExpM1Axis::back_transform(T t) const { + return this->low * std::expm1(t * this->stepsize + M_LN2); +} +template T ExpM1Axis::derive(T x) const { + return 1. / (this->stepsize * (x + this->low)); +} + +template T ExpM1Axis::back_derive(T t) const { return this->low * this->stepsize * std::exp(t * this->stepsize + M_LN2); } +// +// LinAxis +// + template LinAxis::LinAxis(T _low, T _high, T _stepsize) : Axis(_low, _high, _stepsize) {} @@ -60,4 +95,6 @@ template class LinAxis; template class LinAxis; template class ExpAxis; template class ExpAxis; +template class ExpM1Axis; +template class ExpM1Axis; } // namespace cubic_splines diff --git a/tests/TestCubicSplines.cpp b/tests/TestCubicSplines.cpp index d8ac626..fa03516 100644 --- a/tests/TestCubicSplines.cpp +++ b/tests/TestCubicSplines.cpp @@ -92,7 +92,7 @@ TEST(CubicSplines, evaluate_exp_distributed_nodes) { auto func = [](double x) { return std::exp(x); }; auto def = spline_def_t(); def.f = func; - def.axis = std::make_unique>(low, high, N); + def.axis = std::make_unique>(low, high, N); auto spline = cubic_splines::Interpolant(std::move(def), "", ""); std::uniform_real_distribution dis(low, high); for (int i = 0; i < 10'000; ++i) { @@ -108,8 +108,8 @@ TEST(CubicSplines, evaluate_log_func_values_and_axis) { auto def = spline_def_t(); auto func = [](double x) { return x*x*std::log(x); }; def.f = func; - def.f_trafo = std::make_unique>(1, 0); - def.axis = std::make_unique>(low, high, N); + def.f_trafo = std::make_unique>(1, 0); + def.axis = std::make_unique>(low, high, N); auto spline = cubic_splines::Interpolant(std::move(def), "", ""); std::uniform_real_distribution dis(low, high); for (int i = 0; i < 10'000; ++i) { @@ -126,8 +126,8 @@ TEST(CubicSplines, evaluate_log_func_values_and_axis_trivial_integrand) { // this function should be a straight line in the transformed space auto func = [](double x) { return std::pow(x+1, 10); }; def.f = func; - def.f_trafo = std::make_unique>(1, 0); - def.axis = std::make_unique>(low, high, N); + def.f_trafo = std::make_unique>(1, 0); + def.axis = std::make_unique>(low, high, N); auto spline = cubic_splines::Interpolant(std::move(def), "", ""); std::uniform_real_distribution dis(low, high); for (int i = 0; i < 10'000; ++i) { @@ -161,7 +161,7 @@ TEST(CubicSplines, prime_x_trafo) { auto func = [](double x) { return x * x + x + 1; }; auto df_dx = [](double x) { return 2 * x + 1; }; def.f = func; - def.axis = std::make_unique>(low, high, N); + def.axis = std::make_unique>(low, high, N); auto spline = cubic_splines::Interpolant(std::move(def), "", ""); std::uniform_real_distribution dis(low, high); for (int i = 0; i < 10'000; ++i) { @@ -178,7 +178,7 @@ TEST(CubicSplines, prime_y_trafo) { auto func = [](double x) { return x * x + x + 1; }; auto df_dx = [](double x) { return 2 * x + 1; }; def.f = func; - def.f_trafo = std::make_unique>(1, 0); + def.f_trafo = std::make_unique>(1, 0); def.axis = std::make_unique>(low, high, N); auto spline = cubic_splines::Interpolant(std::move(def), "", ""); std::uniform_real_distribution dis(low, high);