Skip to content

GitLab

Commit 39b2cc63 authored by Dipl.-Ing. Jonas Stienen's avatar Dipl.-Ing. Jonas Stienen
Browse files
Options

Removing SPLINE integration (GPLv2) and created ExternalLibs packaged with... 

Removing SPLINE integration (GPLv2) and created ExternalLibs packaged with same name using shared lib. Only required by ITAOps (solely required by RAVEN).
parent 0444f011
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 5 additions and 5591 deletions
CMakeLists.txt
View file @ 39b2cc63
......@@ -11,6 +11,7 @@ vista_use_package( sndfile )
vista_use_package( IPP QUIET )
vista_use_package( PCRE QUIET )
vista_use_package( SimpleIni QUIET )
vista_use_package( SPLINE QUIET )
if( NOT DEFINED ITA_BASE_WITH_FASTMATH_IPP )
set( ITA_BASE_WITH_FASTMATH_IPP OFF CACHE BOOL "Build with IPP implementation of fast math ops" )
......@@ -37,7 +38,7 @@ if( NOT DEFINED ITA_BASE_WITH_CONFIG_OLD_IMPL )
endif( )
if( NOT DEFINED ITA_BASE_WITH_OLD_RAVEN_OPS )
set( ITA_BASE_WITH_OLD_RAVEN_OPS OFF CACHE BOOL "Build with old ITAOps helper functions implementation (legacy code for RAVEN compatibility)" )
set( ITA_BASE_WITH_OLD_RAVEN_OPS OFF CACHE BOOL "Build with old ITAOps helper functions implementation (legacy code for RAVEN compatibility, adds GPL-licensed libs!)" )
endif( )
......@@ -76,7 +77,6 @@ set( ITABaseHeader
"include/ITATypes.h"
"include/ITAUncopyable.h"
"include/ITAWinPCClock.h"
"include/spline.h"
)
set( ITABaseSources
"src/ITAASCIITable.cpp"
......@@ -102,7 +102,6 @@ set( ITABaseSources
"src/ITASimpleConvolution.cpp"
"src/ITAStopWatch.cpp"
"src/ITAWinPCClock.cpp"
"src/spline.cpp"
)
......@@ -153,6 +152,9 @@ if( ITA_BASE_WITH_OLD_RAVEN_OPS )
if( NOT ITA_BASE_WITH_CONFIG_OLD_IMPL )
message( FATAL_ERROR "ITABase old RAVEN ops requires regular expressions and old config implementation. Please activate." )
endif( )
if( NOT VSPLINE_FOUND )
message( FATAL_ERROR "ITABase old RAVEN ops requires SPLINE, which was not found. Please provide." )
endif( )
set( ITABaseHeader "${ITABaseHeader}" "include/ITAOps.h" )
set( ITABaseSources "${ITABaseSources}" "src/ITAOps.cpp" )
endif( )
......
include/spline.h deleted 100644 → 0
View file @ 0444f011
// spline.h GPL!!!
#ifndef INCLUDE_WATCHER_SPLINE
#define INCLUDE_WATCHER_SPLINE
#include <ITABaseDefinitions.h>
float ITA_BASE_API basis_function_b_val ( float tdata[], float tval );
float ITA_BASE_API basis_function_beta_val(float beta1, float beta2, float tdata[], float tval );
ITA_BASE_API float *basis_matrix_b_uni(void);
ITA_BASE_API float *basis_matrix_beta_uni(float beta1, float beta2);
ITA_BASE_API float *basis_matrix_bezier(void);
ITA_BASE_API float *basis_matrix_hermite(void);
ITA_BASE_API float *basis_matrix_overhauser_nonuni(float alpha, float beta);
ITA_BASE_API float *basis_matrix_overhauser_nul(float alpha);
ITA_BASE_API float *basis_matrix_overhauser_nur(float beta);
ITA_BASE_API float *basis_matrix_overhauser_uni(void);
ITA_BASE_API float *basis_matrix_overhauser_uni_l(void);
ITA_BASE_API float *basis_matrix_overhauser_uni_r(void);
ITA_BASE_API float basis_matrix_tmp(int left, int n, float mbasis[], int ndata,
float tdata[], float ydata[], float tval );
ITA_BASE_API void bc_val(int n, float t, float xcon[], float ycon[], float *xval,
float *yval );
ITA_BASE_API float bez_val(int n, float x, float a, float b, float y[]);
ITA_BASE_API float bp_approx(int n, float a, float b, float ydata[], float xval);
ITA_BASE_API float *bp01(int n, float x);
ITA_BASE_API float *bpab(int n, float a, float b, float x);
ITA_BASE_API float d_max(float x, float y);
ITA_BASE_API float d_min(float x, float y);
ITA_BASE_API float d_random(float rlo, float rhi, int *seed);
ITA_BASE_API float d_uniform_01(int *seed);
ITA_BASE_API int d3_fs(float a1[], float a2[], float a3[], int n, float b[], float x[]);
ITA_BASE_API float *d3_mxv(int n, float a[], float x[]);
ITA_BASE_API float *d3_np_fs(int n, float a[], float b[]);
ITA_BASE_API void d3_print(int n, float a[], char *title);
ITA_BASE_API void d3_print_some(int n, float a[], int ilo, int jlo, int ihi, int jhi);
ITA_BASE_API float *d3_random(int n, int *seed);
ITA_BASE_API void data_to_dif(float diftab[], int ntab, float xtab[], float ytab[]);
ITA_BASE_API void dif_val(float diftab[], int ntab, float xtab[], float xval,
float *yval );
ITA_BASE_API void dvec_bracket(int n, float x[], float xval, int *left, int *right);
ITA_BASE_API void dvec_bracket3(int n, float t[], float tval, int *left);
ITA_BASE_API float *dvec_even(int n, float alo, float ahi);
ITA_BASE_API float *dvec_indicator(int n);
ITA_BASE_API void dvec_order_type(int n, float x[], int *order);
ITA_BASE_API void dvec_print(int n, float a[], char *title);
ITA_BASE_API float *dvec_random(int n, float alo, float ahi, int *seed);
ITA_BASE_API void dvec_sort_bubble_a(int n, float a[]);
ITA_BASE_API int i_max(int i1, int i2);
ITA_BASE_API int i_min(int i1, int i2);
ITA_BASE_API void least_set(int ntab, float xtab[], float ytab[], int ndeg,
float ptab[], float b[], float c[], float d[], float *eps, int *ierror );
ITA_BASE_API float least_val(float x, int ndeg, float b[], float c[], float d[]);
ITA_BASE_API void parabola_val2(int ndim, int ndata, float tdata[], float ydata[],
int left, float tval, float yval[] );
ITA_BASE_API int s_len_trim(char* s);
ITA_BASE_API float spline_b_val(int ndata, float tdata[], float ydata[], float tval);
ITA_BASE_API float spline_beta_val(float beta1, float beta2, int ndata, float tdata[],
float ydata[], float tval );
ITA_BASE_API float spline_constant_val(int ndata, float tdata[], float ydata[], float tval);
ITA_BASE_API float* spline_cubic_set(
const int n,
const float t[],
const float y[],
const int ibcbeg,
const float ybcbeg,
const int ibcend,
const float ybcend );
ITA_BASE_API float spline_cubic_val(
const int n,
const float t[],
const float tval,
const float y[],
const float ypp[],
float *ypval,
float *yppval );
ITA_BASE_API float spline_cubic_val(
const int n,
const float t[],
const float tval,
const float y[],
const float ypp[]);
ITA_BASE_API void spline_cubic_val2(int n, float t[], float tval, int *left, float y[],
float ypp[], float *yval, float *ypval, float *yppval );
ITA_BASE_API float *spline_hermite_set(int ndata, float tdata[], float ydata[],
float ypdata[] );
ITA_BASE_API void spline_hermite_val(int ndata, float tdata[], float c[], float tval,
float *sval, float *spval );
ITA_BASE_API float spline_linear_int(int ndata, float tdata[], float ydata[], float a, float b);
ITA_BASE_API void spline_linear_intset(int int_n, float int_x[], float int_v[],
float data_x[], float data_y[] );
ITA_BASE_API void spline_linear_val(int ndata, float tdata[], float ydata[],
float tval, float *yval, float *ypval );
ITA_BASE_API float spline_overhauser_nonuni_val(int ndata, float tdata[],
float ydata[], float tval );
ITA_BASE_API float spline_overhauser_uni_val(int ndata, float tdata[], float ydata[],
float tval );
ITA_BASE_API void spline_overhauser_val(int ndim, int ndata, float tdata[], float ydata[],
float tval, float yval[] );
ITA_BASE_API void spline_quadratic_val(int ndata, float tdata[], float ydata[],
float tval, float *yval, float *ypval );
ITA_BASE_API void timestamp(void);
#endif // INCLUDE_WATCHER_SPLINE
\ No newline at end of file
include/spline_tk.h deleted 100644 → 0
View file @ 0444f011
/*
* spline.h
*
* simple cubic spline interpolation library without external
* dependencies
*
* ---------------------------------------------------------------------
* Copyright (C) 2011, 2014 Tino Kluge (ttk448 at gmail.com)
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* 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 General Public License for more details.
*
* You should have received a copy of the GNU 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>
// unnamed namespace only because the implementation is in this
// header file and we don't want to export symbols to the obj files
namespace
{
namespace tk
{
// band matrix solver
class band_matrix
{
private:
std::vector< std::vector<double> > m_upper; // upper band
std::vector< std::vector<double> > m_lower; // lower band
public:
band_matrix() {}; // constructor
band_matrix(int dim, int n_u, int n_l); // constructor
~band_matrix() {}; // destructor
void resize(int dim, int n_u, int n_l); // init with dim,n_u,n_l
int dim() const; // matrix dimension
int num_upper() const
{
return m_upper.size()-1;
}
int num_lower() const
{
return m_lower.size()-1;
}
// 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)
double& saved_diag(int i);
double saved_diag(int i) const;
void lu_decompose();
std::vector<double> r_solve(const std::vector<double>& b) const;
std::vector<double> l_solve(const std::vector<double>& b) const;
std::vector<double> lu_solve(const std::vector<double>& b,
bool is_lu_decomposed=false);
};
// spline interpolation
class spline
{
public:
enum bd_type {
first_deriv = 1,
second_deriv = 2
};
private:
std::vector<double> m_x,m_y; // x,y coordinates of points
// 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, m_right;
double m_left_value, m_right_value;
bool m_force_linear_extrapolation;
public:
// set default boundary condition to be zero curvature at both ends
spline(): m_left(second_deriv), m_right(second_deriv),
m_left_value(0.0), m_right_value(0.0),
m_force_linear_extrapolation(false)
{
;
}
// 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 std::vector<double>& x,
const std::vector<double>& y, bool cubic_spline=true);
double operator() (double x) const;
};
// ---------------------------------------------------------------------
// implementation part, which could be separated into a cpp file
// ---------------------------------------------------------------------
// band_matrix implementation
// -------------------------
band_matrix::band_matrix(int dim, int n_u, int n_l)
{
resize(dim, n_u, n_l);
}
void band_matrix::resize(int dim, int n_u, int n_l)
{
assert(dim>0);
assert(n_u>=0);
assert(n_l>=0);
m_upper.resize(n_u+1);
m_lower.resize(n_l+1);
for(size_t i=0; i<m_upper.size(); i++) {
m_upper[i].resize(dim);
}
for(size_t i=0; i<m_lower.size(); i++) {
m_lower[i].resize(dim);
}
}
int band_matrix::dim() const
{
if(m_upper.size()>0) {
return m_upper[0].size();
} else {
return 0;
}
}
// 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( (-num_lower()<=k) && (k<=num_upper()) );
// k=0 -> diogonal, k<0 lower left part, k>0 upper right part
if(k>=0) return m_upper[k][i];
else return m_lower[-k][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()) );
assert( (-num_lower()<=k) && (k<=num_upper()) );
// k=0 -> diogonal, k<0 lower left part, k>0 upper right part
if(k>=0) return m_upper[k][i];
else return m_lower[-k][i];
}
// second diag (used in LU decomposition), saved in m_lower
double band_matrix::saved_diag(int i) const
{
assert( (i>=0) && (i<dim()) );
return m_lower[0][i];
}
double & band_matrix::saved_diag(int i)
{
assert( (i>=0) && (i<dim()) );
return m_lower[0][i];
}
// LR-Decomposition of a band matrix
void band_matrix::lu_decompose()
{
int i_max,j_max;
int j_min;
double x;
// preconditioning
// normalize column i so that a_ii=1
for(int i=0; i<this->dim(); i++) {
assert(this->operator()(i,i)!=0.0);
this->saved_diag(i)=1.0/this->operator()(i,i);
j_min=std::max(0,i-this->num_lower());
j_max=std::min(this->dim()-1,i+this->num_upper());
for(int j=j_min; j<=j_max; j++) {
this->operator()(i,j) *= this->saved_diag(i);
}
this->operator()(i,i)=1.0; // prevents rounding errors
}
// Gauss LR-Decomposition
for(int k=0; k<this->dim(); k++) {
i_max=std::min(this->dim()-1,k+this->num_lower()); // num_lower not a mistake!
for(int i=k+1; i<=i_max; i++) {
assert(this->operator()(k,k)!=0.0);
x=-this->operator()(i,k)/this->operator()(k,k);
this->operator()(i,k)=-x; // assembly part of L
j_max=std::min(this->dim()-1,k+this->num_upper());
for(int j=k+1; j<=j_max; j++) {
// assembly part of R
this->operator()(i,j)=this->operator()(i,j)+x*this->operator()(k,j);
}
}
}
}
// solves Ly=b
std::vector<double> band_matrix::l_solve(const std::vector<double>& b) const
{
assert( this->dim()==(int)b.size() );
std::vector<double> x(this->dim());
int j_start;
double sum;
for(int i=0; i<this->dim(); i++) {
sum=0;
j_start=std::max(0,i-this->num_lower());
for(int j=j_start; j<i; j++) sum += this->operator()(i,j)*x[j];
x[i]=(b[i]*this->saved_diag(i)) - sum;
}
return x;
}
// solves Rx=y
std::vector<double> band_matrix::r_solve(const std::vector<double>& b) const
{
assert( this->dim()==(int)b.size() );
std::vector<double> x(this->dim());
int j_stop;
double sum;
for(int i=this->dim()-1; i>=0; i--) {
sum=0;
j_stop=std::min(this->dim()-1,i+this->num_upper());
for(int j=i+1; j<=j_stop; j++) sum += this->operator()(i,j)*x[j];
x[i]=( b[i] - sum ) / this->operator()(i,i);
}
return x;
}
std::vector<double> band_matrix::lu_solve(const std::vector<double>& b,
bool is_lu_decomposed)
{
assert( this->dim()==(int)b.size() );
std::vector<double> x,y;
if(is_lu_decomposed==false) {
this->lu_decompose();
}
y=this->l_solve(b);
x=this->r_solve(y);
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(m_x.size()==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 std::vector<double>& x,
const std::vector<double>& y, bool cubic_spline)
{
assert(x.size()==y.size());
assert(x.size()>2);
m_x=x;
m_y=y;
int n=x.size();
// 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]);
}
if(cubic_spline==true) { // cubic spline interpolation
// setting up the matrix and right hand side of the equation system
// for the parameters b[]
band_matrix A(n,1,1);
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::second_deriv) {
// 2*b[0] = f''
A(0,0)=2.0;
A(0,1)=0.0;
rhs[0]=m_left_value;
} else if(m_left == spline::first_deriv) {
// 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);
} else {
assert(false);
}
if(m_right == spline::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 if(m_right == spline::first_deriv) {
// 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]));
} else {
assert(false);
}
// solve the equation system to obtain the parameters b[]
m_b=A.lu_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]);
}
} else { // linear interpolation
m_a.resize(n);
m_b.resize(n);
m_c.resize(n);
for(int i=0; i<n-1; i++) {
m_a[i]=0.0;
m_b[i]=0.0;
m_c[i]=(m_y[i+1]-m_y[i])/(m_x[i+1]-m_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
{
size_t n=m_x.size();
// find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
std::vector<double>::const_iterator it;
it=std::lower_bound(m_x.begin(),m_x.end(),x);
int idx=std::max( int(it-m_x.begin())-1, 0);
double h=x-m_x[idx];
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;
}
} // namespace tk
} // namespace