lsst::afw::math::GaussianProcess< T > Class Template Reference

Stores values of a function sampled on an image and allows you to interpolate the function to unsampled points. More...

#include <GaussianProcess.h>

Public Member Functions
	GaussianProcess (const GaussianProcess &)=delete
GaussianProcess &	operator= (const GaussianProcess &)=delete
	GaussianProcess (GaussianProcess &&)=delete
GaussianProcess &	operator= (GaussianProcess &&)=delete
	GaussianProcess (ndarray::Array< T, 2, 2 > const &dataIn, ndarray::Array< T, 1, 1 > const &ff, std::shared_ptr< Covariogram< T > > const &covarIn)
	This is the constructor you call if you do not wish to normalize the positions of your data points and you have only one function.
	GaussianProcess (ndarray::Array< T, 2, 2 > const &dataIn, ndarray::Array< T, 1, 1 > const &mn, ndarray::Array< T, 1, 1 > const &mx, ndarray::Array< T, 1, 1 > const &ff, std::shared_ptr< Covariogram< T > > const &covarIn)
	This is the constructor you call if you want the positions of your data points normalized by the span of each dimension and you have only one function.
	GaussianProcess (ndarray::Array< T, 2, 2 > const &dataIn, ndarray::Array< T, 2, 2 > const &ff, std::shared_ptr< Covariogram< T > > const &covarIn)
	this is the constructor to use in the case of a vector of input functions and an unbounded/unnormalized parameter space
	GaussianProcess (ndarray::Array< T, 2, 2 > const &dataIn, ndarray::Array< T, 1, 1 > const &mn, ndarray::Array< T, 1, 1 > const &mx, ndarray::Array< T, 2, 2 > const &ff, std::shared_ptr< Covariogram< T > > const &covarIn)
	this is the constructor to use in the case of a vector of input functions using minima and maxima in parameter space
int	getNPoints () const
	return the number of data points stored in the GaussianProcess
int	getDim () const
	return the dimensionality of data points stored in the GaussianProcess
void	getData (ndarray::Array< T, 2, 2 > pts, ndarray::Array< T, 1, 1 > fn, ndarray::Array< int, 1, 1 > indices) const
	Return a sub-sample the data underlying the Gaussian Process.
void	getData (ndarray::Array< T, 2, 2 > pts, ndarray::Array< T, 2, 2 > fn, ndarray::Array< int, 1, 1 > indices) const
	Return a sub-sample the data underlying the Gaussian Process.
T	interpolate (ndarray::Array< T, 1, 1 > variance, ndarray::Array< T, 1, 1 > const &vin, int numberOfNeighbors) const
	Interpolate the function value at one point using a specified number of nearest neighbors.
void	interpolate (ndarray::Array< T, 1, 1 > mu, ndarray::Array< T, 1, 1 > variance, ndarray::Array< T, 1, 1 > const &vin, int numberOfNeighbors) const
	This is the version of GaussianProcess::interpolate for a vector of functions.
T	selfInterpolate (ndarray::Array< T, 1, 1 > variance, int dex, int numberOfNeighbors) const
	This method will interpolate the function on a data point for purposes of optimizing hyper parameters.
void	selfInterpolate (ndarray::Array< T, 1, 1 > mu, ndarray::Array< T, 1, 1 > variance, int dex, int numberOfNeighbors) const
	The version of selfInterpolate called for a vector of functions.
void	batchInterpolate (ndarray::Array< T, 1, 1 > mu, ndarray::Array< T, 1, 1 > variance, ndarray::Array< T, 2, 2 > const &queries) const
	Interpolate a list of query points using all of the input data (rather than nearest neighbors)
void	batchInterpolate (ndarray::Array< T, 1, 1 > mu, ndarray::Array< T, 2, 2 > const &queries) const
	Interpolate a list of points using all of the data.
void	batchInterpolate (ndarray::Array< T, 2, 2 > mu, ndarray::Array< T, 2, 2 > variance, ndarray::Array< T, 2, 2 > const &queries) const
	This is the version of batchInterpolate (with variances) that is called for a vector of functions.
void	batchInterpolate (ndarray::Array< T, 2, 2 > mu, ndarray::Array< T, 2, 2 > const &queries) const
	This is the version of batchInterpolate (without variances) that is called for a vector of functions.
void	addPoint (ndarray::Array< T, 1, 1 > const &vin, T f)
	Add a point to the pool of data used by GaussianProcess for interpolation.
void	addPoint (ndarray::Array< T, 1, 1 > const &vin, ndarray::Array< T, 1, 1 > const &f)
	This is the version of addPoint that is called for a vector of functions.
void	removePoint (int dex)
	This will remove a point from the data set.
void	setKrigingParameter (T kk)
	Assign a value to the Kriging paramter.
void	setCovariogram (std::shared_ptr< Covariogram< T > > const &covar)
	Assign a different covariogram to this GaussianProcess.
void	setLambda (T lambda)
	set the value of the hyperparameter _lambda
GaussianProcessTimer &	getTimes () const
	Give the user acces to _timer, an object keeping track of the time spent on various processes within interpolate.

Detailed Description

template<typename T>
class lsst::afw::math::GaussianProcess< T >

Stores values of a function sampled on an image and allows you to interpolate the function to unsampled points.

The data will be stored in a KD Tree for easy nearest neighbor searching when interpolating.

The array _function[] will contain the values of the function being interpolated. You can provide a two dimensional array _function[][] if you wish to interpolate a vector of functions. In this case _function[i][j] is the jth function associated with the ith data point. Note: presently, the covariance matrices do not relate elements of _function[i][] to each other, so the variances returned will be identical for all functions evaluated at the same point in parameter space.

_data[i][j] will be the jth component of the ith data point.

_max and _min contain the maximum and minimum values of each dimension in parameter space (if applicable) so that data points can be normalized by _max-_min to keep distances between points reasonable. This is an option specified by calling the relevant constructor.

Definition at line 471 of file GaussianProcess.h.

Constructor & Destructor Documentation

GaussianProcess() [1/6]

template<typename T >

lsst::afw::math::GaussianProcess< T >::GaussianProcess ( const GaussianProcess< T > & )

delete

GaussianProcess() [2/6]

template<typename T >

lsst::afw::math::GaussianProcess< T >::GaussianProcess ( GaussianProcess< T > && )

delete

GaussianProcess() [3/6]

template<typename T >

lsst::afw::math::GaussianProcess< T >::GaussianProcess	(	ndarray::Array< T, 2, 2 > const &	dataIn,
		ndarray::Array< T, 1, 1 > const &	ff,
		std::shared_ptr< Covariogram< T > > const &	covarIn )

This is the constructor you call if you do not wish to normalize the positions of your data points and you have only one function.

Parameters

[in]	dataIn	an ndarray containing the data points; the ith row of datain is the ith data point
[in]	ff	a one-dimensional ndarray containing the values of the scalar function associated with each data point. This is the function you are interpolating
[in]	covarIn	is the input covariogram

Definition at line 751 of file GaussianProcess.cc.

{
    int i;
 
    _covariogram = covarIn;
 
    _npts = dataIn.template getSize<0>();
    _dimensions = dataIn.template getSize<1>();
 
    _room = _npts;
    _roomStep = 5000;
 
    _nFunctions = 1;
    _function = allocate(ndarray::makeVector(_npts, 1));
 
    if (ff.getNumElements() != static_cast<ndarray::Size>(_npts)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You did not pass in the same number of data points as function values\n");
    }
 
    for (i = 0; i < _npts; i++) _function[i][0] = ff[i];
 
    _krigingParameter = T(1.0);
    _lambda = T(1.0e-5);
 
    _useMaxMin = 0;
 
    _kdTree.Initialize(dataIn);
 
    _npts = _kdTree.getNPoints();
}

