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"
#include "lsst/pex/policy.h"
#include "lsst/pex/policy/Policy.h"

namespace lsst {
namespace afw {
namespace math {

class GaussianProcessTimer{

public:

    GaussianProcessTimer();

    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();

#ifndef SWIG
    // No copying
    Covariogram (const Covariogram&) = delete;
    Covariogram& operator=(const Covariogram&) = delete;

    // No moving
    Covariogram (Covariogram&&) = delete;
    Covariogram& operator=(Covariogram&&) = delete;
#endif

    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:
    virtual ~SquaredExpCovariogram();

    explicit SquaredExpCovariogram();

    void setEllSquared(double ellSquared);

    virtual T operator() (ndarray::Array<const T,1,1> const &,
                          ndarray::Array<const T,1,1> const &
                          ) const;

private:
    double _ellSquared;

};

template <typename T>
class NeuralNetCovariogram : public Covariogram<T>{

public:
    virtual ~NeuralNetCovariogram();

    explicit NeuralNetCovariogram();

    void setSigma0(double sigma0);

    void setSigma1(double sigma1);

    virtual T operator() (ndarray::Array<const T,1,1> const &,
                          ndarray::Array<const T,1,1> const &
                          ) const;

private:
    double _sigma0,_sigma1;


};

template <typename T>
class KdTree {
public:

#ifndef SWIG
    // 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;
#endif

    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:

#ifndef SWIG
    // No copying
    GaussianProcess (const GaussianProcess&) = delete;
    GaussianProcess& operator=(const GaussianProcess&) = delete;

    // No moving
    GaussianProcess (GaussianProcess&&) = delete;
    GaussianProcess& operator=(GaussianProcess&&) = delete;
#endif

    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;

};

}}}



#endif //#ifndef LSST_AFW_MATH_GAUSSIAN_PROCESS_H