LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Public Member Functions | 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>

Public Member Functions

 GaussianProcess (const GaussianProcess &)=delete
 
GaussianProcessoperator= (const GaussianProcess &)=delete
 
 GaussianProcess (GaussianProcess &&)=delete
 
GaussianProcessoperator= (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. 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, 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. More...
 
 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 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, 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 More...
 
int getNPoints () const
 return the number of data points stored in the GaussianProcess More...
 
int getDim () const
 return the dimensionality of data points stored in the GaussianProcess More...
 
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. More...
 
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. 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. 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 (std::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...
 

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]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 751 of file GaussianProcess.cc.

754 {
755  int i;
756 
757  _covariogram = covarIn;
758 
759  _npts = dataIn.template getSize<0>();
760  _dimensions = dataIn.template getSize<1>();
761 
762  _room = _npts;
763  _roomStep = 5000;
764 
765  _nFunctions = 1;
766  _function = allocate(ndarray::makeVector(_npts, 1));
767 
768  if (ff.getNumElements() != static_cast<ndarray::Size>(_npts)) {
770  "You did not pass in the same number of data points as function values\n");
771  }
772 
773  for (i = 0; i < _npts; i++) _function[i][0] = ff[i];
774 
775  _krigingParameter = T(1.0);
776  _lambda = T(1.0e-5);
777 
778  _useMaxMin = 0;
779 
780  _kdTree.Initialize(dataIn);
781 
782  _npts = _kdTree.getNPoints();
783 }
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104

◆ 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]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 786 of file GaussianProcess.cc.

788  {
789  int i, j;
790  ndarray::Array<T, 2, 2> normalizedData;
791 
792  _covariogram = covarIn;
793 
794  _npts = dataIn.template getSize<0>();
795  _dimensions = dataIn.template getSize<1>();
796  _room = _npts;
797  _roomStep = 5000;
798 
799  if (ff.getNumElements() != static_cast<ndarray::Size>(_npts)) {
801  "You did not pass in the same number of data points as function values\n");
802  }
803 
804  if (mn.getNumElements() != static_cast<ndarray::Size>(_dimensions) ||
805  mx.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
807  "Your min/max values have different dimensionality than your data points\n");
808  }
809 
810  _krigingParameter = T(1.0);
811 
812  _lambda = T(1.0e-5);
813  _krigingParameter = T(1.0);
814 
815  _max = allocate(ndarray::makeVector(_dimensions));
816  _min = allocate(ndarray::makeVector(_dimensions));
817  _max.deep() = mx;
818  _min.deep() = mn;
819  _useMaxMin = 1;
820  normalizedData = allocate(ndarray::makeVector(_npts, _dimensions));
821  for (i = 0; i < _npts; i++) {
822  for (j = 0; j < _dimensions; j++) {
823  normalizedData[i][j] = (dataIn[i][j] - _min[j]) / (_max[j] - _min[j]);
824  // note the normalization by _max - _min in each dimension
825  }
826  }
827 
828  _kdTree.Initialize(normalizedData);
829 
830  _npts = _kdTree.getNPoints();
831  _nFunctions = 1;
832  _function = allocate(ndarray::makeVector(_npts, 1));
833  for (i = 0; i < _npts; i++) _function[i][0] = ff[i];
834 }

◆ 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]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 837 of file GaussianProcess.cc.

838  {
839  _covariogram = covarIn;
840 
841  _npts = dataIn.template getSize<0>();
842  _dimensions = dataIn.template getSize<1>();
843 
844  _room = _npts;
845  _roomStep = 5000;
846 
847  if (ff.template getSize<0>() != static_cast<ndarray::Size>(_npts)) {
849  "You did not pass in the same number of data points as function values\n");
850  }
851 
852  _nFunctions = ff.template getSize<1>();
853  _function = allocate(ndarray::makeVector(_npts, _nFunctions));
854  _function.deep() = ff;
855 
856  _krigingParameter = T(1.0);
857 
858  _lambda = T(1.0e-5);
859 
860  _useMaxMin = 0;
861 
862  _kdTree.Initialize(dataIn);
863 
864  _npts = _kdTree.getNPoints();
865 }

◆ 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]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 868 of file GaussianProcess.cc.

870  {
871  int i, j;
872  ndarray::Array<T, 2, 2> normalizedData;
873 
874  _covariogram = covarIn;
875 
876  _npts = dataIn.template getSize<0>();
877  _dimensions = dataIn.template getSize<1>();
878 
879  _room = _npts;
880  _roomStep = 5000;
881 
882  if (ff.template getSize<0>() != static_cast<ndarray::Size>(_npts)) {
884  "You did not pass in the same number of data points as function values\n");
885  }
886 
887  if (mn.getNumElements() != static_cast<ndarray::Size>(_dimensions) ||
888  mx.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
890  "Your min/max values have different dimensionality than your data points\n");
891  }
892 
893  _krigingParameter = T(1.0);
894 
895  _lambda = T(1.0e-5);
896  _krigingParameter = T(1.0);
897 
898  _max = allocate(ndarray::makeVector(_dimensions));
899  _min = allocate(ndarray::makeVector(_dimensions));
900  _max.deep() = mx;
901  _min.deep() = mn;
902  _useMaxMin = 1;
903  normalizedData = allocate(ndarray::makeVector(_npts, _dimensions));
904  for (i = 0; i < _npts; i++) {
905  for (j = 0; j < _dimensions; j++) {
906  normalizedData[i][j] = (dataIn[i][j] - _min[j]) / (_max[j] - _min[j]);
907  // note the normalization by _max - _min in each dimension
908  }
909  }
910 
911  _kdTree.Initialize(normalizedData);
912  _npts = _kdTree.getNPoints();
913  _nFunctions = ff.template getSize<1>();
914  _function = allocate(ndarray::makeVector(_npts, _nFunctions));
915  _function.deep() = ff;
916 }

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::RuntimeErrorif 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.