GaussianProcess() [4/6]

template<typename T >

lsst::afw::math::GaussianProcess< T >::GaussianProcess	(	ndarray::Array< T, 2, 2 > const &	dataIn,
		ndarray::Array< T, 1, 1 > const &	mn,
		ndarray::Array< T, 1, 1 > const &	mx,
		ndarray::Array< T, 1, 1 > const &	ff,
		std::shared_ptr< Covariogram< T > > const &	covarIn )

This is the constructor you call if you want the positions of your data points normalized by the span of each dimension and you have only one function.

Parameters

[in]	dataIn	an ndarray containing the data points; the ith row of datain is the ith data point
[in]	mn	a one-dimensional ndarray containing the minimum values of each dimension (for normalizing the positions of data points)
[in]	mx	a one-dimensional ndarray containing the maximum values of each dimension (for normalizing the positions of data points)
[in]	ff	a one-dimensional ndarray containing the values of the scalar function associated with each data point. This is the function you are interpolating
[in]	covarIn	is the input covariogram

Note: the member variable _useMaxMin will allow the code to remember which constructor you invoked

Definition at line 786 of file GaussianProcess.cc.

                                                                                   {
    int i, j;
    ndarray::Array<T, 2, 2> normalizedData;
 
    _covariogram = covarIn;
 
    _npts = dataIn.template getSize<0>();
    _dimensions = dataIn.template getSize<1>();
    _room = _npts;
    _roomStep = 5000;
 
    if (ff.getNumElements() != static_cast<ndarray::Size>(_npts)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You did not pass in the same number of data points as function values\n");
    }
 
    if (mn.getNumElements() != static_cast<ndarray::Size>(_dimensions) ||
        mx.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your min/max values have different dimensionality than your data points\n");
    }
 
    _krigingParameter = T(1.0);
 
    _lambda = T(1.0e-5);
    _krigingParameter = T(1.0);
 
    _max = allocate(ndarray::makeVector(_dimensions));
    _min = allocate(ndarray::makeVector(_dimensions));
    _max.deep() = mx;
    _min.deep() = mn;
    _useMaxMin = 1;
    normalizedData = allocate(ndarray::makeVector(_npts, _dimensions));
    for (i = 0; i < _npts; i++) {
        for (j = 0; j < _dimensions; j++) {
            normalizedData[i][j] = (dataIn[i][j] - _min[j]) / (_max[j] - _min[j]);
            // note the normalization by _max - _min in each dimension
        }
    }
 
    _kdTree.Initialize(normalizedData);
 
    _npts = _kdTree.getNPoints();
    _nFunctions = 1;
    _function = allocate(ndarray::makeVector(_npts, 1));
    for (i = 0; i < _npts; i++) _function[i][0] = ff[i];
}

GaussianProcess() [5/6]

template<typename T >

lsst::afw::math::GaussianProcess< T >::GaussianProcess	(	ndarray::Array< T, 2, 2 > const &	dataIn,
		ndarray::Array< T, 2, 2 > const &	ff,
		std::shared_ptr< Covariogram< T > > const &	covarIn )

this is the constructor to use in the case of a vector of input functions and an unbounded/unnormalized parameter space

Parameters

[in]	dataIn	contains the data points, as in other constructors
[in]	ff	contains the functions. Each row of ff corresponds to a datapoint. Each column corresponds to a function (ff[i][j] is the jth function associated with the ith data point)
[in]	covarIn	is the input covariogram

Definition at line 837 of file GaussianProcess.cc.

                                                                                   {
    _covariogram = covarIn;
 
    _npts = dataIn.template getSize<0>();
    _dimensions = dataIn.template getSize<1>();
 
    _room = _npts;
    _roomStep = 5000;
 
    if (ff.template getSize<0>() != static_cast<ndarray::Size>(_npts)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You did not pass in the same number of data points as function values\n");
    }
 
    _nFunctions = ff.template getSize<1>();
    _function = allocate(ndarray::makeVector(_npts, _nFunctions));
    _function.deep() = ff;
 
    _krigingParameter = T(1.0);
 
    _lambda = T(1.0e-5);
 
    _useMaxMin = 0;
 
    _kdTree.Initialize(dataIn);
 
    _npts = _kdTree.getNPoints();
}

GaussianProcess() [6/6]

template<typename T >

lsst::afw::math::GaussianProcess< T >::GaussianProcess	(	ndarray::Array< T, 2, 2 > const &	dataIn,
		ndarray::Array< T, 1, 1 > const &	mn,
		ndarray::Array< T, 1, 1 > const &	mx,
		ndarray::Array< T, 2, 2 > const &	ff,
		std::shared_ptr< Covariogram< T > > const &	covarIn )

this is the constructor to use in the case of a vector of input functions using minima and maxima in parameter space

Parameters

[in]	dataIn	contains the data points, as in other constructors
[in]	mn	contains the minimum allowed values of the parameters in parameter space
[in]	mx	contains the maximum allowed values of the parameters in parameter space
[in]	ff	contains the functions. Each row of ff corresponds to a datapoint. Each column corresponds to a function (ff[i][j] is the jth function associated with the ith data point)
[in]	covarIn	is the input covariogram

Definition at line 868 of file GaussianProcess.cc.

                                                                                   {
    int i, j;
    ndarray::Array<T, 2, 2> normalizedData;
 
    _covariogram = covarIn;
 
    _npts = dataIn.template getSize<0>();
    _dimensions = dataIn.template getSize<1>();
 
    _room = _npts;
    _roomStep = 5000;
 
    if (ff.template getSize<0>() != static_cast<ndarray::Size>(_npts)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You did not pass in the same number of data points as function values\n");
    }
 
    if (mn.getNumElements() != static_cast<ndarray::Size>(_dimensions) ||
        mx.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your min/max values have different dimensionality than your data points\n");
    }
 
    _krigingParameter = T(1.0);
 
    _lambda = T(1.0e-5);
    _krigingParameter = T(1.0);
 
    _max = allocate(ndarray::makeVector(_dimensions));
    _min = allocate(ndarray::makeVector(_dimensions));
    _max.deep() = mx;
    _min.deep() = mn;
    _useMaxMin = 1;
    normalizedData = allocate(ndarray::makeVector(_npts, _dimensions));
    for (i = 0; i < _npts; i++) {
        for (j = 0; j < _dimensions; j++) {
            normalizedData[i][j] = (dataIn[i][j] - _min[j]) / (_max[j] - _min[j]);
            // note the normalization by _max - _min in each dimension
        }
    }
 
    _kdTree.Initialize(normalizedData);
    _npts = _kdTree.getNPoints();
    _nFunctions = ff.template getSize<1>();
    _function = allocate(ndarray::makeVector(_npts, _nFunctions));
    _function.deep() = ff;
}

Member Function Documentation

addPoint() [1/2]

template<typename T >

void lsst::afw::math::GaussianProcess< T >::addPoint	(	ndarray::Array< T, 1, 1 > const &	vin,
		ndarray::Array< T, 1, 1 > const &	f )

This is the version of addPoint that is called for a vector of functions.

Exceptions

pex::exceptions::RuntimeError

if the tree does not end up properly constructed (the exception is actually thrown by KdTree<T>::addPoint() )

Note: excessive use of addPoint and removePoint can result in an unbalanced KdTree, which will slow down nearest neighbor searches

