USML
data_grid< DATA_TYPE, NUM_DIMS > Class Template Reference
Collaboration diagram for data_grid< DATA_TYPE, NUM_DIMS >:

Detailed Description

template<class DATA_TYPE, unsigned NUM_DIMS>
class usml::types::data_grid< DATA_TYPE, NUM_DIMS >

N-dimensional data set and its associated axes.

Supports interpolation in any number of dimensions.

Parameters
DATA_TYPEType of data to be interpolated. Must support +,-,*,/ with itself and double precision scalars.
NUM_DIMSNumber of dimensions in this grid. Specifying this at compile time allows for some loop un-wrapping.

Constructor & Destructor Documentation

data_grid ( )
inlineprotected

Default constructor for sub-classes.

data_grid ( seq_vector axis[])
inline

Create data grid from its associated axes.

Allocates new memory for the data at each grid point. Initialize all of the interpolation types to GRID_INTERP_LINEAR.

Parameters
axisAxes to use for each dimension of the grid. The seq_vector::clone() routine is used to make a local copy of each axis within the data grid.
data_grid ( const data_grid< DATA_TYPE, NUM_DIMS > &  other,
bool  copy_data 
)
inline

Create data grid from an existing grid.

Allocates new memory for the data at each grid point. Copies interpolation types from the original grid.

Parameters
otherGrid to be copied.
copy_dataCopy both axes and data if true. Copy just axes and structure if false.
~data_grid ( )
inline

Destroys memory area for field data.

Member Function Documentation

const seq_vector* axis ( int unsigned  dim) const
inline

Extract a reference to one of the axes.

DATA_TYPE data ( const unsigned *  index) const
inline

Extract a data value at a specific combination of indices.

Parameters
indexIndex number in each dimension.
void data ( const unsigned *  index,
DATA_TYPE  value 
)
inline

Define a new data value at a specific combination of indices.

Parameters
indexIndex number in each dimension.
valueValue to insert at this location.
bool edge_limit ( int unsigned  dimension) const
inline

Returns the edge_limit flag for the requested dimension.

void edge_limit ( int unsigned  dimension,
bool  flag 
)
inline

Set the edge_limit flag for the requested dimension.

Note: Default is true

Parameters
dimensionDimension number to be modified.
flagLimits locations when true.
DATA_TYPE interp ( int  dim,
const unsigned *  index,
const double *  location,
DATA_TYPE &  deriv,
DATA_TYPE *  deriv_vec 
) const
private

Private recursion engine for multi-dimensional interpolation.

Recursion engine for multi-dimensional interpolation.

The type of interpolation for each dimension is determined using the _interp_type[] field. Interpolation coefficients are computed on the fly to make arbitrary combinations of interpolation types viable.

Parameters
dimIndex of the dimension currently being processed. Recursion starts at dim=NUM_DIMS-1 and reduces to element retrieval when dim=-1.
indexPosition of the corner before the desired field point. Must have the same rank as the data grid.
locationLocation at which field value is desired. Must have the same NUM_DIM as the data grid or higher.
derivDerivative for this iteration.
deriv_vecResults vector for derivative. Derivative not computed if NULL.
Returns
Estimate of the field after interpolation.
enum GRID_INTERP_TYPE interp_type ( int unsigned  dimension) const
inline

Retrieve the type of interpolation for one of the axes.

void interp_type ( int unsigned  dimension,
enum GRID_INTERP_TYPE  type 
)
inline

Define the type of interpolation for one of the axes.

Modifies the axis buffer size as a side effect. Note that linear interpolation requires a minimum of 2 points; pchip requires a minimum of 4 points.

Parameters
dimensionDimension number to be modified.
typeType of interpolation for this dimension.
DATA_TYPE interpolate ( double *  location,
DATA_TYPE *  derivative = NULL 
)
inline

Multi-dimensional interpolation with the derivative calculation.

So many calculations are shared between the determination of an interpolate value and its derivative, that it is computationally efficient to compute them both at the same time.

Limit interpolation to axis domain if _edge_limit turned on for that dimension. Allow extrapolation if _edge_limit turned off.

Parameters
locationLocation at which field value is desired. Must have the same rank as the data grid or higher. WARNING: The contents of the location vector may be modified if edge_limit() is true for any dimension.
derivativeIf this is not null, the first derivative of the field at this point will also be computed.
Returns
Value of the field at this point.
void interpolate ( const matrix< double > &  x,
matrix< double > *  result,
matrix< double > *  dx = NULL 
)
inline

Interpolation 1-D specialization where the arguments, and results, are matrix<double>.

This is used frequently in the WaveQ3D model to interpolate environmental parameters.

Parameters
xFirst dimension of location.
resultInterpolated values at each location (output).
dxFirst dimension of derivative (output).
void interpolate ( const matrix< double > &  x,
const matrix< double > &  y,
matrix< double > *  result,
matrix< double > *  dx = NULL,
matrix< double > *  dy = NULL 
)
inline

Interpolation 2-D specialization where the arguments, and results, are matrix<double>.

This is used frequently in the WaveQ3D model to interpolate environmental parameters.

Parameters
xFirst dimension of location.
ySecond dimension of location.
resultInterpolated values at each location (output).
dxFirst dimension of derivative (output).
dySecond dimension of derivative (output).
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 
)
inline