1907  {
1908  if (vin.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
1910  "You are trying to add a point of the wrong dimensionality to "
1911  "your GaussianProcess.\n");
1912  }
1913 
1914  if (f.template getSize<0>() != static_cast<ndarray::Size>(_nFunctions)) {
1916  "You are not adding the correct number of function values to "
1917  "your GaussianProcess.\n");
1918  }
1919 
1920  int i, j;
1921 
1922  ndarray::Array<T, 1, 1> v;
1923  v = allocate(ndarray::makeVector(_dimensions));
1924 
1925  for (i = 0; i < _dimensions; i++) {
1926  v[i] = vin[i];
1927  if (_useMaxMin == 1) {
1928  v[i] = (v[i] - _min[i]) / (_max[i] - _min[i]);
1929  }
1930  }
1931 
1932  if (_npts == _room) {
1933  ndarray::Array<T, 2, 2> buff;
1934  buff = allocate(ndarray::makeVector(_npts, _nFunctions));
1935  buff.deep() = _function;
1936 
1937  _room += _roomStep;
1938  _function = allocate(ndarray::makeVector(_room, _nFunctions));
1939  for (i = 0; i < _npts; i++) {
1940  for (j = 0; j < _nFunctions; j++) {
1941  _function[i][j] = buff[i][j];
1942  }
1943  }
1944  }
1945  for (i = 0; i < _nFunctions; i++) _function[_npts][i] = f[i];
1946 
1947  _kdTree.addPoint(v);
1948  _npts = _kdTree.getNPoints();
1949 }

◆ addPoint() [2/2]

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::exceptions::RuntimeErrorif you call this when you should have called the version taking a vector of functions (below)
pex::exceptions::RuntimeErrorif 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.

1862  {
1863  int i, j;
1864 
1865  if (_nFunctions != 1) {
1867  "You are calling the wrong addPoint; you need a vector of functions\n");
1868  }
1869 
1870  if (vin.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
1872  "You are trying to add a point of the wrong dimensionality to "
1873  "your GaussianProcess.\n");
1874  }
1875 
1876  ndarray::Array<T, 1, 1> v;
1877  v = allocate(ndarray::makeVector(_dimensions));
1878 
1879  for (i = 0; i < _dimensions; i++) {
1880  v[i] = vin[i];
1881  if (_useMaxMin == 1) {
1882  v[i] = (v[i] - _min[i]) / (_max[i] - _min[i]);
1883  }
1884  }
1885 
1886  if (_npts == _room) {
1887  ndarray::Array<T, 2, 2> buff;
1888  buff = allocate(ndarray::makeVector(_npts, _nFunctions));
1889  buff.deep() = _function;
1890 
1891  _room += _roomStep;
1892  _function = allocate(ndarray::makeVector(_room, _nFunctions));
1893  for (i = 0; i < _npts; i++) {
1894  for (j = 0; j < _nFunctions; j++) {
1895  _function[i][j] = buff[i][j];
1896  }
1897  }
1898  }
1899 
1900  _function[_npts][0] = f;
1901 
1902  _kdTree.addPoint(v);
1903  _npts = _kdTree.getNPoints();
1904 }

◆ 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]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 _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.

1489  {
1490  int i, j;
1491 
1492  ndarray::Size nQueries = queries.template getSize<0>();
1493 
1494  if (_nFunctions != 1) {
1496  "Your mu and variance arrays do not have room for all of the functions "
1497  "as you are trying to interpolate\n");
1498  }
1499 
1500  if (mu.getNumElements() != nQueries || variance.getNumElements() != nQueries) {
1502  "Your mu and variance arrays do not have room for all of the points "
1503  "at which you are trying to interpolate your function.\n");
1504  }
1505 
1506  if (queries.template getSize<1>() != static_cast<ndarray::Size>(_dimensions)) {
1508  "The points you passed to batchInterpolate are of the wrong "
1509  "dimensionality for your Gaussian Process\n");
1510  }
1511 
1512  T fbar;
1513  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> batchCovariance, batchbb, batchxx;
1514  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> queryCovariance;
1515  Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
1516 
1517  ndarray::Array<T, 1, 1> v1;
1518 
1519  _timer.start();
1520 
1521  v1 = allocate(ndarray::makeVector(_dimensions));
1522  batchbb.resize(_npts, 1);
1523  batchxx.resize(_npts, 1);
1524  batchCovariance.resize(_npts, _npts);
1525  queryCovariance.resize(_npts, 1);
1526 
1527  for (i = 0; i < _npts; i++) {
1528  batchCovariance(i, i) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(i)) + _lambda;
1529  for (j = i + 1; j < _npts; j++) {
1530  batchCovariance(i, j) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(j));
1531  batchCovariance(j, i) = batchCovariance(i, j);
1532  }
1533  }
1534  _timer.addToIteration();
1535 
1536  ldlt.compute(batchCovariance);
1537 
1538  fbar = 0.0;
1539  for (i = 0; i < _npts; i++) {
1540  fbar += _function[i][0];
1541  }
1542  fbar = fbar / T(_npts);
1543 
1544  for (i = 0; i < _npts; i++) {
1545  batchbb(i, 0) = _function[i][0] - fbar;
1546  }
1547  batchxx = ldlt.solve(batchbb);
1548  _timer.addToEigen();
1549 
1550  for (ndarray::Size ii = 0; ii < nQueries; ii++) {
1551  for (i = 0; i < _dimensions; i++) v1[i] = queries[ii][i];
1552  if (_useMaxMin == 1) {
1553  for (i = 0; i < _dimensions; i++) v1[i] = (v1[i] - _min[i]) / (_max[i] - _min[i]);
1554  }
1555  mu(ii) = fbar;
1556  for (i = 0; i < _npts; i++) {
1557  mu(ii) += batchxx(i) * (*_covariogram)(v1, _kdTree.getData(i));
1558  }
1559  }
1560  _timer.addToIteration();
1561 
1562  for (ndarray::Size ii = 0; ii < nQueries; ii++) {
1563  for (i = 0; i < _dimensions; i++) v1[i] = queries[ii][i];
1564  if (_useMaxMin == 1) {
1565  for (i = 0; i < _dimensions; i++) v1[i] = (v1[i] - _min[i]) / (_max[i] - _min[i]);
1566  }
1567 
1568  for (i = 0; i < _npts; i++) {
1569  batchbb(i, 0) = (*_covariogram)(v1, _kdTree.getData(i));
1570  queryCovariance(i, 0) = batchbb(i, 0);
1571  }
1572  batchxx = ldlt.solve(batchbb);
1573 
1574  variance[ii] = (*_covariogram)(v1, v1) + _lambda;
1575 
1576  for (i = 0; i < _npts; i++) {
1577  variance[ii] -= queryCovariance(i, 0) * batchxx(i);
1578  }
1579 
1580  variance[ii] = variance[ii] * _krigingParameter;
1581  }
1582 
1583  _timer.addToVariance();
1584  _timer.addToTotal(nQueries);
1585 }
afw::table::Key< afw::table::Array< VariancePixelT > > variance
void addToVariance()
Adds time to _varianceTime.
void addToIteration()
Adds time to _iterationTime.
void start()
Starts the timer for an individual call to an interpolation routine.
void addToEigen()
Adds time to _eigenTime.
void addToTotal(int i)
Adds time to _totalTime and adds counts to _interpolationCount.

