LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Public Member Functions | Private Attributes | List of all members
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>

Inheritance diagram for lsst.afw.math::GaussianProcess< T >:

Public Member Functions

 GaussianProcess (ndarray::Array< T, 2, 2 > const &dataIn, ndarray::Array< T, 1, 1 > const &ff, boost::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. More...
 
 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, boost::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. More...
 
 GaussianProcess (ndarray::Array< T, 2, 2 > const &dataIn, ndarray::Array< T, 2, 2 > const &ff, boost::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 More...
 
 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, boost::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 More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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) More...
 
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. Do not return variances for the interpolation. More...
 
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. More...
 
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. More...
 
void addPoint (ndarray::Array< T, 1, 1 > const &vin, T f)
 Add a point to the pool of data used by GaussianProcess for interpolation. More...
 
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. More...
 
void removePoint (int dex)
 This will remove a point from the data set. More...
 
void setKrigingParameter (T kk)
 Assign a value to the Kriging paramter. More...
 
void setCovariogram (boost::shared_ptr< Covariogram< T > > const &covar)
 Assign a different covariogram to this GaussianProcess. More...
 
void setLambda (T lambda)
 set the value of the hyperparameter _lambda More...
 
GaussianProcessTimergetTimes () const
 Give the user acces to _timer, an object keeping track of the time spent on various processes within interpolate. More...
 

Private Attributes

int _pts
 
int _useMaxMin
 
int _dimensions
 
int _room
 
int _roomStep
 
int _nFunctions
 
_krigingParameter
 
_lambda
 
ndarray::Array< T, 1, 1 > _max
 
ndarray::Array< T, 1, 1 > _min
 
ndarray::Array< T, 2, 2 > _function
 
KdTree< T > _kdTree
 
boost::shared_ptr< Covariogram
< T > > 
_covariogram
 
GaussianProcessTimer _timer
 

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 499 of file GaussianProcess.h.

Constructor & Destructor Documentation

template<typename T >
lsst.afw.math::GaussianProcess< T >::GaussianProcess ( ndarray::Array< T, 2, 2 > const &  dataIn,
ndarray::Array< T, 1, 1 > const &  ff,
boost::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]dataInan ndarray containing the data points; the ith row of datain is the ith data point
[in]ffa one-dimensional ndarray containing the values of the scalar function associated with each data point. This is the function you are interpolating
[in]covarInis the input covariogram

Definition at line 745 of file GaussianProcess.cc.

749 {
750  int i;
751 
752  _covariogram = covarIn;
753 
754  _pts = dataIn.template getSize < 0 > ();
755  _dimensions = dataIn.template getSize < 1 > ();
756 
757  _room = _pts;
758  _roomStep = 5000;
759 
760  _nFunctions = 1;
762 
763  for(i = 0; i < _pts; i++ ) _function[i][0] = ff[i];
764 
765  _krigingParameter = T(1.0);
766  _lambda = T(1.0e-5);
767 
768  _useMaxMin = 0;
769 
770 
771  _kdTree.Initialize(dataIn);
772 
773  _pts = _kdTree.getPoints();
774 
775 }
boost::shared_ptr< Covariogram< T > > _covariogram
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
ndarray::Array< T, 2, 2 > _function
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,
boost::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]dataInan ndarray containing the data points; the ith row of datain is the ith data point
[in]mna one-dimensional ndarray containing the minimum values of each dimension (for normalizing the positions of data points)
[in]mxa one-dimensional ndarray containing the maximum values of each dimension (for normalizing the positions of data points)
[in]ffa one-dimensional ndarray containing the values of the scalar function associated with each data point. This is the function you are interpolating
[in]covarInis the input covariogram

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

Definition at line 778 of file GaussianProcess.cc.

784 {
785 
786  int i,j;
787  ndarray::Array < T,2,2 > normalizedData;
788 
789  _covariogram = covarIn;
790 
791  _pts = dataIn.template getSize < 0 > ();
792  _dimensions = dataIn.template getSize < 1 > ();
793  _room = _pts;
794  _roomStep = 5000;
795 
796  _krigingParameter = T(1.0);
797 
798  _lambda = T(1.0e-5);
799  _krigingParameter = T(1.0);
800 
801  _max = allocate(ndarray::makeVector(_dimensions));
802  _min = allocate(ndarray::makeVector(_dimensions));
803  _max.deep() = mx;
804  _min.deep() = mn;
805  _useMaxMin = 1;
806  normalizedData = allocate(ndarray::makeVector(_pts, _dimensions));
807  for(i = 0; i < _pts; i++ ) {
808  for(j = 0; j < _dimensions; j++ ) {
809  normalizedData[i][j] = (dataIn[i][j] - _min[j])/(_max[j] - _min[j]);
810  //note the normalization by _max - _min in each dimension
811  }
812  }
813 
814  _kdTree.Initialize(normalizedData);
815 
816  _pts = _kdTree.getPoints();
817  _nFunctions = 1;
819  for(i = 0; i < _pts; i++ )_function[i][0] = ff[i];
820 }
ndarray::Array< T, 1, 1 > _max
boost::shared_ptr< Covariogram< T > > _covariogram
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
ndarray::Array< T, 1, 1 > _min
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
Deep const deep() const
Return an ArrayRef view to this.
Definition: ArrayBase.h:178
ndarray::Array< T, 2, 2 > _function
template<typename T >
lsst.afw.math::GaussianProcess< T >::GaussianProcess ( ndarray::Array< T, 2, 2 > const &  dataIn,
ndarray::Array< T, 2, 2 > const &  ff,
boost::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]dataIncontains the data points, as in other constructors
[in]ffcontains 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]covarInis the input covariogram

Definition at line 823 of file GaussianProcess.cc.

