drwnGaussian.h

/******************************************************************************
** DARWIN: A FRAMEWORK FOR MACHINE LEARNING RESEARCH AND DEVELOPMENT
** Distributed under the terms of the BSD license (see the LICENSE file)
** Copyright (c) 2007-2015, Stephen Gould
** All rights reserved.
**
******************************************************************************
** FILENAME:    drwnGaussian.h
** AUTHOR(S):   Stephen Gould <stephen.gould@anu.edu.au>
**
*****************************************************************************/

#pragma once

#include <vector>
#include <limits>

#include "Eigen/Core"
#include "Eigen/Cholesky"
#include "Eigen/LU"

#include "drwnBase.h"
#include "drwnSuffStats.h"

using namespace std;
using namespace Eigen;

class drwnConditionalGaussian;

// drwnGaussian class --------------------------------------------------------

class drwnGaussian : public drwnStdObjIface {
 public:
    static bool AUTO_RIDGE;

 private:
    int _n;                   
    VectorXd _mu;             
    MatrixXd _mSigma;         

    mutable MatrixXd *_invSigma; 
    mutable double _logZ;        
    mutable MatrixXd *_mL;       

 public:
    drwnGaussian(int n = 1);
    drwnGaussian(const VectorXd& mu, double sigma2);
    drwnGaussian(const VectorXd& mu, const MatrixXd& sigma2);
    drwnGaussian(const vector<double>& mu, double sigma2);
    drwnGaussian(const drwnSuffStats& stats);
    drwnGaussian(const drwnGaussian& model);
    ~drwnGaussian();

    // initialization
    void initialize(int n);
    void initialize(const VectorXd& mu, double sigma2);
    void initialize(const VectorXd& mu, const MatrixXd& sigma2);

    // marginalization
    drwnGaussian *marginalize(const vector<int> & indx) const;

    // conditioning
    drwnGaussian reduce(const vector<double>& x, const vector<int>& indx) const;
    drwnGaussian reduce(const map<int, double>& x) const;
    drwnConditionalGaussian conditionOn(const vector<int>& indx) const;

    // i/o
    const char *type() const { return "drwnGaussian"; }
    drwnGaussian *clone() const { return new drwnGaussian(*this); }
    bool save(drwnXMLNode& xml) const;
    bool load(drwnXMLNode& xml);

    // evaluate (log-likelihood)
    void evaluate(const MatrixXd& x, VectorXd& p) const;
    void evaluate(const vector<vector<double> >& x, vector<double>& p) const;
    double evaluateSingle(const VectorXd& x) const;
    double evaluateSingle(const vector<double>& x) const;
    double evaluateSingle(double x) const;

    // sampling
    void sample(VectorXd& x) const;
    void sample(vector<double>& x) const;

    // learn parameters
    void train(const MatrixXd& x, double lambda = 0.0);
    void train(const vector<vector<double> >& x, double lambda = 0.0);
    void train(const vector<double> &x, double lambda = 0.0);
    void train(const drwnSuffStats& stats, double lambda = 0.0);

    // access
    unsigned dimension() const { return _n; }
    const VectorXd& mean() const { return _mu; }
    const MatrixXd& covariance() const { return _mSigma; }
    double logPartitionFunction() const;

    double klDivergence(const drwnGaussian& model) const;
    double klDivergence(const drwnSuffStats& stats) const;

    // standard operators
    drwnGaussian& operator=(const drwnGaussian& model);

 protected:
    inline void guaranteeInvSigma() const;
    void freeCachedParameters();
    void updateCachedParameters();
};

// drwnConditionalGaussian class ---------------------------------------------

class drwnConditionalGaussian : public drwnStdObjIface {
 protected:
    int _n;               
    int _m;               
    VectorXd _mu;         
    MatrixXd _mSigma;     
    MatrixXd _mSigmaGain; 

 public:
    drwnConditionalGaussian(const VectorXd& mu, const MatrixXd& Sigma,
        const MatrixXd& SigmaGain);
    drwnConditionalGaussian(const drwnConditionalGaussian& model);
    ~drwnConditionalGaussian();

    // i/o
    const char *type() const { return "drwnConditionalGaussian"; }
    drwnConditionalGaussian *clone() const { return new drwnConditionalGaussian(*this); }
    bool save(drwnXMLNode& xml) const;
    bool load(drwnXMLNode& xml);

    // construct gaussian given observation
    drwnGaussian reduce(const VectorXd& x);
    drwnGaussian reduce(const vector<double>& x);
};