◆ 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]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 _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.

1697  {
1698  int i, j;
1699 
1700  ndarray::Size nQueries = queries.template getSize<0>();
1701 
1702  if (_nFunctions != 1) {
1704  "Your output array does not have enough room for all of the functions "
1705  "you are trying to interpolate.\n");
1706  }
1707 
1708  if (queries.template getSize<1>() != static_cast<ndarray::Size>(_dimensions)) {
1710  "The points at which you are trying to interpolate your function are "
1711  "of the wrong dimensionality.\n");
1712  }
1713 
1714  if (mu.getNumElements() != nQueries) {
1716  "Your output array does not have enough room for all of the points "
1717  "at which you are trying to interpolate your function.\n");
1718  }
1719 
1720  T fbar;
1721  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> batchCovariance, batchbb, batchxx;
1722  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> queryCovariance;
1723  Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
1724 
1725  ndarray::Array<T, 1, 1> v1;
1726 
1727  _timer.start();
1728 
1729  v1 = allocate(ndarray::makeVector(_dimensions));
1730 
1731  batchbb.resize(_npts, 1);
1732  batchxx.resize(_npts, 1);
1733  batchCovariance.resize(_npts, _npts);
1734  queryCovariance.resize(_npts, 1);
1735 
1736  for (i = 0; i < _npts; i++) {
1737  batchCovariance(i, i) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(i)) + _lambda;
1738  for (j = i + 1; j < _npts; j++) {
1739  batchCovariance(i, j) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(j));
1740  batchCovariance(j, i) = batchCovariance(i, j);
1741  }
1742  }
1743  _timer.addToIteration();
1744 
1745  ldlt.compute(batchCovariance);
1746 
1747  fbar = 0.0;
1748  for (i = 0; i < _npts; i++) {
1749  fbar += _function[i][0];
1750  }
1751  fbar = fbar / T(_npts);
1752 
1753  for (i = 0; i < _npts; i++) {
1754  batchbb(i, 0) = _function[i][0] - fbar;
1755  }
1756  batchxx = ldlt.solve(batchbb);
1757  _timer.addToEigen();
1758 
1759  for (ndarray::Size ii = 0; ii < nQueries; ii++) {
1760  for (i = 0; i < _dimensions; i++) v1[i] = queries[ii][i];
1761  if (_useMaxMin == 1) {
1762  for (i = 0; i < _dimensions; i++) v1[i] = (v1[i] - _min[i]) / (_max[i] - _min[i]);
1763  }
1764 
1765  mu(ii) = fbar;
1766  for (i = 0; i < _npts; i++) {
1767  mu(ii) += batchxx(i) * (*_covariogram)(v1, _kdTree.getData(i));
1768  }
1769  }
1770  _timer.addToIteration();
1771  _timer.addToTotal(nQueries);
1772 }

◆ 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.

1776  {
1777  int i, j, ifn;
1778 
1779  ndarray::Size nQueries = queries.template getSize<0>();
1780 
1781  if (mu.template getSize<0>() != nQueries) {
1783  "Your output array does not have enough room for all of the points "
1784  "at which you want to interpolate your functions.\n");
1785  }
1786 
1787  if (mu.template getSize<1>() != static_cast<ndarray::Size>(_nFunctions)) {
1789  "Your output array does not have enough room for all of the functions "
1790  "you are trying to interpolate.\n");
1791  }
1792 
1793  if (queries.template getSize<1>() != static_cast<ndarray::Size>(_dimensions)) {
1795  "The points at which you are interpolating your functions do not "
1796  "have the correct dimensionality.\n");
1797  }
1798 
1799  T fbar;
1800  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> batchCovariance, batchbb, batchxx;
1801  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> queryCovariance;
1802  Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
1803 
1804  ndarray::Array<T, 1, 1> v1;
1805 
1806  _timer.start();
1807 
1808  v1 = allocate(ndarray::makeVector(_dimensions));
1809 
1810  batchbb.resize(_npts, 1);
1811  batchxx.resize(_npts, 1);
1812  batchCovariance.resize(_npts, _npts);
1813  queryCovariance.resize(_npts, 1);
1814 
1815  for (i = 0; i < _npts; i++) {
1816  batchCovariance(i, i) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(i)) + _lambda;
1817  for (j = i + 1; j < _npts; j++) {
1818  batchCovariance(i, j) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(j));
1819  batchCovariance(j, i) = batchCovariance(i, j);
1820  }
1821  }
1822 
1823  _timer.addToIteration();
1824 
1825  ldlt.compute(batchCovariance);
1826 
1827  _timer.addToEigen();
1828 
1829  for (ifn = 0; ifn < _nFunctions; ifn++) {
1830  fbar = 0.0;
1831  for (i = 0; i < _npts; i++) {
1832  fbar += _function[i][ifn];
1833  }
1834  fbar = fbar / T(_npts);
1835 
1836  _timer.addToIteration();
1837 
1838  for (i = 0; i < _npts; i++) {
1839  batchbb(i, 0) = _function[i][ifn] - fbar;
1840  }
1841  batchxx = ldlt.solve(batchbb);
1842  _timer.addToEigen();
1843 
1844  for (ndarray::Size ii = 0; ii < nQueries; ii++) {
1845  for (i = 0; i < _dimensions; i++) v1[i] = queries[ii][i];
1846  if (_useMaxMin == 1) {
1847  for (i = 0; i < _dimensions; i++) v1[i] = (v1[i] - _min[i]) / (_max[i] - _min[i]);
1848  }
1849 
1850  mu[ii][ifn] = fbar;
1851  for (i = 0; i < _npts; i++) {
1852  mu[ii][ifn] += batchxx(i) * (*_covariogram)(v1, _kdTree.getData(i));
1853  }
1854  }
1855 
1856  } // ifn = 0 through _nFunctions
1857 
1858  _timer.addToTotal(nQueries);
1859 }