826 {
827 
828  _covariogram = covarIn;
829 
830  _pts = dataIn.template getSize < 0 > ();
831  _dimensions = dataIn.template getSize < 1 > ();
832 
833  _room = _pts;
834  _roomStep = 5000;
835 
836  _nFunctions = ff.template getSize < 1 > ();
838  _function.deep() = ff;
839 
840  _krigingParameter = T(1.0);
841 
842  _lambda = T(1.0e-5);
843 
844  _useMaxMin = 0;
845 
846 
847  _kdTree.Initialize(dataIn);
848 
849  _pts = _kdTree.getPoints();
850 }
boost::shared_ptr< Covariogram< T > > _covariogram
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
ndarray::Array< T, 2, 2 > _function
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,
boost::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]dataIncontains the data points, as in other constructors
[in]mncontains the minimum allowed values of the parameters in parameter space
[in]mxcontains the maximum allowed values of the parameters in parameter space
[in]ffcontains 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]covarInis the input covariogram

Definition at line 853 of file GaussianProcess.cc.

859 {
860 
861  int i,j;
862  ndarray::Array < T,2,2 > normalizedData;
863 
864  _covariogram = covarIn;
865 
866  _pts = dataIn.template getSize < 0 > ();
867  _dimensions = dataIn.template getSize < 1 > ();
868 
869  _room = _pts;
870  _roomStep = 5000;
871 
872  _krigingParameter = T(1.0);
873 
874  _lambda = T(1.0e-5);
875  _krigingParameter = T(1.0);
876 
877  _max = allocate(ndarray::makeVector(_dimensions));
878  _min = allocate(ndarray::makeVector(_dimensions));
879  _max.deep() = mx;
880  _min.deep() = mn;
881  _useMaxMin = 1;
882  normalizedData = allocate(ndarray::makeVector(_pts,_dimensions));
883  for(i = 0; i < _pts; i++ ) {
884  for(j = 0; j < _dimensions; j++ ) {
885  normalizedData[i][j] = (dataIn[i][j] - _min[j])/(_max[j] - _min[j]);
886  //note the normalization by _max - _min in each dimension
887  }
888  }
889 
890  _kdTree.Initialize(normalizedData);
891  _pts = _kdTree.getPoints();
892  _nFunctions = ff.template getSize < 1 > ();
894  _function.deep() = ff;
895 }
ndarray::Array< T, 1, 1 > _max
boost::shared_ptr< Covariogram< T > > _covariogram
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
ndarray::Array< T, 1, 1 > _min
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
Deep const deep() const
Return an ArrayRef view to this.
Definition: ArrayBase.h:178
ndarray::Array< T, 2, 2 > _function

Member Function Documentation

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

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

Parameters
[in]vina one-dimensional ndarray storing the point in parameter space that you are adding
[in]fthe value of the function at that point
Exceptions
pex_exceptionsRuntimeError if you call this when you should have called the version taking a vector of functions (below)
pex_exceptionsRuntimeError 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 1697 of file GaussianProcess.cc.

1698 {
1699 
1700  int i,j;
1701 
1702  if(_nFunctions!= 1){
1703 
1704  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
1705  "You are calling the wrong addPoint; you need a vector of functions\n");
1706 
1707  }
1708 
1711 
1712  for(i = 0; i < _dimensions; i++ ){
1713  v[i] = vin[i];
1714  if(_useMaxMin == 1){
1715  v[i] = (v[i] - _min[i])/(_max[i] - _min[i]);
1716  }
1717  }
1718 
1719  if(_pts == _room){
1722  buff.deep() = _function;
1723 
1724  _room += _roomStep;
1726  for(i = 0; i < _pts; i++ ) {
1727 
1728  for(j = 0; j < _nFunctions; j++ ) {
1729  _function[i][j] = buff[i][j];
1730  }
1731 
1732  }
1733  }
1734 
1735  _function[_pts][0] = f;
1736 
1737  _kdTree.addPoint(v);
1738  _pts = _kdTree.getPoints();
1739 
1740 }
ndarray::Array< T, 1, 1 > _max
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
ndarray::Array< T, 1, 1 > _min
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
Deep const deep() const
Return an ArrayRef view to this.
Definition: ArrayBase.h:178
ndarray::Array< T, 2, 2 > _function
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_exceptionsRuntimeError 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 1743 of file GaussianProcess.cc.

1745 {
1746 
1747  int i,j;
1748 
1751 
1752  for(i = 0; i < _dimensions; i++ ) {
1753  v[i] = vin[i];
1754  if(_useMaxMin == 1) {
1755  v[i] = (v[i] - _min[i])/(_max[i] - _min[i]);
1756  }
1757  }
1758 
1759  if(_pts == _room) {
1762  buff.deep() = _function;
1763 
1764  _room += _roomStep;
1766  for(i = 0; i < _pts; i++ ) {
1767  for(j = 0; j < _nFunctions; j++ ) {
1768  _function[i][j] = buff[i][j];
1769  }
1770  }
1771 
1772  }
1773  for(i = 0; i < _nFunctions; i++ )_function[_pts][i] = f[i];
1774 
1775  _kdTree.addPoint(v);
1776  _pts = _kdTree.getPoints();
1777 
1778 
1779 }
ndarray::Array< T, 1, 1 > _max
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
ndarray::Array< T, 1, 1 > _min
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
Deep const deep() const
Return an ArrayRef view to this.
Definition: ArrayBase.h:178
ndarray::Array< T, 2, 2 > _function
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]mua 1-dimensional ndarray where the interpolated function values will be stored
[out]variancea 1-dimensional ndarray where the corresponding variances in the function value will be stored
[in]queriesa 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 _pts X _pts covariance matrix C and solve the problem Cx=b. Be wary of using it in the case where _pts 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 1365 of file GaussianProcess.cc.

