MultiDimensionalInterpolationTempl< T > Class Template Reference

This class interpolates multi-dimensional data sets. More...

#include <MultiDimensionalInterpolation.h>

Public Member Functions
	MultiDimensionalInterpolationTempl (const std::vector< std::vector< Real > > &base_points, const MultiIndex< Real > &data)
	Constructor, Takes a double vector containing the interpolation base points, and a MultiIndex object.
	MultiDimensionalInterpolationTempl ()
virtual	~MultiDimensionalInterpolationTempl ()=default
unsigned int	dim () const
	returns the dimensionality of this interpolation object
void	setData (const std::vector< std::vector< Real > > &base_points, const MultiIndex< Real > &data)
	sets data but also fixes degenerate dimensions in data
void	linearSearch (std::vector< T > &values, MultiIndex< Real >::size_type &indices) const
	linearSearch finds the indices i_k of the base point such that base_point[d][i_k] <= values[d] < base_point[d][i_k + 1] Note that argument is non-const to handle smaller/larger than tabulation range For tabulations < 100 entries, linearSearch is expected to be better than bisection search
T	multiLinearInterpolation (const std::vector< T > &x) const
	This function will take an independent variable input and will return the dependent variable based on multi-linear interpolation.

Protected Member Functions
void	errorCheck ()
	checks consistency of the data
unsigned int	linearSearchHelper (T &x, const std::vector< Real > &vector) const
	linear search helper for a single std::vector; note that argument is non-const to handle smaller/larger than tabulation range

Protected Attributes
unsigned int	_original_dim
	original dimension is to allow checks on user inputs for cases where arrays are sliced
std::vector< bool >	_degenerate_index
	this variable keeps track on which dimension is degenerate and was removed
bool	_degenerate_interpolation = false
	if all dimensions have size one then there is only one value to return this corner case should be supported so the calling routine does not need to check for it

Private Attributes
std::vector< std::vector< Real > >	_base_points
MultiIndex< Real >	_data
bool	_setup_complete = false

Detailed Description

template<typename T>
class MultiDimensionalInterpolationTempl< T >

This class interpolates multi-dimensional data sets.

Definition at line 23 of file MultiDimensionalInterpolation.h.

Constructor & Destructor Documentation

MultiDimensionalInterpolationTempl() [1/2]

template<typename T >

MultiDimensionalInterpolationTempl< T >::MultiDimensionalInterpolationTempl	(	const std::vector< std::vector< Real > > &	base_points,
		const MultiIndex< Real > &	data
	)

Constructor, Takes a double vector containing the interpolation base points, and a MultiIndex object.

The double vector base_points and the MultiIndex object data must

Definition at line 17 of file MultiDimensionalInterpolation.C.

  : _base_points(std::vector<std::vector<Real>>()), _data(MultiIndex<Real>({}))
{
  setData(base_points, data);
}

MultiDimensionalInterpolationTempl() [2/2]

template<typename T >

MultiDimensionalInterpolationTempl< T >::MultiDimensionalInterpolationTempl

(

)

Definition at line 25 of file MultiDimensionalInterpolation.C.

  : _base_points(std::vector<std::vector<Real>>()), _data(MultiIndex<Real>({}))
{
}

~MultiDimensionalInterpolationTempl()

template<typename T >

virtual MultiDimensionalInterpolationTempl< T >::~MultiDimensionalInterpolationTempl ( )

virtualdefault

Member Function Documentation

dim()

template<typename T >

unsigned int MultiDimensionalInterpolationTempl< T >::dim ( ) const

inline

returns the dimensionality of this interpolation object

Definition at line 37 of file MultiDimensionalInterpolation.h.

{ return _original_dim; }

errorCheck()

template<typename T >

void MultiDimensionalInterpolationTempl< T >::errorCheck ( )

protected

checks consistency of the data

Definition at line 91 of file MultiDimensionalInterpolation.C.

{
  if (_base_points.size() != _data.dim())
    throw std::domain_error("Size of the leading dimension of base points must be equal to "
                            "dimensionality of data array");
 
  MultiIndex<Real>::size_type shape = _data.size();
  for (unsigned int j = 0; j < _data.dim(); ++j)
  {
    if (shape[j] != _base_points[j].size())
    {
      std::ostringstream oss;
      oss << "Dimension " << j << " has " << _base_points[j].size()
          << " base points but dimension of data array is " << shape[j];
      throw std::domain_error(oss.str());
    }
 
    for (unsigned int i = 1; i < _base_points[j].size(); ++i)
    {
      if (_base_points[j][i - 1] >= _base_points[j][i])
      {
        std::ostringstream oss;
        oss << "Base point values for dimension " << j << " are not strictly increasing: bp["
            << i - 1 << "]: " << _base_points[j][i - 1] << " bc[" << i
            << "]: " << _base_points[j][i];
        throw std::domain_error(oss.str());
      }
    }
  }
 
  // now the setup is complete
  _setup_complete = true;
}