◆ 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.

1589  {
1590  int i, j, ifn;
1591 
1592  ndarray::Size nQueries = queries.template getSize<0>();
1593 
1594  if (mu.template getSize<0>() != nQueries || variance.template getSize<0>() != nQueries) {
1596  "Your output arrays do not have room for all of the points at which "
1597  "you are interpolating your functions.\n");
1598  }
1599 
1600  if (mu.template getSize<1>() != static_cast<ndarray::Size>(_nFunctions) ||
1601  variance.template getSize<1>() != static_cast<ndarray::Size>(_nFunctions)) {
1603  "Your output arrays do not have room for all of the functions you are "
1604  "interpolating\n");
1605  }
1606 
1607  if (queries.template getSize<1>() != static_cast<ndarray::Size>(_dimensions)) {
1609  "The points at which you are interpolating your functions have the "
1610  "wrong dimensionality.\n");
1611  }
1612 
1613  T fbar;
1614  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> batchCovariance, batchbb, batchxx;
1615  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> queryCovariance;
1616  Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
1617 
1618  ndarray::Array<T, 1, 1> v1;
1619 
1620  _timer.start();
1621 
1622  v1 = allocate(ndarray::makeVector(_dimensions));
1623  batchbb.resize(_npts, 1);
1624  batchxx.resize(_npts, 1);
1625  batchCovariance.resize(_npts, _npts);
1626  queryCovariance.resize(_npts, 1);
1627 
1628  for (i = 0; i < _npts; i++) {
1629  batchCovariance(i, i) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(i)) + _lambda;
1630  for (j = i + 1; j < _npts; j++) {
1631  batchCovariance(i, j) = (*_covariogram)(_kdTree.getData(i), _kdTree.getData(j));
1632  batchCovariance(j, i) = batchCovariance(i, j);
1633  }
1634  }
1635 
1636  _timer.addToIteration();
1637 
1638  ldlt.compute(batchCovariance);
1639 
1640  _timer.addToEigen();
1641 
1642  for (ifn = 0; ifn < _nFunctions; ifn++) {
1643  fbar = 0.0;
1644  for (i = 0; i < _npts; i++) {
1645  fbar += _function[i][ifn];
1646  }
1647  fbar = fbar / T(_npts);
1648  _timer.addToIteration();
1649 
1650  for (i = 0; i < _npts; i++) {
1651  batchbb(i, 0) = _function[i][ifn] - fbar;
1652  }
1653  batchxx = ldlt.solve(batchbb);
1654  _timer.addToEigen();
1655 
1656  for (ndarray::Size ii = 0; ii < nQueries; ii++) {
1657  for (i = 0; i < _dimensions; i++) v1[i] = queries[ii][i];
1658  if (_useMaxMin == 1) {
1659  for (i = 0; i < _dimensions; i++) v1[i] = (v1[i] - _min[i]) / (_max[i] - _min[i]);
1660  }
1661  mu[ii][ifn] = fbar;
1662  for (i = 0; i < _npts; i++) {
1663  mu[ii][ifn] += batchxx(i) * (*_covariogram)(v1, _kdTree.getData(i));
1664  }
1665  }
1666 
1667  } // ifn = 0 to _nFunctions
1668 
1669  _timer.addToIteration();
1670  for (ndarray::Size ii = 0; ii < nQueries; ii++) {
1671  for (i = 0; i < _dimensions; i++) v1[i] = queries[ii][i];
1672  if (_useMaxMin == 1) {
1673  for (i = 0; i < _dimensions; i++) v1[i] = (v1[i] - _min[i]) / (_max[i] - _min[i]);
1674  }
1675 
1676  for (i = 0; i < _npts; i++) {
1677  batchbb(i, 0) = (*_covariogram)(v1, _kdTree.getData(i));
1678  queryCovariance(i, 0) = batchbb(i, 0);
1679  }
1680  batchxx = ldlt.solve(batchbb);
1681 
1682  variance[ii][0] = (*_covariogram)(v1, v1) + _lambda;
1683 
1684  for (i = 0; i < _npts; i++) {
1685  variance[ii][0] -= queryCovariance(i, 0) * batchxx(i);
1686  }
1687 
1688  variance[ii][0] = variance[ii][0] * _krigingParameter;
1689  for (i = 1; i < _nFunctions; i++) variance[ii][i] = variance[ii][0];
1690  }
1691  _timer.addToVariance();
1692  _timer.addToTotal(nQueries);
1693 }

◆ 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]ptswill contain the data points from the Gaussian Process
[out]fnwill contain the function values from the Gaussian Process
[in]indicesis an array of indices indicating the points to return

