Line data Source code
1 : /** \file gslInterpolator.hpp
2 : * \brief Class for managing 1-D interpolation using the GNU Scientific Library
3 : *
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : * \ingroup gen_math_files
7 : *
8 : */
9 :
10 : //***********************************************************************//
11 : // Copyright 2015-2022 Jared R. Males (jaredmales@gmail.com)
12 : //
13 : // This file is part of mxlib.
14 : //
15 : // mxlib is free software: you can redistribute it and/or modify
16 : // it under the terms of the GNU General Public License as published by
17 : // the Free Software Foundation, either version 3 of the License, or
18 : // (at your option) any later version.
19 : //
20 : // mxlib is distributed in the hope that it will be useful,
21 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
22 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 : // GNU General Public License for more details.
24 : //
25 : // You should have received a copy of the GNU General Public License
26 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
27 : //***********************************************************************//
28 :
29 : #ifndef gslInterpolator_hpp
30 : #define gslInterpolator_hpp
31 :
32 : #include <vector>
33 : #include <gsl/gsl_interp.h>
34 : #include <gsl/gsl_errno.h>
35 :
36 : #include "../mxlib.hpp"
37 :
38 : namespace mx
39 : {
40 : namespace math
41 : {
42 :
43 : /// GSL Linear Interpolation
44 : /** \ingroup interpolation
45 : */
46 : template <typename _realT>
47 : struct gsl_interp_linear
48 : {
49 : typedef _realT realT;
50 :
51 10 : static const gsl_interp_type *interpolator()
52 : {
53 10 : return ::gsl_interp_linear;
54 : }
55 : };
56 :
57 : #ifndef MX_OLD_GSL
58 : /// GSL Steffen Interpolation
59 : /** \ingroup interpolation
60 : */
61 : template <typename _realT>
62 : struct gsl_interp_steffen
63 : {
64 : typedef _realT realT;
65 :
66 0 : static const gsl_interp_type *interpolator()
67 : {
68 0 : return ::gsl_interp_steffen;
69 : }
70 : };
71 : #endif
72 :
73 : /// Class to manage interpolation using the GSL interpolation library.
74 : /**
75 : *
76 : * \tparam interpT is the interpolation type, which also specifies the floating point precision.
77 : *
78 : * \ingroup interpolation
79 : */
80 : template <typename interpT>
81 : class gslInterpolator
82 : {
83 : public:
84 : typedef typename interpT::realT realT;
85 :
86 : static_assert( std::is_same<double, typename std::remove_cv<realT>::type>::value,
87 : "GSL Interpolation only works with double" );
88 :
89 : protected:
90 : gsl_interp *m_interp{ nullptr }; ///< the gsl interpolator structure
91 : gsl_interp_accel *m_acc{ nullptr }; ///< the gsl interpolation accelerator structure
92 :
93 : realT *m_xin; ///< the input x data
94 : realT *m_yin; ///< the input y data
95 :
96 : public:
97 : /// Default constructor
98 : gslInterpolator();
99 :
100 : /// Raw pointer constructor
101 : gslInterpolator( realT *xin, ///< [in] the input x data, this pointer is stored and must remain valid
102 : realT *yin, ///< [in] the input y data, this pointer is stored and must remain valid
103 : size_t Nin ///< [in] the size of data vectors
104 : );
105 :
106 : /// Vector constructor
107 : gslInterpolator( std::vector<realT> &xin, ///< [in] the input x data, the pointer to xin.data() is stored and xin
108 : ///< must remain unchanged
109 : std::vector<realT> &yin ///< [in] the input y data, the pointer to yin.data() is stored and yin
110 : ///< must remain unchanged
111 : );
112 :
113 : /// Destructor
114 : /** Deallocates working memory.
115 : */
116 : ~gslInterpolator();
117 :
118 : /// Setup the interpolator for the supplied data pointers
119 : /**
120 : * \throws mx::err::allocerr if allocation fails
121 : * \throws mx::err::liberr if GSL initialization fails
122 : */
123 : void setup( realT *xin, ///< [in] the input x data, this pointer is stored and must remain valid
124 : realT *yin, ///< [in] the input y data, this pointer is stored and must remain valid
125 : size_t Nin ///< [in] the size of data vectors
126 : );
127 :
128 : /// Setup the interpolator for the supplied data vectors
129 : /**
130 : * \throws mx::err::sizeerr if the vectors are not the same size
131 : * \throws mx::err::allocerr if allocation fails
132 : * \throws mx::err::liberr if GSL initialization fails
133 : */
134 : void setup( std::vector<realT>
135 : &xin, ///< [in] the input x data, the pointer to xin.data() is stored and xin must remain unchanged
136 : std::vector<realT>
137 : &yin ///< [in] the input y data, the pointer to yin.data() is stored and yin must remain unchanged
138 : );
139 :
140 : /// Calculate the interpolated function value at a point
141 : /**
142 : * \returns the interpolated value at \p x
143 : */
144 : realT operator()( const realT &x /**< [in] the point at which to interpolate */ );
145 : };
146 :
147 : template <typename interpT>
148 0 : gslInterpolator<interpT>::gslInterpolator()
149 : {
150 0 : }
151 :
152 : template <typename interpT>
153 0 : gslInterpolator<interpT>::gslInterpolator( realT *xin, realT *yin, size_t Nin )
154 : {
155 0 : setup( xin, yin, Nin );
156 0 : }
157 :
158 : template <typename interpT>
159 10 : gslInterpolator<interpT>::gslInterpolator( std::vector<realT> &xin, std::vector<realT> &yin )
160 : {
161 10 : setup( xin.data(), yin.data(), xin.size() );
162 10 : }
163 :
164 : template <typename interpT>
165 10 : gslInterpolator<interpT>::~gslInterpolator()
166 : {
167 10 : if( m_interp != nullptr )
168 10 : gsl_interp_free( m_interp );
169 10 : if( m_acc != nullptr )
170 10 : gsl_interp_accel_free( m_acc );
171 10 : }
172 :
173 : template <typename interpT>
174 10 : void gslInterpolator<interpT>::setup( realT *xin, realT *yin, size_t Nin )
175 : {
176 10 : if( m_interp != nullptr )
177 0 : gsl_interp_free( m_interp );
178 10 : if( m_acc != nullptr )
179 0 : gsl_interp_accel_free( m_acc );
180 :
181 10 : m_interp = gsl_interp_alloc( interpT::interpolator(), Nin );
182 10 : if( !m_interp )
183 : {
184 0 : throw(mx::exception(error_t::allocerr, "gsl_interp_alloc failed"));
185 : }
186 :
187 10 : m_acc = gsl_interp_accel_alloc();
188 10 : if( !m_acc )
189 : {
190 0 : throw(mx::exception(error_t::allocerr, "gsl_interp_accel_alloc failed" ));
191 : }
192 :
193 10 : int errv = gsl_interp_init( m_interp, xin, yin, Nin );
194 10 : if( errv != 0 )
195 : {
196 0 : throw( mx::exception( error_t::liberr, std::string( "gsl_interp_init failed: " ) + gsl_strerror( errv ) ) );
197 : }
198 :
199 10 : errv = gsl_interp_accel_reset( m_acc );
200 10 : if( errv != 0 )
201 : {
202 0 : throw(
203 0 : mx::exception( error_t::liberr, std::string( "gsl_interp_accel_reset failed: " ) + gsl_strerror( errv ) ) );
204 : }
205 10 : m_xin = xin;
206 10 : m_yin = yin;
207 10 : }
208 :
209 : template <typename interpT>
210 0 : void gslInterpolator<interpT>::setup( std::vector<realT> &xin, std::vector<realT> &yin )
211 : {
212 0 : if( xin.size() != yin.size() )
213 : {
214 0 : throw(mx::exception(error_t::sizeerr,"input vectors must be the same size" ));
215 : }
216 0 : setup( xin.data(), yin.data(), xin.size() );
217 0 : }
218 :
219 : template <typename interpT>
220 543 : typename interpT::realT gslInterpolator<interpT>::operator()( const realT &x )
221 : {
222 : realT y;
223 543 : int errv = gsl_interp_eval_e( m_interp, m_xin, m_yin, x, m_acc, &y );
224 543 : if( errv != 0 && errv != GSL_EDOM )
225 : {
226 0 : throw(mx::exception(error_t::liberr,"gsl_interp_eval_e failed" ));
227 : }
228 543 : return y;
229 : }
230 :
231 : extern template class gslInterpolator<gsl_interp_linear<double>>;
232 : extern template class gslInterpolator<gsl_interp_steffen<double>>;
233 :
234 : } // namespace math
235 : } // namespace mx
236 :
237 : #endif // gslInterpolator_hpp
|