linearSearch()

template<typename T >

void MultiDimensionalInterpolationTempl< T >::linearSearch	(	std::vector< T > &	values,
		MultiIndex< Real >::size_type &	indices
	)		const

linearSearch finds the indices i_k of the base point such that base_point[d][i_k] <= values[d] < base_point[d][i_k + 1] Note that argument is non-const to handle smaller/larger than tabulation range For tabulations < 100 entries, linearSearch is expected to be better than bisection search

Definition at line 127 of file MultiDimensionalInterpolation.C.

{
  assert(values.size() == _data.dim());
 
  indices.resize(_data.dim());
  for (unsigned int d = 0; d < _data.dim(); ++d)
    indices[d] = linearSearchHelper(values[d], _base_points[d]);
}

linearSearchHelper()

template<typename T >

unsigned int MultiDimensionalInterpolationTempl< T >::linearSearchHelper ( T &  x, const std::vector< Real > &  vector  ) const

protected

linear search helper for a single std::vector; note that argument is non-const to handle smaller/larger than tabulation range

Definition at line 139 of file MultiDimensionalInterpolation.C.

{
// check order of vector only in debug
#ifndef NDEBUG
  for (unsigned int i = 1; i < vector.size(); ++i)
  {
    if (vector[i - 1] >= vector[i])
    {
      std::ostringstream oss;
      oss << "In linearSearchHelper vector is not strictly increasing: vector[" << i - 1
          << "]: " << vector[i - 1] << " vector[" << i << "]: " << vector[i];
      throw std::domain_error(oss.str());
    }
  }
#endif
 
  // smaller/larger cases for tabulation
  if (x < vector[0])
  {
    x = vector[0];
    return 0;
  }
  else if (x >= vector.back())
  {
    // this requires explanation because it looks like a bug (I swear it isn't, it's a feature)
    // if x >= largest base point value and the index i = vector.size() - 1 is returned, then during
    // the interpolation, i + 1 will be accessed, causing a segementation fault. Logic in the
    // interpolation is also bad because it may slow it down. By returning i = vector.size() - 2,
    // the last two values in the base point array are used, but x is reset to be equal to the last
    // value. This will give the right interpolation behavior.
    x = vector.back();
    return vector.size() - 2;
  }
 
  for (unsigned int j = 0; j < vector.size() - 2; ++j)
    if (x >= vector[j] && x < vector[j + 1])
      return j;
  return vector.size() - 2;
}

multiLinearInterpolation()

template<typename T >

T MultiDimensionalInterpolationTempl< T >::multiLinearInterpolation

(

const std::vector< T > &

x

)

const

This function will take an independent variable input and will return the dependent variable based on multi-linear interpolation.

Multi-linear interpolation uses all 2^d values in the base points surrounding x.

Definition at line 182 of file MultiDimensionalInterpolation.C.

{
  // ensure that object has been set up properly
  if (!_setup_complete)
    throw std::domain_error(
        "Object has been constructed properly. This happens when the empty constructor is called "
        "and setData is not called before interpolation.");
 
  // assert does not seem to be enough to check data consistency here
  // because this function is exposed to user and will be frequently used
  if (x.size() != _original_dim)
  {
    std::ostringstream oss;
    oss << "In sample the parameter x has size " << x.size() << " but data has dimension "
        << _data.dim();
    throw std::domain_error(oss.str());
  }
 
  // in this case there is only one entry that resides at flat index 0
  if (_degenerate_interpolation)
    return _data[0];
 
  // make a copy of x because linearSearch needs to be able to modify it
  // also adjust for degenerate indices that were removed
  std::vector<T> y(_data.dim());
  unsigned int l = 0;
  for (unsigned int j = 0; j < _original_dim; ++j)
    if (!_degenerate_index[j])
    {
      y[l] = x[j];
      ++l;
    }
 
  // obtain the indices using linearSearch & get flat index in _data array
  MultiIndex<Real>::size_type indices;
  linearSearch(y, indices);
 
  // pre-compute the volume of the hypercube and weights used for the interpolation
  Real volume = 1;
  std::vector<std::vector<T>> weights(_data.dim());
  for (unsigned int j = 0; j < _data.dim(); ++j)
  {
    volume *= _base_points[j][indices[j] + 1] - _base_points[j][indices[j]];
 
    // now compute weight for each dimension, note that
    // weight[j][0] = bp[j][high] - y[j]
    // weight[j][0] = y[j] - bp[j][low]
    // because in the interpolation the value on the "left" is weighted with the
    // distance of y to the "right"
    weights[j].resize(2);
    weights[j][0] = _base_points[j][indices[j] + 1] - y[j];
    weights[j][1] = y[j] - _base_points[j][indices[j]];
  }
 
  // get flat index and stride
  unsigned int flat_index = _data.flatIndex(indices);
  MultiIndex<Real>::size_type stride = _data.stride();
 
  // this is the interpolation routine, evaluation can be expensive if dim is large
 
  // first we retrieve the data around the point y
  // we use a 2^dim MultiIndex hypercube for generating the permulations
  MultiIndex<Real>::size_type shape(_data.dim(), 2);
  MultiIndex<Real> hypercube = MultiIndex<Real>(shape);
  T interpolated_value = 0;
  for (const auto & p : hypercube)
  {
    // this loop computes the flat index for the point in the
    // hypercube and accumulates the total weight for this point
    // the selection of the dimensional weight is based on the swapping
    // done in the pre-computing loop before
    unsigned int index = flat_index;
    T w = 1;
    for (unsigned int j = 0; j < p.first.size(); ++j)
    {
      index += stride[j] * p.first[j];
      w *= weights[j][p.first[j]];
    }
 
    // add the contribution of this base point to the interpolated value
    interpolated_value += w * _data[index];
  }
 
  return interpolated_value / volume;
}

