LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
GaussianProcess.h
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 #ifndef LSST_AFW_MATH_GAUSSIAN_PROCESS_H
26 #define LSST_AFW_MATH_GAUSSIAN_PROCESS_H
27 
28 #include <Eigen/Dense>
29 #include <stdexcept>
30 
31 #include "ndarray/eigen.h"
32 #include <memory>
33 
34 #include "lsst/daf/base/DateTime.h"
35 #include "lsst/pex/exceptions.h"
36 
37 namespace lsst {
38 namespace afw {
39 namespace math {
40 
61 public:
63 
68  ~GaussianProcessTimer() = default;
69 
74  void reset();
75 
79  void start();
80 
84  void addToEigen();
85 
89  void addToVariance();
90 
94  void addToSearch();
95 
99  void addToIteration();
100 
106  void addToTotal(int i);
107 
111  void display();
112 
113 private:
114  double _before, _beginning;
115  double _eigenTime, _iterationTime, _searchTime, _varianceTime, _totalTime;
116  int _interpolationCount;
117 };
118 
125 template <typename T>
126 class Covariogram {
127 public:
128  virtual ~Covariogram();
129 
130  // No copying
131  Covariogram(const Covariogram &) = delete;
132  Covariogram &operator=(const Covariogram &) = delete;
133 
134  // No moving
135  Covariogram(Covariogram &&) = delete;
137 
141  explicit Covariogram(){};
142 
150  virtual T operator()(ndarray::Array<const T, 1, 1> const &p1,
151  ndarray::Array<const T, 1, 1> const &p2) const;
152 };
153 
163 template <typename T>
165 public:
167 
168  explicit SquaredExpCovariogram();
169 
174  void setEllSquared(double ellSquared);
175 
176  T operator()(ndarray::Array<const T, 1, 1> const &, ndarray::Array<const T, 1, 1> const &) const override;
177 
178 private:
179  double _ellSquared;
180 };
181 
192 template <typename T>
194 public:
196 
197  explicit NeuralNetCovariogram();
198 
202  void setSigma0(double sigma0);
203 
207  void setSigma1(double sigma1);
208 
209  T operator()(ndarray::Array<const T, 1, 1> const &, ndarray::Array<const T, 1, 1> const &) const override;
210 
211 private:
212  double _sigma0, _sigma1;
213 };
214 
223 template <typename T>
224 class KdTree {
225 public:
226  // No copying
227  KdTree(const KdTree &) = delete;
228  KdTree &operator=(const KdTree &) = delete;
229 
230  // No moving
231  KdTree(KdTree &&) = delete;
232  KdTree &operator=(KdTree &&) = delete;
233 
234  // Add default constructor (which is no longer generated)
235  KdTree() = default;
236 
245  void Initialize(ndarray::Array<T, 2, 2> const &dt);
246 
262  void findNeighbors(ndarray::Array<int, 1, 1> neighdex, ndarray::Array<double, 1, 1> dd,
263  ndarray::Array<const T, 1, 1> const &v, int n_nn) const;
264 
272  T getData(int ipt, int idim) const;
273 
295  ndarray::Array<T, 1, 1> getData(int ipt) const;
296 
304  void addPoint(ndarray::Array<const T, 1, 1> const &v);
305 
314  void removePoint(int dex);
315 
319  int getNPoints() const;
320 
328  void getTreeNode(ndarray::Array<int, 1, 1> const &v, int dex) const;
329 
330 private:
331  ndarray::Array<int, 2, 2> _tree;
332  ndarray::Array<int, 1, 1> _inn;
333  ndarray::Array<T, 2, 2> _data;
334 
335  enum { DIMENSION, LT, GEQ, PARENT };
336 
337  //_tree stores the relationships between data points
338  //_tree[i][DIMENSION] is the dimension on which the ith node segregates its daughters
339  //
340  //_tree[i][LT] is the branch of the tree down which the daughters' DIMENSIONth component
341  // is less than the parent's
342  //
343  //_tree[i][GEQ] is the branch of the tree down which the daughters' DIMENSIONth component is
344  // greather than or equal to the parent's
345  //
346  //_tree[i][PARENT] is the parent node of the ith node
347 
348  //_data actually stores the data points
349 
350  int _npts, _dimensions, _room, _roomStep, _masterParent;
351  mutable int _neighborsFound, _neighborsWanted;
352 
353  //_room denotes the capacity of _data and _tree. It will usually be larger
354  // than _npts so that we do not have to reallocate
355  //_tree and _data every time we add a new point to the tree
356 
357  mutable ndarray::Array<double, 1, 1> _neighborDistances;
358  mutable ndarray::Array<int, 1, 1> _neighborCandidates;
359 
372  void _organize(ndarray::Array<int, 1, 1> const &use, int ct, int parent, int dir);
373 
379  int _findNode(ndarray::Array<const T, 1, 1> const &v) const;
380 
397  void _lookForNeighbors(ndarray::Array<const T, 1, 1> const &v, int consider, int from) const;
398 
402  int _testTree() const;
403 
417  int _walkUpTree(int target, int dir, int root) const;
418 
427  void _count(int where, int *ct) const;
428 
434  void _descend(int root);
435 
441  void _reassign(int target);
442 
446  double _distance(ndarray::Array<const T, 1, 1> const &p1, ndarray::Array<const T, 1, 1> const &p2) const;
447 };
448 
470 template <typename T>
472 public:
473  // No copying
474  GaussianProcess(const GaussianProcess &) = delete;
476 
477  // No moving
480 
493  GaussianProcess(ndarray::Array<T, 2, 2> const &dataIn, ndarray::Array<T, 1, 1> const &ff,
494  std::shared_ptr<Covariogram<T> > const &covarIn);
495 
515  GaussianProcess(ndarray::Array<T, 2, 2> const &dataIn, ndarray::Array<T, 1, 1> const &mn,
516  ndarray::Array<T, 1, 1> const &mx, ndarray::Array<T, 1, 1> const &ff,
517  std::shared_ptr<Covariogram<T> > const &covarIn);
530  GaussianProcess(ndarray::Array<T, 2, 2> const &dataIn, ndarray::Array<T, 2, 2> const &ff,
531  std::shared_ptr<Covariogram<T> > const &covarIn);
548  GaussianProcess(ndarray::Array<T, 2, 2> const &dataIn, ndarray::Array<T, 1, 1> const &mn,
549  ndarray::Array<T, 1, 1> const &mx, ndarray::Array<T, 2, 2> const &ff,
550  std::shared_ptr<Covariogram<T> > const &covarIn);
551 
555  int getNPoints() const;
556 
561  int getDim() const;
562 
572  void getData(ndarray::Array<T, 2, 2> pts, ndarray::Array<T, 1, 1> fn,
573  ndarray::Array<int, 1, 1> indices) const;
574 
584  void getData(ndarray::Array<T, 2, 2> pts, ndarray::Array<T, 2, 2> fn,
585  ndarray::Array<int, 1, 1> indices) const;
586 
604  T interpolate(ndarray::Array<T, 1, 1> variance, ndarray::Array<T, 1, 1> const &vin,
605  int numberOfNeighbors) const;
621  void interpolate(ndarray::Array<T, 1, 1> mu, ndarray::Array<T, 1, 1> variance,
622  ndarray::Array<T, 1, 1> const &vin, int numberOfNeighbors) const;
623 
643  T selfInterpolate(ndarray::Array<T, 1, 1> variance, int dex, int numberOfNeighbors) const;
644 
659  void selfInterpolate(ndarray::Array<T, 1, 1> mu, ndarray::Array<T, 1, 1> variance, int dex,
660  int numberOfNeighbors) const;
661 
683  void batchInterpolate(ndarray::Array<T, 1, 1> mu, ndarray::Array<T, 1, 1> variance,
684  ndarray::Array<T, 2, 2> const &queries) const;
685 
703  void batchInterpolate(ndarray::Array<T, 1, 1> mu, ndarray::Array<T, 2, 2> const &queries) const;
704 
709  void batchInterpolate(ndarray::Array<T, 2, 2> mu, ndarray::Array<T, 2, 2> variance,
710  ndarray::Array<T, 2, 2> const &queries) const;
711 
716  void batchInterpolate(ndarray::Array<T, 2, 2> mu, ndarray::Array<T, 2, 2> const &queries) const;
717 
734  void addPoint(ndarray::Array<T, 1, 1> const &vin, T f);
735 
745  void addPoint(ndarray::Array<T, 1, 1> const &vin, ndarray::Array<T, 1, 1> const &f);
746 
758  void removePoint(int dex);
759 
766  void setKrigingParameter(T kk);
767 
774  void setCovariogram(std::shared_ptr<Covariogram<T> > const &covar);
775 
787  void setLambda(T lambda);
788 
804 
805 private:
806  int _npts, _useMaxMin, _dimensions, _room, _roomStep, _nFunctions;
807 
808  T _krigingParameter, _lambda;
809 
810  ndarray::Array<T, 1, 1> _max, _min;
811  ndarray::Array<T, 2, 2> _function;
812 
813  KdTree<T> _kdTree;
814 
815  std::shared_ptr<Covariogram<T> > _covariogram;
816  mutable GaussianProcessTimer _timer;
817 };
818 } // namespace math
819 } // namespace afw
820 } // namespace lsst
821 
822 #endif //#ifndef LSST_AFW_MATH_GAUSSIAN_PROCESS_H
Key< Flag > const & target
Interface for DateTime class.
table::Key< double > sigma1
afw::table::Key< afw::table::Array< VariancePixelT > > variance
table::Key< int > from
The parent class of covariogram functions for use in Gaussian Process interpolation.
Covariogram()
construct a Covariogram assigning default values to the hyper parameters
Covariogram & operator=(Covariogram &&)=delete
Covariogram(const Covariogram &)=delete
virtual T operator()(ndarray::Array< const T, 1, 1 > const &p1, ndarray::Array< const T, 1, 1 > const &p2) const
Actually evaluate the covariogram function relating two points you want to interpolate from.
Covariogram(Covariogram &&)=delete
Covariogram & operator=(const Covariogram &)=delete
Stores values of a function sampled on an image and allows you to interpolate the function to unsampl...
GaussianProcessTimer & getTimes() const
Give the user acces to _timer, an object keeping track of the time spent on various processes within ...
void setLambda(T lambda)
set the value of the hyperparameter _lambda
int getDim() const
return the dimensionality of data points stored in the GaussianProcess
GaussianProcess & operator=(GaussianProcess &&)=delete
void setKrigingParameter(T kk)
Assign a value to the Kriging paramter.
void addPoint(ndarray::Array< T, 1, 1 > const &vin, T f)
Add a point to the pool of data used by GaussianProcess for interpolation.
void 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)
GaussianProcess(const GaussianProcess &)=delete
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.
T selfInterpolate(ndarray::Array< T, 1, 1 > variance, int dex, int numberOfNeighbors) const
This method will interpolate the function on a data point for purposes of optimizing hyper parameters...
void removePoint(int dex)
This will remove a point from the data set.
void setCovariogram(std::shared_ptr< Covariogram< T > > const &covar)
Assign a different covariogram to this GaussianProcess.
void getData(ndarray::Array< T, 2, 2 > pts, ndarray::Array< T, 1, 1 > fn, ndarray::Array< int, 1, 1 > indices) const
Return a sub-sample the data underlying the Gaussian Process.
GaussianProcess & operator=(const GaussianProcess &)=delete
GaussianProcess(GaussianProcess &&)=delete
int getNPoints() const
return the number of data points stored in the GaussianProcess
This is a structure for keeping track of how long the interpolation methods spend on different parts ...
void addToSearch()
Adds time to _searchTime.
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 display()
Displays the current values of all times and _interpolationCount.
GaussianProcessTimer(GaussianProcessTimer const &)=default
GaussianProcessTimer & operator=(GaussianProcessTimer &&)=default
void reset()
Resets all of the data members of GaussianProcessTimer to zero.
void addToEigen()
Adds time to _eigenTime.
GaussianProcessTimer & operator=(GaussianProcessTimer const &)=default
void addToTotal(int i)
Adds time to _totalTime and adds counts to _interpolationCount.
GaussianProcessTimer(GaussianProcessTimer &&)=default
The data for GaussianProcess is stored in a KD tree to facilitate nearest-neighbor searches.
KdTree(const KdTree &)=delete
KdTree & operator=(const KdTree &)=delete
void removePoint(int dex)
Remove a point from the tree.
void findNeighbors(ndarray::Array< int, 1, 1 > neighdex, ndarray::Array< double, 1, 1 > dd, ndarray::Array< const T, 1, 1 > const &v, int n_nn) const
Find the nearest neighbors of a point.
void addPoint(ndarray::Array< const T, 1, 1 > const &v)
Add a point to the tree.
void Initialize(ndarray::Array< T, 2, 2 > const &dt)
Build a KD Tree to store the data for GaussianProcess.
void getTreeNode(ndarray::Array< int, 1, 1 > const &v, int dex) const
Return the _tree information for a given data point.
int getNPoints() const
return the number of data points stored in the tree
KdTree & operator=(KdTree &&)=delete
T getData(int ipt, int idim) const
Return one element of one node on the tree.
KdTree(KdTree &&)=delete
a Covariogram that recreates a neural network with one hidden layer and infinite units in that layer
void setSigma0(double sigma0)
set the _sigma0 hyper parameter
T operator()(ndarray::Array< const T, 1, 1 > const &, ndarray::Array< const T, 1, 1 > const &) const override
Actually evaluate the covariogram function relating two points you want to interpolate from.
void setSigma1(double sigma1)
set the _sigma1 hyper parameter
void setEllSquared(double ellSquared)
set the _ellSquared hyper parameter (the square of the characteristic length scale of the covariogram...
T operator()(ndarray::Array< const T, 1, 1 > const &, ndarray::Array< const T, 1, 1 > const &) const override
Actually evaluate the covariogram function relating two points you want to interpolate from.
A base class for image defects.