Definition at line 1907 of file GaussianProcess.cc.

                                                                                                  {
    if (vin.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You are trying to add a point of the wrong dimensionality to "
                          "your GaussianProcess.\n");
    }
 
    if (f.template getSize<0>() != static_cast<ndarray::Size>(_nFunctions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You are not adding the correct number of function values to "
                          "your GaussianProcess.\n");
    }
 
    int i, j;
 
    ndarray::Array<T, 1, 1> v;
    v = allocate(ndarray::makeVector(_dimensions));
 
    for (i = 0; i < _dimensions; i++) {
        v[i] = vin[i];
        if (_useMaxMin == 1) {
            v[i] = (v[i] - _min[i]) / (_max[i] - _min[i]);
        }
    }
 
    if (_npts == _room) {
        ndarray::Array<T, 2, 2> buff;
        buff = allocate(ndarray::makeVector(_npts, _nFunctions));
        buff.deep() = _function;
 
        _room += _roomStep;
        _function = allocate(ndarray::makeVector(_room, _nFunctions));
        for (i = 0; i < _npts; i++) {
            for (j = 0; j < _nFunctions; j++) {
                _function[i][j] = buff[i][j];
            }
        }
    }
    for (i = 0; i < _nFunctions; i++) _function[_npts][i] = f[i];
 
    _kdTree.addPoint(v);
    _npts = _kdTree.getNPoints();
}

addPoint() [2/2]

template<typename T >

void lsst::afw::math::GaussianProcess< T >::addPoint	(	ndarray::Array< T, 1, 1 > const &	vin,
		T	f )

Add a point to the pool of data used by GaussianProcess for interpolation.

Parameters

[in]	vin	a one-dimensional ndarray storing the point in parameter space that you are adding
[in]	f	the value of the function at that point

Exceptions

pex::exceptions::RuntimeError	if you call this when you should have called the version taking a vector of functions (below)
pex::exceptions::RuntimeError	if the tree does not end up properly constructed (the exception is actually thrown by KdTree<T>::addPoint() )

Note: excessive use of addPoint and removePoint can result in an unbalanced KdTree, which will slow down nearest neighbor searches

Definition at line 1862 of file GaussianProcess.cc.

                                                                       {
    int i, j;
 
    if (_nFunctions != 1) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You are calling the wrong addPoint; you need a vector of functions\n");
    }
 
    if (vin.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You are trying to add a point of the wrong dimensionality to "
                          "your GaussianProcess.\n");
    }
 
    ndarray::Array<T, 1, 1> v;
    v = allocate(ndarray::makeVector(_dimensions));
 
    for (i = 0; i < _dimensions; i++) {
        v[i] = vin[i];
        if (_useMaxMin == 1) {
            v[i] = (v[i] - _min[i]) / (_max[i] - _min[i]);
        }
    }
 
    if (_npts == _room) {
        ndarray::Array<T, 2, 2> buff;
        buff = allocate(ndarray::makeVector(_npts, _nFunctions));
        buff.deep() = _function;
 
        _room += _roomStep;
        _function = allocate(ndarray::makeVector(_room, _nFunctions));
        for (i = 0; i < _npts; i++) {
            for (j = 0; j < _nFunctions; j++) {
                _function[i][j] = buff[i][j];
            }
        }
    }
 
    _function[_npts][0] = f;
 
    _kdTree.addPoint(v);
    _npts = _kdTree.getNPoints();
}

batchInterpolate() [1/4]

template<typename T >

void lsst::afw::math::GaussianProcess< T >::batchInterpolate	(	ndarray::Array< T, 1, 1 >	mu,
		ndarray::Array< T, 1, 1 >	variance,
		ndarray::Array< T, 2, 2 > const &	queries ) const

Interpolate a list of query points using all of the input data (rather than nearest neighbors)

Parameters

[out]	mu	a 1-dimensional ndarray where the interpolated function values will be stored
[out]	variance	a 1-dimensional ndarray where the corresponding variances in the function value will be stored
[in]	queries	a 2-dimensional ndarray containing the points to be interpolated. queries[i][j] is the jth component of the ith point

This method will attempt to construct a _npts X _npts covariance matrix C and solve the problem Cx=b. Be wary of using it in the case where _npts is very large.

This version of the method will also return variances for all of the query points. That is a very time consuming calculation relative to just returning estimates for the function. Consider calling the version of this method that does not calculate variances (below). The difference in time spent is an order of magnitude for 189 data points and 1,000,000 interpolations.

Definition at line 1488 of file GaussianProcess.cc.

                                                                                      {
    int i, j;
 
    ndarray::Size nQueries = queries.template getSize<0>();
 
    if (_nFunctions != 1) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your mu and variance arrays do not have room for all of the functions "
                          "as you are trying to interpolate\n");
    }
 
    if (mu.getNumElements() != nQueries || variance.getNumElements() != nQueries) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your mu and variance arrays do not have room for all of the points "
                          "at which you are trying to interpolate your function.\n");
    }
 
    if (queries.template getSize<1>() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "The points you passed to batchInterpolate are of the wrong "
                          "dimensionality for your Gaussian Process\n");
    }
 
    T fbar;
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> batchCovariance, batchbb, batchxx;
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> queryCovariance;
    Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
 
    ndarray::Array<T, 1, 1> v1;
 
    _timer.start();
 
    v1 = allocate(ndarray::makeVector(_dimensions));
    batchbb.resize(_npts, 1);
    batchxx.resize(_npts, 1);
    batchCovariance.resize(_npts, _npts);
    queryCovariance.resize(_npts, 1);
 
    for (i = 0; i < _npts; i++) {
        batchCovariance(i, i) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(i)) + _lambda;
        for (j = i + 1; j < _npts; j++) {
            batchCovariance(i, j) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(j));
            batchCovariance(j, i) = batchCovariance(i, j);
        }
    }
    _timer.addToIteration();
 
    ldlt.compute(batchCovariance);
 
    fbar = 0.0;
    for (i = 0; i < _npts; i++) {
        fbar += _function[i][0];
    }
    fbar = fbar / T(_npts);
 
    for (i = 0; i < _npts; i++) {
        batchbb(i, 0) = _function[i][0] - fbar;
    }
    batchxx = ldlt.solve(batchbb);
    _timer.addToEigen();
 
    for (ndarray::Size ii = 0; ii < nQueries; ii++) {
        for (i = 0; i < _dimensions; i++) v1[i] = queries[ii][i];
        if (_useMaxMin == 1) {
            for (i = 0; i < _dimensions; i++) v1[i] = (v1[i] - _min[i]) / (_max[i] - _min[i]);
        }
        mu(ii) = fbar;
        for (i = 0; i < _npts; i++) {
            mu(ii) += batchxx(i) * (*_covariogram)(v1, _kdTree.getData(i));
        }
    }
    _timer.addToIteration();
 
    for (ndarray::Size ii = 0; ii < nQueries; ii++) {
        for (i = 0; i < _dimensions; i++) v1[i] = queries[ii][i];
        if (_useMaxMin == 1) {
            for (i = 0; i < _dimensions; i++) v1[i] = (v1[i] - _min[i]) / (_max[i] - _min[i]);
        }
 
        for (i = 0; i < _npts; i++) {
            batchbb(i, 0) = (*_covariogram)(v1, _kdTree.getData(i));
            queryCovariance(i, 0) = batchbb(i, 0);
        }
        batchxx = ldlt.solve(batchbb);
 
        variance[ii] = (*_covariogram)(v1, v1) + _lambda;
 
        for (i = 0; i < _npts; i++) {
            variance[ii] -= queryCovariance(i, 0) * batchxx(i);
        }
 
        variance[ii] = variance[ii] * _krigingParameter;
    }
 
    _timer.addToVariance();
    _timer.addToTotal(nQueries);
}