1368 {
1369 
1370  int i,j,ii,nQueries;
1371 
1372  T fbar;
1373  Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > batchCovariance,batchbb,batchxx;
1374  Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > queryCovariance;
1375  Eigen::LDLT < Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > > ldlt;
1376 
1378 
1379  _timer.start();
1380 
1381  nQueries = queries.template getSize < 0 > ();
1382 
1383 
1385  batchbb.resize(_pts, 1);
1386  batchxx.resize(_pts, 1);
1387  batchCovariance.resize(_pts, _pts);
1388  queryCovariance.resize(_pts, 1);
1389 
1390 
1391  for(i = 0; i < _pts; i++ ) {
1392 
1393  batchCovariance(i, i) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(i)) + _lambda;
1394  for(j = i + 1; j < _pts; j++ ) {
1395  batchCovariance(i, j) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(j));
1396  batchCovariance(j, i) = batchCovariance(i, j);
1397  }
1398  }
1400 
1401  ldlt.compute(batchCovariance);
1402 
1403  fbar = 0.0;
1404  for(i = 0; i < _pts; i++ ) {
1405  fbar += _function[i][0];
1406  }
1407  fbar = fbar/T(_pts);
1408 
1409  for(i = 0; i < _pts; i++ ){
1410  batchbb(i, 0) = _function[i][0] - fbar;
1411  }
1412  batchxx = ldlt.solve(batchbb);
1413  _timer.addToEigen();
1414 
1415 
1416  for(ii = 0; ii < nQueries; ii++ ) {
1417  for(i = 0; i < _dimensions; i++ )v1[i] = queries[ii][i];
1418  if(_useMaxMin == 1) {
1419  for(i = 0; i < _dimensions; i++ )v1[i] = (v1[i] - _min[i])/(_max[i] - _min[i]);
1420  }
1421  mu(ii) = fbar;
1422  for(i = 0; i < _pts; i++ ){
1423  mu(ii) += batchxx(i)*(*_covariogram)(v1, _kdTree.getData(i));
1424  }
1425  }
1427 
1428  for(ii = 0; ii < nQueries; ii++ ) {
1429  for(i = 0; i < _dimensions; i++ )v1[i] = queries[ii][i];
1430  if(_useMaxMin == 1) {
1431  for(i = 0; i < _dimensions; i++ )v1[i] = (v1[i] - _min[i])/(_max[i] - _min[i]);
1432  }
1433 
1434  for(i = 0; i < _pts; i++ ) {
1435  batchbb(i, 0) = (*_covariogram)(v1, _kdTree.getData(i));
1436  queryCovariance(i, 0) = batchbb(i, 0);
1437  }
1438  batchxx = ldlt.solve(batchbb);
1439 
1440  variance[ii] = (*_covariogram)(v1, v1) + _lambda;
1441 
1442  for(i = 0; i < _pts; i++ ){
1443  variance[ii] -= queryCovariance(i, 0)*batchxx(i);
1444  }
1445 
1446  variance[ii] = variance[ii]*_krigingParameter;
1447 
1448  }
1449 
1451  _timer.addToTotal(nQueries);
1452 }
ndarray::Array< T, 1, 1 > _max
void addToEigen()
Adds time to _eigenTime.
void addToIteration()
Adds time to _iterationTime.
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
ndarray::Array< T, 1, 1 > _min
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
ndarray::Array< T, 2, 2 > _function
void start()
Starts the timer for an individual call to an interpolation routine.
void addToVariance()
Adds time to _varianceTime.
void addToTotal(int i)
Adds time to _totalTime and adds counts to _interpolationCount.
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]mua 1-dimensional ndarray where the interpolated function values will be stored
[in]queriesa 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 _pts X _pts covariance matrix C and solve the problem Cx=b. Be wary of using it in the case where _pts 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 1552 of file GaussianProcess.cc.

1554 {
1555 
1556  int i,j,ii,nQueries;
1557 
1558  T fbar;
1559  Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > batchCovariance,batchbb,batchxx;
1560  Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > queryCovariance;
1561  Eigen::LDLT < Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > > ldlt;
1562 
1564 
1565  _timer.start();
1566 
1567  nQueries = queries.template getSize < 0 > ();
1568 
1570 
1571  batchbb.resize(_pts, 1);
1572  batchxx.resize(_pts, 1);
1573  batchCovariance.resize(_pts, _pts);
1574  queryCovariance.resize(_pts, 1);
1575 
1576 
1577  for(i = 0; i < _pts; i++ ) {
1578  batchCovariance(i, i) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(i)) + _lambda;
1579  for(j = i + 1; j < _pts; j++ ) {
1580  batchCovariance(i, j) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(j));
1581  batchCovariance(j, i) = batchCovariance(i, j);
1582  }
1583  }
1585 
1586  ldlt.compute(batchCovariance);
1587 
1588  fbar = 0.0;
1589  for(i = 0; i < _pts; i++ ) {
1590  fbar += _function[i][0];
1591  }
1592  fbar = fbar/T(_pts);
1593 
1594  for(i = 0; i < _pts; i++ ) {
1595  batchbb(i, 0) = _function[i][0] - fbar;
1596  }
1597  batchxx = ldlt.solve(batchbb);
1598  _timer.addToEigen();
1599 
1600  for(ii = 0; ii < nQueries; ii++ ) {
1601 
1602  for(i = 0; i < _dimensions; i++ ) v1[i] = queries[ii][i];
1603  if(_useMaxMin == 1) {
1604  for(i = 0; i < _dimensions; i++ )v1[i] = (v1[i] - _min[i])/(_max[i] - _min[i]);
1605  }
1606 
1607  mu(ii) = fbar;
1608  for(i = 0; i < _pts; i++ ) {
1609  mu(ii) += batchxx(i)*(*_covariogram)(v1, _kdTree.getData(i));
1610  }
1611  }
1613  _timer.addToTotal(nQueries);
1614 
1615 }
ndarray::Array< T, 1, 1 > _max
void addToEigen()
Adds time to _eigenTime.
void addToIteration()
Adds time to _iterationTime.
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
ndarray::Array< T, 1, 1 > _min
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
ndarray::Array< T, 2, 2 > _function
void start()
Starts the timer for an individual call to an interpolation routine.
void addToTotal(int i)
Adds time to _totalTime and adds counts to _interpolationCount.
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 1455 of file GaussianProcess.cc.

