GaussianProcess.h

Go to the documentation of this file.

// -*- LSST-C++ -*-
 
/*
 * LSST Data Management System
 * Copyright 2008, 2009, 2010 LSST Corporation.
 *
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the LSST License Statement and
 * the GNU General Public License along with this program.  If not,
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */
 
#ifndef LSST_AFW_MATH_GAUSSIAN_PROCESS_H
#define LSST_AFW_MATH_GAUSSIAN_PROCESS_H
 
#include <Eigen/Dense>
#include <stdexcept>
 
#include "ndarray/eigen.h"
#include <memory>
 
#include "lsst/daf/base/DateTime.h"
#include "lsst/pex/exceptions.h"
 
namespace lsst {
namespace afw {
namespace math {
 
class GaussianProcessTimer {
public:
    GaussianProcessTimer();
 
    GaussianProcessTimer(GaussianProcessTimer const &) = default;
    GaussianProcessTimer(GaussianProcessTimer &&) = default;
    GaussianProcessTimer &operator=(GaussianProcessTimer const &) = default;
    GaussianProcessTimer &operator=(GaussianProcessTimer &&) = default;
    ~GaussianProcessTimer() = default;
 
    void reset();
 
    void start();
 
    void addToEigen();
 
    void addToVariance();
 
    void addToSearch();
 
    void addToIteration();
 
    void addToTotal(int i);
 
    void display();
 
private:
    double _before, _beginning;
    double _eigenTime, _iterationTime, _searchTime, _varianceTime, _totalTime;
    int _interpolationCount;
};
 
template <typename T>
class Covariogram {
public:
    virtual ~Covariogram();
 
    // No copying
    Covariogram(const Covariogram &) = delete;
    Covariogram &operator=(const Covariogram &) = delete;
 
    // No moving
    Covariogram(Covariogram &&) = delete;
    Covariogram &operator=(Covariogram &&) = delete;
 
    explicit Covariogram(){};
 
    virtual T operator()(ndarray::Array<const T, 1, 1> const &p1,
                         ndarray::Array<const T, 1, 1> const &p2) const;
};
 
template <typename T>
class SquaredExpCovariogram : public Covariogram<T> {
public:
    ~SquaredExpCovariogram() override;
 
    explicit SquaredExpCovariogram();
 
    void setEllSquared(double ellSquared);
 
    T operator()(ndarray::Array<const T, 1, 1> const &, ndarray::Array<const T, 1, 1> const &) const override;
 
private:
    double _ellSquared;
};
 
template <typename T>
class NeuralNetCovariogram : public Covariogram<T> {
public:
    ~NeuralNetCovariogram() override;
 
    explicit NeuralNetCovariogram();
 
    void setSigma0(double sigma0);
 
    void setSigma1(double sigma1);
 
    T operator()(ndarray::Array<const T, 1, 1> const &, ndarray::Array<const T, 1, 1> const &) const override;
 
private:
    double _sigma0, _sigma1;
};
 
template <typename T>
class KdTree {
public:
    // No copying
    KdTree(const KdTree &) = delete;
    KdTree &operator=(const KdTree &) = delete;
 
    // No moving
    KdTree(KdTree &&) = delete;
    KdTree &operator=(KdTree &&) = delete;
 
    // Add default constructor (which is no longer generated)
    KdTree() = default;
 
    void Initialize(ndarray::Array<T, 2, 2> const &dt);
 
    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;
 
    T getData(int ipt, int idim) const;
 
    ndarray::Array<T, 1, 1> getData(int ipt) const;
 
    void addPoint(ndarray::Array<const T, 1, 1> const &v);
 
    void removePoint(int dex);
 
    int getNPoints() const;
 
    void getTreeNode(ndarray::Array<int, 1, 1> const &v, int dex) const;
 
private:
    ndarray::Array<int, 2, 2> _tree;
    ndarray::Array<int, 1, 1> _inn;
    ndarray::Array<T, 2, 2> _data;
 
    enum { DIMENSION, LT, GEQ, PARENT };
 
    //_tree stores the relationships between data points
    //_tree[i][DIMENSION] is the dimension on which the ith node segregates its daughters
    //
    //_tree[i][LT] is the branch of the tree down which the daughters' DIMENSIONth component
    // is less than the parent's
    //
    //_tree[i][GEQ] is the branch of the tree down which the daughters' DIMENSIONth component is
    // greather than or equal to the parent's
    //
    //_tree[i][PARENT] is the parent node of the ith node
 
    //_data actually stores the data points
 
    int _npts, _dimensions, _room, _roomStep, _masterParent;
    mutable int _neighborsFound, _neighborsWanted;
 