setData()

template<typename T >

void MultiDimensionalInterpolationTempl< T >::setData	(	const std::vector< std::vector< Real > > &	base_points,
		const MultiIndex< Real > &	data
	)

sets data but also fixes degenerate dimensions in data

Definition at line 32 of file MultiDimensionalInterpolation.C.

{
  _degenerate_interpolation = false;
  _setup_complete = false;
 
  // stash away original dimension of data because we might modify it
  _original_dim = data.dim();
  _degenerate_index.resize(_original_dim, false);
 
  // need to check this here to make sure the next loop is safe
  if (data.dim() != base_points.size())
  {
    std::ostringstream oss;
    oss << "Leading dimension of base_points " << base_points.size() << " and data dimensionality "
        << data.dim() << " are inconsistent.";
    throw std::domain_error(oss.str());
  }
 
  // first a corner case is addressed; if the size of the MultiIndex in one
  // dimension is 1, then finding the interpolation "stencil" becomes much
  // harder. It is easier to slice the data array and keep track of which values
  // to discard in the interpolation routines
  MultiIndex<Real>::size_type slice_dimensions;
  MultiIndex<Real>::size_type slice_indices;
  MultiIndex<Real>::size_type size = data.size();
  std::vector<std::vector<Real>> modified_base_points;
  for (unsigned int j = 0; j < data.dim(); ++j)
    if (size[j] == 1)
    {
      slice_dimensions.push_back(j);
      slice_indices.push_back(0);
      _degenerate_index[j] = true;
    }
    else
      modified_base_points.push_back(base_points[j]);
 
  // if all dimensions have size 1, then there is only one value to return
  // this is a weird corner case that we should support to make life easier
  // for whoever uses this object
  // In this case, we can just use the passed base_points and data objects
  // and short-cut the computation in the interpolation function
  if (modified_base_points.size() == 0)
  {
    _degenerate_interpolation = true;
    _base_points = base_points;
    _data = data;
  }
  else
  {
    // set the data
    _base_points = modified_base_points;
    _data = data.slice(slice_dimensions, slice_indices);
  }
  errorCheck();
}

Member Data Documentation

_base_points

template<typename T >

std::vector<std::vector<Real> > MultiDimensionalInterpolationTempl< T >::_base_points

private

Definition at line 81 of file MultiDimensionalInterpolation.h.

_data

template<typename T >

MultiIndex<Real> MultiDimensionalInterpolationTempl< T >::_data

private

Definition at line 82 of file MultiDimensionalInterpolation.h.

_degenerate_index

template<typename T >

std::vector<bool> MultiDimensionalInterpolationTempl< T >::_degenerate_index

protected

this variable keeps track on which dimension is degenerate and was removed

Definition at line 71 of file MultiDimensionalInterpolation.h.

_degenerate_interpolation

template<typename T >

bool MultiDimensionalInterpolationTempl< T >::_degenerate_interpolation = false

protected

if all dimensions have size one then there is only one value to return this corner case should be supported so the calling routine does not need to check for it

Definition at line 78 of file MultiDimensionalInterpolation.h.

_original_dim

template<typename T >

unsigned int MultiDimensionalInterpolationTempl< T >::_original_dim

protected

original dimension is to allow checks on user inputs for cases where arrays are sliced

Definition at line 68 of file MultiDimensionalInterpolation.h.

Referenced by MultiDimensionalInterpolationTempl< T >::dim().

_setup_complete

template<typename T >

bool MultiDimensionalInterpolationTempl< T >::_setup_complete = false

private

Definition at line 83 of file MultiDimensionalInterpolation.h.

---

The documentation for this class was generated from the following files:

• include/utils/MultiDimensionalInterpolation.h
• src/utils/MultiDimensionalInterpolation.C