1458 {
1459 
1460  int i,j,ii,nQueries,ifn;
1461 
1462  T fbar;
1463  Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > batchCovariance,batchbb,batchxx;
1464  Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > queryCovariance;
1465  Eigen::LDLT < Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > > ldlt;
1466 
1468 
1469  _timer.start();
1470 
1471  nQueries = queries.template getSize < 0 > ();
1472 
1473 
1474 
1476  batchbb.resize(_pts, 1);
1477  batchxx.resize(_pts, 1);
1478  batchCovariance.resize(_pts, _pts);
1479  queryCovariance.resize(_pts, 1);
1480 
1481  for(i = 0; i < _pts; i++ ){
1482  batchCovariance(i, i) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(i)) + _lambda;
1483  for(j = i + 1; j < _pts; j++ ) {
1484  batchCovariance(i, j) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(j));
1485  batchCovariance(j, i) = batchCovariance(i, j);
1486  }
1487  }
1488 
1490 
1491  ldlt.compute(batchCovariance);
1492 
1493  _timer.addToEigen();
1494 
1495  for(ifn = 0; ifn < _nFunctions; ifn++ ) {
1496 
1497  fbar = 0.0;
1498  for(i = 0; i < _pts; i++ ){
1499  fbar += _function[i][ifn];
1500  }
1501  fbar = fbar/T(_pts);
1503 
1504  for(i = 0; i < _pts; i++ ){
1505  batchbb(i,0) = _function[i][ifn] - fbar;
1506  }
1507  batchxx = ldlt.solve(batchbb);
1508  _timer.addToEigen();
1509 
1510 
1511  for(ii = 0; ii < nQueries; ii++ ){
1512  for(i = 0; i < _dimensions; i++ ) v1[i] = queries[ii][i];
1513  if(_useMaxMin == 1) {
1514  for(i = 0; i < _dimensions; i++ ) v1[i] = (v1[i] - _min[i])/(_max[i] - _min[i]);
1515  }
1516  mu[ii][ifn] = fbar;
1517  for(i = 0; i < _pts; i++ ){
1518  mu[ii][ifn] += batchxx(i)*(*_covariogram)(v1, _kdTree.getData(i));
1519  }
1520  }
1521 
1522  }//ifn = 0 to _nFunctions
1523 
1525  for(ii = 0; ii < nQueries; ii++ ){
1526  for(i = 0; i < _dimensions; i++ ) v1[i] = queries[ii][i];
1527  if(_useMaxMin == 1){
1528  for(i = 0; i < _dimensions; i++ ) v1[i] = (v1[i] - _min[i])/(_max[i] - _min[i]);
1529  }
1530 
1531  for(i = 0;i < _pts;i++ ) {
1532  batchbb(i,0) = (*_covariogram)(v1,_kdTree.getData(i));
1533  queryCovariance(i,0) = batchbb(i,0);
1534  }
1535  batchxx = ldlt.solve(batchbb);
1536 
1537  variance[ii][0] = (*_covariogram)(v1, v1) + _lambda;
1538 
1539  for(i = 0; i < _pts; i++ ) {
1540  variance[ii][0] -= queryCovariance(i, 0)*batchxx(i);
1541  }
1542 
1543  variance[ii][0] = variance[ii][0]*_krigingParameter;
1544  for(i = 1; i < _nFunctions; i++ ) variance[ii][i] = variance[ii][0];
1545 
1546  }
1548  _timer.addToTotal(nQueries);
1549 }
ndarray::Array< T, 1, 1 > _max
void addToEigen()
Adds time to _eigenTime.
void addToIteration()
Adds time to _iterationTime.
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
ndarray::Array< T, 1, 1 > _min
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
ndarray::Array< T, 2, 2 > _function
void start()
Starts the timer for an individual call to an interpolation routine.
void addToVariance()
Adds time to _varianceTime.
void addToTotal(int i)
Adds time to _totalTime and adds counts to _interpolationCount.
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 1618 of file GaussianProcess.cc.

1620 {
1621 
1622  int i,j,ii,nQueries,ifn;
1623 
1624 
1625  T fbar;
1626  Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > batchCovariance,batchbb,batchxx;
1627  Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > queryCovariance;
1628  Eigen::LDLT < Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > > ldlt;
1629 
1631 
1632  _timer.start();
1633 
1634  nQueries = queries.template getSize < 0 > ();
1635 
1636 
1638 
1639  batchbb.resize(_pts, 1);
1640  batchxx.resize(_pts, 1);
1641  batchCovariance.resize(_pts, _pts);
1642  queryCovariance.resize(_pts, 1);
1643 
1644 
1645  for(i = 0; i < _pts; i++ ){
1646  batchCovariance(i, i) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(i)) + _lambda;
1647  for(j = i + 1; j < _pts; j++ ){
1648  batchCovariance(i, j) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(j));
1649  batchCovariance(j, i) = batchCovariance(i, j);
1650  }
1651  }
1652 
1654 
1655  ldlt.compute(batchCovariance);
1656 
1657  _timer.addToEigen();
1658 
1659  for(ifn = 0; ifn < _nFunctions; ifn++ ){
1660 
1661  fbar = 0.0;
1662  for(i = 0; i < _pts; i++ ){
1663  fbar += _function[i][ifn];
1664  }
1665  fbar = fbar/T(_pts);
1666 
1668 
1669  for(i = 0; i < _pts; i++ ){
1670  batchbb(i, 0) = _function[i][ifn] - fbar;
1671  }
1672  batchxx = ldlt.solve(batchbb);
1673  _timer.addToEigen();
1674 
1675 
1676  for(ii = 0; ii < nQueries; ii++ ){
1677  for(i = 0; i < _dimensions; i++ )v1[i] = queries[ii][i];
1678  if(_useMaxMin == 1){
1679  for(i = 0;i < _dimensions;i++ )v1[i] = (v1[i] - _min[i])/(_max[i] - _min[i]);
1680  }
1681 
1682  mu[ii][ifn] = fbar;
1683  for(i = 0; i < _pts; i++ ){
1684  mu[ii][ifn] += batchxx(i)*(*_covariogram)(v1, _kdTree.getData(i));
1685  }
1686  }
1687 
1688 
1689  }//ifn = 0 through _nFunctions
1690 
1691  _timer.addToTotal(nQueries);
1692 
1693 }
ndarray::Array< T, 1, 1 > _max
void addToEigen()
Adds time to _eigenTime.
void addToIteration()
Adds time to _iterationTime.
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
ndarray::Array< T, 1, 1 > _min
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
ndarray::Array< T, 2, 2 > _function
void start()
Starts the timer for an individual call to an interpolation routine.
void addToTotal(int i)
Adds time to _totalTime and adds counts to _interpolationCount.
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 1815 of file GaussianProcess.cc.