Interpolation 3-D specialization where the arguments, and results, are matrix<double>.

This is used frequently in the WaveQ3D model to interpolate environmental parameters.

Parameters
xFirst dimension of location.
ySecond dimension of location.
zThird dimension of location.
resultInterpolated values at each location (output).
dxFirst dimension of derivative (output).
dySecond dimension of derivative (output).
dzThird dimension of derivative (output).
DATA_TYPE linear ( int  dim,
const unsigned *  index,
const double *  location,
DATA_TYPE &  deriv,
DATA_TYPE *  deriv_vec 
) const
inlineprivate

Perform a linear interpolation on this dimension.

Parameters
dimIndex of the dimension currently being processed.
indexPosition of the corner before the desired field point. Must have the same rank as the data grid.
locationLocation at which field value is desired. Must have the same rank as the data grid or higher.
derivDerivative for this iteration. Constant across the interval for linear interpolation.
deriv_vecResults vector for derivative. Derivative not computed if NULL.
Returns
Estimate of the field after interpolation.
DATA_TYPE nearest ( int  dim,
const unsigned *  index,
const double *  location,
DATA_TYPE &  deriv,
DATA_TYPE *  deriv_vec 
) const
inlineprivate

Perform a nearest neighbor interpolation on this dimension.

Parameters
dimIndex of the dimension currently being processed.
indexPosition of the corner before the desired field point. Must have the same NUM_DIM as the data grid.
locationLocation at which field value is desired. Must have the same NUM_DIM as the data grid or higher.
derivDerivative for this iteration. Always zero for nearest neighbor interpolation.
deriv_vecResults vector for derivative. Derivative not computed if NULL.
Returns
Estimate of the field after interpolation.
DATA_TYPE pchip ( int  dim,
const unsigned *  index,
const double *  location,
DATA_TYPE &  deriv,
DATA_TYPE *  deriv_vec 
) const
inlineprivate

Interpolate this dimension using the Piecewise Cubic Hermite Interpolation Polynomial (PCHIP) algorithm from Matlab.

Matlab uses shape preserving, "visually pleasing" version of the cubic interpolant that does not suffer from the overshooting issues prevalent in the cubic spline. Although the the first derivative of the PCHIP result is guaranteed to be continuous, the second derivative has no such guarantee.

This algorithm differs from the Matlab implementation in that is simultaneously interpolates the function value for the current dimension, and interpolates the derivative for the previous dimension.

When using a gridded data set, it is recommended that edge_limit be set to TRUE for any dimensional axis that uses the PCHIP interpolation. This is because of PCHIP allowing for extreme values when extrapolating data.

References:

Cleve Moler, Numerical Computing in Matlab, Chapter 3 Interpolation, http://www.mathworks.com/moler/chapters.html, accessed 5/15/2012.

F. N. Fritsch and R. E. Carlson, Monotone Piecewise Cubic Interpolation, SIAM Journal on Numerical Analysis, 17 (1980), pp. 238-246.

D. Kahaner, C. Moler, and S. Nash, Numerical Methods and Software, Prentice{Hall, Englewood Cli®s, NJ, 1989.

The basic algorithm assumes that the interpolation location is in the interval [ x[k], x[k+1] ), where "k" is known as the "interval index". The result is then calculated from four unevenly spaced points, and their forward (one-sided) derivatives.

                y0 = y[k-1]             h0 = x[k]-x[k-1]        deriv0 = (y1-y0)/h0
                y1 = y[k]               h1 = x[k+1]-x[k]        deriv1 = (y2-y1)/h1
                y2 = y[k+1]             h2 = x[k+2]-x[k+1]      deriv2 = (y3-y2)/h2
                y3 = y[k+2]             s = x - x[k]
such that
     p(x) = y[k+1]  * ( 3 h1 s^2 - 2 s^3 ) / h1^3
            y[k]    * ( h1^3 - 3 h1 s^2 - 2 s^3 ) / h1^3
            slope[k+1] * ( s1^2 (s-h) ) / h1^2
            slope[k]   * ( s (s-h[k])^2 ) / h1^2
where:
        slope[k] = weighted harmonic average of deriv0, deriv1, deriv2 terms

At the end-points, y'[0] and y'[N-1] must be estimated. This implementation uses Matlab's non-centered, shape-preserving, three-point formula for the end-point slope.

Parameters
dimIndex of the dimension currently being processed.
indexPosition of the corner before the desired field point. Must have the same rank as the data grid.
locationLocation at which field value is desired. Must have the same rank as the data grid or higher.
derivDerivative for this iteration. Constant across the interval for linear interpolation.
deriv_vecResults vector for derivative. Derivative not computed if NULL.
Returns
Estimate of the field after interpolation.

Member Data Documentation

seq_vector* _axis[NUM_DIMS]
protected

Axis associated with each dimension of the data grid.

DATA_TYPE* _data
protected

Multi-dimensional data stored as a linear array in column major order.

This format is used to support an N-dimensional data set with any number of dimensions. This memory is created in the constructor and deleted in the destructor.

bool _edge_limit[NUM_DIMS]
private

Limits locations to values inside axis when true.

enum GRID_INTERP_TYPE _interp_type[NUM_DIMS]
private

Defines the type of interpolation for each axis.

unsigned _offset[NUM_DIMS]
protected

Used during interpolation to hold the axis offsets.