batchInterpolate() [2/4]

template<typename T >

void lsst::afw::math::GaussianProcess< T >::batchInterpolate	(	ndarray::Array< T, 1, 1 >	mu,
		ndarray::Array< T, 2, 2 > const &	queries ) const

Interpolate a list of points using all of the data.

Do not return variances for the interpolation.

Parameters

[out]	mu	a 1-dimensional ndarray where the interpolated function values will be stored
[in]	queries	a 2-dimensional ndarray containing the points to be interpolated. queries[i][j] is the jth component of the ith point

This method will attempt to construct a _npts X _npts covariance matrix C and solve the problem Cx=b. Be wary of using it in the case where _npts is very large.

This version of the method does not return variances. It is an order of magnitude faster than the version of the method that does return variances (timing done on a case with 189 data points and 1 million query points).

Definition at line 1696 of file GaussianProcess.cc.

                                                                                      {
    int i, j;
 
    ndarray::Size nQueries = queries.template getSize<0>();
 
    if (_nFunctions != 1) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your output array does not have enough room for all of the functions "
                          "you are trying to interpolate.\n");
    }
 
    if (queries.template getSize<1>() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "The points at which you are trying to interpolate your function are "
                          "of the wrong dimensionality.\n");
    }
 
    if (mu.getNumElements() != nQueries) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your output array does not have enough room for all of the points "
                          "at which you are trying to interpolate your function.\n");
    }
 
    T fbar;
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> batchCovariance, batchbb, batchxx;
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> queryCovariance;
    Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
 
    ndarray::Array<T, 1, 1> v1;
 
    _timer.start();
 
    v1 = allocate(ndarray::makeVector(_dimensions));
 
    batchbb.resize(_npts, 1);
    batchxx.resize(_npts, 1);
    batchCovariance.resize(_npts, _npts);
    queryCovariance.resize(_npts, 1);
 
    for (i = 0; i < _npts; i++) {
        batchCovariance(i, i) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(i)) + _lambda;
        for (j = i + 1; j < _npts; j++) {
            batchCovariance(i, j) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(j));
            batchCovariance(j, i) = batchCovariance(i, j);
        }
    }
    _timer.addToIteration();
 
    ldlt.compute(batchCovariance);
 
    fbar = 0.0;
    for (i = 0; i < _npts; i++) {
        fbar += _function[i][0];
    }
    fbar = fbar / T(_npts);
 
    for (i = 0; i < _npts; i++) {
        batchbb(i, 0) = _function[i][0] - fbar;
    }
    batchxx = ldlt.solve(batchbb);
    _timer.addToEigen();
 
    for (ndarray::Size ii = 0; ii < nQueries; ii++) {
        for (i = 0; i < _dimensions; i++) v1[i] = queries[ii][i];
        if (_useMaxMin == 1) {
            for (i = 0; i < _dimensions; i++) v1[i] = (v1[i] - _min[i]) / (_max[i] - _min[i]);
        }
 
        mu(ii) = fbar;
        for (i = 0; i < _npts; i++) {
            mu(ii) += batchxx(i) * (*_covariogram)(v1, _kdTree.getData(i));
        }
    }
    _timer.addToIteration();
    _timer.addToTotal(nQueries);
}

batchInterpolate() [3/4]

template<typename T >

void lsst::afw::math::GaussianProcess< T >::batchInterpolate	(	ndarray::Array< T, 2, 2 >	mu,
		ndarray::Array< T, 2, 2 > const &	queries ) const

This is the version of batchInterpolate (without variances) that is called for a vector of functions.

Definition at line 1775 of file GaussianProcess.cc.

                                                                                      {
    int i, j, ifn;
 
    ndarray::Size nQueries = queries.template getSize<0>();
 
    if (mu.template getSize<0>() != nQueries) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your output array does not have enough room for all of the points "
                          "at which you want to interpolate your functions.\n");
    }
 
    if (mu.template getSize<1>() != static_cast<ndarray::Size>(_nFunctions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your output array does not have enough room for all of the functions "
                          "you are trying to interpolate.\n");
    }
 
    if (queries.template getSize<1>() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "The points at which you are interpolating your functions do not "
                          "have the correct dimensionality.\n");
    }
 
    T fbar;
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> batchCovariance, batchbb, batchxx;
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> queryCovariance;
    Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
 
    ndarray::Array<T, 1, 1> v1;
 
    _timer.start();
 
    v1 = allocate(ndarray::makeVector(_dimensions));
 
    batchbb.resize(_npts, 1);
    batchxx.resize(_npts, 1);
    batchCovariance.resize(_npts, _npts);
    queryCovariance.resize(_npts, 1);
 
    for (i = 0; i < _npts; i++) {
        batchCovariance(i, i) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(i)) + _lambda;
        for (j = i + 1; j < _npts; j++) {
            batchCovariance(i, j) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(j));
            batchCovariance(j, i) = batchCovariance(i, j);
        }
    }
 
    _timer.addToIteration();
 
    ldlt.compute(batchCovariance);
 
    _timer.addToEigen();
 
    for (ifn = 0; ifn < _nFunctions; ifn++) {
        fbar = 0.0;
        for (i = 0; i < _npts; i++) {
            fbar += _function[i][ifn];
        }
        fbar = fbar / T(_npts);
 
        _timer.addToIteration();
 
        for (i = 0; i < _npts; i++) {
            batchbb(i, 0) = _function[i][ifn] - fbar;
        }
        batchxx = ldlt.solve(batchbb);
        _timer.addToEigen();
 
        for (ndarray::Size ii = 0; ii < nQueries; ii++) {
            for (i = 0; i < _dimensions; i++) v1[i] = queries[ii][i];
            if (_useMaxMin == 1) {
                for (i = 0; i < _dimensions; i++) v1[i] = (v1[i] - _min[i]) / (_max[i] - _min[i]);
            }
 
            mu[ii][ifn] = fbar;
            for (i = 0; i < _npts; i++) {
                mu[ii][ifn] += batchxx(i) * (*_covariogram)(v1, _kdTree.getData(i));
            }
        }
 
    }  // ifn = 0 through _nFunctions
 
    _timer.addToTotal(nQueries);
}

batchInterpolate() [4/4]

template<typename T >

void lsst::afw::math::GaussianProcess< T >::batchInterpolate	(	ndarray::Array< T, 2, 2 >	mu,
		ndarray::Array< T, 2, 2 >	variance,
		ndarray::Array< T, 2, 2 > const &	queries ) const

This is the version of batchInterpolate (with variances) that is called for a vector of functions.

