Differentiate between exp and expm1 (#15)

* need to differentiate between exp and exp1m because it provides different usecases.

src/CubicInterpolation/Axis.h

@@ -81,7 +81,7 @@ template <typename T> std::ostream &operator<<(std::ostream &out, const Axis<T>
 /**
  * @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 <typename T> class ExpAxis : public Axis<T> {
 
@@ -109,6 +109,32 @@ template <typename T> class ExpAxis : public Axis<T> {
   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 <typename T> class ExpM1Axis : public Axis<T> {
+
+  void print(std::ostream &os) const {
+    os << "ExpM1Axis(";
+    Axis<T>::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.

src/detail/Axis.cxx

@@ -11,29 +11,64 @@ template <typename T> unsigned int Axis<T>::required_nodes() const {
   return std::round(transform(high) + 1u);
 }
 
+//
+// Exp Axis
+//
+
 template <typename T>
 ExpAxis<T>::ExpAxis(T _low, T _high, T _stepsize) : Axis<T>(_low, _high, _stepsize) {}
 
 template <typename T>
 ExpAxis<T>::ExpAxis(T _low, T _high, size_t _nodes) : Axis<T>(_low, _high, 0.f) {
-  this->stepsize = (std::log1p(_high / _low) - M_LN2) / static_cast<T>(_nodes - 1);
+  this->stepsize = std::log(_high / _low) / static_cast<T>(_nodes - 1);
 }
 
 template <typename T> T ExpAxis<T>::transform(T x) const {
-  return (std::log1p(x / this->low) - M_LN2) / this->stepsize;
+  return std::log(x / this->low) / this->stepsize;
 }
 
 template <typename T> T ExpAxis<T>::back_transform(T t) const {
-  return this->low * std::expm1(t * this->stepsize + M_LN2);
+  return this->low * std::exp(t * this->stepsize);
 }
 template <typename T> T ExpAxis<T>::derive(T x) const {
-  return 1. / (this->stepsize * (x + this->low));
+  return 1. / (x * this->stepsize);
 }
 
 template <typename T> T ExpAxis<T>::back_derive(T t) const {
+  return this->low * this->stepsize * std::exp(t * this->stepsize);
+}
+
+//
+// ExpM1 Axis
+//
+
+template <typename T>
+ExpM1Axis<T>::ExpM1Axis(T _low, T _high, T _stepsize) : Axis<T>(_low, _high, _stepsize) {}
+
+template <typename T>
+ExpM1Axis<T>::ExpM1Axis(T _low, T _high, size_t _nodes) : Axis<T>(_low, _high, 0.f) {
+  this->stepsize = (std::log1p(_high / _low) - M_LN2) / static_cast<T>(_nodes - 1);
+}
+
+template <typename T> T ExpM1Axis<T>::transform(T x) const {
+  return (std::log1p(x / this->low) - M_LN2) / this->stepsize;
+}
+
+template <typename T> T ExpM1Axis<T>::back_transform(T t) const {
+  return this->low * std::expm1(t * this->stepsize + M_LN2);
+}
+template <typename T> T ExpM1Axis<T>::derive(T x) const {
+  return 1. / (this->stepsize * (x + this->low));
+}
+
+template <typename T> T ExpM1Axis<T>::back_derive(T t) const {
   return this->low * this->stepsize * std::exp(t * this->stepsize + M_LN2);
 }
 
+//
+// LinAxis
+//
+
 template <typename T>
 LinAxis<T>::LinAxis(T _low, T _high, T _stepsize) : Axis<T>(_low, _high, _stepsize) {}
 
@@ -60,4 +95,6 @@ template class LinAxis<double>;
 template class LinAxis<float>;
 template class ExpAxis<double>;
 template class ExpAxis<float>;
+template class ExpM1Axis<double>;
+template class ExpM1Axis<float>;
 } // namespace cubic_splines

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<cubic_splines::ExpAxis<double>>(low, high, N);
+  def.axis = std::make_unique<cubic_splines::ExpM1Axis<double>>(low, high, N);
   auto spline = cubic_splines::Interpolant<spline_t>(std::move(def), "", "");
   std::uniform_real_distribution<double> 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<cubic_splines::ExpAxis<double>>(1, 0);
-  def.axis = std::make_unique<cubic_splines::ExpAxis<double>>(low, high, N);
+  def.f_trafo = std::make_unique<cubic_splines::ExpM1Axis<double>>(1, 0);
+  def.axis = std::make_unique<cubic_splines::ExpM1Axis<double>>(low, high, N);
   auto spline = cubic_splines::Interpolant<spline_t>(std::move(def), "", "");
   std::uniform_real_distribution<double> 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<cubic_splines::ExpAxis<double>>(1, 0);
-  def.axis = std::make_unique<cubic_splines::ExpAxis<double>>(low, high, N);
+  def.f_trafo = std::make_unique<cubic_splines::ExpM1Axis<double>>(1, 0);
+  def.axis = std::make_unique<cubic_splines::ExpM1Axis<double>>(low, high, N);
   auto spline = cubic_splines::Interpolant<spline_t>(std::move(def), "", "");
   std::uniform_real_distribution<double> 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<cubic_splines::ExpAxis<double>>(low, high, N);
+  def.axis = std::make_unique<cubic_splines::ExpM1Axis<double>>(low, high, N);
   auto spline = cubic_splines::Interpolant<spline_t>(std::move(def), "", "");
   std::uniform_real_distribution<double> 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<cubic_splines::ExpAxis<double>>(1, 0);
+  def.f_trafo = std::make_unique<cubic_splines::ExpM1Axis<double>>(1, 0);
   def.axis = std::make_unique<cubic_splines::LinAxis<double>>(low, high, N);
   auto spline = cubic_splines::Interpolant<spline_t>(std::move(def), "", "");
   std::uniform_real_distribution<double> dis(low, high);