Definition at line 929 of file GaussianProcess.cc.

930  {
931  if (_nFunctions != 1) {
933  "Your function value array does not have enough room for all of the functions "
934  "in your GaussianProcess.\n");
935  }
936 
937  if (pts.template getSize<1>() != static_cast<ndarray::Size>(_dimensions)) {
939  "Your pts array is constructed for points of the wrong dimensionality.\n");
940  }
941 
942  if (pts.template getSize<0>() != indices.getNumElements()) {
944  "You did not put enough room in your pts array to fit all the points "
945  "you asked for in your indices array.\n");
946  }
947 
948  if (fn.template getSize<0>() != indices.getNumElements()) {
950  "You did not provide enough room in your function value array "
951  "for all of the points you requested in your indices array.\n");
952  }
953 
954  for (ndarray::Size i = 0; i < indices.template getSize<0>(); i++) {
955  pts[i] = _kdTree.getData(indices[i]); // do this first in case one of the indices is invalid.
956  // _kdTree.getData() will raise an exception in that case
957  fn[i] = _function[indices[i]][0];
958  if (_useMaxMin == 1) {
959  for (int j = 0; j < _dimensions; j++) {
960  pts[i][j] *= (_max[j] - _min[j]);
961  pts[i][j] += _min[j];
962  }
963  }
964  }
965 }

◆ 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]ptswill contain the data points from the Gaussian Process
[out]fnwill contain the function values from the Gaussian Process
[in]indicesis an array of indices indicating the points to return

Definition at line 968 of file GaussianProcess.cc.

969  {
970  if (fn.template getSize<1>() != static_cast<ndarray::Size>(_nFunctions)) {
972  "Your function value array does not have enough room for all of the functions "
973  "in your GaussianProcess.\n");
974  }
975 
976  if (pts.template getSize<1>() != static_cast<ndarray::Size>(_dimensions)) {
978  "Your pts array is constructed for points of the wrong dimensionality.\n");
979  }
980 
981  if (pts.template getSize<0>() != indices.getNumElements()) {
983  "You did not put enough room in your pts array to fit all the points "
984  "you asked for in your indices array.\n");
985  }
986 
987  if (fn.template getSize<0>() != indices.getNumElements()) {
989  "You did not provide enough room in your function value array "
990  "for all of the points you requested in your indices array.\n");
991  }
992 
993  for (ndarray::Size i = 0; i < indices.template getSize<0>(); i++) {
994  pts[i] = _kdTree.getData(indices[i]); // do this first in case one of the indices is invalid.
995  // _kdTree.getData() will raise an exception in that case
996  for (int j = 0; j < _nFunctions; j++) {
997  fn[i][j] = _function[indices[i]][j];
998  }
999  if (_useMaxMin == 1) {
1000  for (int j = 0; j < _dimensions; j++) {
1001  pts[i][j] *= (_max[j] - _min[j]);
1002  pts[i][j] += _min[j];
1003  }
1004  }
1005  }
1006 }

◆ getDim()

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

return the dimensionality of data points stored in the GaussianProcess

Definition at line 924 of file GaussianProcess.cc.

924  {
925  return _dimensions;
926 }

◆ getNPoints()

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

return the number of data points stored in the GaussianProcess

Definition at line 919 of file GaussianProcess.cc.

919  {
920  return _npts;
921 }

◆ getTimes()

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

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.

1981  {
1982  return _timer;
1983 }

◆ 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]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 1131 of file GaussianProcess.cc.

1132  {
1133  if (numberOfNeighbors <= 0) {
1135  "Asked for zero or negative number of neighbors\n");
1136  }
1137 
1138  if (numberOfNeighbors > _kdTree.getNPoints()) {
1140  "Asked for more neighbors than you have data points\n");
1141  }
1142 
1143  if (vin.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
1145  "You are interpolating at a point with different dimensionality than you data\n");
1146  }
1147 
1148  if (mu.getNumElements() != static_cast<ndarray::Size>(_nFunctions) ||
1149  variance.getNumElements() != static_cast<ndarray::Size>(_nFunctions)) {
1151  "Your mu and/or var arrays are improperly sized for the number of functions "
1152  "you are interpolating\n");
1153  }
1154 
1155  int i, j, ii;
1156  T fbar;
1157 
1158  ndarray::Array<T, 1, 1> covarianceTestPoint;
1159  ndarray::Array<int, 1, 1> neighbors;
1160  ndarray::Array<double, 1, 1> neighborDistances, vv;
1161 
1162  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> covariance, bb, xx;
1163  Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
1164 
1165  _timer.start();
1166 
1167  bb.resize(numberOfNeighbors, 1);
1168  xx.resize(numberOfNeighbors, 1);
1169  covariance.resize(numberOfNeighbors, numberOfNeighbors);
1170  covarianceTestPoint = allocate(ndarray::makeVector(numberOfNeighbors));
1171  neighbors = allocate(ndarray::makeVector(numberOfNeighbors));
1172  neighborDistances = allocate(ndarray::makeVector(numberOfNeighbors));
1173 
1174  vv = allocate(ndarray::makeVector(_dimensions));
1175 
1176  if (_useMaxMin == 1) {
1177  // if you constructed this Gaussian process with minimum and maximum
1178  // values for the dimensions of your parameter space,
1179  // the point you are interpolating must be scaled to match the data so
1180  // that the selected nearest neighbors are appropriate
1181 
1182  for (i = 0; i < _dimensions; i++) vv[i] = (vin[i] - _min[i]) / (_max[i] - _min[i]);
1183  } else {
1184  vv = vin;
1185  }
1186 
1187  _kdTree.findNeighbors(neighbors, neighborDistances, vv, numberOfNeighbors);
1188 
1189  _timer.addToSearch();
1190 
1191  for (i = 0; i < numberOfNeighbors; i++) {
1192  covarianceTestPoint[i] = (*_covariogram)(vv, _kdTree.getData(neighbors[i]));
1193 
1194  covariance(i, i) =
1195  (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[i])) + _lambda;
1196 
1197  for (j = i + 1; j < numberOfNeighbors; j++) {
1198  covariance(i, j) = (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[j]));
1199  covariance(j, i) = covariance(i, j);
1200  }
1201  }
1202 
1203  _timer.addToIteration();
1204 
1205  // use Eigen's ldlt solver in place of matrix inversion (for speed purposes)
1206  ldlt.compute(covariance);
1207 
1208  for (ii = 0; ii < _nFunctions; ii++) {
1209  fbar = 0.0;
1210  for (i = 0; i < numberOfNeighbors; i++) fbar += _function[neighbors[i]][ii];
1211  fbar = fbar / double(numberOfNeighbors);
1212 
1213  for (i = 0; i < numberOfNeighbors; i++) bb(i, 0) = _function[neighbors[i]][ii] - fbar;
1214  xx = ldlt.solve(bb);
1215 
1216  mu[ii] = fbar;
1217 
1218  for (i = 0; i < numberOfNeighbors; i++) {
1219  mu[ii] += covarianceTestPoint[i] * xx(i, 0);
1220  }
1221 
1222  } // ii = 0 through _nFunctions
1223 
1224  _timer.addToEigen();
1225 
1226  variance[0] = (*_covariogram)(vv, vv) + _lambda;
1227 
1228  for (i = 0; i < numberOfNeighbors; i++) bb(i) = covarianceTestPoint[i];
1229 
1230  xx = ldlt.solve(bb);
1231 
1232  for (i = 0; i < numberOfNeighbors; i++) {
1233  variance[0] -= covarianceTestPoint[i] * xx(i, 0);
1234  }
1235  variance[0] = variance[0] * _krigingParameter;
1236 
1237  for (i = 1; i < _nFunctions; i++) variance[i] = variance[0];
1238 
1239  _timer.addToVariance();
1240  _timer.addToTotal(1);
1241 }
void addToSearch()
Adds time to _searchTime.
MatrixQ covariance
Definition: simpleShape.cc:152