Definition at line 1588 of file GaussianProcess.cc.

                                                                                      {
    int i, j, ifn;
 
    ndarray::Size nQueries = queries.template getSize<0>();
 
    if (mu.template getSize<0>() != nQueries || variance.template getSize<0>() != nQueries) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your output arrays do not have room for all of the points at which "
                          "you are interpolating your functions.\n");
    }
 
    if (mu.template getSize<1>() != static_cast<ndarray::Size>(_nFunctions) ||
        variance.template getSize<1>() != static_cast<ndarray::Size>(_nFunctions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your output arrays do not have room for all of the functions you are "
                          "interpolating\n");
    }
 
    if (queries.template getSize<1>() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "The points at which you are interpolating your functions have the "
                          "wrong dimensionality.\n");
    }
 
    T fbar;
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> batchCovariance, batchbb, batchxx;
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> queryCovariance;
    Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
 
    ndarray::Array<T, 1, 1> v1;
 
    _timer.start();
 
    v1 = allocate(ndarray::makeVector(_dimensions));
    batchbb.resize(_npts, 1);
    batchxx.resize(_npts, 1);
    batchCovariance.resize(_npts, _npts);
    queryCovariance.resize(_npts, 1);
 
    for (i = 0; i < _npts; i++) {
        batchCovariance(i, i) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(i)) + _lambda;
        for (j = i + 1; j < _npts; j++) {
            batchCovariance(i, j) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(j));
            batchCovariance(j, i) = batchCovariance(i, j);
        }
    }
 
    _timer.addToIteration();
 
    ldlt.compute(batchCovariance);
 
    _timer.addToEigen();
 
    for (ifn = 0; ifn < _nFunctions; ifn++) {
        fbar = 0.0;
        for (i = 0; i < _npts; i++) {
            fbar += _function[i][ifn];
        }
        fbar = fbar / T(_npts);
        _timer.addToIteration();
 
        for (i = 0; i < _npts; i++) {
            batchbb(i, 0) = _function[i][ifn] - fbar;
        }
        batchxx = ldlt.solve(batchbb);
        _timer.addToEigen();
 
        for (ndarray::Size ii = 0; ii < nQueries; ii++) {
            for (i = 0; i < _dimensions; i++) v1[i] = queries[ii][i];
            if (_useMaxMin == 1) {
                for (i = 0; i < _dimensions; i++) v1[i] = (v1[i] - _min[i]) / (_max[i] - _min[i]);
            }
            mu[ii][ifn] = fbar;
            for (i = 0; i < _npts; i++) {
                mu[ii][ifn] += batchxx(i) * (*_covariogram)(v1, _kdTree.getData(i));
            }
        }
 
    }  // ifn = 0 to _nFunctions
 
    _timer.addToIteration();
    for (ndarray::Size ii = 0; ii < nQueries; ii++) {
        for (i = 0; i < _dimensions; i++) v1[i] = queries[ii][i];
        if (_useMaxMin == 1) {
            for (i = 0; i < _dimensions; i++) v1[i] = (v1[i] - _min[i]) / (_max[i] - _min[i]);
        }
 
        for (i = 0; i < _npts; i++) {
            batchbb(i, 0) = (*_covariogram)(v1, _kdTree.getData(i));
            queryCovariance(i, 0) = batchbb(i, 0);
        }
        batchxx = ldlt.solve(batchbb);
 
        variance[ii][0] = (*_covariogram)(v1, v1) + _lambda;
 
        for (i = 0; i < _npts; i++) {
            variance[ii][0] -= queryCovariance(i, 0) * batchxx(i);
        }
 
        variance[ii][0] = variance[ii][0] * _krigingParameter;
        for (i = 1; i < _nFunctions; i++) variance[ii][i] = variance[ii][0];
    }
    _timer.addToVariance();
    _timer.addToTotal(nQueries);
}

getData() [1/2]

template<typename T >

void lsst::afw::math::GaussianProcess< T >::getData	(	ndarray::Array< T, 2, 2 >	pts,
		ndarray::Array< T, 1, 1 >	fn,
		ndarray::Array< int, 1, 1 >	indices ) const

Return a sub-sample the data underlying the Gaussian Process.

Parameters

[out]	pts	will contain the data points from the Gaussian Process
[out]	fn	will contain the function values from the Gaussian Process
[in]	indices	is an array of indices indicating the points to return

Definition at line 929 of file GaussianProcess.cc.

                                                                        {
    if (_nFunctions != 1) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your function value array does not have enough room for all of the functions "
                          "in your GaussianProcess.\n");
    }
 
    if (pts.template getSize<1>() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your pts array is constructed for points of the wrong dimensionality.\n");
    }
 
    if (pts.template getSize<0>() != indices.getNumElements()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You did not put enough room in your pts array to fit all the points "
                          "you asked for in your indices array.\n");
    }
 
    if (fn.template getSize<0>() != indices.getNumElements()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You did not provide enough room in your function value array "
                          "for all of the points you requested in your indices array.\n");
    }
 
    for (ndarray::Size i = 0; i < indices.template getSize<0>(); i++) {
        pts[i] = _kdTree.getData(indices[i]);  // do this first in case one of the indices is invalid.
                                               // _kdTree.getData() will raise an exception in that case
        fn[i] = _function[indices[i]][0];
        if (_useMaxMin == 1) {
            for (int j = 0; j < _dimensions; j++) {
                pts[i][j] *= (_max[j] - _min[j]);
                pts[i][j] += _min[j];
            }
        }
    }
}

getData() [2/2]

template<typename T >

void lsst::afw::math::GaussianProcess< T >::getData	(	ndarray::Array< T, 2, 2 >	pts,
		ndarray::Array< T, 2, 2 >	fn,
		ndarray::Array< int, 1, 1 >	indices ) const

Return a sub-sample the data underlying the Gaussian Process.

Parameters

[out]	pts	will contain the data points from the Gaussian Process
[out]	fn	will contain the function values from the Gaussian Process
[in]	indices	is an array of indices indicating the points to return

Definition at line 968 of file GaussianProcess.cc.

                                                                        {
    if (fn.template getSize<1>() != static_cast<ndarray::Size>(_nFunctions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your function value array does not have enough room for all of the functions "
                          "in your GaussianProcess.\n");
    }
 
    if (pts.template getSize<1>() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your pts array is constructed for points of the wrong dimensionality.\n");
    }
 
    if (pts.template getSize<0>() != indices.getNumElements()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You did not put enough room in your pts array to fit all the points "
                          "you asked for in your indices array.\n");
    }
 
    if (fn.template getSize<0>() != indices.getNumElements()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You did not provide enough room in your function value array "
                          "for all of the points you requested in your indices array.\n");
    }
 
    for (ndarray::Size i = 0; i < indices.template getSize<0>(); i++) {
        pts[i] = _kdTree.getData(indices[i]);  // do this first in case one of the indices is invalid.
                                               // _kdTree.getData() will raise an exception in that case
        for (int j = 0; j < _nFunctions; j++) {
            fn[i][j] = _function[indices[i]][j];
        }
        if (_useMaxMin == 1) {
            for (int j = 0; j < _dimensions; j++) {
                pts[i][j] *= (_max[j] - _min[j]);
                pts[i][j] += _min[j];
            }
        }
    }
}

getDim()

template<typename T >

int lsst::afw::math::GaussianProcess< T >::getDim

(

)

const

return the dimensionality of data points stored in the GaussianProcess

Definition at line 924 of file GaussianProcess.cc.

                                     {
    return _dimensions;
}

getNPoints()

template<typename T >

int lsst::afw::math::GaussianProcess< T >::getNPoints

(

)

const

return the number of data points stored in the GaussianProcess

Definition at line 919 of file GaussianProcess.cc.

                                         {
    return _npts;
}

getTimes()

template<typename T >

GaussianProcessTimer & lsst::afw::math::GaussianProcess< T >::getTimes

(

)

const

Give the user acces to _timer, an object keeping track of the time spent on various processes within interpolate.

This will return a GaussianProcessTimer object. The user can, for example, see how much time has been spent on Eigen's linear algebra package (see the comments on the GaussianProcessTimer class) using code like

gg=GaussianProcess(....)

ticktock=gg.getTimes()

ticktock.display()

Definition at line 1981 of file GaussianProcess.cc.

                                                         {
    return _timer;
}