1816 {
1817  return _timer;
1818 }
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]variancea one-dimensional ndarray. The value of the variance predicted by the Gaussian process will be stored in the zeroth element
[in]vina one-dimensional ndarray representing the point at which you want to interpolate the function
[in]numberOfNeighborsthe 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 899 of file GaussianProcess.cc.

902 {
903 
904  if(numberOfNeighbors <= 0){
905  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
906  "Asked for zero or negative number of neighbors\n");
907  }
908 
909  if(numberOfNeighbors > _kdTree.getPoints()){
910  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
911  "Asked for more neighbors than you have data points\n");
912  }
913 
914  int i,j;
915  T fbar,mu;
916 
917  ndarray::Array < T,1,1 > covarianceTestPoint;
918  ndarray::Array < int,1,1 > neighbors;
919  ndarray::Array < double,1,1 > neighborDistances,vv;
920 
921  Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > covariance,bb,xx;
922  Eigen::LDLT < Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > > ldlt;
923 
924 
925  _timer.start();
926 
927  bb.resize(numberOfNeighbors, 1);
928  xx.resize(numberOfNeighbors, 1);
929 
930  covariance.resize(numberOfNeighbors,numberOfNeighbors);
931 
932  covarianceTestPoint = allocate(ndarray::makeVector(numberOfNeighbors));
933 
934  neighbors = allocate(ndarray::makeVector(numberOfNeighbors));
935 
936  neighborDistances = allocate(ndarray::makeVector(numberOfNeighbors));
937 
939  if(_useMaxMin == 1){
940  //if you constructed this Gaussian process with minimum and maximum
941  //values for the dimensions of your parameter space,
942  //the point you are interpolating must be scaled to match the data so
943  //that the selected nearest neighbors are appropriate
944 
945  for(i = 0; i < _dimensions; i++ ) vv[i] = (vin[i] - _min[i])/(_max[i] - _min[i]);
946  }
947  else{
948  vv = vin;
949  }
950 
951 
952  _kdTree.findNeighbors(neighbors, neighborDistances, vv, numberOfNeighbors);
953 
955 
956  fbar = 0.0;
957  for(i = 0; i < numberOfNeighbors; i++ )fbar += _function[neighbors[i]][0];
958  fbar = fbar/double(numberOfNeighbors);
959 
960  for(i = 0; i < numberOfNeighbors; i++ ){
961  covarianceTestPoint[i] = (*_covariogram)(vv, _kdTree.getData(neighbors[i]));
962 
963  covariance(i,i) = (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[i]))
964  + _lambda;
965 
966  for(j = i + 1; j < numberOfNeighbors; j++ ){
967  covariance(i,j) = (*_covariogram)(_kdTree.getData(neighbors[i]),
968  _kdTree.getData(neighbors[j]));
969  covariance(j,i) = covariance(i, j);
970  }
971  }
972 
974 
975  //use Eigen's ldlt solver in place of matrix inversion (for speed purposes)
976  ldlt.compute(covariance);
977 
978  for(i = 0; i < numberOfNeighbors; i++) bb(i,0) = _function[neighbors[i]][0] - fbar;
979 
980  xx = ldlt.solve(bb);
981  _timer.addToEigen();
982 
983  mu = fbar;
984 
985  for(i = 0; i < numberOfNeighbors; i++ ){
986  mu += covarianceTestPoint[i]*xx(i, 0);
987  }
988 
990 
991  variance(0) = (*_covariogram)(vv, vv) + _lambda;
992 
993  for(i = 0; i < numberOfNeighbors; i++ ) bb(i) = covarianceTestPoint[i];
994 
995  xx = ldlt.solve(bb);
996 
997  for(i = 0; i < numberOfNeighbors; i++ ){
998  variance(0) -= covarianceTestPoint[i]*xx(i, 0);
999  }
1000 
1001  variance(0) = variance(0)*_krigingParameter;
1002 
1004  _timer.addToTotal(1);
1005 
1006  return mu;
1007 }
ndarray::Array< T, 1, 1 > _max
void addToEigen()
Adds time to _eigenTime.
void addToSearch()
Adds time to _searchTime.
void addToIteration()
Adds time to _iterationTime.
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
ndarray::Array< T, 1, 1 > _min
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
ndarray::Array< T, 2, 2 > _function
void start()
Starts the timer for an individual call to an interpolation routine.
void addToVariance()
Adds time to _varianceTime.
void addToTotal(int i)
Adds time to _totalTime and adds counts to _interpolationCount.
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]muwill store the vector of interpolated function values
[out]variancewill store the vector of interpolated variances on mu
[in]vinthe point at which you wish to interpolate the functions
[in]numberOfNeighborsis 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 1010 of file GaussianProcess.cc.