    //_room denotes the capacity of _data and _tree.  It will usually be larger
    // than _npts so that we do not have to reallocate
    //_tree and _data every time we add a new point to the tree
 
    mutable ndarray::Array<double, 1, 1> _neighborDistances;
    mutable ndarray::Array<int, 1, 1> _neighborCandidates;
 
    void _organize(ndarray::Array<int, 1, 1> const &use, int ct, int parent, int dir);
 
    int _findNode(ndarray::Array<const T, 1, 1> const &v) const;
 
    void _lookForNeighbors(ndarray::Array<const T, 1, 1> const &v, int consider, int from) const;
 
    int _testTree() const;
 
    int _walkUpTree(int target, int dir, int root) const;
 
    void _count(int where, int *ct) const;
 
    void _descend(int root);
 
    void _reassign(int target);
 
    double _distance(ndarray::Array<const T, 1, 1> const &p1, ndarray::Array<const T, 1, 1> const &p2) const;
};
 
template <typename T>
class GaussianProcess {
public:
    // No copying
    GaussianProcess(const GaussianProcess &) = delete;
    GaussianProcess &operator=(const GaussianProcess &) = delete;
 
    // No moving
    GaussianProcess(GaussianProcess &&) = delete;
    GaussianProcess &operator=(GaussianProcess &&) = delete;
 
    GaussianProcess(ndarray::Array<T, 2, 2> const &dataIn, ndarray::Array<T, 1, 1> const &ff,
                    std::shared_ptr<Covariogram<T> > const &covarIn);
 
    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);
    GaussianProcess(ndarray::Array<T, 2, 2> const &dataIn, ndarray::Array<T, 2, 2> const &ff,
                    std::shared_ptr<Covariogram<T> > const &covarIn);
    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);
 
    int getNPoints() const;
 
    int getDim() const;
 
    void getData(ndarray::Array<T, 2, 2> pts, ndarray::Array<T, 1, 1> fn,
                 ndarray::Array<int, 1, 1> indices) const;
 
    void getData(ndarray::Array<T, 2, 2> pts, ndarray::Array<T, 2, 2> fn,
                 ndarray::Array<int, 1, 1> indices) const;
 
    T interpolate(ndarray::Array<T, 1, 1> variance, ndarray::Array<T, 1, 1> const &vin,
                  int numberOfNeighbors) const;
    void interpolate(ndarray::Array<T, 1, 1> mu, ndarray::Array<T, 1, 1> variance,
                     ndarray::Array<T, 1, 1> const &vin, int numberOfNeighbors) const;
 
    T selfInterpolate(ndarray::Array<T, 1, 1> variance, int dex, int numberOfNeighbors) const;
 
    void selfInterpolate(ndarray::Array<T, 1, 1> mu, ndarray::Array<T, 1, 1> variance, int dex,
                         int numberOfNeighbors) const;
 
    void batchInterpolate(ndarray::Array<T, 1, 1> mu, ndarray::Array<T, 1, 1> variance,
                          ndarray::Array<T, 2, 2> const &queries) const;
 
    void batchInterpolate(ndarray::Array<T, 1, 1> mu, ndarray::Array<T, 2, 2> const &queries) const;
 
    void batchInterpolate(ndarray::Array<T, 2, 2> mu, ndarray::Array<T, 2, 2> variance,
                          ndarray::Array<T, 2, 2> const &queries) const;
 
    void batchInterpolate(ndarray::Array<T, 2, 2> mu, ndarray::Array<T, 2, 2> const &queries) const;
 
    void addPoint(ndarray::Array<T, 1, 1> const &vin, T f);
 
    void addPoint(ndarray::Array<T, 1, 1> const &vin, ndarray::Array<T, 1, 1> const &f);
 
    void removePoint(int dex);
 
    void setKrigingParameter(T kk);
 
    void setCovariogram(std::shared_ptr<Covariogram<T> > const &covar);
 
    void setLambda(T lambda);
 
    GaussianProcessTimer &getTimes() const;
 
private:
    int _npts, _useMaxMin, _dimensions, _room, _roomStep, _nFunctions;
 
    T _krigingParameter, _lambda;
 
    ndarray::Array<T, 1, 1> _max, _min;
    ndarray::Array<T, 2, 2> _function;
 
    KdTree<T> _kdTree;
 
    std::shared_ptr<Covariogram<T> > _covariogram;
    mutable GaussianProcessTimer _timer;
};
}  // namespace math
}  // namespace afw
}  // namespace lsst
 
#endif  //#ifndef LSST_AFW_MATH_GAUSSIAN_PROCESS_H