StatisticalModel.h

/*
 * This file is part of the statismo library.
 *
 * Author: Marcel Luethi (marcel.luethi@unibas.ch)
 *
 * Copyright (c) 2011 University of Basel
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * Redistributions in binary form must reproduce the above copyright
 * notice, this list of conditions and the following disclaimer in the
 * documentation and/or other materials provided with the distribution.
 *
 * Neither the name of the project's author nor the names of its
 * contributors may be used to endorse or promote products derived from
 * this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
 * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 */


#ifndef STATISTICALMODEL_H_
#define STATISTICALMODEL_H_

#include <limits>
#include <vector>

#include "CommonTypes.h"
#include "Config.h"
#include "DataManager.h"
#include "ModelInfo.h"
#include "Representer.h"

namespace H5 {
class Group;
}

namespace statismo {

/*
template <typename Representer>
class PointValuePair {
public:
  typedef typename Representer::PointType PointType;
  typedef typename Representer::ValueType ValueType;

  PointValuePair(const PointType& pt, const ValueType& val) : point(pt), value(val) {}

  PointType point;
  ValueType value;
};
*/

template <typename T>
class StatisticalModel {
  public:
    typedef Representer<T> RepresenterType ;
    typedef typename RepresenterType::DatasetPointerType DatasetPointerType;
    typedef typename RepresenterType::DatasetConstPointerType DatasetConstPointerType;
    typedef typename RepresenterType::ValueType ValueType;
    typedef typename RepresenterType::PointType PointType;


    typedef Domain<PointType> DomainType;

    typedef unsigned PointIdType;


    //typedef  PointValuePair<Representer>  PointValuePairType;
    typedef std::pair<PointType, ValueType> PointValuePairType;
    typedef std::pair<unsigned, ValueType> PointIdValuePairType;
    typedef std::list<PointValuePairType> PointValueListType;
    typedef std::list<PointIdValuePairType> PointIdValueListType;




    virtual ~StatisticalModel();



    /*
    * Factory method that creates a new Model.
    *
    * \warning The use of this constructor is discouraged. If possible, use a ModelBuilder to create
    * a new model or call Load to load an existing model
    *
    * \param representer the represener
    * \param m the mean
    * \param orthonormalPCABasis An orthonormal matrix with the principal Axes.
    * \param pcaVariance The Variance for each principal Axis
    * \param noiseVariance The variance of the (N(0,noiseVariance)) noise on each point
    */
    static StatisticalModel* Create(const RepresenterType* representer,
                                    const VectorType& m,
                                    const MatrixType& orthonormalPCABasis,
                                    const VectorType& pcaVariance,
                                    double noiseVariance) {
        return new StatisticalModel(representer, m, orthonormalPCABasis, pcaVariance, noiseVariance);
    }



    static StatisticalModel* Load(Representer<T>* representer, const std::string& filename,  unsigned maxNumberOfPCAComponents = std::numeric_limits<unsigned>::max());


    static StatisticalModel* Load(Representer<T>* representer, const H5::Group& modelroot, unsigned maxNumberOfPCAComponents = std::numeric_limits<unsigned>::max());



    void Delete() const {
        delete this;
    }



    void Save(const std::string& filename) const;

    void Save(const H5::Group& modelRoot) const;




    unsigned int GetNumberOfPrincipalComponents() const;

    const ModelInfo& GetModelInfo() const;



    DatasetPointerType DatasetToSample(DatasetConstPointerType dataset) const;

    ValueType EvaluateSampleAtPoint(DatasetConstPointerType sample, unsigned ptId) const ;


    ValueType EvaluateSampleAtPoint(DatasetConstPointerType sample, const PointType& pt) const;


    DatasetPointerType DrawMean() const;

    DatasetPointerType DrawSample(const VectorType& coefficients, bool addNoise = false) const ;

    DatasetPointerType DrawSample(bool addNoise = false) const ;



    DatasetPointerType DrawPCABasisSample(unsigned componentNumber) const;



    ValueType DrawMeanAtPoint( const PointType& point) const;

    ValueType DrawMeanAtPoint( unsigned pointId) const;

    ValueType DrawSampleAtPoint(const VectorType& coefficients, const PointType& point, bool addNoise = false) const;

    ValueType DrawSampleAtPoint(const VectorType& coefficients, unsigned pointId, bool addNoise = false) const;


    MatrixType GetJacobian(const PointType& pt) const ;

    MatrixType GetCovarianceAtPoint(const PointType& pt1, const PointType& pt2) const;

    MatrixType GetCovarianceAtPoint(unsigned ptId1, unsigned ptId2) const;



    MatrixType GetCovarianceMatrix() const;

    double ComputeProbabilityOfDataset(DatasetConstPointerType dataset) const ;

    double ComputeLogProbabilityOfDataset(DatasetConstPointerType dataset) const ;


    double ComputeProbabilityOfCoefficients(const VectorType& coefficients) const ;

    double ComputeLogProbabilityOfCoefficients(const VectorType& coefficients) const ;


    double ComputeMahalanobisDistanceForDataset(DatasetConstPointerType dataset) const;

    VectorType ComputeCoefficientsForDataset(DatasetConstPointerType dataset) const;


    VectorType ComputeCoefficientsForSample(DatasetConstPointerType sample) const;


    VectorType ComputeCoefficientsForPointValues(const PointValueListType&  pointValues, double pointValueNoiseVariance=0.0) const;

    //RB: I had to modify the method name, to avoid prototype collisions when the PointType corresponds to unsigned (= type of the point id)
    VectorType ComputeCoefficientsForPointIDValues(const PointIdValueListType&  pointValues, double pointValueNoiseVariance=0.0) const;

    VectorType RobustlyComputeCoefficientsForDataset(DatasetConstPointerType dataset, unsigned nIterations=100, unsigned nu=6, double sigma2=1) const;


    float GetNoiseVariance() const;

    const VectorType& GetPCAVarianceVector() const;

    const VectorType& GetMeanVector() const;

    const MatrixType& GetPCABasisMatrix() const;

    MatrixType GetOrthonormalPCABasisMatrix() const;

    VectorType DrawSampleVector(const VectorType& coefficients, bool addNoise = false) const ;



    void SetModelInfo(const ModelInfo& modelInfo);

    VectorType ComputeCoefficientsForSampleVector(const VectorType& sample) const;


    const RepresenterType* GetRepresenter() const {
        return m_representer;
    }


    const DomainType& GetDomain() const {
        return m_representer->GetDomain();
    }


  private:
    // computes the M Matrix for the PPCA Method (see Bishop, PRML, Chapter 12)
    void CheckAndUpdateCachedParameters() const;



    StatisticalModel(const RepresenterType* representer, const VectorType& m, const MatrixType& orthonormalPCABasis, const VectorType& pcaVariance, double noiseVariance);

    StatisticalModel(const RepresenterType* representer);

    // to prevent use
    StatisticalModel(const StatisticalModel& rhs);
    StatisticalModel& operator=(const StatisticalModel& rhs);

    const RepresenterType* m_representer;

    VectorType m_mean;
    MatrixType m_pcaBasisMatrix;
    VectorType m_pcaVariance;
    float m_noiseVariance;


    // caching
    mutable bool m_cachedValuesValid;

    //the matrix M^{-1} in Bishops PRML book. This is roughly the Latent Covariance matrix (but not exactly)
    mutable MatrixType m_MInverseMatrix;

    ModelInfo m_modelInfo;
    bool m_modelLoaded;
};

} // namespace statismo

#include "StatisticalModel.hxx"

#endif /* STATISTICALMODEL_H_ */