1014 {
1015 
1016 
1017  if(numberOfNeighbors <= 0){
1018  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
1019  "Asked for zero or negative number of neighbors\n");
1020  }
1021 
1022  if(numberOfNeighbors > _kdTree.getPoints()){
1023  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
1024  "Asked for more neighbors than you have data points\n");
1025  }
1026 
1027 
1028  int i,j,ii;
1029  T fbar;
1030 
1031  ndarray::Array < T,1,1 > covarianceTestPoint;
1032  ndarray::Array < int,1,1 > neighbors;
1033  ndarray::Array < double,1,1 > neighborDistances,vv;
1034 
1035  Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > covariance,bb,xx;
1036  Eigen::LDLT < Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > > ldlt;
1037 
1038  _timer.start();
1039 
1040 
1041  bb.resize(numberOfNeighbors,1);
1042  xx.resize(numberOfNeighbors,1);
1043  covariance.resize(numberOfNeighbors,numberOfNeighbors);
1044  covarianceTestPoint = allocate(ndarray::makeVector(numberOfNeighbors));
1045  neighbors = allocate(ndarray::makeVector(numberOfNeighbors));
1046  neighborDistances = allocate(ndarray::makeVector(numberOfNeighbors));
1047 
1049 
1050 
1051  if(_useMaxMin == 1) {
1052  //if you constructed this Gaussian process with minimum and maximum
1053  //values for the dimensions of your parameter space,
1054  //the point you are interpolating must be scaled to match the data so
1055  //that the selected nearest neighbors are appropriate
1056 
1057  for(i = 0; i < _dimensions; i++ )vv[i] = (vin[i] - _min[i])/(_max[i] - _min[i]);
1058  }
1059  else {
1060  vv = vin;
1061  }
1062 
1063 
1064  _kdTree.findNeighbors(neighbors, neighborDistances, vv, numberOfNeighbors);
1065 
1066  _timer.addToSearch();
1067 
1068  for(i = 0; i < numberOfNeighbors; i++ ) {
1069  covarianceTestPoint[i] = (*_covariogram)(vv, _kdTree.getData(neighbors[i]));
1070 
1071  covariance(i,i) = (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[i]))
1072  + _lambda;
1073 
1074  for(j = i + 1; j < numberOfNeighbors; j++ ){
1075  covariance(i,j) = (*_covariogram)(_kdTree.getData(neighbors[i]),
1076  _kdTree.getData(neighbors[j]));
1077  covariance(j,i) = covariance(i, j);
1078  }
1079  }
1080 
1082 
1083  //use Eigen's ldlt solver in place of matrix inversion (for speed purposes)
1084  ldlt.compute(covariance);
1085 
1086  for(ii = 0; ii < _nFunctions; ii++ ) {
1087 
1088  fbar = 0.0;
1089  for(i = 0; i < numberOfNeighbors; i++ )fbar += _function[neighbors[i]][ii];
1090  fbar = fbar/double(numberOfNeighbors);
1091 
1092  for(i = 0; i < numberOfNeighbors; i++ )bb(i,0) = _function[neighbors[i]][ii] - fbar;
1093  xx = ldlt.solve(bb);
1094 
1095  mu[ii] = fbar;
1096 
1097  for(i = 0; i < numberOfNeighbors; i++ ) {
1098  mu[ii] += covarianceTestPoint[i]*xx(i, 0);
1099  }
1100 
1101  }//ii = 0 through _nFunctions
1102 
1103  _timer.addToEigen();
1104 
1105  variance[0] = (*_covariogram)(vv, vv) + _lambda;
1106 
1107  for(i = 0; i < numberOfNeighbors; i++ )bb(i) = covarianceTestPoint[i];
1108 
1109  xx = ldlt.solve(bb);
1110 
1111  for(i = 0; i < numberOfNeighbors; i++ ) {
1112  variance[0] -= covarianceTestPoint[i]*xx(i, 0);
1113  }
1114  variance[0] = variance[0]*_krigingParameter;
1115 
1116 
1117  for(i = 1; i < _nFunctions; i++ )variance[i] = variance[0];
1118 
1120  _timer.addToTotal(1);
1121 }
ndarray::Array< T, 1, 1 > _max
void addToEigen()
Adds time to _eigenTime.
void addToSearch()
Adds time to _searchTime.
void addToIteration()
Adds time to _iterationTime.
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
ndarray::Array< T, 1, 1 > _min
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
ndarray::Array< T, 2, 2 > _function
void start()
Starts the timer for an individual call to an interpolation routine.
void addToVariance()
Adds time to _varianceTime.
void addToTotal(int i)
Adds time to _totalTime and adds counts to _interpolationCount.
template<typename T >
void lsst.afw.math::GaussianProcess< T >::removePoint ( int  dex)

This will remove a point from the data set.

Parameters
[in]dexthe index of the point you want to remove from your data set
Exceptions
pex_exceptionsRuntimeError 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 1782 of file GaussianProcess.cc.

1783 {
1784 
1785  int i,j;
1786 
1787  _kdTree.removePoint(dex);
1788 
1789  for(i = dex; i < _pts; i++ ) {
1790  for(j = 0; j < _nFunctions; j++ ) {
1791  _function[i][j] = _function[i + 1][j];
1792  }
1793  }
1794  _pts = _kdTree.getPoints();
1795 }
ndarray::Array< T, 2, 2 > _function
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]variancea one-dimensional ndarray. The value of the variance predicted by the Gaussian process will be stored in the zeroth element
[in]dexthe index of the point you wish to self interpolate
[in]numberOfNeighborsthe number of nearest neighbors to be used in the interpolation
Exceptions
pex_exceptionsRuntimeError 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 1125 of file GaussianProcess.cc.