◆ 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]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 1009 of file GaussianProcess.cc.

1010  {
1011  if (_nFunctions > 1) {
1013  "You need to call the version of GaussianProcess.interpolate() "
1014  "that accepts mu and variance arrays (which it populates with results). "
1015  "You are interpolating more than one function.");
1016  }
1017 
1018  if (numberOfNeighbors <= 0) {
1020  "Asked for zero or negative number of neighbors\n");
1021  }
1022 
1023  if (numberOfNeighbors > _kdTree.getNPoints()) {
1025  "Asked for more neighbors than you have data points\n");
1026  }
1027 
1028  if (variance.getNumElements() != static_cast<ndarray::Size>(_nFunctions)) {
1030  "Your variance array is the incorrect size for the number "
1031  "of functions you are trying to interpolate\n");
1032  }
1033 
1034  if (vin.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
1036  "You are interpolating at a point with different dimensionality than you data\n");
1037  }
1038 
1039  int i, j;
1040  T fbar, mu;
1041 
1042  ndarray::Array<T, 1, 1> covarianceTestPoint;
1043  ndarray::Array<int, 1, 1> neighbors;
1044  ndarray::Array<double, 1, 1> neighborDistances, vv;
1045 
1046  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> covariance, bb, xx;
1047  Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
1048 
1049  _timer.start();
1050 
1051  bb.resize(numberOfNeighbors, 1);
1052  xx.resize(numberOfNeighbors, 1);
1053 
1054  covariance.resize(numberOfNeighbors, numberOfNeighbors);
1055 
1056  covarianceTestPoint = allocate(ndarray::makeVector(numberOfNeighbors));
1057 
1058  neighbors = allocate(ndarray::makeVector(numberOfNeighbors));
1059 
1060  neighborDistances = allocate(ndarray::makeVector(numberOfNeighbors));
1061 
1062  vv = allocate(ndarray::makeVector(_dimensions));
1063  if (_useMaxMin == 1) {
1064  // if you constructed this Gaussian process with minimum and maximum
1065  // values for the dimensions of your parameter space,
1066  // the point you are interpolating must be scaled to match the data so
1067  // that the selected nearest neighbors are appropriate
1068 
1069  for (i = 0; i < _dimensions; i++) vv[i] = (vin[i] - _min[i]) / (_max[i] - _min[i]);
1070  } else {
1071  vv = vin;
1072  }
1073 
1074  _kdTree.findNeighbors(neighbors, neighborDistances, vv, numberOfNeighbors);
1075 
1076  _timer.addToSearch();
1077 
1078  fbar = 0.0;
1079  for (i = 0; i < numberOfNeighbors; i++) fbar += _function[neighbors[i]][0];
1080  fbar = fbar / double(numberOfNeighbors);
1081 
1082  for (i = 0; i < numberOfNeighbors; i++) {
1083  covarianceTestPoint[i] = (*_covariogram)(vv, _kdTree.getData(neighbors[i]));
1084 
1085  covariance(i, i) =
1086  (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[i])) + _lambda;
1087 
1088  for (j = i + 1; j < numberOfNeighbors; j++) {
1089  covariance(i, j) = (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[j]));
1090  covariance(j, i) = covariance(i, j);
1091  }
1092  }
1093 
1094  _timer.addToIteration();
1095 
1096  // use Eigen's ldlt solver in place of matrix inversion (for speed purposes)
1097  ldlt.compute(covariance);
1098 
1099  for (i = 0; i < numberOfNeighbors; i++) bb(i, 0) = _function[neighbors[i]][0] - fbar;
1100 
1101  xx = ldlt.solve(bb);
1102  _timer.addToEigen();
1103 
1104  mu = fbar;
1105 
1106  for (i = 0; i < numberOfNeighbors; i++) {
1107  mu += covarianceTestPoint[i] * xx(i, 0);
1108  }
1109 
1110  _timer.addToIteration();
1111 
1112  variance(0) = (*_covariogram)(vv, vv) + _lambda;
1113 
1114  for (i = 0; i < numberOfNeighbors; i++) bb(i) = covarianceTestPoint[i];
1115 
1116  xx = ldlt.solve(bb);
1117 
1118  for (i = 0; i < numberOfNeighbors; i++) {
1119  variance(0) -= covarianceTestPoint[i] * xx(i, 0);
1120  }
1121 
1122  variance(0) = variance(0) * _krigingParameter;
1123 
1124  _timer.addToVariance();
1125  _timer.addToTotal(1);
1126 
1127  return mu;
1128 }

