forked from plumed/plumed2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
CubicInterpolation.h
109 lines (91 loc) · 4.16 KB
/
CubicInterpolation.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Copyright (c) 2012-2016 The plumed team
(see the PEOPLE file at the root of the distribution for a list of names)
See http://www.plumed.org for more information.
This file is part of plumed, version 2.
plumed is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
plumed 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 Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with plumed. If not, see <http://www.gnu.org/licenses/>.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
#ifndef __PLUMED_tools_CubicInterpolation_h
#define __PLUMED_tools_CubicInterpolation_h
#include <vector>
#include "Matrix.h"
#include "core/Value.h"
namespace PLMD {
// Abstract base class for cubic interpolation
class CInterpolation {
private:
unsigned bold;
Matrix<double> splinepoints;
unsigned search1( const unsigned& kk, const double& x, const unsigned& jold ) const ;
protected:
std::vector<double> lb, ub;
std::vector<unsigned> np, stride;
unsigned findBox( const std::vector<double>& pos );
double getPointSpacing( const unsigned dir, const unsigned k ) const ;
double getCrossTermDenominator( const unsigned i, const unsigned j ) const ;
public:
CInterpolation( const std::vector<unsigned>& dd, const std::vector<double>& fmin, const std::vector<double>& fmax );
virtual ~CInterpolation();
unsigned getNumberOfSplinePoints() const ;
void getNumbersOfPoints( std::vector<unsigned>& nspline ) const ;
void getSplinePoint( const unsigned nn, std::vector<double>& pp ) const ;
void getGridBoundaries( std::vector<double>& gmin, std::vector<double>& gmax ) const ;
virtual void set_table( const std::vector<Value>& ff )=0;
virtual double get_fdf( const std::vector<double>& pos )=0;
};
inline
unsigned CInterpolation::getNumberOfSplinePoints() const {
return splinepoints.nrows();
}
inline
void CInterpolation::getSplinePoint( const unsigned nn, std::vector<double>& pp ) const {
plumed_dbg_assert( nn<splinepoints.nrows() && pp.size()==np.size() );
for(unsigned i=0;i<np.size();++i) pp[i]=splinepoints(nn,i);
}
inline
double CInterpolation::getPointSpacing( const unsigned dir, const unsigned k ) const {
unsigned i=k*stride[dir];
return splinepoints(i+stride[dir], dir) - splinepoints(i, dir);
}
inline
double CInterpolation::getCrossTermDenominator( const unsigned i, const unsigned j ) const {
plumed_dbg_assert( splinepoints.ncols()==2 );
unsigned iplus, iminus; iplus=(i+1)*stride[0]; iminus=(i-1)*stride[0];
return ( splinepoints(iplus,0) - splinepoints(iminus,0) ) * ( splinepoints(iplus+j+1,1) - splinepoints(iplus+j-1,1) );
}
inline
void CInterpolation::getGridBoundaries( std::vector<double>& gmin, std::vector<double>& gmax ) const {
getSplinePoint( 0, gmin ); getSplinePoint( splinepoints.nrows()-1, gmax );
}
class InterpolateCubic : public CInterpolation {
private:
std::vector<double> clist;
public:
InterpolateCubic( const std::vector<unsigned>& dd, const std::vector<double>& fmin, const std::vector<double>& fmax );
void set_table( const std::vector<Value>& ff );
double get_fdf( const std::vector<double>& pos );
};
class InterpolateBicubic : public CInterpolation {
private:
Matrix<int> wt;
std::vector<double> t1, t2;
Matrix<double> dcross;
std::vector<double> clist;
void IBicCoeff( const std::vector<double>& y, const std::vector<double>& dy1, const std::vector<double>& dy2,
const std::vector<double>& d2y12, const double& d1, const double& d2, Matrix<double>& c );
public:
InterpolateBicubic( const std::vector<unsigned>& dd, const std::vector<double>& fmin, const std::vector<double>& fmax );
void set_table( const std::vector<Value>& ff );
double get_fdf( const std::vector<double>& pos );
};
}
#endif