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 int npt; ///< No. of grid points
65 std::vector<T> a; ///< (1+4)*npt vector of x and polynomial coefficients
66 /// Below variables meaningful if user does specifies uniform spacing
67 double h; ///< Grid spacing.
68 double rh; ///< 1/h
69 /// Below variables meaningful if user provides grid points
70 std::vector<double> pts_; /// Grid points. Only stored if not evenly spaced. Empty otherwise.
71
72 // Cubic interp thru 4 points ... not good for noisy data
73 static void cubic_fit(const double *x, const T *f, T *a) {
74
75 const auto base0 = f[0] / ((x[0] - x[1]) * (x[0] - x[2]) * (x[0] - x[3]));
76 const auto base1 = f[1] / ((x[1] - x[0]) * (x[1] - x[2]) * (x[1] - x[3]));
77 const auto base2 = f[2] / ((x[2] - x[0]) * (x[2] - x[1]) * (x[2] - x[3]));
78 const auto base3 = f[3] / ((x[3] - x[0]) * (x[3] - x[1]) * (x[3] - x[2]));
79 a[3] = base0 + base1 + base2 + base3;
80 auto temp0 = -base0 * (x[1] + x[2] + x[3]);
81 auto temp1 = -base1 * (x[0] + x[2] + x[3]);
82 auto temp2 = -base2 * (x[0] + x[1] + x[3]);
83 auto temp3 = -base3 * (x[0] + x[1] + x[2]);
84 a[2] = temp0 + temp1 + temp2 + temp3;
85 temp0 = base0 * (x[1] * x[2] + x[2] * x[3] + x[3] * x[1]);
86 temp1 = base1 * (x[0] * x[2] + x[2] * x[3] + x[3] * x[0]);
87 temp2 = base2 * (x[0] * x[1] + x[1] * x[3] + x[3] * x[0]);
88 temp3 = base3 * (x[0] * x[1] + x[1] * x[2] + x[2] * x[0]);
89 a[1] = temp0 + temp1 + temp2 + temp3;
90 temp0 = -base0 * (x[1] * x[2] * x[3]);
91 temp1 = -base1 * (x[0] * x[2] * x[3]);
92 temp2 = -base2 * (x[0] * x[1] * x[3]);
93 temp3 = -base3 * (x[0] * x[1] * x[2]);
94 a[0] = temp0 + temp1 + temp2 + temp3;
95 }
96
97 // Use the x- and y-points to make the interpolation
98 void make_interpolation(const std::vector<double> &x, const std::vector<T> &p,
99 const int npts_per_task = std::numeric_limits<int>::max() - 1,
100 World* world_ptr = nullptr) {
101 const bool use_threads = npts_per_task < npt && world_ptr;
102
103 // Generate interior polynomial coeffs
104 const auto iend = npt - 2;
105 for (int i = 1; i < iend; i += npts_per_task) {
106 auto task = [istart = i, this, &x, &p, npts_per_task, iend]() {
107 const auto ifence = std::min(istart + npts_per_task, iend);
108 for (int i = istart; i < ifence; ++i) {
109 // Center x points for numerical stability
110 double mid = (x[i] + x[i + 1]) * 0.5;
111 double y[4] = {x[i - 1] - mid, x[i] - mid, x[i + 1] - mid,
112 x[i + 2] - mid};
113 this->a[i * 5] = mid;
114 cubic_fit(y, &p[i - 1], &this->a[i * 5 + 1]);
115 }
116 };
117 if (use_threads)
118 world_ptr->taskq.add(task);
119 else
120 task();
121 }
122 if (use_threads)
123 world_ptr->taskq.fence();
124
125 // Fixup end points
126 for (int j = 0; j < 5; ++j) {
127 a[j] = a[5 + j];
128 a[5 * npt - 5 + j] = a[5 * npt - 10 + j] = a[5 * npt - 15 + j];
129 }
130 }
131
132 /// constructs the interpolation table using optional tasking
133 /// \param world_ptr pointer to the World object whose local taskq to use for tasking; if null, will compute serially
134 /// \param lo the lower bound of the interpolation interval
135 /// \param up the upper bound of the interpolation interval
136 /// \param npt the number of interpolation points
137 /// \param[in] f a `T(T)` callable; should be reentrant if \p npt is less than \p min_npts_per_task
138 /// \param[in] min_npts_per_task if \p npt is greater than this and there is more than 1 thread will use tasks
139 /// \warning computes data locally even if \p world has more than 1 rank
140 template <typename functionT>
142 World* world_ptr, double lo, double hi, int npt, const functionT &f,
143 const int min_npts_per_task = std::numeric_limits<double>::max())
144 : lo(lo), hi(hi), h((hi - lo) / (npt - 1)), rh(1.0 / h), npt(npt),
145 a(npt * 5) {
146
147 // Evaluate the function to be interpolated
148 std::vector<T> p(npt);
149 std::vector<double> x(npt);
150 const int nthreads = 1 + ThreadPool::size();
151 const auto npts_per_task = std::max(min_npts_per_task, (npt + nthreads - 1)/ nthreads);
152 const auto use_threads = world_ptr && nthreads && npts_per_task < npt;
153 for (int i = 0; i < npt; i += npts_per_task) {
154 auto task = [istart = i, npts_per_task, npt, &x, &f, &p, this]() {
155 const auto ifence = std::min(istart + npts_per_task, npt);
156 for (int i = istart; i < ifence; ++i) {
157 x[i] = this->lo + i * this->h;
158 p[i] = f(x[i]);
159 }
160 };
161 if (use_threads) {
162 world_ptr->taskq.add(task);
163 }
164 else
165 task();
166 }
167 if (use_threads)
168 world_ptr->taskq.fence();
169
171 }
172
173public:
174 static int min_npts_per_task_default; //!< minimum # of points per task
175
176 CubicInterpolationTable() : lo(0.0), hi(-1.0), h(0.0), rh(0.0), npt(0) {}
177
178 /// constructs the interpolation table 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 template <typename functionT>
186 double lo, double hi, int npt, const functionT &f)
188
189 /// constructs the interpolation table using optional tasking
190 /// \param world the World object whose local taskq to use for tasking
191 /// \param lo the lower bound of the interpolation interval
192 /// \param up the upper bound of the interpolation interval
193 /// \param npt the number of interpolation points
194 /// \param[in] f a `T(T)` callable; should be reentrant if \p npt is less than \p min_npts_per_task
195 /// \param[in] min_npts_per_task if \p npt is greater than this and there is more than 1 thread will use tasks
196 /// \warning computes data locally even if \p world has more than 1 rank
197 template <typename functionT>
199 World& world, double lo, double hi, int npt, const functionT &f,
202
203 CubicInterpolationTable(std::vector<T> x, std::vector<T> y)
204 : npt(static_cast<int>(x.size())), a(npt * 5) {
205
206 if (x.size() != y.size()) {
207 throw std::runtime_error("Sizes of input and output arrays do not equal.");
208 }
209 // Out of paranoia, re-sort the data.
210 std::vector<int> indices(x.size());
211 std::vector<T> sorted_y(x.size());
212 std::iota(indices.begin(), indices.end(), 0);
213 std::stable_sort(indices.begin(), indices.end(), [&x](int i, int j) { return x[i] < x[j]; });
214 for (size_t i = 0; i < x.size(); i++) {
215 sorted_y[i] = y[indices[i]];
216 }
217 std::stable_sort(x.begin(), x.end());
218 lo = x[0];
219 hi = x[x.size() - 1];
220 pts_ = x;
221
223 }
224
225 CubicInterpolationTable(double lo, double hi, int npt,
226 const std::vector<T> &y)
227 : lo(lo), hi(hi), h((hi - lo) / (npt - 1)), rh(1.0 / h), npt(npt),
228 a(npt * 5) {
229
230 if ((int)y.size() < npt)
231 throw "Insufficient y-points";
232
233 std::vector<double> x(npt);
234 for (int i = 0; i < npt; ++i)
235 x[i] = lo + i * h;
236
237 make_interpolation(x, y);
238 }
239
240 T operator()(double y) const {
241 int i;
242
243 if (!pts_.empty()) {
244 i = std::lower_bound(pts_.begin(), pts_.end(), y) - pts_.begin();
245 if (y < lo || y > hi) throw std::runtime_error("Out of range point");
246 } else {
247 i = static_cast<int>((y - lo) * rh);
248 if (i < 0 || i >= npt) throw std::runtime_error("Out of range point");
249 }
250 i *= 5;
251 T y1 = y - a[i];
252 T yy = y1 * y1;
253 return (a[i + 1] + y1 * a[i + 2]) + yy * (a[i + 3] + y1 * a[i + 4]);
254 }
255
256 double get_lo() const { return lo; }
257
258 double get_hi() const { return hi; }
259
260 template <typename functionT> double err(const functionT &f) const {
261 if (! pts_.empty()) {
262 throw "This error test is only meaningful if f generated the interpolation.";
263 }
264 double maxabserr = 0.0;
265 double h7 = h / 7.0;
266 for (int i = 0; i < 7 * npt; ++i) {
267 double x = lo + h7 * i;
268 T fit = (*this)(x);
269 T exact = f(x);
270 maxabserr = std::max(fabs(fit - exact), maxabserr);
271 }
272 return maxabserr;
273 }
274
276};
277
278template <typename T>
280
281} // namespace madness
282
283#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
Below variables meaningful if user does specifies uniform spacing.
Definition interpolation_1d.h:67
int npt
No. of grid points.
Definition interpolation_1d.h:64
CubicInterpolationTable(std::vector< T > x, std::vector< T > y)
Definition interpolation_1d.h:203
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:198
double get_hi() const
Definition interpolation_1d.h:258
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:141
double rh
Definition interpolation_1d.h:68
std::vector< T > a
Definition interpolation_1d.h:65
CubicInterpolationTable(double lo, double hi, int npt, const std::vector< T > &y)
Definition interpolation_1d.h:225
double hi
Interpolation is in range [lo,hi].
Definition interpolation_1d.h:63
CubicInterpolationTable()
Definition interpolation_1d.h:176
virtual ~CubicInterpolationTable()
Definition interpolation_1d.h:275
static int min_npts_per_task_default
minimum # of points per task
Definition interpolation_1d.h:174
static void cubic_fit(const double *x, const T *f, T *a)
Grid points. Only stored if not evenly spaced. Empty otherwise.
Definition interpolation_1d.h:73
double get_lo() const
Definition interpolation_1d.h:256
void make_interpolation(const std::vector< double > &x, const std::vector< T > &p, const int npts_per_task=std::numeric_limits< int >::max() - 1, World *world_ptr=nullptr)
Definition interpolation_1d.h:98
CubicInterpolationTable(double lo, double hi, int npt, const functionT &f)
Definition interpolation_1d.h:185
double lo
Interpolation is in range [lo,hi].
Definition interpolation_1d.h:62
T operator()(double y) const
Definition interpolation_1d.h:240
std::vector< double > pts_
Below variables meaningful if user provides grid points.
Definition interpolation_1d.h:70
double err(const functionT &f) const
Definition interpolation_1d.h:260
static std::size_t size()
Returns the number of threads in the pool.
Definition thread.h:1413
A parallel world class.
Definition world.h:132
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:2448
static XNonlinearSolver< std::vector< Function< T, NDIM > >, T, vector_function_allocator< T, NDIM > > nonlinear_vector_solver(World &world, const long nvec)
Definition nonlinsol.h:284
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