◆ 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]dexthe index of the point you want to remove from your data set
Exceptions
pex::exceptions::RuntimeErrorif 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.

1952  {
1953  int i, j;
1954 
1955  _kdTree.removePoint(dex);
1956 
1957  for (i = dex; i < _npts; i++) {
1958  for (j = 0; j < _nFunctions; j++) {
1959  _function[i][j] = _function[i + 1][j];
1960  }
1961  }
1962  _npts = _kdTree.getNPoints();
1963 }

◆ 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]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::exceptions::RuntimeErrorif the nearest neighbor search does not find the data point itself as the nearest neighbor

Definition at line 1367 of file GaussianProcess.cc.

1368  {
1369  if (mu.getNumElements() != static_cast<ndarray::Size>(_nFunctions) ||
1370  variance.getNumElements() != static_cast<ndarray::Size>(_nFunctions)) {
1372  "Your mu and/or var arrays are improperly sized for the number of functions "
1373  "you are interpolating\n");
1374  }
1375 
1376  if (numberOfNeighbors <= 0) {
1378  "Asked for zero or negative number of neighbors\n");
1379  }
1380 
1381  if (numberOfNeighbors + 1 > _kdTree.getNPoints()) {
1383  "Asked for more neighbors than you have data points\n");
1384  }
1385 
1386  if (dex < 0 || dex >= _kdTree.getNPoints()) {
1388  "Asked to self interpolate on a point that does not exist\n");
1389  }
1390 
1391  int i, j, ii;
1392  T fbar;
1393 
1394  ndarray::Array<T, 1, 1> covarianceTestPoint;
1395  ndarray::Array<int, 1, 1> selfNeighbors;
1396  ndarray::Array<double, 1, 1> selfDistances;
1397  ndarray::Array<int, 1, 1> neighbors;
1398  ndarray::Array<double, 1, 1> neighborDistances;
1399 
1400  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> covariance, bb, xx;
1401  Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
1402 
1403  _timer.start();
1404 
1405  bb.resize(numberOfNeighbors, 1);
1406  xx.resize(numberOfNeighbors, 1);
1407  covariance.resize(numberOfNeighbors, numberOfNeighbors);
1408  covarianceTestPoint = allocate(ndarray::makeVector(numberOfNeighbors));
1409  neighbors = allocate(ndarray::makeVector(numberOfNeighbors));
1410  neighborDistances = allocate(ndarray::makeVector(numberOfNeighbors));
1411 
1412  selfNeighbors = allocate(ndarray::makeVector(numberOfNeighbors + 1));
1413  selfDistances = allocate(ndarray::makeVector(numberOfNeighbors + 1));
1414 
1415  // we don't use _useMaxMin because the data has already been normalized
1416 
1417  _kdTree.findNeighbors(selfNeighbors, selfDistances, _kdTree.getData(dex), numberOfNeighbors + 1);
1418 
1419  _timer.addToSearch();
1420 
1421  if (selfNeighbors[0] != dex) {
1423  "Nearest neighbor search in selfInterpolate did not find self\n");
1424  }
1425 
1426  // SelfNeighbors[0] will be the point itself (it is its own nearest neighbor)
1427  // We discard that for the interpolation calculation
1428  //
1429  // If you do not wish to do this, simply call the usual ::interpolate() method instead of
1430  //::selfInterpolate()
1431  for (i = 0; i < numberOfNeighbors; i++) {
1432  neighbors[i] = selfNeighbors[i + 1];
1433  neighborDistances[i] = selfDistances[i + 1];
1434  }
1435 
1436  for (i = 0; i < numberOfNeighbors; i++) {
1437  covarianceTestPoint[i] = (*_covariogram)(_kdTree.getData(dex), _kdTree.getData(neighbors[i]));
1438 
1439  covariance(i, i) =
1440  (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[i])) + _lambda;
1441 
1442  for (j = i + 1; j < numberOfNeighbors; j++) {
1443  covariance(i, j) = (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[j]));
1444  covariance(j, i) = covariance(i, j);
1445  }
1446  }
1447  _timer.addToIteration();
1448 
1449  // use Eigen's ldlt solver in place of matrix inversion (for speed purposes)
1450  ldlt.compute(covariance);
1451 
1452  for (ii = 0; ii < _nFunctions; ii++) {
1453  fbar = 0.0;
1454  for (i = 0; i < numberOfNeighbors; i++) fbar += _function[neighbors[i]][ii];
1455  fbar = fbar / double(numberOfNeighbors);
1456 
1457  for (i = 0; i < numberOfNeighbors; i++) bb(i, 0) = _function[neighbors[i]][ii] - fbar;
1458  xx = ldlt.solve(bb);
1459 
1460  mu[ii] = fbar;
1461 
1462  for (i = 0; i < numberOfNeighbors; i++) {
1463  mu[ii] += covarianceTestPoint[i] * xx(i, 0);
1464  }
1465  } // ii = 0 through _nFunctions
1466 
1467  _timer.addToEigen();
1468 
1469  variance[0] = (*_covariogram)(_kdTree.getData(dex), _kdTree.getData(dex)) + _lambda;
1470 
1471  for (i = 0; i < numberOfNeighbors; i++) bb(i) = covarianceTestPoint[i];
1472 
1473  xx = ldlt.solve(bb);
1474 
1475  for (i = 0; i < numberOfNeighbors; i++) {
1476  variance[0] -= covarianceTestPoint[i] * xx(i, 0);
1477  }
1478 
1479  variance[0] = variance[0] * _krigingParameter;
1480 
1481  for (i = 1; i < _nFunctions; i++) variance[i] = variance[0];
1482 
1483  _timer.addToVariance();
1484  _timer.addToTotal(1);
1485 }

