alternative spline

catima.h

@@ -27,11 +27,10 @@
 #include "catima/structures.h"
 #include "catima/calculations.h"
 #include "catima/material_database.h"
+#include "catima/storage.h"
 
 namespace catima{
     
-    class Interpolator;
-    
     /**
       * calculate dEdx for projectile-Material combination
       * @param p - Projectile

spline.cpp

@@ -0,0 +1,258 @@
+/*
+ * This is modification of Tino Kluge tk spline
+ * https://github.com/ttk592/spline/
+ *
+ *  the modification is in LU caclulation, 
+ * optimized for tridiagonal matrices
+ */
+#include "spline.h"
+
+
+namespace catima{
+
+band_matrix::band_matrix(int dim)
+{
+    resize(dim);
+}
+void band_matrix::resize(int dim)
+{
+    assert(dim>0);
+    a.resize(dim);
+    d.resize(dim);
+    c.resize(dim);
+}
+int band_matrix::dim() const
+{
+    return d.size();
+}
+
+
+// defines the new operator (), so that we can access the elements
+// by A(i,j), index going from i=0,...,dim()-1
+double & band_matrix::operator () (int i, int j)
+{
+    int k=j-i;       // what band is the entry
+    assert( (i>=0) && (i<dim()) && (j>=0) && (j<dim()) );
+    assert(k<2 && k>-2);
+    if(k>0)return c[i];
+    else if(k==0) return d[i];
+    else return a[i];
+}
+double band_matrix::operator () (int i, int j) const
+{
+    int k=j-i;       // what band is the entry
+    assert( (i>=0) && (i<dim()) && (j>=0) && (j<dim()) );
+    if(k>0)return c[i];
+    else if(k==0) return d[i];
+    else return a[i];
+}
+
+
+
+std::vector<double> band_matrix::trig_solve(const std::vector<double>& b) const
+{
+    assert( this->dim()==(int)b.size() );
+    std::vector<double> x(this->dim());
+    std::vector<double> g(this->dim());
+    int j_stop;
+    double sum;
+    
+    assert(d[0]!=0.0);
+    x[0] = b[0]/d[0];
+    double bet = d[0];
+    for(int j=1;j<this->dim();j++){
+        g[j] = c[j-1]/bet;
+        bet = d[j] - (a[j]*g[j]);
+        assert(bet != 0.0);
+        x[j] = (b[j]-a[j]*x[j-1])/bet;
+    }
+    for(int j=this->dim()-2;j>=0;j--){
+        x[j] -= g[j+1]*x[j+1];
+    }
+    
+    return x;
+}
+
+
+// spline implementation
+// -----------------------
+
+void spline::set_boundary(spline::bd_type left, double left_value,
+                          spline::bd_type right, double right_value,
+                          bool force_linear_extrapolation)
+{
+    assert(n==0);          // set_points() must not have happened yet
+    m_left=left;
+    m_right=right;
+    m_left_value=left_value;
+    m_right_value=right_value;
+    m_force_linear_extrapolation=force_linear_extrapolation;
+}
+
+
+void spline::set_points(const double *x,
+                        const double *y,
+                        const size_t num
+                        )
+{
+    assert(num>2);
+    m_x=x;
+    m_y=y;
+    n=num;
+    // TODO: maybe sort x and y, rather than returning an error
+    for(int i=0; i<n-1; i++) {
+        assert(m_x[i]<m_x[i+1]);
+    }
+
+    
+        // setting up the matrix and right hand side of the equation system
+        // for the parameters b[]
+        band_matrix A(n);
+        std::vector<double>  rhs(n);
+        for(int i=1; i<n-1; i++) {
+            A(i,i-1)=1.0/3.0*(x[i]-x[i-1]);
+            A(i,i)=2.0/3.0*(x[i+1]-x[i-1]);
+            A(i,i+1)=1.0/3.0*(x[i+1]-x[i]);
+            rhs[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
+        }
+        // boundary conditions
+        if(m_left == spline::bd_type::second_deriv) {
+            // 2*b[0] = f''
+            A(0,0)=2.0;
+            A(0,1)=0.0;
+            rhs[0]=m_left_value;
+        } else{
+            // c[0] = f', needs to be re-expressed in terms of b:
+            // (2b[0]+b[1])(x[1]-x[0]) = 3 ((y[1]-y[0])/(x[1]-x[0]) - f')
+            A(0,0)=2.0*(x[1]-x[0]);
+            A(0,1)=1.0*(x[1]-x[0]);
+            rhs[0]=3.0*((y[1]-y[0])/(x[1]-x[0])-m_left_value);
+        } 
+
+        if(m_right == spline::bd_type::second_deriv) {
+            // 2*b[n-1] = f''
+            A(n-1,n-1)=2.0;
+            A(n-1,n-2)=0.0;
+            rhs[n-1]=m_right_value;
+        } else{
+            // c[n-1] = f', needs to be re-expressed in terms of b:
+            // (b[n-2]+2b[n-1])(x[n-1]-x[n-2])
+            // = 3 (f' - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]))
+            A(n-1,n-1)=2.0*(x[n-1]-x[n-2]);
+            A(n-1,n-2)=1.0*(x[n-1]-x[n-2]);
+            rhs[n-1]=3.0*(m_right_value-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
+        }
+
+        // solve the equation system to obtain the parameters b[]
+        //m_b=A.lu_solve(rhs);
+        m_b=A.trig_solve(rhs);
+
+        // calculate parameters a[] and c[] based on b[]
+        m_a.resize(n);
+        m_c.resize(n);
+        for(int i=0; i<n-1; i++) {
+            m_a[i]=1.0/3.0*(m_b[i+1]-m_b[i])/(x[i+1]-x[i]);
+            m_c[i]=(y[i+1]-y[i])/(x[i+1]-x[i])
+                   - 1.0/3.0*(2.0*m_b[i]+m_b[i+1])*(x[i+1]-x[i]);
+        }
+    
+    // for left extrapolation coefficients
+    m_b0 = (m_force_linear_extrapolation==false) ? m_b[0] : 0.0;
+    m_c0 = m_c[0];
+
+    // for the right extrapolation coefficients
+    // f_{n-1}(x) = b*(x-x_{n-1})^2 + c*(x-x_{n-1}) + y_{n-1}
+    double h=x[n-1]-x[n-2];
+    // m_b[n-1] is determined by the boundary condition
+    m_a[n-1]=0.0;
+    m_c[n-1]=3.0*m_a[n-2]*h*h+2.0*m_b[n-2]*h+m_c[n-2];   // = f'_{n-2}(x_{n-1})
+    if(m_force_linear_extrapolation==true)
+        m_b[n-1]=0.0;
+}
+
+double spline::operator() (double x) const
+{
+    assert(n>2);
+    // find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
+    auto it=std::lower_bound(m_x,m_x+n,x);
+    //int idx=std::max( int(it-m_x)-1, 0);
+    if(it!=m_x)it--;
+    int idx = std::distance(m_x,it);
+    double mx = *it;
+    double h=x-mx;
+    double interpol;
+    if(x<m_x[0]) {
+        // extrapolation to the left
+        interpol=(m_b0*h + m_c0)*h + m_y[0];
+    } else if(x>m_x[n-1]) {
+        // extrapolation to the right
+        interpol=(m_b[n-1]*h + m_c[n-1])*h + m_y[n-1];
+    } else {
+        // interpolation
+        interpol=((m_a[idx]*h + m_b[idx])*h + m_c[idx])*h + m_y[idx];
+    }
+    return interpol;
+}
+
+double spline::deriv(int order, double x) const
+{
+    assert(order>0);
+    // find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
+    auto it=std::lower_bound(m_x,m_x+n,x);
+    //int idx=std::max( int(it-m_x)-1, 0);
+    if(it!=m_x)it--;
+    int idx = std::distance(m_x,it);
+    double mx = *it;
+    double h=x-mx;
+    double interpol;
+    if(x<m_x[0]) {
+        // extrapolation to the left
+        switch(order) {
+        case 1:
+            interpol=2.0*m_b0*h + m_c0;
+            break;
+        case 2:
+            interpol=2.0*m_b0*h;
+            break;
+        default:
+            interpol=0.0;
+            break;
+        }
+    } else if(x>m_x[n-1]) {
+        // extrapolation to the right
+        switch(order) {
+        case 1:
+            interpol=2.0*m_b[n-1]*h + m_c[n-1];
+            break;
+        case 2:
+            interpol=2.0*m_b[n-1];
+            break;
+        default:
+            interpol=0.0;
+            break;
+        }
+    } else {
+        // interpolation
+        switch(order) {
+        case 1:
+            interpol=(3.0*m_a[idx]*h + 2.0*m_b[idx])*h + m_c[idx];
+            break;
+        case 2:
+            interpol=6.0*m_a[idx]*h + 2.0*m_b[idx];
+            break;
+        case 3:
+            interpol=6.0*m_a[idx];
+            break;
+        default:
+            interpol=0.0;
+            break;
+        }
+    }
+    return interpol;
+}
+
+
+
+} // namespace tk
+
+

spline.h

@@ -0,0 +1,93 @@
+/*
+ * 
+ * This is modification of Tino Kluge tk spline
+ * calculation is optimized for tridiagonal matrices
+ * 
+ *  Copyright(C) 2017
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Affero General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU Affero General Public License for more details.
+
+ *  You should have received a copy of the GNU Affero General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef TK_SPLINE_H
+#define TK_SPLINE_H
+
+#include <cstdio>
+#include <cassert>
+#include <vector>
+#include <algorithm>
+
+namespace catima
+{
+
+// band matrix solver
+class band_matrix
+{
+private:
+    std::vector<double> a;
+    std::vector<double> d;
+    std::vector<double> c;
+    std::vector<double> save;
+public:
+    band_matrix() {};                             // constructor
+    band_matrix(int dim);       // constructor
+    ~band_matrix() {};                            // destructor
+    void resize(int dim);      // init with dim,n_u,n_l
+    int dim() const;                             // matrix dimension
+
+    // access operator
+    double & operator () (int i, int j);            // write
+    double   operator () (int i, int j) const;      // read
+    // we can store an additional diogonal (in m_lower)
+    std::vector<double> trig_solve(const std::vector<double>& b) const;
+};
+
+
+// spline interpolation
+class spline
+{
+public:
+    enum class bd_type {
+        first_deriv = 1,
+        second_deriv = 2
+    };
+
+private:
+    const double *m_x, *m_y;            // x,y coordinates of points
+    size_t n=0;
+    // interpolation parameters
+    // f(x) = a*(x-x_i)^3 + b*(x-x_i)^2 + c*(x-x_i) + y_i
+    std::vector<double> m_a,m_b,m_c;        // spline coefficients
+    double  m_b0, m_c0;                     // for left extrapol
+    bd_type m_left = bd_type::second_deriv;
+    bd_type m_right = bd_type::second_deriv;
+    double  m_left_value = 0.0;
+    double  m_right_value = 0.0;
+    bool    m_force_linear_extrapolation = false;
+
+public:
+    // set default boundary condition to be zero curvature at both ends
+    spline(){}
+
+    // optional, but if called it has to come be before set_points()
+    void set_boundary(bd_type left, double left_value,
+                      bd_type right, double right_value,
+                      bool force_linear_extrapolation=false);
+    void set_points(const double *x,
+                    const double *y,
+                    const size_t num);
+    double operator() (double x) const;
+    double deriv(int order, double x) const;
+};
+
+} // namespace tk
+
+#endif /* TK_SPLINE_H */

storage.cpp

@@ -1,6 +1,22 @@
+/*
+ *  Copyright(C) 2017
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Affero General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU Affero General Public License for more details.
+
+ *  You should have received a copy of the GNU Affero General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
 #include <math.h>
 #include <iostream>
 #include "storage.h"
+#include "catima/catima.h"
 namespace catima {
     Data _storage;
     EnergyTable<max_datapoints> energy_table(logEmin,logEmax);
@@ -58,7 +74,7 @@ DataPoint& Data::Get(const Projectile &p, const Material &t, const Config &c){
     }
 
 //////////// Interpolator ////////////////////////////////
-Interpolator::Interpolator(const double *x, const double *y, int num,interpolation_t type){
+InterpolatorGSL::InterpolatorGSL(const double *x, const double *y, int num,interpolation_t type){
     acc = gsl_interp_accel_alloc ();
 
     if(type==cspline)
@@ -71,7 +87,7 @@ Interpolator::Interpolator(const double *x, const double *y, int num,interpolati
     max= x[num-1];
 
 }
-Interpolator::Interpolator(const std::vector<double>& x, const std::vector<double>& y,interpolation_t type){
+InterpolatorGSL::InterpolatorGSL(const std::vector<double>& x, const std::vector<double>& y,interpolation_t type){
     //Interpolator(x.data(),y.data(),x.size());
     acc = gsl_interp_accel_alloc ();
     if(type==cspline)
@@ -84,21 +100,44 @@ Interpolator::Interpolator(const std::vector<double>& x, const std::vector<doubl
     max= x[x.size()-1];
 }
 
-Interpolator::~Interpolator(){
+InterpolatorGSL::~InterpolatorGSL(){
     gsl_interp_accel_free (acc);
     gsl_spline_free (spline);
 }
 
-double Interpolator::eval(double x){
+double InterpolatorGSL::eval(double x){
     if(x<min)x=min;
     if(x>max)x=max;
     return gsl_spline_eval(spline, x, acc);
 }
 
-double Interpolator::derivative(double x){
+double InterpolatorGSL::derivative(double x){
     if(x<min)x=min;
     if(x>max)x=max;
     return gsl_spline_eval_deriv (spline, x, acc);
 }
 
+
+//////////// Interpolator2 ////////////////////////////////
+#ifdef BUILTIN_SPLINE
+Interpolator2::Interpolator2(const double *x, const double *y, int num){
+    ss.set_points(x,y,num);
+    min= x[0];
+    max= x[num-1];
+
+}
+
+double Interpolator2::eval(double x){
+    if(x<min)x=min;
+    if(x>max)x=max;
+    return ss(x);
+}
+
+double Interpolator2::derivative(double x){
+    if(x<min)x=min;
+    if(x>max)x=max;
+    return ss.deriv(1,x);
+}
+#endif
+
 }

storage.h

@@ -19,10 +19,16 @@
 
 #include <vector>
 #include <iterator>
+#include <cmath>
 //#include <unordered_set>
 #include <gsl/gsl_spline.h>
 #include "catima/constants.h"
-#include "catima/catima.h"
+#include "catima/structures.h"
+#include "catima/config.h"
+
+#ifdef BUILTIN_SPLINE
+#include "catima/spline.h"
+#endif
 
 
 namespace catima{
@@ -49,19 +55,19 @@ namespace catima{
 
 	template<int N>
 	int EnergyTable_index(const EnergyTable<N> &table, double val){
-        double lxval = log(val)/M_LN10;
 		if(val<table.values[0] || val>table.values[table.num-1])return -1;
-		int i = (int)lxval/table.step;
+		double lxval = (log(val/table.values[0])/M_LN10);
+		int i = (int)std::floor(lxval/table.step);
 		return i;
 	}
 	
 	template<int N>
 	double EnergyTable_interpolate(const EnergyTable<N> &table, double xval, double *y){
 	    double r;
-	    double lxval = log(xval)/M_LN10;
 	    if(xval<table.values[0] || xval>table.values[table.num-1])return 0.0;
 	    if(xval==table.values[table.num-1])return y[table.num-1];
-	    int i = (int)(lxval/table.step);
+	    double lxval = (log(xval/table.values[0])/M_LN10);
+	    int i = (int)std::floor(lxval/table.step);
 	    double linstep = table.values[i+1] - table.values[i];
 	    double x = 1.0 - ((xval - table.values[i])/linstep);
 	    r = (x*y[i]) + ((1-x)*y[i+1]);
@@ -113,11 +119,11 @@ namespace catima{
     };
     
 /// Interpolation class, to store interpolated values
-	class Interpolator{
+	class InterpolatorGSL{
         public:
-        Interpolator(const double *x, const double *y, int num,interpolation_t type=cspline);
+        InterpolatorGSL(const double *x, const double *y, int num,interpolation_t type=cspline);
-        Interpolator(const std::vector<double>& x, const std::vector<double>& y,interpolation_t type=cspline);
+        InterpolatorGSL(const std::vector<double>& x, const std::vector<double>& y,interpolation_t type=cspline);
-        ~Interpolator();
+        ~InterpolatorGSL();
         double operator()(double x){return eval(x);};
         double eval(double x);
         double derivative(double x);
@@ -131,6 +137,22 @@ namespace catima{
         gsl_spline *spline;
     };
     
+#ifdef BUILTIN_SPLINE
+	class Interpolator2{
+        public:
+        Interpolator2(const double *x, const double *y, int num);
+        double operator()(double x){return eval(x);};
+        double eval(double x);
+        double derivative(double x);
+        double get_min(){return min;};
+        double get_max(){return max;};
+
+        private:
+        double min=0;
+        double max=0;
+        spline ss;
+    };
+#endif    
     extern Data _storage;
 	
 	inline DataPoint& get_data(const Projectile &p, const Material &t, const Config &c=default_config){
@@ -138,6 +160,9 @@ namespace catima{
 	}
 
     bool operator==(const DataPoint &a, const DataPoint &b);
+    
+    using InterpolatorLinear = InterpolatorGSL;
+    using Interpolator = InterpolatorGSL;
 }
 
 #endif

tests/test_storage.cpp

@@ -121,6 +121,19 @@ const lest::test specification[] =
       EXPECT(catima::energy_table.values[0]==exp(M_LN10*(catima::logEmin)));
       EXPECT(catima::energy_table.values[1]==exp(M_LN10*(catima::logEmin+step)));
       EXPECT(catima::energy_table.values[catima::max_datapoints-1]==approx(exp(M_LN10*(catima::logEmax))).epsilon(1e-6));
+    },
+    CASE("indexing"){
+      double step = catima::energy_table.step;
+      double val, dif;
+      
+      EXPECT(EnergyTable_index(catima::energy_table, 0.0)==-1);
+
+      for(int i:{5,10,100,498}){
+          val = catima::energy_table.values[i];
+          dif = catima::energy_table.values[i+1] - val;
+          EXPECT(EnergyTable_index(catima::energy_table, val)==i);
+          EXPECT(EnergyTable_index(catima::energy_table, val+0.5*dif)==i);
+      }
     }
 };