MMSD_FluryGautschiGaussianLaw.h

Go to the documentation of this file.

#ifndef MMSD_FluryGautschiGaussianLaw_H
#define MMSD_FluryGautschiGaussianLaw_H


#include "MMSD_GaussianLaw.h"


DEFINE_SPTR(MMSD_FluryGautschiGaussianLaw);
class MMSD_FluryGautschiGaussianLaw : public virtual MMSD_GaussianLaw {
    SP_OBJECT(MMSD_FluryGautschiGaussianLaw);
    // ATTRIBUTES
public:
  
    int mMaxIterationsNumber;
  
private:

    MMSD_IntegerFullMatrix mCombinations;
    MMSD_DoubleFullMatrix mBestMatrix,mMatrix,mYM;
    MMSD_DoubleFullMatrix mQ;
  
    
    
protected: 
    // METHODS
  
    // CONSTRUCTORS 
  
    MMSD_FluryGautschiGaussianLaw(void);
  
  
  
    // DESTRUCTORS 
  
  
    virtual ~MMSD_FluryGautschiGaussianLaw(void);   
  
  
  
public:
    //----------------
    // NEW constructors
    // ---------------

    inline static SP::MMSD_FluryGautschiGaussianLaw New() {
        return New(6,1.e-16);
    }
    inline static SP::MMSD_FluryGautschiGaussianLaw New(const int& nFGIterations,const double& minEigenValue) {
        SP::MMSD_FluryGautschiGaussianLaw p(new MMSD_FluryGautschiGaussianLaw(),
                                            MMSD_FluryGautschiGaussianLaw::Delete());
        p->setThis(p);
        p->setMaxIterationsNumber(nFGIterations);
        p->setMinEigenValue(minEigenValue);
        return p;
    }

    // -----------------
    // SET & GET Methods
    // -----------------

    inline void setMaxIterationsNumber(const int& n) {
        mMaxIterationsNumber=n;
    };
    inline int getMaxIterationsNumber() const {
        return mMaxIterationsNumber;
    };


protected:
   
   
public:
    virtual void restore(const MMSD_DoubleFullMatrix& Y) {
        MMSD_GaussianLaw::restore(Y);

        int nProperties=Y.getColumnsNumber();
        if (nProperties>1) mCombinations.computeCombinations(2,nProperties);
        
        
    }
    
    virtual void initialize(const int& clusterIndex,
                            const int& samplesNumber,
                            const MMSD_IntegerVector& clustersIndex,
                            const MMSD_DoubleFullMatrix& Y) {
        MMSD_GaussianLaw::initialize(clusterIndex,samplesNumber,clustersIndex,Y);
        
      
        //compute the combinations
        int nProperties=Y.getColumnsNumber();
        if (nProperties>1) mCombinations.computeCombinations(2,nProperties);
        
        
     
    }
    
    // -----------------------------
    // Parameters estimation methods
    // ------------------------------

private:
    void updateCovarianceDecompositionMatrix(const MMSD_DoubleFullMatrix& Y,
                                             const MMSD_DoubleVector& Diag,
                                             const MMSD_DoubleVector& mean,
                                             const MMSD_DoubleFullMatrix& W,
                                             const MMSD_DoubleVector& probs,
                                             MMSD_IntegerFullMatrix& combinations,
                                             MMSD_DoubleFullMatrix& workOptimalMatrix,
                                             MMSD_DoubleFullMatrix& workMatrix,
                                             MMSD_DoubleFullMatrix& workYMatrix,
                                             MMSD_DoubleFullMatrix& Q,
                                             MMSD_DoubleVector& work,
                                             MMSD_DoubleFullMatrix& P,
                                             MMSD_DoubleFullMatrix& YP,
                                             const int& maxIterationsNumber);
    
    static void updateCovarianceDecompositionMatrix(const MMSD_DoubleFullMatrix& Y,
                                                    MMSD_IntegerFullMatrix& combinations,
                                                    MMSD_DoubleFullMatrix& P,
                                                    MMSD_DoubleFullMatrix& YP,
                                                    const MMSD_DoubleVector& Diag,
                                                    const MMSD_DoubleFullMatrix& W,
                                                    const MMSD_DoubleVector& mean,
                                                    MMSD_DoubleFullMatrix& Q,
                                                    MMSD_DoubleVector& work,
                                                    const MMSD_DoubleVector& probs); 
    
   
    static void  computeGStepForFluryGautschiAlgorithm(const int& C0,const int& C1,
                                                       const MMSD_DoubleFullMatrix& YP,
                                                       const MMSD_DoubleFullMatrix& P,
                                                       const MMSD_DoubleFullMatrix& W,
                                                       const MMSD_DoubleVector& D,
                                                       const MMSD_DoubleVector& probs,
                                                       const MMSD_DoubleVector& mean,
                                                       MMSD_DoubleVector& work,
                                                       MMSD_DoubleFullMatrix& Q);
    
    static double computeQuality(const MMSD_DoubleFullMatrix& YP,
                                 const MMSD_DoubleFullMatrix& P,
                                 const MMSD_DoubleVector& D,
                                 const MMSD_DoubleVector& mean,
                                 const MMSD_DoubleFullMatrix& W,
                                 const MMSD_DoubleVector& probs);
    
 
    
protected:
 
    virtual void updateCovarianceDecompositionMatrix(const MMSD_DoubleFullMatrix& Y,
                                                     const MMSD_DoubleVector& Diag,
                                                     const MMSD_DoubleVector& mean,
                                                     const MMSD_DoubleFullMatrix& W,
                                                     const MMSD_DoubleVector& probs,
                                                     MMSD_DoubleFullMatrix& P,
                                                     MMSD_DoubleFullMatrix& YP) {
     
        updateCovarianceDecompositionMatrix(Y,Diag,mean,
                                            W,probs,
                                            mCombinations,mBestMatrix,mMatrix,mYM,mQ,getWork(),
                                            P,YP,getMaxIterationsNumber());
        
    }


    // ----------------------
    // String representation
    // ----------------------
public:
    virtual tString toString() const {
        tString ret=MMSD_GaussianLaw::toString();
        return ret;
    }
    
};

#endif