interpolate() [1/2]

template<typename T >

void lsst::afw::math::GaussianProcess< T >::interpolate	(	ndarray::Array< T, 1, 1 >	mu,
		ndarray::Array< T, 1, 1 >	variance,
		ndarray::Array< T, 1, 1 > const &	vin,
		int	numberOfNeighbors ) const

This is the version of GaussianProcess::interpolate for a vector of functions.

Parameters

[out]	mu	will store the vector of interpolated function values
[out]	variance	will store the vector of interpolated variances on mu
[in]	vin	the point at which you wish to interpolate the functions
[in]	numberOfNeighbors	is the number of nearest neighbor points to use in the interpolation

Note: Because the variance currently only depends on the covariance function and the covariance function currently does not include any terms relating different elements of mu to each other, all of the elements of variance will be identical

Definition at line 1131 of file GaussianProcess.cc.

                                                                                                    {
    if (numberOfNeighbors <= 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for zero or negative number of neighbors\n");
    }
 
    if (numberOfNeighbors > _kdTree.getNPoints()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for more neighbors than you have data points\n");
    }
 
    if (vin.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You are interpolating at a point with different dimensionality than you data\n");
    }
 
    if (mu.getNumElements() != static_cast<ndarray::Size>(_nFunctions) ||
        variance.getNumElements() != static_cast<ndarray::Size>(_nFunctions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your mu and/or var arrays are improperly sized for the number of functions "
                          "you are interpolating\n");
    }
 
    int i, j, ii;
    T fbar;
 
    ndarray::Array<T, 1, 1> covarianceTestPoint;
    ndarray::Array<int, 1, 1> neighbors;
    ndarray::Array<double, 1, 1> neighborDistances, vv;
 
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> covariance, bb, xx;
    Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
 
    _timer.start();
 
    bb.resize(numberOfNeighbors, 1);
    xx.resize(numberOfNeighbors, 1);
    covariance.resize(numberOfNeighbors, numberOfNeighbors);
    covarianceTestPoint = allocate(ndarray::makeVector(numberOfNeighbors));
    neighbors = allocate(ndarray::makeVector(numberOfNeighbors));
    neighborDistances = allocate(ndarray::makeVector(numberOfNeighbors));
 
    vv = allocate(ndarray::makeVector(_dimensions));
 
    if (_useMaxMin == 1) {
        // if you constructed this Gaussian process with minimum and maximum
        // values for the dimensions of your parameter space,
        // the point you are interpolating must be scaled to match the data so
        // that the selected nearest neighbors are appropriate
 
        for (i = 0; i < _dimensions; i++) vv[i] = (vin[i] - _min[i]) / (_max[i] - _min[i]);
    } else {
        vv = vin;
    }
 
    _kdTree.findNeighbors(neighbors, neighborDistances, vv, numberOfNeighbors);
 
    _timer.addToSearch();
 
    for (i = 0; i < numberOfNeighbors; i++) {
        covarianceTestPoint[i] = (*_covariogram)(vv, _kdTree.getData(neighbors[i]));
 
        covariance(i, i) =
                (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[i])) + _lambda;
 
        for (j = i + 1; j < numberOfNeighbors; j++) {
            covariance(i, j) = (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[j]));
            covariance(j, i) = covariance(i, j);
        }
    }
 
    _timer.addToIteration();
 
    // use Eigen's ldlt solver in place of matrix inversion (for speed purposes)
    ldlt.compute(covariance);
 
    for (ii = 0; ii < _nFunctions; ii++) {
        fbar = 0.0;
        for (i = 0; i < numberOfNeighbors; i++) fbar += _function[neighbors[i]][ii];
        fbar = fbar / double(numberOfNeighbors);
 
        for (i = 0; i < numberOfNeighbors; i++) bb(i, 0) = _function[neighbors[i]][ii] - fbar;
        xx = ldlt.solve(bb);
 
        mu[ii] = fbar;
 
        for (i = 0; i < numberOfNeighbors; i++) {
            mu[ii] += covarianceTestPoint[i] * xx(i, 0);
        }
 
    }  // ii = 0 through _nFunctions
 
    _timer.addToEigen();
 
    variance[0] = (*_covariogram)(vv, vv) + _lambda;
 
    for (i = 0; i < numberOfNeighbors; i++) bb(i) = covarianceTestPoint[i];
 
    xx = ldlt.solve(bb);
 
    for (i = 0; i < numberOfNeighbors; i++) {
        variance[0] -= covarianceTestPoint[i] * xx(i, 0);
    }
    variance[0] = variance[0] * _krigingParameter;
 
    for (i = 1; i < _nFunctions; i++) variance[i] = variance[0];
 
    _timer.addToVariance();
    _timer.addToTotal(1);
}

interpolate() [2/2]

template<typename T >

T lsst::afw::math::GaussianProcess< T >::interpolate	(	ndarray::Array< T, 1, 1 >	variance,
		ndarray::Array< T, 1, 1 > const &	vin,
		int	numberOfNeighbors ) const

Interpolate the function value at one point using a specified number of nearest neighbors.

Parameters

[out]	variance	a one-dimensional ndarray. The value of the variance predicted by the Gaussian process will be stored in the zeroth element
[in]	vin	a one-dimensional ndarray representing the point at which you want to interpolate the function
[in]	numberOfNeighbors	the number of nearest neighbors to be used in the interpolation

the interpolated value of the function will be returned at the end of this method

Note: if you used a normalized parameter space, you should not normalize vin before inputting. The code will remember that you want a normalized parameter space, and will apply the normalization when you call interpolate

Definition at line 1009 of file GaussianProcess.cc.

                                                               {
    if (_nFunctions > 1) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You need to call the version of GaussianProcess.interpolate() "
                          "that accepts mu and variance arrays (which it populates with results). "
                          "You are interpolating more than one function.");
    }
 
    if (numberOfNeighbors <= 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for zero or negative number of neighbors\n");
    }
 
    if (numberOfNeighbors > _kdTree.getNPoints()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for more neighbors than you have data points\n");
    }
 
    if (variance.getNumElements() != static_cast<ndarray::Size>(_nFunctions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your variance array is the incorrect size for the number "
                          "of functions you are trying to interpolate\n");
    }
 
    if (vin.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You are interpolating at a point with different dimensionality than you data\n");
    }
 
    int i, j;
    T fbar, mu;
 
    ndarray::Array<T, 1, 1> covarianceTestPoint;
    ndarray::Array<int, 1, 1> neighbors;
    ndarray::Array<double, 1, 1> neighborDistances, vv;
 
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> covariance, bb, xx;
    Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
 
    _timer.start();
 
    bb.resize(numberOfNeighbors, 1);
    xx.resize(numberOfNeighbors, 1);
 
    covariance.resize(numberOfNeighbors, numberOfNeighbors);
 
    covarianceTestPoint = allocate(ndarray::makeVector(numberOfNeighbors));
 
    neighbors = allocate(ndarray::makeVector(numberOfNeighbors));
 
    neighborDistances = allocate(ndarray::makeVector(numberOfNeighbors));
 
    vv = allocate(ndarray::makeVector(_dimensions));
    if (_useMaxMin == 1) {
        // if you constructed this Gaussian process with minimum and maximum
        // values for the dimensions of your parameter space,
        // the point you are interpolating must be scaled to match the data so
        // that the selected nearest neighbors are appropriate
 
        for (i = 0; i < _dimensions; i++) vv[i] = (vin[i] - _min[i]) / (_max[i] - _min[i]);
    } else {
        vv = vin;
    }
 
    _kdTree.findNeighbors(neighbors, neighborDistances, vv, numberOfNeighbors);
 
    _timer.addToSearch();
 
    fbar = 0.0;
    for (i = 0; i < numberOfNeighbors; i++) fbar += _function[neighbors[i]][0];
    fbar = fbar / double(numberOfNeighbors);
 
    for (i = 0; i < numberOfNeighbors; i++) {
        covarianceTestPoint[i] = (*_covariogram)(vv, _kdTree.getData(neighbors[i]));
 
        covariance(i, i) =
                (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[i])) + _lambda;
 
        for (j = i + 1; j < numberOfNeighbors; j++) {
            covariance(i, j) = (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[j]));
            covariance(j, i) = covariance(i, j);
        }
    }
 
    _timer.addToIteration();
 
    // use Eigen's ldlt solver in place of matrix inversion (for speed purposes)
    ldlt.compute(covariance);
 
    for (i = 0; i < numberOfNeighbors; i++) bb(i, 0) = _function[neighbors[i]][0] - fbar;
 
    xx = ldlt.solve(bb);
    _timer.addToEigen();
 
    mu = fbar;
 
    for (i = 0; i < numberOfNeighbors; i++) {
        mu += covarianceTestPoint[i] * xx(i, 0);
    }
 
    _timer.addToIteration();
 
    variance(0) = (*_covariogram)(vv, vv) + _lambda;
 
    for (i = 0; i < numberOfNeighbors; i++) bb(i) = covarianceTestPoint[i];
 
    xx = ldlt.solve(bb);
 
    for (i = 0; i < numberOfNeighbors; i++) {
        variance(0) -= covarianceTestPoint[i] * xx(i, 0);
    }
 
    variance(0) = variance(0) * _krigingParameter;
 
    _timer.addToVariance();
    _timer.addToTotal(1);
 
    return mu;
}

