USML
data_grid.h
1 
5 #ifndef USML_TYPES_DATA_GRID_H
6 #define USML_TYPES_DATA_GRID_H
7 
8 #include <string.h>
9 #include <usml/types/wvector.h>
10 #include <usml/types/seq_vector.h>
11 
12 using namespace usml::ublas;
13 
14 namespace usml {
15 namespace types {
16 
19 
22 {
25 };
26 
38 template<unsigned Dim> inline size_t data_grid_compute_offset(
39  seq_vector *axis[], const unsigned* index)
40 {
41  return index[Dim] + axis[Dim]->size() * data_grid_compute_offset<Dim - 1> (
42  axis, index);
43 }
44 
49 template<> inline size_t data_grid_compute_offset<0> (seq_vector *axis[],
50  const unsigned* index)
51 {
52  return index[0];
53 }
54 
64 template<class DATA_TYPE, unsigned NUM_DIMS> class data_grid
65 {
66 
67  //*************************************************************************
68  // Axis and Data properties
69 
70 protected:
71 
73  seq_vector* _axis[NUM_DIMS];
74 
81  DATA_TYPE *_data;
82 
84  unsigned _offset[NUM_DIMS];
85 
86 
87 public:
88 
92  inline const seq_vector* axis(int unsigned dim) const
93  {
94  return _axis[dim];
95  }
96 
102  inline DATA_TYPE data(const unsigned* index) const
103  {
104  const size_t offset = data_grid_compute_offset<NUM_DIMS - 1> (
105  (seq_vector**) _axis, index);
106  return _data[offset];
107  }
108 
115  inline void data(const unsigned* index, DATA_TYPE value)
116  {
117  const size_t offset = data_grid_compute_offset<NUM_DIMS - 1> (_axis,
118  index);
119  _data[offset] = value;
120  }
121 
122  //*************************************************************************
123  // interp_type property
124 
125 private:
126 
128  enum GRID_INTERP_TYPE _interp_type[NUM_DIMS];
129 
130 public:
131 
135  inline enum GRID_INTERP_TYPE interp_type(int unsigned dimension) const
136  {
137  return _interp_type[dimension];
138  }
139 
149  inline void interp_type(int unsigned dimension, enum GRID_INTERP_TYPE type)
150  {
151  const unsigned size = _axis[dimension]->size() ;
152  if ( type > GRID_INTERP_NEAREST && size < 2 ) {
153  type = GRID_INTERP_NEAREST ;
154  } else if ( type > GRID_INTERP_LINEAR && size < 4 ) {
155  type = GRID_INTERP_LINEAR ;
156  } else {
157  _interp_type[dimension] = type;
158  }
159  }
160 
161  //*************************************************************************
162  // edge_limit property
163 
164 private:
165 
167  bool _edge_limit[NUM_DIMS];
168 
169 public:
170 
174  inline bool edge_limit(int unsigned dimension) const
175  {
176  return _edge_limit[dimension];
177  }
178 
186  inline void edge_limit(int unsigned dimension, bool flag)
187  {
188  _edge_limit[dimension] = flag;
189  }
190 
191  //*************************************************************************
192  // interpolation methods
193 
194 private:
195 
214  DATA_TYPE interp(int dim, const unsigned* index, const double* location,
215  DATA_TYPE& deriv, DATA_TYPE* deriv_vec) const;
216  // forward reference needed for recursion
217 
232  DATA_TYPE nearest(int dim, const unsigned* index, const double* location,
233  DATA_TYPE& deriv, DATA_TYPE* deriv_vec) const
234  {
235  DATA_TYPE result, da;
236 
237  // compute field value in this dimension
238 
239  const unsigned k = index[dim];
240  seq_vector* ax = _axis[dim];
241  const double u = (location[dim] - (*ax)(k)) / ax->increment(k);
242  if (u < 0.5) {
243  result = interp(dim - 1, index, location, da, deriv_vec);
244  } else {
245  unsigned next[NUM_DIMS];
246  memcpy(next, index, NUM_DIMS * sizeof(unsigned));
247  ++next[dim];
248  result = interp(dim - 1, next, location, da, deriv_vec);
249  }
250 
251  // compute derivative in this dimension
252 
253  if (deriv_vec) {
254  deriv = 0.0;
255  deriv_vec[dim] = deriv;
256  if (dim > 0) deriv_vec[dim - 1] = da;
257  }
258 
259  // use results for dim+1 iteration
260 
261  return result;
262  }
263 
278  DATA_TYPE linear(int dim, const unsigned* index, const double* location,
279  DATA_TYPE& deriv, DATA_TYPE* deriv_vec) const
280  {
281  DATA_TYPE result, da, db;
282 
283  // build interpolation coefficients
284 
285  unsigned next[NUM_DIMS];
286  memcpy(next, index, NUM_DIMS * sizeof(unsigned));
287  ++next[dim];
288 
289  const DATA_TYPE a = interp(dim - 1, index, location, da, deriv_vec);
290  const DATA_TYPE b = interp(dim - 1, next, location, db, deriv_vec);
291  const unsigned k = index[dim];
292  seq_vector* ax = _axis[dim];
293 
294  // compute field value in this dimension
295 
296  const DATA_TYPE h = (DATA_TYPE) ax->increment(k) ;
297  const DATA_TYPE u = (location[dim] - (*ax)(k)) / h ;
298  result = a * (1.0 - u) + b * u;
299 
300  // compute derivative in this dimension and prior dimension
301 
302  if (deriv_vec) {
303  deriv = (b - a) / h ;
304  deriv_vec[dim] = deriv;
305  if (dim > 0) {
306  deriv_vec[dim - 1] = da * (1.0 - u) + db * u;
307  }
308  }
309 
310  // use results for dim+1 iteration
311 
312  return result;
313  }
314 
378  DATA_TYPE pchip(int dim, const unsigned* index, const double* location,
379  DATA_TYPE& deriv, DATA_TYPE* deriv_vec) const
380  {
381  DATA_TYPE result ;
382  DATA_TYPE y0, y1, y2, y3 ; // dim-1 values at k-1, k, k+1, k+2
383  DATA_TYPE dy0, dy1, dy2, dy3 ; // dim-1 derivs at k-1, k, k+1, k+2
384  const unsigned kmin = 1u ; // at endpt if k-1 < 0
385  const unsigned kmax = _axis[dim]->size()-3u ; // at endpt if k+2 > N-1
386 
387  // interpolate in dim-1 dimension to find values and derivs at k, k-1
388 
389  const unsigned k = index[dim];
390  seq_vector* ax = _axis[dim];
391  y1 = interp( dim-1, index, location, dy1, deriv_vec );
392 
393  if ( k >= kmin ) {
394  unsigned prev[NUM_DIMS];
395  memcpy(prev, index, NUM_DIMS * sizeof(unsigned));
396  --prev[dim];
397  y0 = interp( dim-1, prev, location, dy0, deriv_vec );
398  } else { // use harmless values at left end-point
399  y0 = y1 ;
400  dy0 = dy1 ;
401  }
402 
403  // interpolate in dim-1 dimension to find values and derivs at k+1, k+2
404 
405  unsigned next[NUM_DIMS];
406  memcpy(next, index, NUM_DIMS * sizeof(unsigned));
407  ++next[dim];
408  y2 = interp( dim-1, next, location, dy2, deriv_vec );
409 
410  if ( k <= kmax ) {
411  unsigned last[NUM_DIMS];
412  memcpy(last, next, NUM_DIMS * sizeof(unsigned));
413  ++last[dim];
414  y3 = interp(dim - 1, last, location, dy3, deriv_vec);
415  } else { // use harmless values at right end-point
416  y3 = y2 ;
417  dy3 = dy2 ;
418  }
419 
420  // compute difference values used frequently in computation
421 
422  const DATA_TYPE h0 = (DATA_TYPE) ax->increment(k - 1); // interval from k-1 to k
423  const DATA_TYPE h1 = (DATA_TYPE) ax->increment(k); // interval from k to k+1
424  const DATA_TYPE h2 = (DATA_TYPE) ax->increment(k + 1); // interval from k+1 to k+2
425  const DATA_TYPE h1_2 = h1 * h1; // k to k+1 interval squared
426  const DATA_TYPE h1_3 = h1_2 * h1; // k to k+1 interval cubed
427 
428  const DATA_TYPE s = location[dim]-(*ax)(k); // local variable
429  const DATA_TYPE s_2 = s * s, s_3 = s_2 * s; // s squared and cubed
430  const DATA_TYPE sh_minus = s - h1;
431  const DATA_TYPE sh_term = 3.0 * h1 * s_2 - 2.0 * s_3;
432 
433  // compute first divided differences (forward derivative)
434  // for both the values, and their derivatives
435 
436  const DATA_TYPE deriv0 = (y1 - y0) / h0; // fwd deriv from k-1 to k
437  const DATA_TYPE deriv1 = (y2 - y1) / h1; // fwd deriv from k to k+1
438  const DATA_TYPE deriv2 = (y3 - y2) / h2; // fwd deriv from k+1 to k+2
439 
440  DATA_TYPE dderiv0=0.0, dderiv1=0.0, dderiv2=0.0 ;
441  if (deriv_vec) { // fwd deriv of dim-1 derivatives
442  dderiv0 = (dy1 - dy0) / h0;
443  dderiv1 = (dy2 - dy1) / h1;
444  dderiv2 = (dy3 - dy2) / h2;
445  }
446 
447  //*************
448  // compute weighted harmonic mean of slopes around index k
449  // for both the values, and their derivatives
450  // set it zero at local maxima or minima
451  // deriv0 * deriv1 condition guards against division by zero
452 
453  DATA_TYPE slope1=0.0, dslope1=0.0;
454 
455  // when not at an end-point, slope1 is the harmonic, weighted
456  // average of deriv0 and deriv1.
457 
458  if ( k >= kmin ) {
459  const DATA_TYPE w0 = 2.0 * h1 + h0;
460  const DATA_TYPE w1 = h1 + 2.0 * h0;
461  if ( deriv0 * deriv1 > 0.0 ) {
462  slope1 = (w0 + w1) / ( w0 / deriv0 + w1 / deriv1 );
463  }
464  if ( deriv_vec != NULL && dderiv0 * dderiv1 > 0.0 ) {
465  dslope1 = (w0 + w1) / ( w0 / dderiv0 + w1 / dderiv1 );
466  }
467 
468  // at left end-point, use Matlab end-point formula with slope limits
469  // note that the deriv0 value is bogus values when this is true
470 
471  } else {
472  slope1 = ( (2.0+h1+h2) * deriv1 - h1 * deriv2 ) / (h1+h2) ;
473  if ( slope1 * deriv1 < 0.0 ) {
474  slope1 = 0.0 ;
475  } else if ( (deriv1*deriv2 < 0.0) && (abs(slope1) > abs(3.0*deriv1)) ) {
476  slope1 = 3.0*deriv1 ;
477  }
478  if ( deriv_vec ) {
479  dslope1 = ( (2.0+h1+h2) * dderiv1 - h1 * dderiv2 ) / (h1+h2) ;
480  if ( dslope1 * dderiv1 < 0.0 ) {
481  dslope1 = 0.0 ;
482  } else if ( (dderiv1*dderiv2 < 0.0) && (abs(dslope1) > abs(3.0*dderiv1)) ) {
483  dslope1 = 3.0*dderiv1 ;
484  }
485  }
486  }
487 
488  //*************
489  // compute weighted harmonic mean of slopes around index k+1
490  // for both the values, and their derivatives
491  // set it zero at local maxima or minima
492  // deriv1 * deriv2 condition guards against division by zero
493 
494  DATA_TYPE slope2=0.0, dslope2=0.0;
495 
496  // when not at an end-point, slope2 is the harmonic, weighted
497  // average of deriv1 and deriv2.
498 
499  if ( k <= kmax ) {
500  const DATA_TYPE w1 = 2.0 * h1 + h0;
501  const DATA_TYPE w2 = h1 + 2.0 * h0;
502  if ( deriv1 * deriv2 > 0.0 ) {
503  slope2 = (w1 + w2) / ( w1 / deriv1 + w2 / deriv2 );
504  }
505  if ( deriv_vec != NULL && dderiv1 * dderiv2 > 0.0 ) {
506  dslope2 = (w1 + w2) / ( w1 / dderiv1 + w2 / dderiv2 );
507  }
508 
509  // at right end-point, use Matlab end-point formula with slope limits
510  // note that the deriv2 value is bogus values when this is true
511 
512  } else { // otherwise, compute harmonic weighted average
513  slope2 = ( (2.0+h1+h2) * deriv1 - h1 * deriv0 ) / (h1+h0) ;
514  if ( slope2 * deriv1 < 0.0 ) {
515  slope2 = 0.0 ;
516  } else if ( (deriv1*deriv0 < 0.0) && (abs(slope2) > abs(3.0*deriv1)) ) {
517  slope2 = 3.0*deriv1 ;
518  }
519  if ( deriv_vec ) {
520  dslope2 = ( (2.0+h1+h2) * dderiv1 - h1 * dderiv0 ) / (h1+h0) ;
521  if ( dslope2 * dderiv1 < 0.0 ) {
522  dslope2 = 0.0 ;
523  } else if ( (dderiv1*dderiv0 < 0.0) && (abs(dslope2) > abs(3.0*dderiv1)) ) {
524  dslope2 = 3.0*dderiv1 ;
525  }
526  }
527  }
528 
529  // compute interpolation value in this dimension
530 
531  result = y2 * sh_term / h1_3
532  + y1 * (h1_3 - sh_term) / h1_3
533  + slope2 * s_2 * sh_minus / h1_2
534  + slope1 * s * sh_minus * sh_minus / h1_2;
535 
536  // compute derivative in this dimension
537  // assume linear change of slope across interval
538 
539  if (deriv_vec) {
540  const DATA_TYPE u = s / h1;
541  deriv = slope1 * (1.0 - u) + slope2 * u;
542  deriv_vec[dim] = deriv;
543  if (dim > 0) {
544  deriv_vec[dim-1] = dy2 * sh_term / h1_3
545  + dy1 * (h1_3 - sh_term) / h1_3
546  + dslope2 * s_2 * sh_minus / h1_2
547  + dslope1 * s * sh_minus * sh_minus / h1_2;
548  }
549  }
550 
551  // use results for dim+1 iteration
552 
553  return result;
554  }
555 
556 public:
557 
576  DATA_TYPE interpolate(double* location, DATA_TYPE* derivative = NULL)
577  {
578  // find the "interval index" in each dimension
579 
580  for (unsigned dim = 0; dim < NUM_DIMS; ++dim) {
581 
582  // limit interpolation to axis domain if _edge_limit turned on
583 
584  if ( _edge_limit[dim] ) {
585  double a = *(_axis[dim]->begin()) ;
586  double b = *(_axis[dim]->rbegin()) ;
587  double inc = _axis[dim]->increment(0);
588  if ( inc < 0) { // a > b
589  if ( location[dim] >= a ) { //left of the axis
590  location[dim] = a ;
591  _offset[dim] = 0 ;
592  } else if ( location[dim] <= b ) { //right of the axis
593  location[dim] = b ;
594  _offset[dim] = _axis[dim]->size()-2 ;
595  } else {
596  _offset[dim] = _axis[dim]->find_index(location[dim]); //somewhere in-between the endpoints of the axis
597  }
598  }
599  if (inc > 0 ) { // a < b
600  if ( location[dim] <= a ) { //left of the axis
601  location[dim] = a ;
602  _offset[dim] = 0 ;
603  } else if ( location[dim] >= b ) { //right of the axis
604  location[dim] = b ;
605  _offset[dim] = _axis[dim]->size()-2 ;
606  } else {
607  _offset[dim] = _axis[dim]->find_index(location[dim]); //somewhere in-between the endpoints of the axis
608  }
609  }
610 
611  // allow extrapolation if _edge_limit turned off
612 
613  } else {
614  _offset[dim] = _axis[dim]->find_index(location[dim]);
615  }
616  }
617 
618  // compute interpolation results for value and derivative
619 
620  DATA_TYPE dresult;
621  return interp(NUM_DIMS - 1, _offset, location, dresult, derivative);
622  }
623 
633  void interpolate(const matrix<double>& x, matrix<double>* result, matrix<
634  double>* dx = NULL)
635  {
636  double location[1];
637  DATA_TYPE derivative[1];
638  for (unsigned n = 0; n < x.size1(); ++n) {
639  for (unsigned m = 0; m < x.size2(); ++m) {
640  location[0] = x(n, m);
641  if (dx == NULL) {
642  (*result)(n, m) = (double) interpolate(location);
643  } else {
644  (*result)(n, m)
645  = (double) interpolate(location, derivative);
646  (*dx)(n, m) = (double) derivative[0];
647  }
648  }
649  }
650  }
651 
663  void interpolate(const matrix<double>& x, const matrix<double>& y, matrix<
664  double>* result, matrix<double>* dx = NULL, matrix<double>* dy =
665  NULL)
666  {
667  double location[2];
668  DATA_TYPE derivative[2];
669  for (unsigned n = 0; n < x.size1(); ++n) {
670  for (unsigned m = 0; m < x.size2(); ++m) {
671  location[0] = x(n, m);
672  location[1] = y(n, m);
673  if (dx == NULL || dy == NULL) {
674  (*result)(n, m) = (double) interpolate(location);
675  } else {
676  (*result)(n, m)
677  = (double) interpolate(location, derivative);
678  (*dx)(n, m) = (double) derivative[0];
679  (*dy)(n, m) = (double) derivative[1];
680  }
681  }
682  }
683  }
684 
698  void interpolate(const matrix<double>& x, const matrix<double>& y,
699  const matrix<double>& z, matrix<double>* result,
700  matrix<double>* dx = NULL, matrix<double>* dy = NULL,
701  matrix<double>* dz = NULL)
702  {
703  double location[3];
704  DATA_TYPE derivative[3];
705  for (unsigned n = 0; n < x.size1(); ++n) {
706  for (unsigned m = 0; m < x.size2(); ++m) {
707  location[0] = x(n, m);
708  location[1] = y(n, m);
709  location[2] = z(n, m);
710  if (dx == NULL || dy == NULL || dz == NULL) {
711  (*result)(n, m) = (double) interpolate(location);
712  } else {
713  (*result)(n, m)
714  = (double) interpolate(location, derivative);
715  (*dx)(n, m) = (double) derivative[0];
716  (*dy)(n, m) = (double) derivative[1];
717  (*dz)(n, m) = (double) derivative[2];
718  }
719  }
720  }
721  }
722 
723  //*************************************************************************
724  // initialization
725 
726 protected:
727 
732 
733  for (unsigned n = 0; n < NUM_DIMS; ++n) {
734  _axis[n] = NULL;
735  }
736  _data = NULL;
737 
738  memset(_interp_type, GRID_INTERP_LINEAR, NUM_DIMS * sizeof(enum GRID_INTERP_TYPE));
739  memset(_edge_limit, true, NUM_DIMS * sizeof(bool));
740  }
741 
742 public:
743 
754  {
755  size_t N = 1;
756  for (unsigned n = 0; n < NUM_DIMS; ++n) {
757  _axis[n] = axis[n]->clone();
758  N *= _axis[n]->size();
759  interp_type( n, GRID_INTERP_LINEAR ) ;
760  }
761  _data = new DATA_TYPE[N];
762  memset(_data, 0, N * sizeof(DATA_TYPE));
763  memset(_edge_limit, true, NUM_DIMS * sizeof(bool));
764  }
765 
775  data_grid(const data_grid& other, bool copy_data)
776  {
777  size_t N = 1;
778  for (unsigned n = 0; n < NUM_DIMS; ++n) {
779  _axis[n] = other._axis[n]->clone();
780  N *= _axis[n]->size();
781  }
782  _data = new DATA_TYPE[N];
783  if (copy_data) {
784  memcpy(_data, other._data, N * sizeof(DATA_TYPE));
785  } else {
786  memset(_data, 0, N * sizeof(DATA_TYPE));
787  }
788  memcpy(_interp_type, other._interp_type, NUM_DIMS
789  * sizeof(enum GRID_INTERP_TYPE));
790  memset(_edge_limit, true, NUM_DIMS * sizeof(bool));
791  }
792 
797  {
798  for (unsigned n = 0; n < NUM_DIMS; ++n) {
799  if (_axis[n] != NULL) {
800  delete _axis[n];
801  }
802  }
803  if (_data != NULL) {
804  delete[] _data;
805  }
806  }
807 
808 }; // end data_grid class
809 
810 //*************************************************************************
811 // implementation of forward methods
812 
816 template<class DATA_TYPE, unsigned NUM_DIMS>
818  const unsigned* index, const double* location, DATA_TYPE& deriv,
819  DATA_TYPE* deriv_vec) const
820 {
821  DATA_TYPE result;
822 
823  if (dim < 0) {
824  result = data(index);
825  // terminates recursion
826 
827  } else {
828  if (_interp_type[dim] == GRID_INTERP_LINEAR) {
829  result = linear(dim, index, location, deriv, deriv_vec);
830  } else if (_interp_type[dim] == GRID_INTERP_PCHIP) {
831  result = pchip(dim, index, location, deriv, deriv_vec);
832  } else {
833  result = nearest(dim, index, location, deriv, deriv_vec);
834  }
835  }
836  return result;
837 }
838 
839 } // end of namespace types
840 } // end of namespace usml
841 
842 #endif
seq_vector * _axis[NUM_DIMS]
Axis associated with each dimension of the data grid.
Definition: data_grid.h:73
const seq_vector * axis(int unsigned dim) const
Extract a reference to one of the axes.
Definition: data_grid.h:92
value_type increment(size_type index) const
Retrieves the increment between two elements in this sequence.
Definition: seq_vector.h:119
data_grid(seq_vector *axis[])
Create data grid from its associated axes.
Definition: data_grid.h:753
DATA_TYPE pchip(int dim, const unsigned *index, const double *location, DATA_TYPE &deriv, DATA_TYPE *deriv_vec) const
Interpolate this dimension using the Piecewise Cubic Hermite Interpolation Polynomial (PCHIP) algorit...
Definition: data_grid.h:378
data_grid(const data_grid &other, bool copy_data)
Create data grid from an existing grid.
Definition: data_grid.h:775
Definition: data_grid.h:23
virtual seq_vector * clone() const =0
Create a copy using a reference to the base class.
size_t data_grid_compute_offset< 0 >(seq_vector *axis[], const unsigned *index)
Definition: data_grid.h:49
DATA_TYPE interpolate(double *location, DATA_TYPE *derivative=NULL)
Multi-dimensional interpolation with the derivative calculation.
Definition: data_grid.h:576
Definition: data_grid.h:23
BOOST_UBLAS_INLINE matrix_unary1_traits< E, scalar_abs< typename E::value_type > >::result_type abs(const matrix_expression< E > &e)
Magnitude of a complex matrix.
Definition: matrix_math.h:257
DATA_TYPE * _data
Multi-dimensional data stored as a linear array in column major order.
Definition: data_grid.h:81
data_grid()
Default constructor for sub-classes.
Definition: data_grid.h:731
~data_grid()
Destroys memory area for field data.
Definition: data_grid.h:796
size_t data_grid_compute_offset(seq_vector *axis[], const unsigned *index)
Definition: data_grid.h:38
void interp_type(int unsigned dimension, enum GRID_INTERP_TYPE type)
Define the type of interpolation for one of the axes.
Definition: data_grid.h:149
void interpolate(const matrix< double > &x, const matrix< double > &y, const matrix< double > &z, matrix< double > *result, matrix< double > *dx=NULL, matrix< double > *dy=NULL, matrix< double > *dz=NULL)
Interpolation 3-D specialization where the arguments, and results, are matrix<double>.
Definition: data_grid.h:698
void interpolate(const matrix< double > &x, matrix< double > *result, matrix< double > *dx=NULL)
Interpolation 1-D specialization where the arguments, and results, are matrix<double>.
Definition: data_grid.h:633
bool edge_limit(int unsigned dimension) const
Returns the edge_limit flag for the requested dimension.
Definition: data_grid.h:174
DATA_TYPE linear(int dim, const unsigned *index, const double *location, DATA_TYPE &deriv, DATA_TYPE *deriv_vec) const
Perform a linear interpolation on this dimension.
Definition: data_grid.h:278
size_type size() const
Returns the number of elements in this sequence.
Definition: seq_vector.h:138
void interpolate(const matrix< double > &x, const matrix< double > &y, matrix< double > *result, matrix< double > *dx=NULL, matrix< double > *dy=NULL)
Interpolation 2-D specialization where the arguments, and results, are matrix<double>.
Definition: data_grid.h:663
enum GRID_INTERP_TYPE _interp_type[NUM_DIMS]
Defines the type of interpolation for each axis.
Definition: data_grid.h:128
N-dimensional data set and its associated axes.
Definition: data_grid.h:64
GRID_INTERP_TYPE
Type of interpolation used for each axis.
Definition: data_grid.h:21
void edge_limit(int unsigned dimension, bool flag)
Set the edge_limit flag for the requested dimension.
Definition: data_grid.h:186
void data(const unsigned *index, DATA_TYPE value)
Define a new data value at a specific combination of indices.
Definition: data_grid.h:115
Definition: data_grid.h:24
DATA_TYPE nearest(int dim, const unsigned *index, const double *location, DATA_TYPE &deriv, DATA_TYPE *deriv_vec) const
Perform a nearest neighbor interpolation on this dimension.
Definition: data_grid.h:232
A read-only, monotonic sequence of values.
Definition: seq_vector.h:36
DATA_TYPE data(const unsigned *index) const
Extract a data value at a specific combination of indices.
Definition: data_grid.h:102