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/Citizen.h"
#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 lsst::daf::base::Citizen {
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() : lsst::daf::base::Citizen(typeid(this)){};

    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