operator=() [1/2]

template<typename T >

GaussianProcess & lsst::afw::math::GaussianProcess< T >::operator= ( const GaussianProcess< T > & )

delete

operator=() [2/2]

template<typename T >

GaussianProcess & lsst::afw::math::GaussianProcess< T >::operator= ( GaussianProcess< T > && )

delete

removePoint()

template<typename T >

void lsst::afw::math::GaussianProcess< T >::removePoint

(

int

dex

)

This will remove a point from the data set.

Parameters

[in]

dex

the index of the point you want to remove from your data set

Exceptions

pex::exceptions::RuntimeError

if the tree does not end up properly constructed (the exception is actually thrown by KdTree<T>::removePoint() )

Note: excessive use of addPoint and removePoint can result in an unbalanced KdTree, which will slow down nearest neighbor searches

Definition at line 1952 of file GaussianProcess.cc.

                                            {
    int i, j;
 
    _kdTree.removePoint(dex);
 
    for (i = dex; i < _npts; i++) {
        for (j = 0; j < _nFunctions; j++) {
            _function[i][j] = _function[i + 1][j];
        }
    }
    _npts = _kdTree.getNPoints();
}

selfInterpolate() [1/2]

template<typename T >

void lsst::afw::math::GaussianProcess< T >::selfInterpolate	(	ndarray::Array< T, 1, 1 >	mu,
		ndarray::Array< T, 1, 1 >	variance,
		int	dex,
		int	numberOfNeighbors ) const

The version of selfInterpolate called for a vector of functions.

Parameters

[out]	mu	this is where the interpolated function values will be stored
[out]	variance	the variance on mu will be stored here
[in]	dex	the index of the point you wish to interpolate
[in]	numberOfNeighbors	the number of nearest neighbors to use in the interpolation

Exceptions

pex::exceptions::RuntimeError

if the nearest neighbor search does not find the data point itself as the nearest neighbor

Definition at line 1367 of file GaussianProcess.cc.

                                                                               {
    if (mu.getNumElements() != static_cast<ndarray::Size>(_nFunctions) ||
        variance.getNumElements() != static_cast<ndarray::Size>(_nFunctions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your mu and/or var arrays are improperly sized for the number of functions "
                          "you are interpolating\n");
    }
 
    if (numberOfNeighbors <= 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for zero or negative number of neighbors\n");
    }
 
    if (numberOfNeighbors + 1 > _kdTree.getNPoints()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for more neighbors than you have data points\n");
    }
 
    if (dex < 0 || dex >= _kdTree.getNPoints()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked to self interpolate on a point that does not exist\n");
    }
 
    int i, j, ii;
    T fbar;
 
    ndarray::Array<T, 1, 1> covarianceTestPoint;
    ndarray::Array<int, 1, 1> selfNeighbors;
    ndarray::Array<double, 1, 1> selfDistances;
    ndarray::Array<int, 1, 1> neighbors;
    ndarray::Array<double, 1, 1> neighborDistances;
 
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> covariance, bb, xx;
    Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
 
    _timer.start();
 
    bb.resize(numberOfNeighbors, 1);
    xx.resize(numberOfNeighbors, 1);
    covariance.resize(numberOfNeighbors, numberOfNeighbors);
    covarianceTestPoint = allocate(ndarray::makeVector(numberOfNeighbors));
    neighbors = allocate(ndarray::makeVector(numberOfNeighbors));
    neighborDistances = allocate(ndarray::makeVector(numberOfNeighbors));
 
    selfNeighbors = allocate(ndarray::makeVector(numberOfNeighbors + 1));
    selfDistances = allocate(ndarray::makeVector(numberOfNeighbors + 1));
 
    // we don't use _useMaxMin because the data has already been normalized
 
    _kdTree.findNeighbors(selfNeighbors, selfDistances, _kdTree.getData(dex), numberOfNeighbors + 1);
 
    _timer.addToSearch();
 
    if (selfNeighbors[0] != dex) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Nearest neighbor search in selfInterpolate did not find self\n");
    }
 
    // SelfNeighbors[0] will be the point itself (it is its own nearest neighbor)
    // We discard that for the interpolation calculation
    //
    // If you do not wish to do this, simply call the usual ::interpolate() method instead of
    //::selfInterpolate()
    for (i = 0; i < numberOfNeighbors; i++) {
        neighbors[i] = selfNeighbors[i + 1];
        neighborDistances[i] = selfDistances[i + 1];
    }
 
    for (i = 0; i < numberOfNeighbors; i++) {
        covarianceTestPoint[i] = (*_covariogram)(_kdTree.getData(dex), _kdTree.getData(neighbors[i]));
 
        covariance(i, i) =
                (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[i])) + _lambda;
 
        for (j = i + 1; j < numberOfNeighbors; j++) {
            covariance(i, j) = (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[j]));
            covariance(j, i) = covariance(i, j);
        }
    }
    _timer.addToIteration();
 
    // use Eigen's ldlt solver in place of matrix inversion (for speed purposes)
    ldlt.compute(covariance);
 
    for (ii = 0; ii < _nFunctions; ii++) {
        fbar = 0.0;
        for (i = 0; i < numberOfNeighbors; i++) fbar += _function[neighbors[i]][ii];
        fbar = fbar / double(numberOfNeighbors);
 
        for (i = 0; i < numberOfNeighbors; i++) bb(i, 0) = _function[neighbors[i]][ii] - fbar;
        xx = ldlt.solve(bb);
 
        mu[ii] = fbar;
 
        for (i = 0; i < numberOfNeighbors; i++) {
            mu[ii] += covarianceTestPoint[i] * xx(i, 0);
        }
    }  // ii = 0 through _nFunctions
 
    _timer.addToEigen();
 
    variance[0] = (*_covariogram)(_kdTree.getData(dex), _kdTree.getData(dex)) + _lambda;
 
    for (i = 0; i < numberOfNeighbors; i++) bb(i) = covarianceTestPoint[i];
 
    xx = ldlt.solve(bb);
 
    for (i = 0; i < numberOfNeighbors; i++) {
        variance[0] -= covarianceTestPoint[i] * xx(i, 0);
    }
 
    variance[0] = variance[0] * _krigingParameter;
 
    for (i = 1; i < _nFunctions; i++) variance[i] = variance[0];
 
    _timer.addToVariance();
    _timer.addToTotal(1);
}

