MADNESS 0.10.1
interpolation_1d.h
Go to the documentation of this file.
1/*
2 This file is part of MADNESS.
3
4 Copyright (C) 2007,2010 Oak Ridge National Laboratory
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19
20 For more information please contact:
21
22 Robert J. Harrison
23 Oak Ridge National Laboratory
24 One Bethel Valley Road
25 P.O. Box 2008, MS-6367
26
27 email: harrisonrj@ornl.gov
28 tel: 865-241-3937
29 fax: 865-572-0680
30
31 $Id$
32*/
33#ifndef MADNESS_MISC_INTERPOLATION_1D_H__INCLUDED
34#define MADNESS_MISC_INTERPOLATION_1D_H__INCLUDED
35
36#include <iostream>
37#include <cmath>
38#include <vector>
39
40#include "../world/world_task_queue.h"
41
42namespace madness {
43
44/*!
45 \file misc/interpolation_1d.h
46 \brief Provides 1D cubic interpolation class
47 \ingroup misc
48 */
49
50/// An class for 1-D data interpolation based on cubic polynomials.
51
52/// \ingroup misc
53/// Needs to be passed the endpoints of the interpolation: [lo,hi] and the
54/// number of grid points.
55///
56/// Two methods for generating the interpolation are presently supported:
57/// 1) Pass in a std::vector containing the y-points.
58/// 2) Pass in some object that provides an appropriate () operator, perhaps
59/// a function pointer.
60template <typename T> class CubicInterpolationTable {
61protected:
62 double lo; ///< Interpolation is in range [lo,hi]
63 double hi; ///< Interpolation is in range [lo,hi]
64 double h; ///< Grid spacing
65 double rh; ///< 1/h
66 int npt; ///< No. of grid points
67 std::vector<T> a; ///< (1+4)*npt vector of x and polynomial coefficients
68
69 // Cubic interp thru 4 points ... not good for noisy data
70 static void cubic_fit(const double *x, const T *f, T *a) {
71 double x0_2 = x[0] * x[0], x1_2 = x[1] * x[1], x2_2 = x[2] * x[2],
72 x3_2 = x[3] * x[3];
73 double x0_3 = x[0] * x[0] * x[0], x1_3 = x[1] * x[1] * x[1],
74 x2_3 = x[2] * x[2] * x[2], x3_3 = x[3] * x[3] * x[3];
75
76 a[0] = -(-x0_3 * x2_2 * x[3] * f[1] + x0_3 * x[2] * x3_2 * f[1] -
77 x0_3 * f[3] * x[2] * x1_2 + x0_3 * x[3] * f[2] * x1_2 +
78 x0_3 * f[3] * x2_2 * x[1] - x0_3 * x3_2 * f[2] * x[1] +
79 x0_2 * x1_3 * f[3] * x[2] - x0_2 * x1_3 * f[2] * x[3] +
80 x0_2 * x3_3 * f[2] * x[1] + x0_2 * f[1] * x2_3 * x[3] -
81 x0_2 * f[1] * x3_3 * x[2] - x0_2 * f[3] * x2_3 * x[1] +
82 x[0] * x3_2 * f[2] * x1_3 - x[0] * f[3] * x2_2 * x1_3 +
83 x[0] * x1_2 * f[3] * x2_3 - x[0] * x1_2 * f[2] * x3_3 -
84 x[0] * f[1] * x3_2 * x2_3 + x[0] * f[1] * x2_2 * x3_3 -
85 f[0] * x2_3 * x1_2 * x[3] + f[0] * x2_2 * x1_3 * x[3] +
86 f[0] * x3_2 * x2_3 * x[1] - f[0] * x3_3 * x2_2 * x[1] +
87 f[0] * x3_3 * x1_2 * x[2] - f[0] * x3_2 * x1_3 * x[2]) /
88 (-x2_2 * x[0] * x3_3 + x2_2 * x[0] * x1_3 - x0_2 * x[3] * x2_3 +
89 x0_2 * x3_3 * x[2] - x0_2 * x[1] * x3_3 + x0_2 * x[1] * x2_3 +
90 x0_2 * x1_3 * x[3] - x0_2 * x1_3 * x[2] + x3_2 * x[0] * x2_3 -
91 x3_2 * x[0] * x1_3 + x[3] * x2_3 * x1_2 - x3_2 * x2_3 * x[1] +
92 x3_3 * x2_2 * x[1] - x3_3 * x[2] * x1_2 + x[0] * x3_3 * x1_2 -
93 x[0] * x2_3 * x1_2 - x0_3 * x3_2 * x[2] - x0_3 * x[3] * x1_2 +
94 x0_3 * x3_2 * x[1] + x[2] * x3_2 * x1_3 - x2_2 * x[3] * x1_3 +
95 x0_3 * x2_2 * x[3] - x0_3 * x2_2 * x[1] + x0_3 * x[2] * x1_2);
96 a[1] = (-x2_3 * x1_2 * f[0] + x3_2 * x2_3 * f[0] + x2_3 * x0_2 * f[1] +
97 x1_2 * f[3] * x2_3 - x2_3 * x0_2 * f[3] - f[1] * x3_2 * x2_3 -
98 f[3] * x2_2 * x1_3 - x3_3 * x2_2 * f[0] + f[1] * x2_2 * x3_3 +
99 x2_2 * x1_3 * f[0] - f[1] * x2_2 * x0_3 + f[3] * x2_2 * x0_3 -
100 x1_3 * x3_2 * f[0] - x0_2 * x1_3 * f[2] - f[3] * x0_3 * x1_2 +
101 f[1] * x3_2 * x0_3 + x1_2 * f[2] * x0_3 + x3_3 * f[0] * x1_2 -
102 x3_2 * f[2] * x0_3 - f[1] * x3_3 * x0_2 + x0_2 * x3_3 * f[2] -
103 x1_2 * f[2] * x3_3 + x3_2 * f[2] * x1_3 + x0_2 * x1_3 * f[3]) /
104 (-x[3] + x[2]) /
105 (-x2_2 * x0_2 * x[3] - x2_2 * x[1] * x3_2 + x2_2 * x1_2 * x[3] +
106 x2_2 * x[0] * x3_2 - x2_2 * x[0] * x1_2 + x2_2 * x0_2 * x[1] +
107 x[2] * x[0] * x1_3 + x[2] * x0_3 * x[3] - x[2] * x0_3 * x[1] -
108 x[2] * x1_3 * x[3] - x[2] * x0_2 * x3_2 + x[2] * x1_2 * x3_2 -
109 x[2] * x[3] * x[0] * x1_2 + x[2] * x[3] * x0_2 * x[1] +
110 x0_3 * x1_2 - x0_2 * x1_3 + x[3] * x[0] * x1_3 -
111 x[3] * x0_3 * x[1] - x3_2 * x[0] * x1_2 + x3_2 * x0_2 * x[1]);
112 a[2] = -(-x1_3 * f[3] * x[2] + x1_3 * f[2] * x[3] + x1_3 * x[0] * f[3] +
113 x1_3 * f[0] * x[2] - x1_3 * x[0] * f[2] - x1_3 * f[0] * x[3] +
114 f[3] * x2_3 * x[1] - f[3] * x0_3 * x[1] - x[1] * x2_3 * f[0] +
115 x[1] * f[2] * x0_3 + x3_3 * f[0] * x[1] - x3_3 * f[2] * x[1] +
116 f[1] * x[3] * x0_3 - f[1] * x[0] * x3_3 - x3_3 * f[0] * x[2] +
117 x3_3 * x[0] * f[2] - f[2] * x0_3 * x[3] + x2_3 * f[0] * x[3] +
118 f[1] * x3_3 * x[2] - f[1] * x2_3 * x[3] + x2_3 * x[0] * f[1] -
119 x2_3 * x[0] * f[3] + f[3] * x0_3 * x[2] - x[2] * f[1] * x0_3) /
120 (x[3] * x2_2 - x3_2 * x[2] + x[1] * x3_2 - x[1] * x2_2 -
121 x1_2 * x[3] + x1_2 * x[2]) /
122 (-x[2] * x[1] * x[3] + x[1] * x[2] * x[0] - x0_2 * x[1] +
123 x[1] * x[3] * x[0] + x[2] * x[3] * x[0] + x0_3 - x0_2 * x[2] -
124 x0_2 * x[3]);
125 a[3] = (x[0] * f[3] * x1_2 - x0_2 * x[3] * f[2] + x2_2 * x[0] * f[1] +
126 x0_2 * f[3] * x[2] - x3_2 * f[2] * x[1] - f[0] * x3_2 * x[2] -
127 f[3] * x[2] * x1_2 - x2_2 * x[0] * f[3] - f[0] * x2_2 * x[1] +
128 f[3] * x2_2 * x[1] + x0_2 * f[1] * x[3] + x[2] * x3_2 * f[1] -
129 x0_2 * f[1] * x[2] + f[0] * x[2] * x1_2 + x[3] * f[2] * x1_2 +
130 f[0] * x3_2 * x[1] + x3_2 * x[0] * f[2] - x[0] * f[2] * x1_2 -
131 f[0] * x[3] * x1_2 - x0_2 * x[1] * f[3] + x0_2 * x[1] * f[2] +
132 f[0] * x2_2 * x[3] - x2_2 * x[3] * f[1] - x3_2 * x[0] * f[1]) /
133 (-x2_2 * x[0] * x3_3 + x2_2 * x[0] * x1_3 - x0_2 * x[3] * x2_3 +
134 x0_2 * x3_3 * x[2] - x0_2 * x[1] * x3_3 + x0_2 * x[1] * x2_3 +
135 x0_2 * x1_3 * x[3] - x0_2 * x1_3 * x[2] + x3_2 * x[0] * x2_3 -
136 x3_2 * x[0] * x1_3 + x[3] * x2_3 * x1_2 - x3_2 * x2_3 * x[1] +
137 x3_3 * x2_2 * x[1] - x3_3 * x[2] * x1_2 + x[0] * x3_3 * x1_2 -
138 x[0] * x2_3 * x1_2 - x0_3 * x3_2 * x[2] - x0_3 * x[3] * x1_2 +
139 x0_3 * x3_2 * x[1] + x[2] * x3_2 * x1_3 - x2_2 * x[3] * x1_3 +
140 x0_3 * x2_2 * x[3] - x0_3 * x2_2 * x[1] + x0_3 * x[2] * x1_2);
141 }
142
143 // Use the x- and y-points to make the interpolation
144 void make_interpolation(const std::vector<double> &x, const std::vector<T> &p,
145 const int npts_per_task = std::numeric_limits<int>::max(),
146 World* world_ptr = nullptr) {
147 const bool use_threads = npts_per_task < npt && world_ptr;
148
149 // Generate interior polynomial coeffs
150 const auto iend = npt - 2;
151 for (int i = 1; i < iend; i += npts_per_task) {
152 auto task = [istart = i, this, &x, &p, npts_per_task, iend]() {
153 const auto ifence = std::min(istart + npts_per_task, iend);
154 for (int i = istart; i < ifence; ++i) {
155 double mid = (x[i] + x[i + 1]) * 0.5;
156 double y[4] = {x[i - 1] - mid, x[i] - mid, x[i + 1] - mid,
157 x[i + 2] - mid};
158 this->a[i * 5] = mid;
159 cubic_fit(y, &p[i - 1], &this->a[i * 5 + 1]);
160 }
161 };
162 if (use_threads)
163 world_ptr->taskq.add(task);
164 else
165 task();
166 }
167 if (use_threads)
168 world_ptr->taskq.fence();
169
170 // Fixup end points
171 for (int j = 0; j < 5; ++j) {
172 a[j] = a[5 + j];
173 a[5 * npt - 5 + j] = a[5 * npt - 10 + j] = a[5 * npt - 15 + j];
174 }
175 }
176
177 /// constructs the interpolation table using optional tasking
178 /// \param world_ptr pointer to the World object whose local taskq to use for tasking; if null, will compute serially
179 /// \param lo the lower bound of the interpolation interval
180 /// \param up the upper bound of the interpolation interval
181 /// \param npt the number of interpolation points
182 /// \param[in] f a `T(T)` callable; should be reentrant if \p npt is less than \p min_npts_per_task
183 /// \param[in] min_npts_per_task if \p npt is greater than this and there is more than 1 thread will use tasks
184 /// \warning computes data locally even if \p world has more than 1 rank
185 template <typename functionT>
187 World* world_ptr, double lo, double hi, int npt, const functionT &f,
188 const int min_npts_per_task = std::numeric_limits<double>::max())
189 : lo(lo), hi(hi), h((hi - lo) / (npt - 1)), rh(1.0 / h), npt(npt),
190 a(npt * 5) {
191
192 // Evaluate the function to be interpolated
193 std::vector<T> p(npt);
194 std::vector<double> x(npt);
195 const int nthreads = 1 + ThreadPool::size();
196 const auto npts_per_task = std::max(min_npts_per_task, (npt + nthreads - 1)/ nthreads);
197 const auto use_threads = world_ptr && nthreads && npts_per_task < npt;
198 for (int i = 0; i < npt; i += npts_per_task) {
199 auto task = [istart = i, npts_per_task, npt, &x, &f, &p, this]() {
200 const auto ifence = std::min(istart + npts_per_task, npt);
201 for (int i = istart; i < ifence; ++i) {
202 x[i] = this->lo + i * this->h;
203 p[i] = f(x[i]);
204 }
205 };
206 if (use_threads) {
207 world_ptr->taskq.add(task);
208 }
209 else
210 task();
211 }
212 if (use_threads)
213 world_ptr->taskq.fence();
214
215 make_interpolation(x, p, npts_per_task, world_ptr);
216 }
217
218public:
219 static int min_npts_per_task_default; //!< minimum # of points per task
220
221 CubicInterpolationTable() : lo(0.0), hi(-1.0), h(0.0), rh(0.0), npt(0) {}
222
223 /// constructs the interpolation table serially
224 /// \param lo the lower bound of the interpolation interval
225 /// \param up the upper bound of the interpolation interval
226 /// \param npt the number of interpolation points
227 /// \param[in] f a `T(T)` callable; should be reentrant if \p npt is less than \p min_npts_per_task
228 /// \param[in] min_npts_per_task if \p npt is greater than this and there is more than 1 thread will use tasks
229 template <typename functionT>
231 double lo, double hi, int npt, const functionT &f)
232 : CubicInterpolationTable(nullptr, lo, hi, npt, f) {}
233
234 /// constructs the interpolation table using optional tasking
235 /// \param world the World object whose local taskq to use for tasking
236 /// \param lo the lower bound of the interpolation interval
237 /// \param up the upper bound of the interpolation interval
238 /// \param npt the number of interpolation points
239 /// \param[in] f a `T(T)` callable; should be reentrant if \p npt is less than \p min_npts_per_task
240 /// \param[in] min_npts_per_task if \p npt is greater than this and there is more than 1 thread will use tasks
241 /// \warning computes data locally even if \p world has more than 1 rank
242 template <typename functionT>
244 World& world, double lo, double hi, int npt, const functionT &f,
245 const int min_npts_per_task = min_npts_per_task_default)
246 : CubicInterpolationTable(&world, lo, hi, npt, f, min_npts_per_task) {}
247
248 CubicInterpolationTable(double lo, double hi, int npt,
249 const std::vector<T> &y)
250 : lo(lo), hi(hi), h((hi - lo) / (npt - 1)), rh(1.0 / h), npt(npt),
251 a(npt * 5) {
252
253 if ((int)y.size() < npt)
254 throw "Insufficient y-points";
255
256 std::vector<double> x(npt);
257 for (int i = 0; i < npt; ++i)
258 x[i] = lo + i * h;
259
260 make_interpolation(x, y);
261 }
262
263 T operator()(double y) const {
264 T y1;
265 int i = int((y - lo) * rh);
266 if (i < 0 || i >= npt)
267 throw "Out of range point";
268 i *= 5;
269 y1 = y - a[i];
270 T yy = y1 * y1;
271 return (a[i + 1] + y1 * a[i + 2]) + yy * (a[i + 3] + y1 * a[i + 4]);
272 }
273
274 template <typename functionT> double err(const functionT &f) const {
275 double maxabserr = 0.0;
276 double h7 = h / 7.0;
277 for (int i = 0; i < 7 * npt; ++i) {
278 double x = lo + h7 * i;
279 T fit = (*this)(x);
280 T exact = f(x);
281 maxabserr = std::max(fabs(fit - exact), maxabserr);
282 }
283 return maxabserr;
284 }
285
287};
288
289template <typename T>
291
292} // namespace madness
293
294#endif // MADNESS_MISC_INTERPOLATION_1D_H__INCLUDED
An class for 1-D data interpolation based on cubic polynomials.
Definition interpolation_1d.h:60
double h
Grid spacing.
Definition interpolation_1d.h:64
int npt
No. of grid points.
Definition interpolation_1d.h:66
CubicInterpolationTable(World &world, double lo, double hi, int npt, const functionT &f, const int min_npts_per_task=min_npts_per_task_default)
Definition interpolation_1d.h:243
CubicInterpolationTable(World *world_ptr, double lo, double hi, int npt, const functionT &f, const int min_npts_per_task=std::numeric_limits< double >::max())
Definition interpolation_1d.h:186
double rh
1/h
Definition interpolation_1d.h:65
std::vector< T > a
(1+4)*npt vector of x and polynomial coefficients
Definition interpolation_1d.h:67
CubicInterpolationTable(double lo, double hi, int npt, const std::vector< T > &y)
Definition interpolation_1d.h:248
double hi
Interpolation is in range [lo,hi].
Definition interpolation_1d.h:63
CubicInterpolationTable()
Definition interpolation_1d.h:221
void make_interpolation(const std::vector< double > &x, const std::vector< T > &p, const int npts_per_task=std::numeric_limits< int >::max(), World *world_ptr=nullptr)
Definition interpolation_1d.h:144
virtual ~CubicInterpolationTable()
Definition interpolation_1d.h:286
static int min_npts_per_task_default
minimum # of points per task
Definition interpolation_1d.h:219
static void cubic_fit(const double *x, const T *f, T *a)
Definition interpolation_1d.h:70
CubicInterpolationTable(double lo, double hi, int npt, const functionT &f)
Definition interpolation_1d.h:230
double lo
Interpolation is in range [lo,hi].
Definition interpolation_1d.h:62
T operator()(double y) const
Definition interpolation_1d.h:263
double err(const functionT &f) const
Definition interpolation_1d.h:274
static std::size_t size()
Returns the number of threads in the pool.
Definition thread.h:1413
void add(TaskInterface *t)
Add a new local task, taking ownership of the pointer.
Definition world_task_queue.h:466
void fence()
Returns after all local tasks have completed.
Definition world_task_queue.h:1384
A parallel world class.
Definition world.h:132
WorldTaskQueue & taskq
Task queue.
Definition world.h:204
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition derivatives.cc:72
auto T(World &world, response_space &f) -> response_space
Definition global_functions.cc:34
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
NDIM & f
Definition mra.h:2416
int task(int i)
Definition test_runtime.cpp:4
double(* exact)(double, double, double)
Definition testfuns.cc:6
std::vector< double > fit(size_t m, size_t n, const std::vector< double > N, const std::vector< double > &f)
Definition testfuns.cc:36