1128 {
1129 
1130  if(numberOfNeighbors <= 0){
1131  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
1132  "Asked for zero or negative number of neighbors\n");
1133  }
1134 
1135  if(numberOfNeighbors > _kdTree.getPoints()){
1136  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
1137  "Asked for more neighbors than you have data points\n");
1138  }
1139 
1140  int i,j;
1141  T fbar,mu;
1142 
1143  ndarray::Array < T,1,1 > covarianceTestPoint;
1144  ndarray::Array < int,1,1 > selfNeighbors;
1145  ndarray::Array < double,1,1 > selfDistances;
1146  ndarray::Array < int,1,1 > neighbors;
1147  ndarray::Array < double,1,1 > neighborDistances;
1148 
1149  Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > covariance,bb,xx;
1150  Eigen::LDLT < Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > > ldlt;
1151 
1152  _timer.start();
1153 
1154  bb.resize(numberOfNeighbors, 1);
1155  xx.resize(numberOfNeighbors, 1);
1156  covariance.resize(numberOfNeighbors,numberOfNeighbors);
1157  covarianceTestPoint = allocate(ndarray::makeVector(numberOfNeighbors));
1158  neighbors = allocate(ndarray::makeVector(numberOfNeighbors));
1159  neighborDistances = allocate(ndarray::makeVector(numberOfNeighbors));
1160 
1161  selfNeighbors = allocate(ndarray::makeVector(numberOfNeighbors + 1));
1162  selfDistances = allocate(ndarray::makeVector(numberOfNeighbors + 1));
1163 
1164  //we don't use _useMaxMin because the data has already been normalized
1165 
1166 
1167  _kdTree.findNeighbors(selfNeighbors, selfDistances, _kdTree.getData(dex),
1168  numberOfNeighbors + 1);
1169 
1170  _timer.addToSearch();
1171 
1172  if(selfNeighbors[0]!= dex) {
1173  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
1174  "Nearest neighbor search in selfInterpolate did not find self\n");
1175  }
1176 
1177  //SelfNeighbors[0] will be the point itself (it is its own nearest neighbor)
1178  //We discard that for the interpolation calculation
1179  //
1180  //If you do not wish to do this, simply call the usual ::interpolate() method instead of
1181  //::selfInterpolate()
1182  for(i = 0; i < numberOfNeighbors; i++ ){
1183  neighbors[i] = selfNeighbors[i + 1];
1184  neighborDistances[i] = selfDistances[i + 1];
1185  }
1186 
1187  fbar = 0.0;
1188  for(i = 0; i < numberOfNeighbors; i++ ) fbar += _function[neighbors[i]][0];
1189  fbar = fbar/double(numberOfNeighbors);
1190 
1191  for(i = 0; i < numberOfNeighbors; i++ ){
1192  covarianceTestPoint[i] = (*_covariogram)(_kdTree.getData(dex),
1193  _kdTree.getData(neighbors[i]));
1194 
1195  covariance(i, i) = (*_covariogram)(_kdTree.getData(neighbors[i]),
1196  _kdTree.getData(neighbors[i])) + _lambda;
1197 
1198  for(j = i + 1; j < numberOfNeighbors; j++ ) {
1199  covariance(i, j) = (*_covariogram)(_kdTree.getData(neighbors[i]),
1200  _kdTree.getData(neighbors[j]));
1201  covariance(j, i) = covariance(i ,j);
1202  }
1203  }
1205 
1206  //use Eigen's ldlt solver in place of matrix inversion (for speed purposes)
1207  ldlt.compute(covariance);
1208 
1209 
1210  for(i = 0; i < numberOfNeighbors;i++ ) bb(i, 0) = _function[neighbors[i]][0] - fbar;
1211  xx = ldlt.solve(bb);
1212  _timer.addToEigen();
1213 
1214  mu = fbar;
1215 
1216  for(i = 0; i < numberOfNeighbors; i++ ) {
1217  mu += covarianceTestPoint[i]*xx(i, 0);
1218  }
1219 
1220 
1221  variance(0) = (*_covariogram)(_kdTree.getData(dex), _kdTree.getData(dex)) + _lambda;
1222 
1223  for(i = 0; i < numberOfNeighbors; i++ )bb(i) = covarianceTestPoint[i];
1224 
1225  xx = ldlt.solve(bb);
1226 
1227  for(i = 0; i < numberOfNeighbors; i++ ){
1228  variance(0) -= covarianceTestPoint[i]*xx(i,0);
1229  }
1230 
1231  variance(0) = variance(0)*_krigingParameter;
1233  _timer.addToTotal(1);
1234 
1235  return mu;
1236 }
void addToEigen()
Adds time to _eigenTime.
void addToSearch()
Adds time to _searchTime.
void addToIteration()
Adds time to _iterationTime.
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
ndarray::Array< T, 2, 2 > _function
void start()
Starts the timer for an individual call to an interpolation routine.
void addToVariance()
Adds time to _varianceTime.
void addToTotal(int i)
Adds time to _totalTime and adds counts to _interpolationCount.
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]muthis is where the interpolated function values will be stored
[out]variancethe variance on mu will be stored here
[in]dexthe index of the point you wish to interpolate
[in]numberOfNeighborsthe number of nearest neighbors to use in the interpolation
Exceptions
pex_exceptionsRuntimeError if the nearest neighbor search does not find the data point itself as the nearest neighbor

Definition at line 1239 of file GaussianProcess.cc.