selfInterpolate() [2/2]

template<typename T >

T lsst::afw::math::GaussianProcess< T >::selfInterpolate	(	ndarray::Array< T, 1, 1 >	variance,
		int	dex,
		int	numberOfNeighbors ) const

This method will interpolate the function on a data point for purposes of optimizing hyper parameters.

Parameters

[out]	variance	a one-dimensional ndarray. The value of the variance predicted by the Gaussian process will be stored in the zeroth element
[in]	dex	the index of the point you wish to self interpolate
[in]	numberOfNeighbors	the number of nearest neighbors to be used in the interpolation

Exceptions

pex::exceptions::RuntimeError

if the nearest neighbor search does not find the data point itself as the nearest neighbor

The interpolated value of the function will be returned at the end of this method

This method ignores the point on which you are interpolating when requesting nearest neighbors

Definition at line 1244 of file GaussianProcess.cc.

                                                                   {
    if (_nFunctions > 1) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You need to call the version of GaussianProcess.selfInterpolate() "
                          "that accepts mu and variance arrays (which it populates with results). "
                          "You are interpolating more than one function.");
    }
 
    if (variance.getNumElements() != static_cast<ndarray::Size>(_nFunctions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Your variance array is the incorrect size for the number "
                          "of functions you are trying to interpolate\n");
    }
 
    if (numberOfNeighbors <= 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for zero or negative number of neighbors\n");
    }
 
    if (numberOfNeighbors + 1 > _kdTree.getNPoints()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for more neighbors than you have data points\n");
    }
 
    if (dex < 0 || dex >= _kdTree.getNPoints()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked to self interpolate on a point that does not exist\n");
    }
 
    int i, j;
    T fbar, mu;
 
    ndarray::Array<T, 1, 1> covarianceTestPoint;
    ndarray::Array<int, 1, 1> selfNeighbors;
    ndarray::Array<double, 1, 1> selfDistances;
    ndarray::Array<int, 1, 1> neighbors;
    ndarray::Array<double, 1, 1> neighborDistances;
 
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> covariance, bb, xx;
    Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
 
    _timer.start();
 
    bb.resize(numberOfNeighbors, 1);
    xx.resize(numberOfNeighbors, 1);
    covariance.resize(numberOfNeighbors, numberOfNeighbors);
    covarianceTestPoint = allocate(ndarray::makeVector(numberOfNeighbors));
    neighbors = allocate(ndarray::makeVector(numberOfNeighbors));
    neighborDistances = allocate(ndarray::makeVector(numberOfNeighbors));
 
    selfNeighbors = allocate(ndarray::makeVector(numberOfNeighbors + 1));
    selfDistances = allocate(ndarray::makeVector(numberOfNeighbors + 1));
 
    // we don't use _useMaxMin because the data has already been normalized
 
    _kdTree.findNeighbors(selfNeighbors, selfDistances, _kdTree.getData(dex), numberOfNeighbors + 1);
 
    _timer.addToSearch();
 
    if (selfNeighbors[0] != dex) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Nearest neighbor search in selfInterpolate did not find self\n");
    }
 
    // SelfNeighbors[0] will be the point itself (it is its own nearest neighbor)
    // We discard that for the interpolation calculation
    //
    // If you do not wish to do this, simply call the usual ::interpolate() method instead of
    //::selfInterpolate()
    for (i = 0; i < numberOfNeighbors; i++) {
        neighbors[i] = selfNeighbors[i + 1];
        neighborDistances[i] = selfDistances[i + 1];
    }
 
    fbar = 0.0;
    for (i = 0; i < numberOfNeighbors; i++) fbar += _function[neighbors[i]][0];
    fbar = fbar / double(numberOfNeighbors);
 
    for (i = 0; i < numberOfNeighbors; i++) {
        covarianceTestPoint[i] = (*_covariogram)(_kdTree.getData(dex), _kdTree.getData(neighbors[i]));
 
        covariance(i, i) =
                (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[i])) + _lambda;
 
        for (j = i + 1; j < numberOfNeighbors; j++) {
            covariance(i, j) = (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[j]));
            covariance(j, i) = covariance(i, j);
        }
    }
    _timer.addToIteration();
 
    // use Eigen's ldlt solver in place of matrix inversion (for speed purposes)
    ldlt.compute(covariance);
 
    for (i = 0; i < numberOfNeighbors; i++) bb(i, 0) = _function[neighbors[i]][0] - fbar;
    xx = ldlt.solve(bb);
    _timer.addToEigen();
 
    mu = fbar;
 
    for (i = 0; i < numberOfNeighbors; i++) {
        mu += covarianceTestPoint[i] * xx(i, 0);
    }
 
    variance(0) = (*_covariogram)(_kdTree.getData(dex), _kdTree.getData(dex)) + _lambda;
 
    for (i = 0; i < numberOfNeighbors; i++) bb(i) = covarianceTestPoint[i];
 
    xx = ldlt.solve(bb);
 
    for (i = 0; i < numberOfNeighbors; i++) {
        variance(0) -= covarianceTestPoint[i] * xx(i, 0);
    }
 
    variance(0) = variance(0) * _krigingParameter;
    _timer.addToVariance();
    _timer.addToTotal(1);
 
    return mu;
}

setCovariogram()

template<typename T >

void lsst::afw::math::GaussianProcess< T >::setCovariogram

(

std::shared_ptr< Covariogram< T > > const &

covar

)

Assign a different covariogram to this GaussianProcess.

Parameters

[in]

covar

the Covariogram object that you wish to assign

Definition at line 1971 of file GaussianProcess.cc.

                                                                                   {
    _covariogram = covar;
}

setKrigingParameter()

template<typename T >

void lsst::afw::math::GaussianProcess< T >::setKrigingParameter

(

T

kk

)

Assign a value to the Kriging paramter.

Parameters

[in]

kk

the value assigned to the Kriging parameters

Definition at line 1966 of file GaussianProcess.cc.

                                                 {
    _krigingParameter = kk;
}

setLambda()

template<typename T >

void lsst::afw::math::GaussianProcess< T >::setLambda

(

T

lambda

)

set the value of the hyperparameter _lambda

Parameters

[in]

lambda

the value you want assigned to _lambda

_lambda is a parameter meant to represent the characteristic variance of the function you are interpolating. Currently, it is a scalar such that all data points must have the same characteristic variance. Future iterations of the code may want to promote _lambda to an array so that different data points can have different variances.

Definition at line 1976 of file GaussianProcess.cc.

                                           {
    _lambda = lambda;
}

---

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

• /j/snowflake/release/lsstsw/stack/lsst-scipipe-9.0.0/Linux64/afw/g9726552aa6+10f999ec6a/include/lsst/afw/math/GaussianProcess.h
• /j/snowflake/release/lsstsw/stack/lsst-scipipe-9.0.0/Linux64/afw/g9726552aa6+10f999ec6a/src/math/GaussianProcess.cc