◆ 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]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::exceptions::RuntimeErrorif 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.

1245  {
1246  if (_nFunctions > 1) {
1248  "You need to call the version of GaussianProcess.selfInterpolate() "
1249  "that accepts mu and variance arrays (which it populates with results). "
1250  "You are interpolating more than one function.");
1251  }
1252 
1253  if (variance.getNumElements() != static_cast<ndarray::Size>(_nFunctions)) {
1255  "Your variance array is the incorrect size for the number "
1256  "of functions you are trying to interpolate\n");
1257  }
1258 
1259  if (numberOfNeighbors <= 0) {
1261  "Asked for zero or negative number of neighbors\n");
1262  }
1263 
1264  if (numberOfNeighbors + 1 > _kdTree.getNPoints()) {
1266  "Asked for more neighbors than you have data points\n");
1267  }
1268 
1269  if (dex < 0 || dex >= _kdTree.getNPoints()) {
1271  "Asked to self interpolate on a point that does not exist\n");
1272  }
1273 
1274  int i, j;
1275  T fbar, mu;
1276 
1277  ndarray::Array<T, 1, 1> covarianceTestPoint;
1278  ndarray::Array<int, 1, 1> selfNeighbors;
1279  ndarray::Array<double, 1, 1> selfDistances;
1280  ndarray::Array<int, 1, 1> neighbors;
1281  ndarray::Array<double, 1, 1> neighborDistances;
1282 
1283  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> covariance, bb, xx;
1284  Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > ldlt;
1285 
1286  _timer.start();
1287 
1288  bb.resize(numberOfNeighbors, 1);
1289  xx.resize(numberOfNeighbors, 1);
1290  covariance.resize(numberOfNeighbors, numberOfNeighbors);
1291  covarianceTestPoint = allocate(ndarray::makeVector(numberOfNeighbors));
1292  neighbors = allocate(ndarray::makeVector(numberOfNeighbors));
1293  neighborDistances = allocate(ndarray::makeVector(numberOfNeighbors));
1294 
1295  selfNeighbors = allocate(ndarray::makeVector(numberOfNeighbors + 1));
1296  selfDistances = allocate(ndarray::makeVector(numberOfNeighbors + 1));
1297 
1298  // we don't use _useMaxMin because the data has already been normalized
1299 
1300  _kdTree.findNeighbors(selfNeighbors, selfDistances, _kdTree.getData(dex), numberOfNeighbors + 1);
1301 
1302  _timer.addToSearch();
1303 
1304  if (selfNeighbors[0] != dex) {
1306  "Nearest neighbor search in selfInterpolate did not find self\n");
1307  }
1308 
1309  // SelfNeighbors[0] will be the point itself (it is its own nearest neighbor)
1310  // We discard that for the interpolation calculation
1311  //
1312  // If you do not wish to do this, simply call the usual ::interpolate() method instead of
1313  //::selfInterpolate()
1314  for (i = 0; i < numberOfNeighbors; i++) {
1315  neighbors[i] = selfNeighbors[i + 1];
1316  neighborDistances[i] = selfDistances[i + 1];
1317  }
1318 
1319  fbar = 0.0;
1320  for (i = 0; i < numberOfNeighbors; i++) fbar += _function[neighbors[i]][0];
1321  fbar = fbar / double(numberOfNeighbors);
1322 
1323  for (i = 0; i < numberOfNeighbors; i++) {
1324  covarianceTestPoint[i] = (*_covariogram)(_kdTree.getData(dex), _kdTree.getData(neighbors[i]));
1325 
1326  covariance(i, i) =
1327  (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[i])) + _lambda;
1328 
1329  for (j = i + 1; j < numberOfNeighbors; j++) {
1330  covariance(i, j) = (*_covariogram)(_kdTree.getData(neighbors[i]), _kdTree.getData(neighbors[j]));
1331  covariance(j, i) = covariance(i, j);
1332  }
1333  }
1334  _timer.addToIteration();
1335 
1336  // use Eigen's ldlt solver in place of matrix inversion (for speed purposes)
1337  ldlt.compute(covariance);
1338 
1339  for (i = 0; i < numberOfNeighbors; i++) bb(i, 0) = _function[neighbors[i]][0] - fbar;
1340  xx = ldlt.solve(bb);
1341  _timer.addToEigen();
1342 
1343  mu = fbar;
1344 
1345  for (i = 0; i < numberOfNeighbors; i++) {
1346  mu += covarianceTestPoint[i] * xx(i, 0);
1347  }
1348 
1349  variance(0) = (*_covariogram)(_kdTree.getData(dex), _kdTree.getData(dex)) + _lambda;
1350 
1351  for (i = 0; i < numberOfNeighbors; i++) bb(i) = covarianceTestPoint[i];
1352 
1353  xx = ldlt.solve(bb);
1354 
1355  for (i = 0; i < numberOfNeighbors; i++) {
1356  variance(0) -= covarianceTestPoint[i] * xx(i, 0);
1357  }
1358 
1359  variance(0) = variance(0) * _krigingParameter;
1360  _timer.addToVariance();
1361  _timer.addToTotal(1);
1362 
1363  return mu;
1364 }

◆ 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]covarthe Covariogram object that you wish to assign

Definition at line 1971 of file GaussianProcess.cc.

1971  {
1972  _covariogram = covar;
1973 }

◆ setKrigingParameter()

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 1966 of file GaussianProcess.cc.

1966  {
1967  _krigingParameter = kk;
1968 }

◆ setLambda()

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 1976 of file GaussianProcess.cc.

1976  {
1977  _lambda = lambda;
1978 }

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