1242  {
1243 
1244  if(numberOfNeighbors <= 0){
1245  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
1246  "Asked for zero or negative number of neighbors\n");
1247  }
1248 
1249  if(numberOfNeighbors + 1 > _kdTree.getPoints()){
1250  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
1251  "Asked for more neighbors than you have data points\n");
1252  }
1253 
1254  if(dex < 0 || dex >=_kdTree.getPoints()){
1255  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
1256  "Asked to self interpolate on a point that does not exist\n");
1257  }
1258 
1259  int i,j,ii;
1260  T fbar;
1261 
1262  ndarray::Array < T,1,1 > covarianceTestPoint;
1263  ndarray::Array < int,1,1 > selfNeighbors;
1264  ndarray::Array < double,1,1 > selfDistances;
1265  ndarray::Array < int,1,1 > neighbors;
1266  ndarray::Array < double,1,1 > neighborDistances;
1267 
1268  Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > covariance,bb,xx;
1269  Eigen::LDLT < Eigen::Matrix < T,Eigen::Dynamic,Eigen::Dynamic > > ldlt;
1270 
1271  _timer.start();
1272 
1273  bb.resize(numberOfNeighbors, 1);
1274  xx.resize(numberOfNeighbors, 1);
1275  covariance.resize(numberOfNeighbors,numberOfNeighbors);
1276  covarianceTestPoint = allocate(ndarray::makeVector(numberOfNeighbors));
1277  neighbors = allocate(ndarray::makeVector(numberOfNeighbors));
1278  neighborDistances = allocate(ndarray::makeVector(numberOfNeighbors));
1279 
1280  selfNeighbors = allocate(ndarray::makeVector(numberOfNeighbors + 1));
1281  selfDistances = allocate(ndarray::makeVector(numberOfNeighbors + 1));
1282 
1283  //we don't use _useMaxMin because the data has already been normalized
1284 
1285 
1286  _kdTree.findNeighbors(selfNeighbors, selfDistances, _kdTree.getData(dex),
1287  numberOfNeighbors + 1);
1288 
1289  _timer.addToSearch();
1290 
1291  if(selfNeighbors[0]!= dex) {
1292 
1293  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
1294  "Nearest neighbor search in selfInterpolate did not find self\n");
1295  }
1296 
1297  //SelfNeighbors[0] will be the point itself (it is its own nearest neighbor)
1298  //We discard that for the interpolation calculation
1299  //
1300  //If you do not wish to do this, simply call the usual ::interpolate() method instead of
1301  //::selfInterpolate()
1302  for(i = 0; i < numberOfNeighbors; i++ ) {
1303  neighbors[i] = selfNeighbors[i + 1];
1304  neighborDistances[i] = selfDistances[i + 1];
1305  }
1306 
1307 
1308 
1309  for(i = 0; i < numberOfNeighbors; i++ ) {
1310  covarianceTestPoint[i] = (*_covariogram)(_kdTree.getData(dex),_kdTree.getData(neighbors[i]));
1311 
1312  covariance(i, i) = (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[i]))
1313  + _lambda;
1314 
1315  for(j = i + 1; j < numberOfNeighbors; j++ ) {
1316  covariance(i, j) = (*_covariogram)(_kdTree.getData(neighbors[i]),
1317  _kdTree.getData(neighbors[j]));
1318  covariance(j, i) = covariance(i, j);
1319  }
1320  }
1322 
1323 
1324  //use Eigen's ldlt solver in place of matrix inversion (for speed purposes)
1325  ldlt.compute(covariance);
1326 
1327  for(ii = 0; ii < _nFunctions; ii++ ) {
1328 
1329  fbar = 0.0;
1330  for(i = 0; i < numberOfNeighbors; i++ )fbar += _function[neighbors[i]][ii];
1331  fbar = fbar/double(numberOfNeighbors);
1332 
1333  for(i = 0; i < numberOfNeighbors; i++ )bb(i,0) = _function[neighbors[i]][ii] - fbar;
1334  xx = ldlt.solve(bb);
1335 
1336  mu[ii] = fbar;
1337 
1338  for(i = 0; i < numberOfNeighbors; i++ ){
1339  mu[ii] += covarianceTestPoint[i]*xx(i,0);
1340  }
1341  }//ii = 0 through _nFunctions
1342 
1343  _timer.addToEigen();
1344 
1345  variance[0] = (*_covariogram)(_kdTree.getData(dex), _kdTree.getData(dex)) + _lambda;
1346 
1347  for(i = 0; i < numberOfNeighbors; i++ )bb(i) = covarianceTestPoint[i];
1348 
1349  xx = ldlt.solve(bb);
1350 
1351  for(i = 0; i < numberOfNeighbors; i++ ) {
1352  variance[0] -= covarianceTestPoint[i]*xx(i,0);
1353  }
1354 
1355  variance[0] = variance[0]*_krigingParameter;
1356 
1357  for(i = 1; i < _nFunctions; i++ )variance[i] = variance[0];
1358 
1360  _timer.addToTotal(1);
1361 }
void addToEigen()
Adds time to _eigenTime.
void addToSearch()
Adds time to _searchTime.
void addToIteration()
Adds time to _iterationTime.
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
ndarray::Array< T, 2, 2 > _function
void start()
Starts the timer for an individual call to an interpolation routine.
void addToVariance()
Adds time to _varianceTime.
void addToTotal(int i)
Adds time to _totalTime and adds counts to _interpolationCount.
template<typename T >
void lsst.afw.math::GaussianProcess< T >::setCovariogram ( boost::shared_ptr< Covariogram< T > > const &  covar)

Assign a different covariogram to this GaussianProcess.

Parameters
[in]covarthe Covariogram object that you wish to assign

Definition at line 1804 of file GaussianProcess.cc.

1804  {
1805  _covariogram = covar;
1806 }
boost::shared_ptr< Covariogram< T > > _covariogram
template<typename T >
void lsst.afw.math::GaussianProcess< T >::setKrigingParameter ( kk)

Assign a value to the Kriging paramter.

Parameters
[in]kkthe value assigned to the Kriging parameters

Definition at line 1798 of file GaussianProcess.cc.

1799 {
1800  _krigingParameter = kk;
1801 }
template<typename T >
void lsst.afw.math::GaussianProcess< T >::setLambda ( lambda)

set the value of the hyperparameter _lambda

Parameters
[in]lambdathe 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 1809 of file GaussianProcess.cc.

1809  {
1810  _lambda = lambda;
1811 }

Member Data Documentation

template<typename T >
boost::shared_ptr< Covariogram<T> > lsst.afw.math::GaussianProcess< T >::_covariogram
private

Definition at line 825 of file GaussianProcess.h.

template<typename T >
int lsst.afw.math::GaussianProcess< T >::_dimensions
private

Definition at line 816 of file GaussianProcess.h.

template<typename T >
ndarray::Array<T,2,2> lsst.afw.math::GaussianProcess< T >::_function
private

Definition at line 821 of file GaussianProcess.h.

template<typename T >
KdTree<T> lsst.afw.math::GaussianProcess< T >::_kdTree
private

Definition at line 823 of file GaussianProcess.h.

template<typename T >
T lsst.afw.math::GaussianProcess< T >::_krigingParameter
private

Definition at line 818 of file GaussianProcess.h.

template<typename T >
T lsst.afw.math::GaussianProcess< T >::_lambda
private

Definition at line 818 of file GaussianProcess.h.

template<typename T >
ndarray::Array<T,1,1> lsst.afw.math::GaussianProcess< T >::_max
private

Definition at line 820 of file GaussianProcess.h.

template<typename T >
ndarray::Array<T,1,1> lsst.afw.math::GaussianProcess< T >::_min
private

Definition at line 820 of file GaussianProcess.h.

template<typename T >
int lsst.afw.math::GaussianProcess< T >::_nFunctions
private

Definition at line 816 of file GaussianProcess.h.

template<typename T >
int lsst.afw.math::GaussianProcess< T >::_pts
private

Definition at line 816 of file GaussianProcess.h.

template<typename T >
int lsst.afw.math::GaussianProcess< T >::_room
private

Definition at line 816 of file GaussianProcess.h.

template<typename T >
int lsst.afw.math::GaussianProcess< T >::_roomStep
private

Definition at line 816 of file GaussianProcess.h.

template<typename T >
GaussianProcessTimer lsst.afw.math::GaussianProcess< T >::_timer
mutableprivate

Definition at line 826 of file GaussianProcess.h.

template<typename T >
int lsst.afw.math::GaussianProcess< T >::_useMaxMin
private

Definition at line 816 of file GaussianProcess.h.


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