33#ifndef MADNESS_MISC_INTERPOLATION_1D_H__INCLUDED
34#define MADNESS_MISC_INTERPOLATION_1D_H__INCLUDED
40#include "../world/world_task_queue.h"
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]));
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]);
99 const int npts_per_task = std::numeric_limits<int>::max() - 1,
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,
113 this->a[i * 5] =
mid;
126 for (
int j = 0; j < 5; ++j) {
128 a[5 *
npt - 5 + j] =
a[5 *
npt - 10 + j] =
a[5 *
npt - 15 + j];
140 template <
typename functionT>
148 std::vector<T>
p(
npt);
149 std::vector<double> x(
npt);
157 x[i] = this->lo + i * this->
h;
184 template <
typename functionT>
197 template <
typename functionT>
206 if (x.size() != y.size()) {
207 throw std::runtime_error(
"Sizes of input and output arrays do not equal.");
210 std::vector<int>
indices(x.size());
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++) {
217 std::stable_sort(x.begin(), x.end());
219 hi = x[x.size() - 1];
226 const std::vector<T> &y)
230 if ((
int)y.size() <
npt)
231 throw "Insufficient y-points";
233 std::vector<double> x(
npt);
234 for (
int i = 0; i <
npt; ++i)
244 i = std::lower_bound(
pts_.begin(),
pts_.end(), y) -
pts_.begin();
247 i =
static_cast<int>((y -
lo) *
rh);
253 return (
a[i + 1] +
y1 *
a[i + 2]) + yy * (
a[i + 3] +
y1 *
a[i + 4]);
261 if (!
pts_.empty()) {
262 throw "This error test is only meaningful if f generated the interpolation.";
266 for (
int i = 0; i < 7 *
npt; ++i) {
267 double x =
lo +
h7 * i;
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