setDiagonalStore(Inv_store);
-
- delete [] Inv_store;
-}
-
-
-double* XEMDiagMatrix::getDiagonalStore(){
- return(_store);
-}
-
-
-double* XEMDiagMatrix::getSymmetricStore(){
- throw wrongMatrixType;
-}
-
-double* XEMDiagMatrix::getGeneralStore(){
- throw wrongMatrixType;
-}
-
-double XEMDiagMatrix::getSphericalStore(){
- throw wrongMatrixType;
-}
-
-double XEMDiagMatrix::norme(double * xMoinsMean){
- int64_t p;
- double termesDiag = 0.0;
- double xMoinsMean_p;
-
- for(p=0 ; p<_s_pbDimension ; p++){
- xMoinsMean_p = xMoinsMean[p];
- termesDiag += xMoinsMean_p * xMoinsMean_p * _store[p];
- }
- return termesDiag;
-
-}
-
-void XEMDiagMatrix::compute_as__multi_O_S_O(double multi, XEMGeneralMatrix* & O, XEMDiagMatrix *& S){
- throw nonImplementedMethod;
-}
-double XEMDiagMatrix::trace_this_O_Sm1_O(XEMGeneralMatrix* & O, XEMDiagMatrix* & S){
- throw nonImplementedMethod;
-}
-double XEMDiagMatrix::compute_trace_W_C(XEMMatrix * C){
- throw nonImplementedMethod;
-}
-void XEMDiagMatrix::computeShape_as__diag_Ot_this_O(XEMDiagMatrix* & Shape, XEMGeneralMatrix* & Ori, double diviseur){
- throw nonImplementedMethod;
-}
-
-
-
-double XEMDiagMatrix::putSphericalValueInStore(double & store){
- store = 0.0;
- int64_t p;
- for(p=0; p<_s_pbDimension; p++){
- store += _store[p];
- }
- store /= _s_pbDimension;
- return(store);
-}
-
-
-double XEMDiagMatrix::addSphericalValueInStore(double & store){
- int64_t p;
- for(p=0; p<_s_pbDimension; p++){
- store += _store[p];
- }
- store /= _s_pbDimension;
- return(store);
-}
-
-
-
-double* XEMDiagMatrix::putDiagonalValueInStore(double * store){
- for (int64_t p=0 ; p<_s_pbDimension; p++){
- store[p] = _store[p];
- }
- return(store);
-}
-
-double* XEMDiagMatrix::addDiagonalValueInStore(double * store){
- for (int64_t p=0 ; p<_s_pbDimension; p++){
- store[p] += _store[p];
- }
- return(store);
-}
-
-double* XEMDiagMatrix::addSymmetricValueInStore(double * store){
- // return the store of of a symmetric matrix with this on the diag
- //int64_t dimStore = _s_pbDimension*(_s_pbDimension+1)/2;
- // double * store = new double[dimStore];
-
- int64_t p,q,r;
- for(p=0,r=0; p<_s_pbDimension; p++,r++){
- for(q=0; qputDiagonalValueInStore(_store);
-
- int64_t p;
- for(p=0; p<_s_pbDimension; p++){
- _store[p] /= d;
- }
-}
-
-
-void XEMDiagMatrix::equalToMatrixMultiplyByDouble(XEMMatrix* D, double d){
- D->putDiagonalValueInStore(_store);
- int64_t p;
- for(p=0; p<_s_pbDimension; p++){
- _store[p] *= d;
- }
-}
-
-
-
-// add : cik * xMoinsMean * xMoinsMean' to this
-void XEMDiagMatrix::add(double * xMoinsMean, double cik){
-
- int64_t p;
- double xMoinsMean_p;
-
- for(p=0; p<_s_pbDimension ; p++){
- xMoinsMean_p = xMoinsMean[p];
- _store[p] += cik * xMoinsMean_p * xMoinsMean_p;
- }//end for p
-
-}
-
-// add : diag( cik * xMoinsMean * xMoinsMean' ) to this
-/*void XEMDiagMatrix::addDiag(double * xMoinsMean, double cik){
-
- int64_t p;
- double xMoinsMean_p;
-
- for(p=0; p<_s_pbDimension ; p++){
- xMoinsMean_p = xMoinsMean[p];
- _store[p] += cik * xMoinsMean_p * xMoinsMean_p;
- }//end for p
-
- }*/
-
-
-
-// set the value of (d x Identity) to this
-void XEMDiagMatrix::operator=(const double& d){
- int64_t p;
-
- for(p=0; p<_s_pbDimension; p++){
- _store[p] = d;
- }
-}
-
-// divide each element by d
-void XEMDiagMatrix::operator/=(const double& d){
- int64_t p;
- for(p=0; p<_s_pbDimension; p++){
- _store[p] /= d;
- }
-}
-
-
-// multiply each element by d
-void XEMDiagMatrix::operator*=(const double& d){
- int64_t p;
- for(p=0; p<_s_pbDimension; p++){
- _store[p] *= d;
- }
-}
-
-
-//add M to this
-void XEMDiagMatrix::operator+=(XEMMatrix* M){
- M -> addDiagonalValueInStore(_store);
-}
-
-
-void XEMDiagMatrix::operator=(XEMMatrix* M){
- M -> putDiagonalValueInStore(_store);
-}
-
-
-
-
-
-void XEMDiagMatrix::input(ifstream & fi){
- int64_t p,q;
- double garbage;
-
- for (p=0; p<_s_pbDimension; p++){
- // useless because all are 0
- for (q=0; q> garbage;
- }
-
- // here i==j so we are in the diagonal
- fi >> _store[p];
-
- // useless because all are 0
- for(q=p+1;q<_s_pbDimension;q++){
- fi >> garbage;
- }
- }
-
-}
-
-double XEMDiagMatrix::detDiag(XEMErrorType errorType){
- return determinant(errorType);
-}
-
-void XEMDiagMatrix ::sortDiagMatrix(){
- int64_t max;
- for (int64_t i=0; i<_s_pbDimension; i++){
- max = i;
- for (int64_t j=i+1; j<_s_pbDimension; j++){
- if (_store[j] > _store[max]){
- max = j;
- }
- }
- if (max != i ){ // swich
- double tmp = _store[i];
- _store[i] = _store[max];
- _store[max] = tmp;
- }
- }
-
- /*for (int64_t i=1; i<= _s_pbDimension;i++){
- //search the max eigenvalue
- max = i;
- for (int64_t j=i; j<_s_pbDimension; j++){
- if (_store[j-1] > _store[max-1])
- max = j;
- if (max != i){
- // switch
- double tmp = _store[max-1];
- _store[max-1] = _store[i-1];
- _store[i-1] = tmp;
- }
- }
- }*/
-
-}
-
-
-double** XEMDiagMatrix::storeToArray() const{
-
- int64_t i,j;
- double** newStore = new double*[_s_pbDimension];
- for (i=0; i<_s_pbDimension ; ++i){
- newStore[i] = new double[_s_pbDimension];
- }
- for (i=0; i<_s_pbDimension ; ++i){
-
- for (j=0; j<_s_pbDimension ; ++j){
- if (i == j){
- newStore[i][j] = _store[i];
- }
- else {
- newStore[i][j] = 0;
- }
- }
- }
-
- return newStore;
-
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/lib/src/mixmod/MIXMOD/XEMDiagMatrix.h b/lib/src/mixmod/MIXMOD/XEMDiagMatrix.h
deleted file mode 100644
index db71b91..0000000
--- a/lib/src/mixmod/MIXMOD/XEMDiagMatrix.h
+++ /dev/null
@@ -1,169 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMDiagMatrix.h description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#ifndef XEMDIAGMATRIX_H
-#define XEMDIAGMATRIX_H
-
-#include "XEMUtil.h"
-#include "XEMMatrix.h"
-#include "XEMGeneralMatrix.h"
-
-/**
- @brief class XEMDiagMatrix
- @author F Langrognet & A Echenim
-*/
-
-class XEMDiagMatrix : public XEMMatrix{
-
- public:
-
- /// Default constructor
- XEMDiagMatrix();
-
- /// contructor : d*Id
- /// default value = Id
- XEMDiagMatrix(int64_t pbDimension, double d=1.0);
-
- XEMDiagMatrix(XEMDiagMatrix * A);
-
-
- /// Desctructor
- virtual ~XEMDiagMatrix();
-
- /// compute determinant of diagonal matrix
- double determinant(XEMErrorType errorType);
- /// return store of diagonal matrix
- double * getStore();
- /// compute inverse of diagonal matrix
- void inverse(XEMMatrix * & A);
-
- void compute_product_Lk_Wk(XEMMatrix* Wk, double L);
-
- /// compute (x - mean)' this (x - mean)
- double norme(double * xMoinsMean);
-
- /// (this) will be A / d
- void equalToMatrixDividedByDouble(XEMMatrix * A, double d);
-
- /// this = matrix * d
- void equalToMatrixMultiplyByDouble(XEMMatrix*D, double d);
-
- ///compute singular vector decomposition
- void computeSVD(XEMDiagMatrix* & S, XEMGeneralMatrix* & O);
-
- /// compute trace of general matrix
- double computeTrace();
-
- /// add : cik * xMoinsMean * xMoinsMean' to this
- void add(double * xMoinsMean, double cik);
-
- // add : diag( cik * xMoinsMean * xMoinsMean' ) to this
- //void addDiag(double * xMoinsMean, double cik);
-
- /// set the value of (d x Identity) to this
- void operator=(const double& d);
- /// this = this / (d * Identity)
- void operator/=(const double& d);
- /// this = this * (d * Identity)
- void operator*=(const double& d);
- /// this = this + matrix
- void operator+=(XEMMatrix* M);
- /// this = matrix
- void operator=(XEMMatrix* M);
-
- /// Return store of a spherical matrix in a diagonal one
- double putSphericalValueInStore(double & store);
- /// Add store of a spherical matrix in a diagonal one
- double addSphericalValueInStore(double & store);
-
- double getSphericalStore();
-
- /// Return store of a diagonal matrix
- double* putDiagonalValueInStore(double * store);
- /// Add store of a diagonal matrix in a diagonal one
- double* addDiagonalValueInStore(double * store);
-
- double* getDiagonalStore();
-
- /// Return store of a diagonal matrix
- double* putSymmetricValueInStore(double * store);
- /// Add store of a diagonal matrix in a diagonal one
- double* addSymmetricValueInStore(double * store);
-
- double* getSymmetricStore();
-
- /// Return store of a diagonal matrix
- double* putGeneralValueInStore(double * store);
- /// Add store of a diagonal matrix in a diagonal one
- double* addGeneralValueInStore(double * store);
-
- double* getGeneralStore();
-
- /// read general matrix in an input file
- void input(ifstream & fi);
-
- ///set store
- void setSymmetricStore(double * store);
- void setGeneralStore(double * store);
- void setDiagonalStore(double * store);
- void setSphericalStore(double store);
- double** storeToArray() const;
-
- /// gives : det(diag(this))
- double detDiag(XEMErrorType errorType);
-
- void compute_as__multi_O_S_O(double multi, XEMGeneralMatrix* & O, XEMDiagMatrix *& S);
- double trace_this_O_Sm1_O(XEMGeneralMatrix* & O, XEMDiagMatrix* & S);
- double compute_trace_W_C(XEMMatrix * C);
- void computeShape_as__diag_Ot_this_O(XEMDiagMatrix* & Shape, XEMGeneralMatrix* & Ori, double diviseur = 1.0);
-
- ///sort diagonal matrix in decreasing order
- void sortDiagMatrix();
- protected :
-
- double * _store;
-
-};
-
-inline double * XEMDiagMatrix::getStore(){
- return _store;
-}
-
-inline void XEMDiagMatrix::setSymmetricStore(double * store){
- throw wrongMatrixType;
-}
-
-inline void XEMDiagMatrix::setGeneralStore(double * store){
- throw wrongMatrixType;
-}
-
-inline void XEMDiagMatrix::setDiagonalStore(double * store){
- //_store = store;
- recopyTab(store,_store,_s_pbDimension);
-}
-
-inline void XEMDiagMatrix::setSphericalStore(double store){
- throw wrongMatrixType;
-}
-#endif
diff --git a/lib/src/mixmod/MIXMOD/XEMEMAlgo.cpp b/lib/src/mixmod/MIXMOD/XEMEMAlgo.cpp
deleted file mode 100644
index 23368df..0000000
--- a/lib/src/mixmod/MIXMOD/XEMEMAlgo.cpp
+++ /dev/null
@@ -1,109 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMEMAlgo.cpp description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#include "XEMEMAlgo.h"
-#include "XEMUtil.h"
-#include
-#include
-
-//-----------
-//Constructor
-//-----------
-XEMEMAlgo::XEMEMAlgo(){
-}
-
-XEMEMAlgo::XEMEMAlgo(XEMAlgoStopName algoStopName, double epsilon, int64_t nbIteration)
- :XEMAlgo( algoStopName, epsilon, nbIteration){
-}
-
-
-/// Copy constructor
-XEMEMAlgo::XEMEMAlgo(const XEMEMAlgo & emAlgo):XEMAlgo(emAlgo){
-}
-
-
-//----------
-//Destructor
-//----------
-XEMEMAlgo::~XEMEMAlgo(){
-}
-
-
-
-// clone
-//------
-XEMAlgo * XEMEMAlgo::clone(){
- return (new XEMEMAlgo(*this));
-}
-
-//---
-//run
-//---
-void XEMEMAlgo::run(XEMModel *& model){
-
- /*std::ostringstream numero;
- std::string nomfic;
-
- std::string basestring ="parameter";
- std::string suffix = ".txt" ;
-
- ofstream likeli("likelihoods.txt", ios::out);
- */
-
- _indexIteration = 1;
- model->setAlgoName(EM);
-
-#if DEBUG > 0
- cout<<"Debut Algo EM :"<editDebugInformation();
-#endif
-
- model->Estep(); // E Step
-
-
- model->Mstep(); // M Step
-
-#if DEBUG > 0
- cout<<"Apres la 1ere iteration de EM :"<editDebugInformation();
-#endif
- _indexIteration++;
-
- while (continueAgain()){
- model->Estep(); // E Step
- model->Mstep(); // M Step
-#if DEBUG > 0
- cout<<"Apres la "<<_indexIteration<<" eme iteration de EM :"<editDebugInformation();
-#endif
- _indexIteration++;
- _xml_old = _xml;
- _xml = model->getLogLikelihood(true); // true : to compute fik
-
- }
-
- model->Estep(); // E step to update Tik an fik
-
-}
-
diff --git a/lib/src/mixmod/MIXMOD/XEMEMAlgo.h b/lib/src/mixmod/MIXMOD/XEMEMAlgo.h
deleted file mode 100644
index 336087c..0000000
--- a/lib/src/mixmod/MIXMOD/XEMEMAlgo.h
+++ /dev/null
@@ -1,64 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMEMAlgo.h description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#ifndef XEMEMALGO_H
-#define XEMEMALGO_H
-
-#include "XEMAlgo.h"
-
-/**
- @brief Derived class of XEMAlgo for EM Algorithm(s)
- @author F Langrognet & A Echenim
-*/
-
-class XEMEMAlgo : public XEMAlgo{
-
- public:
-
- /// Default constructor
- XEMEMAlgo();
-
- /// Copy constructor
- XEMEMAlgo(const XEMEMAlgo & emAlgo);
-
- /// Constructor
- XEMEMAlgo(XEMAlgoStopName algoStopName, double epsilon, int64_t nbIteration);
-
- /// Destructor
- ~XEMEMAlgo();
-
- /// clone
- XEMAlgo * clone();
-
- /// Run method
- virtual void run(XEMModel *& model);
-
- const XEMAlgoName getAlgoName() const;
-};
-
-inline const XEMAlgoName XEMEMAlgo::getAlgoName() const{
- return EM;
-}
-
-#endif
diff --git a/lib/src/mixmod/MIXMOD/XEMError.cpp b/lib/src/mixmod/MIXMOD/XEMError.cpp
deleted file mode 100644
index fd84e97..0000000
--- a/lib/src/mixmod/MIXMOD/XEMError.cpp
+++ /dev/null
@@ -1,65 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMError.cpp description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#include "XEMError.h"
-
-//-----------
-//Constructor
-//-----------
-XEMError::XEMError(){
-}
-
-
-
-//-----------
-//Constructor
-//-----------
-XEMError::XEMError(XEMErrorType errorType){
- _errorType = errorType;
-}
-
-
-
-//----------
-//Destructor
-//----------
-XEMError::~XEMError(){
-}
-
-
-
-//---
-//Run
-//---
-void XEMError::run(ostream & flux){
-
-
- flux << "MIXMOD ERROR ("<<_errorType<<") :"<.
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#ifndef XEMERROR_H
-#define XEMERROR_H
-
-#include "XEMUtil.h"
-
-/**
- @brief Base class for Error(s)
- @author F Langrognet & A Echenim
-*/
-
-class XEMError{
-
- public:
-
- /// Default constructor
- XEMError();
-
- /// Constructor
- XEMError(XEMErrorType errorType);
-
- /// Destructor
- ~XEMError();
-
-
- /// Run method
- void run(ostream & flux=std::cout);
-
-
- private :
-
- /// Type of error
- XEMErrorType _errorType;
-};
-
-#endif
diff --git a/lib/src/mixmod/MIXMOD/XEMEstimation.cpp b/lib/src/mixmod/MIXMOD/XEMEstimation.cpp
deleted file mode 100644
index 1e8b942..0000000
--- a/lib/src/mixmod/MIXMOD/XEMEstimation.cpp
+++ /dev/null
@@ -1,169 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMEstimation.cpp description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#include "XEMEstimation.h"
-#include "XEMStrategy.h"
-#include "XEMClusteringStrategy.h"
-#include "XEMBICCriterion.h"
-#include "XEMICLCriterion.h"
-#include "XEMNECCriterion.h"
-#include "XEMCVCriterion.h"
-#include "XEMDCVCriterion.h"
-
-
-//------------
-// Constructor
-//------------
-XEMEstimation::XEMEstimation(){
- _strategy = NULL;
- _cStrategy = NULL;
- _modelType = NULL;
- _errorType = noError;
- _model = NULL;
-}
-
-
-//------------------
-// copy constructor
-//-----------------
-XEMEstimation::XEMEstimation(const XEMEstimation & estimation){
- throw internalMixmodError;
-}
-
-
-//------------
-// Constructor
-//------------
-XEMEstimation::XEMEstimation(XEMClusteringStrategy *& strategy, XEMModelType * modelType, int64_t nbCluster, XEMData *& data ,vector & criterionName, XEMPartition *& knownPartition, vector corresp){
- _modelType = modelType;
- _nbCluster = nbCluster;
- _cStrategy = strategy;
- _strategy = NULL;
- _errorType = noError;
- _correspondenceOriginDataToReduceData = corresp;
- _model = new XEMModel(modelType, nbCluster, data, knownPartition);
-
- int64_t nbCriterion = criterionName.size();
- _criterionOutput.resize(nbCriterion);
- _criterion.resize(nbCriterion);
- for (int64_t i=0; i corresp){
- _modelType = modelType;
- _nbCluster = nbCluster;
- _strategy = strategy;
- _cStrategy = NULL;
- _errorType = noError;
- _correspondenceOriginDataToReduceData = corresp;
- _model = new XEMModel(modelType, nbCluster, data, knownPartition);
-}
-
-
-
-//-----------
-// Destructor
-//-----------
-XEMEstimation::~XEMEstimation(){
- // if (_strategy){
- // delete _strategy;
- // _strategy = NULL;
- // }
- if (_model){
- delete _model;
- _model = NULL;
- }
-}
-
-
-
-//---
-//run
-//---
-void XEMEstimation::run(){
- try{
- // compute strategy
- //-----------------
- if (_strategy){
- _strategy->run(_model);
- }
- else{
- _cStrategy->run(_model);
- }
- }
- catch(XEMErrorType & e){
- setErrorType(e);
- for (int64_t i=0; i<_criterion.size(); i++){
- _criterionOutput[i].setCriterionName(_criterion[i]->getCriterionName());
- _criterionOutput[i].setError(errorEstimationStrategyRun);
- _criterionOutput[i].setValue(0);
- }
- return;
- }
-
- // compute criterion
- //------------------
- double value;
- XEMErrorType error;
- for (int64_t i=0; i<_criterion.size(); i++){
- _criterion[i]->run(_model, value, error);
- _criterionOutput[i].setCriterionName(_criterion[i]->getCriterionName());
- _criterionOutput[i].setError(error);
- _criterionOutput[i].setValue(value);
- }
-
-}
-
-//---------
-//set error
-//---------
-void XEMEstimation::setErrorType(XEMErrorType errorType){
- _errorType = errorType;
-}
-
-
-
-//--------------
-//Friend Methods
-//--------------
-ostream & operator << (ostream & fo, XEMEstimation & estimation){
- return fo;
-}
-
diff --git a/lib/src/mixmod/MIXMOD/XEMEstimation.h b/lib/src/mixmod/MIXMOD/XEMEstimation.h
deleted file mode 100644
index 62bddba..0000000
--- a/lib/src/mixmod/MIXMOD/XEMEstimation.h
+++ /dev/null
@@ -1,215 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMEstimation.h description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#ifndef XEMESTIMATION_H
-#define XEMESTIMATION_H
-
-
-#include "XEMUtil.h"
-#include "XEMModel.h"
-#include "XEMCriterionOutput.h"
-#include "XEMClusteringStrategy.h"
-
-class XEMCriterion;
-
-/**
- @brief Base class for Estimation(s)
- @author F Langrognet & A Echenim
-*/
-
-class XEMEstimation{
-
- public:
-
- /// Default constructor
- XEMEstimation();
-
- /// copy constructor
- XEMEstimation(const XEMEstimation & estimation);
-
-
- /// Constructor
- XEMEstimation(XEMClusteringStrategy *& strategy, XEMModelType * modelType, int64_t nbCluster, XEMData *& data, vector & criterionName, XEMPartition *& knownPartition, vector corresp);
-
-
- /// Constructor
- XEMEstimation(XEMStrategy *& strategy, XEMModelType * modelType, int64_t nbCluster, XEMData *& data, XEMPartition *& knownPartition, vector corresp);
-
- /// Destructor
- virtual ~XEMEstimation();
-
-
- /** @brief Selector
- @return The best Model of the current strategy
- */
- XEMModel * getModel();
-
-
- /// get Parameter
- XEMParameter * getParameter() const;
-
- /// get model type
- XEMModelType * getModelType();
-
- /// get nbCluster
- int64_t getNbCluster();
-
- /** @brief Set the error type estimation
- @param errorType Type of error to set
- */
- void setErrorType(XEMErrorType errorType);
-
- /** @brief Selector
- @return The type of the error
- */
- XEMErrorType getErrorType();
-
- /// get the strategy
- XEMClusteringStrategy * getClusteringStrategy();
- /// get the strategy
- XEMStrategy * getStrategy();
-
- //getcorrespondenceOriginDataToReduceData
- const vector & getcorrespondenceOriginDataToReduceData()const;
-
- /// Run method
- void run();
-
- ///
- int64_t getCriterionSize() const;
-
- //XEMCriterionName getCriterionName(int64_t index) const;
-
- /// getCriterionOutput
- vector getCriterionOutput() ;
-
- /// getCriterionOutput(index)
- XEMCriterionOutput & getCriterionOutput(int64_t index) ;
-
- // getData
- XEMData * getData();
-
-
- /// getNbSample
- int64_t getNbSample() const;
-
- /** @brief Friend method
- @return Operator << overloaded to write Estimation in output files
- */
- friend ostream & operator << (ostream & fo, XEMEstimation & estimation);
-
- private :
-
- /// Current strategy
- XEMClusteringStrategy * _cStrategy;
-
- /// Current strategy
- XEMStrategy * _strategy;
-
- /// Number of cluster
- int64_t _nbCluster;
-
- /// Current model type
- XEMModelType * _modelType;
-
- /// Current model error
- XEMErrorType _errorType;
-
- /// model
- XEMModel * _model;
-
- vector _correspondenceOriginDataToReduceData;
-
- vector _criterion;
- vector _criterionOutput;
-
-};
-
-
-
-//---------------
-// inline methods
-//---------------
-
-inline XEMModel * XEMEstimation::getModel(){
- return (_model);
-}
-
-inline XEMData * XEMEstimation::getData(){
- return (_model->getData());
-}
-
-inline XEMModelType * XEMEstimation::getModelType(){
- return _modelType;
-}
-
-inline int64_t XEMEstimation::getNbCluster(){
- return _nbCluster;
-}
-
-
-inline int64_t XEMEstimation::getNbSample() const{
- return (_model->getNbSample());
-}
-
-inline XEMClusteringStrategy * XEMEstimation::getClusteringStrategy(){
- return _cStrategy;
-}
-inline XEMStrategy * XEMEstimation::getStrategy(){
- return _strategy;
-}
-
-inline XEMErrorType XEMEstimation::getErrorType(){
- return (_errorType);
-}
-
-inline const vector & XEMEstimation::getcorrespondenceOriginDataToReduceData()const{
- return _correspondenceOriginDataToReduceData;
-}
-
-inline XEMParameter * XEMEstimation::getParameter() const{
- if (_model){
- return _model->getParameter();
- }
- else{
- throw nullPointerError;
- }
-}
-
-inline int64_t XEMEstimation::getCriterionSize() const{
- return _criterion.size();
-}
-
-inline vector XEMEstimation::getCriterionOutput(){
- return _criterionOutput;
-}
-
-inline XEMCriterionOutput & XEMEstimation::getCriterionOutput(int64_t index){
- if (index>=0 && index<_criterionOutput.size())
- return _criterionOutput[index];
- throw internalMixmodError;
-}
-
-
-#endif
diff --git a/lib/src/mixmod/MIXMOD/XEMGaussianData.cpp b/lib/src/mixmod/MIXMOD/XEMGaussianData.cpp
deleted file mode 100644
index 6ba7336..0000000
--- a/lib/src/mixmod/MIXMOD/XEMGaussianData.cpp
+++ /dev/null
@@ -1,355 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMGaussianData.cpp description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#include "XEMGaussianData.h"
-
-
-//-----------
-//Constructor
-//-----------
-XEMGaussianData::XEMGaussianData(){
-}
-
-
-
-//------------
-// Constructor
-//------------
-XEMGaussianData::XEMGaussianData(const XEMGaussianData & iData):XEMData(iData){
- int64_t i;
- XEMSample ** matrix = iData._matrix;
-
- _matrix = new XEMSample*[_nbSample];
- _yStore = new double* [_nbSample];
-
- for (i=0; i<_nbSample; i++){
- _matrix[i] = new XEMGaussianSample((XEMGaussianSample*)matrix[i]);
- _yStore[i] = ((XEMGaussianSample *)_matrix[i])->getTabValue() ;
- }
- _Inv2PiPow = iData.getInv2PiPow();
- _pbDimensionLog2Pi = iData.getPbDimensionLog2Pi();
- _halfPbDimensionLog2Pi = _pbDimensionLog2Pi / 2.0 ;
- __tmpTabOfSizePbDimension = new double[_pbDimension];
- _deleteSamples = true;
-}
-
-
-
-//------------
-// Constructor
-//------------
-XEMGaussianData::XEMGaussianData(int64_t nbSample, int64_t pbDimension):XEMData(nbSample,pbDimension){
- int64_t i;
- _Inv2PiPow = 1.0 / pow(2.0 * XEMPI, pbDimension/2.0 );
- _pbDimensionLog2Pi = pbDimension * log(2.0 * XEMPI);
- _halfPbDimensionLog2Pi = _pbDimensionLog2Pi / 2.0;
- __tmpTabOfSizePbDimension = new double[_pbDimension];
-
-
- _matrix = new XEMSample*[_nbSample];
- _yStore = new double * [_nbSample];
-
- for (i=0; i<_nbSample; i++){
- _weight[i] = 1.0;
- _matrix[i] = new XEMGaussianSample(_pbDimension);
- _yStore[i] = ((XEMGaussianSample *)(_matrix[i]))->getTabValue();
- }
- _weightTotal = _nbSample;
-
-}
-
-
-
-//------------
-// Constructor
-//------------
-XEMGaussianData::XEMGaussianData(int64_t nbSample, int64_t pbDimension, double ** matrix):XEMData(nbSample,pbDimension){
- int64_t i;
-
- if (matrix == NULL) throw internalMixmodError;
-
- _Inv2PiPow = 1.0 / pow(2.0 * XEMPI, pbDimension/2.0 );
- _pbDimensionLog2Pi = pbDimension * log(2.0 * XEMPI);
- _halfPbDimensionLog2Pi = _pbDimensionLog2Pi / 2.0;
- __tmpTabOfSizePbDimension = new double[_pbDimension];
-
- _matrix = new XEMSample*[_nbSample];
- _yStore = new double * [_nbSample];
-
- for (i=0; i<_nbSample; i++){
- _weight[i] = 1.0;
- _matrix[i] = new XEMGaussianSample(_pbDimension, matrix[i]);
- _yStore[i] = ((XEMGaussianSample *)(_matrix[i]))->getTabValue();
- }
- _weightTotal = _nbSample;
-}
-
-
-
-//------------
-// Constructor
-//------------
-XEMGaussianData::XEMGaussianData(int64_t nbSample, int64_t pbDimension, const string & dataFileName):XEMData(nbSample,pbDimension){
- int64_t i;
-
- _Inv2PiPow = 1.0 / pow(2.0 * XEMPI, pbDimension/2.0 );
- _pbDimensionLog2Pi = pbDimension * log(2.0 * XEMPI);
- _halfPbDimensionLog2Pi = _pbDimensionLog2Pi / 2.0;
- __tmpTabOfSizePbDimension = new double[_pbDimension];
-
-
- _matrix = new XEMSample*[_nbSample];
- _yStore = new double * [_nbSample];
-
- for (i=0; i<_nbSample; i++){
- _matrix[i] = new XEMGaussianSample(_pbDimension);
- _yStore[i] = ((XEMGaussianSample *)(_matrix[i]))->getTabValue();
- }
-
- ifstream dataStream((dataFileName).c_str(), ios::in);
- if (! dataStream.is_open()){
- throw wrongDataFileName;
- }
- input(dataStream);
- dataStream.close();
- _deleteSamples = true;
- _fileNameData = dataFileName;
-}
-
-
-
-//------------
-// Constructor for dataReduce
-//------------
-XEMGaussianData::XEMGaussianData(int64_t nbSample, int64_t pbDimension, double weightTotal, XEMSample **& matrix, double * weight):XEMData(nbSample,pbDimension,weightTotal,weight){
-
- // 1/ (2 * pi)^(d/2)
- _Inv2PiPow = 1.0 / pow(2.0 * XEMPI, pbDimension/2.0 );
- _pbDimensionLog2Pi = pbDimension * log(2.0 * XEMPI);
- _halfPbDimensionLog2Pi = _pbDimensionLog2Pi / 2.0;
- _pbDimensionLog2Pi = pbDimension * log(2.0 * XEMPI);
- __tmpTabOfSizePbDimension = new double[_pbDimension];
-
- _matrix = matrix;
- _yStore = new double *[nbSample];
- int64_t i;
-
- for (i=0; i<_nbSample; i++){
- _yStore[i] = ((XEMGaussianSample *)(_matrix[i]))->getTabValue();
- }
- _deleteSamples = true;
-}
-
-
-
-//------------
-// Constructor (used in DCV context)
-//------------
-XEMGaussianData::XEMGaussianData(int64_t nbSample, int64_t pbDimension, XEMData * originalData, XEMCVBlock & block):XEMData(nbSample, pbDimension){
-
- XEMGaussianData * origData = (XEMGaussianData *)(originalData);
- XEMSample ** origMatrix = origData->_matrix ;
-
- // 1/ (2 * pi)^(d/2)
- _Inv2PiPow = 1.0 / pow(2.0 * XEMPI, pbDimension/2.0 );
- _pbDimensionLog2Pi = pbDimension * log(2.0 * XEMPI);
- _halfPbDimensionLog2Pi = _pbDimensionLog2Pi / 2.0;
- __tmpTabOfSizePbDimension = new double[_pbDimension];
- _deleteSamples = false;
-
-
- _weightTotal = block._weightTotal;
- _matrix = new XEMSample*[_nbSample] ;
- for (int64_t i=0; i<_nbSample; i++){
- _matrix[i] = origMatrix[block._tabWeightedIndividual[i].val];
- //cout<<"ind : "<getTabValue();
- }
-}
-
-
-
-
-//----------
-//Destructor
-//----------
-XEMGaussianData::~XEMGaussianData(){
- int64_t i;
- if (_matrix){
-
- if(_deleteSamples){
- for (i=0; i<_nbSample; i++){
- delete _matrix[i];
- _matrix[i] = NULL;
- }
- }
-
- delete[] _matrix;
- _matrix = NULL;
- }
- if(_yStore){
- delete[] _yStore;
- _yStore = NULL;
- }
- if(__tmpTabOfSizePbDimension){
- delete[] __tmpTabOfSizePbDimension;
- __tmpTabOfSizePbDimension = NULL;
- }
-
-}
-
-
-
-//-----------
-// Clone data
-//-----------
-XEMData * XEMGaussianData::clone() const{
- XEMGaussianData * newData = new XEMGaussianData(*this);
- return(newData);
-}
-
-
-
-//------------------
-// Clone data matrix
-//------------------
-XEMSample ** XEMGaussianData::cloneMatrix(){
- int64_t i;
-
- XEMSample ** newMatrix = new XEMSample*[_nbSample];
- for (i=0; i<_nbSample; i++){
- newMatrix[i] = new XEMGaussianSample((XEMGaussianSample*)_matrix[i]);
- }
-
- return newMatrix;
-}
-
-
-
-//------
-// input
-//------
-void XEMGaussianData:: input(ifstream & fi){
-
- int64_t j;
- int64_t i;
- double ** p_y;
- double * p_y_i;
-
- p_y = _yStore;
- for(i=0; i<_nbSample; i++){
- p_y_i = *p_y;
- for(j=0; j<_pbDimension; j++){
- if (fi.eof()){
- throw endDataFileReach;
- }
- fi >> p_y_i[j];
- }
- _weight[i] = 1.0;
- p_y++;
- }
- _weightTotal = _nbSample;
-
-}
-
-
-
-//------
-// input
-//------
-void XEMGaussianData::input(const XEMDataDescription & dataDescription){
- _fileNameData = dataDescription.getFileName();
-
- _weightTotal = 0;
- ifstream fi((_fileNameData).c_str(), ios::in);
- if (! fi.is_open()){
- throw wrongDataFileName;
- }
- int64_t j;
- int64_t i;
- double tmp;
- double ** p_y;
- double * p_y_i;
-
- p_y = _yStore;
- int64_t g=0;
- for(i=0; i<_nbSample; i++){
- _weight[i] = 1.0;
- p_y_i = *p_y;
- g=0;
- for(j=0; j> p_y_i[g];
- g++;
- //cout<>_weight[i];
- }
- else{
- if (typeid(*(dataDescription.getColumnDescription(j)))==typeid(XEMIndividualColumnDescription)){
- string stringTmp;
- fi>>stringTmp;
- //cout<>tmp; // on avance !
- }
- }
- }
- }// end j
- p_y++;
- _weightTotal += _weight[i];
- }// end i
-}
-
-
-
-// output
-//-------
-void XEMGaussianData:: output(ostream & fo){
- cout<<"Sample size: "<<_nbSample << endl;
- cout<<" Dimension: "<<_pbDimension << endl ;
-
- editTab(_yStore,_nbSample,_pbDimension," ","",fo);
-
-}
-
-
-bool XEMGaussianData::verify()const{
- bool res = XEMData::verify();
-
- // others verif ?
- return res;
-}
diff --git a/lib/src/mixmod/MIXMOD/XEMGaussianData.h b/lib/src/mixmod/MIXMOD/XEMGaussianData.h
deleted file mode 100644
index c9b9e57..0000000
--- a/lib/src/mixmod/MIXMOD/XEMGaussianData.h
+++ /dev/null
@@ -1,153 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMGaussianData.h description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#ifndef XEMGAUSSIANDATA_H
-#define XEMGAUSSIANDATA_H
-
-#include "XEMUtil.h"
-#include "XEMData.h"
-#include "XEMGaussianSample.h"
-#include "XEMOldInput.h"
-
-/**
- @brief Base class for Gaussian Data
- @author F Langrognet & A Echenim
-*/
-
-class XEMGaussianData : public XEMData{
-
- public:
-
- /// Default constructor
- XEMGaussianData();
-
- /// Constructor
- XEMGaussianData(const XEMGaussianData & iData);
-
- /// Constructor
- XEMGaussianData(int64_t nbSample, int64_t pbDimension, const string & dataFileName);
-
- /// Constructor (without fill matrix)
- XEMGaussianData(int64_t nbSample, int64_t pbDimension);
-
- /// Constructor (with matrix)
- XEMGaussianData(int64_t nbSample, int64_t pbDimension, double ** matrix);
-
- /// Constructor
- XEMGaussianData(int64_t nbSample, int64_t pbDimension, double weightTotal, XEMSample **& matrix, double * weight);
-
- /// Constructor
- // used in DCV context
- XEMGaussianData(int64_t nbSample, int64_t pbDimension, XEMData * originalData, XEMCVBlock & block);
-
- /// Desctructor
- virtual ~XEMGaussianData();
-
-
- /** @brief Selector
- @return A copy of data
- */
- XEMData * clone() const;
-
- /** @brief Copy
- @return A copy data matrix
- */
- XEMSample ** cloneMatrix();
-
-
- //TODO a enlever XEMInput
- /** @brief Read data from gaussian data file
- @fi Gaussian Data file to read
- */
- void input(ifstream & fi);
-
- /** @brief Read data from XEMDataDescription
- */
- void input(const XEMDataDescription & dataDescription);
-
- /** @brief Write gaussian data in output file
- @f0 Output file to write into
- */
- void output(ostream & fo);
-
- bool verify()const;
-
- /// pointer to stored values
- double ** _yStore ;
-
- double ** getYStore();
-
- double getInv2PiPow() const;
-
- double getHalfPbDimensionLog2Pi() const;
-
- double getPbDimensionLog2Pi() const ;
-
- double * getTmpTabOfSizePbDimension() const;
-
-
- protected :
-
- /// 1/ (2 * pi)^(d/2)
- double _Inv2PiPow ;
-
- /// 0.5 * p * log(2 * PI)
- double _halfPbDimensionLog2Pi;
-
- double _pbDimensionLog2Pi;
-
-
- // tableau de double longueur _pbDimension
- // utilisé pour ne pas avoir à faire sans arret allocation / destruction
- // utilisé par XEMnorme, getLogLikeLihoodOne
- double * __tmpTabOfSizePbDimension;
-
- bool _deleteSamples ;
-
-};
-
-
-inline double ** XEMGaussianData::getYStore(){
- return _yStore;
-}
-
-inline double XEMGaussianData::getInv2PiPow() const{
- return _Inv2PiPow;
-}
-
-inline double XEMGaussianData::getHalfPbDimensionLog2Pi()const{
- return _halfPbDimensionLog2Pi;
-}
-
-inline double XEMGaussianData::getPbDimensionLog2Pi()const{
- return _pbDimensionLog2Pi;
-}
-
-
-inline double* XEMGaussianData::getTmpTabOfSizePbDimension()const{
- return __tmpTabOfSizePbDimension;
-}
-
-
-#endif
diff --git a/lib/src/mixmod/MIXMOD/XEMGaussianDiagParameter.cpp b/lib/src/mixmod/MIXMOD/XEMGaussianDiagParameter.cpp
deleted file mode 100644
index 142d19a..0000000
--- a/lib/src/mixmod/MIXMOD/XEMGaussianDiagParameter.cpp
+++ /dev/null
@@ -1,434 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMGaussianDiagParameter.cpp description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
- ***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#include "XEMGaussianDiagParameter.h"
-#include "XEMGaussianGeneralParameter.h"
-#include "XEMGaussianEDDAParameter.h"
-#include "XEMGaussianData.h"
-#include "XEMModel.h"
-
-
-/****************/
-/* Constructors */
-/****************/
-
-//----------------------------------------------------
-//----------------------------------------------------
-XEMGaussianDiagParameter::XEMGaussianDiagParameter(){
- throw wrongConstructorType;
-}
-
-
-
-//-------------------------------------------------------------------------------------
-// constructor called by XEMModel
-//-------------------------------------------------------------------------------------
-XEMGaussianDiagParameter::XEMGaussianDiagParameter(XEMModel * iModel, XEMModelType * iModelType): XEMGaussianEDDAParameter(iModel, iModelType){
- int64_t k;
- _tabLambda = new double [_nbCluster];
- _tabShape = new XEMDiagMatrix*[_nbCluster];
- _W = new XEMDiagMatrix(_pbDimension); //Id
- for (k=0; k<_nbCluster; k++){
- _tabLambda[k] = 1.0;
- _tabShape[k] = new XEMDiagMatrix(_pbDimension); //Id
-
- _tabSigma[k] = new XEMDiagMatrix(_pbDimension); //Id
- _tabInvSigma[k] = new XEMDiagMatrix(_pbDimension); //Id
- _tabWk[k] = new XEMDiagMatrix(_pbDimension); //Id
- }
-
-}
-
-
-
-//---------------------------------------------------------------------
-// copy Constructor
-//---------------------------------------------------------------------
-XEMGaussianDiagParameter::XEMGaussianDiagParameter(const XEMGaussianDiagParameter * iParameter):XEMGaussianEDDAParameter(iParameter){
- int64_t k;
- _tabLambda = copyTab(iParameter->getTabLambda(), _nbCluster);
- _tabShape = new XEMDiagMatrix*[_nbCluster];
- _W = new XEMDiagMatrix(_pbDimension);
- (* _W) = iParameter->getW();
-
- XEMMatrix ** iTabSigma = iParameter->getTabSigma();
- XEMMatrix ** iTabInvSigma = iParameter->getTabInvSigma();
- XEMMatrix ** iTabWk = iParameter->getTabWk();
- XEMDiagMatrix ** iTabShape = iParameter->getTabShape();
-
- for (k=0; k<_nbCluster; k++){
- _tabSigma[k] = new XEMDiagMatrix(_pbDimension);
- (* _tabSigma[k]) = iTabSigma[k];
- _tabInvSigma[k] = new XEMDiagMatrix(_pbDimension);
- (* _tabInvSigma[k]) = iTabInvSigma[k];
- _tabWk[k] = new XEMDiagMatrix(_pbDimension);
- (* _tabWk[k]) = iTabWk[k];
- _tabShape[k] = new XEMDiagMatrix(_pbDimension);
- (* _tabShape[k]) = iTabShape[k];
- }
-}
-
-
-
-/**************/
-/* Destructor */
-/**************/
-XEMGaussianDiagParameter::~XEMGaussianDiagParameter(){
- int64_t k;
-
- if (_tabLambda){
- delete[] _tabLambda;
- _tabLambda = NULL;
- }
-
- if(_tabShape){
- for(k=0; k<_nbCluster; k++){
- delete _tabShape[k];
- _tabShape[k] = NULL;
- }
- delete[] _tabShape;
- _tabShape = NULL;
- }
-
- if (_tabInvSigma){
- for(k=0; k<_nbCluster; k++){
- delete _tabInvSigma[k];
- _tabInvSigma[k] = NULL;
- }
- }
-
- if (_tabSigma){
- for(k=0; k<_nbCluster; k++){
- delete _tabSigma[k];
- _tabSigma[k] = NULL;
- }
- }
-}
-
-
-
-//------------------------
-// reset to default values
-//------------------------
-void XEMGaussianDiagParameter::reset(){
- int64_t k;
- for (k=0; k<_nbCluster; k++){
- _tabLambda[k] = 1.0;
- *(_tabShape[k]) = 1.0;
- }
- XEMGaussianEDDAParameter::reset();
-}
-
-
-
-
-/*********/
-/* clone */
-/*********/
-XEMParameter * XEMGaussianDiagParameter::clone() const{
- XEMGaussianDiagParameter * newParam = new XEMGaussianDiagParameter(this);
- return(newParam);
-}
-
-
-/************/
-/* initUSER */
-/************/
-void XEMGaussianDiagParameter::initUSER(XEMParameter * iParam){
- XEMGaussianEDDAParameter::initUSER(iParam);
- updateTabInvSigmaAndDet();
-}
-
-
-
-/*******************/
-/* computeTabSigma */
-/*******************/
-void XEMGaussianDiagParameter::computeTabSigma(){
-
- /* Initialization */
- XEMGaussianData * data= (XEMGaussianData*)(_model->getData());
- double * tabNk = _model->getTabNk();
- int64_t k;
- XEMDiagMatrix * B = new XEMDiagMatrix(_pbDimension);
- XEMDiagMatrix * Bk = new XEMDiagMatrix(_pbDimension);
- double detB = 0.0; // Determinant of matrix B
- int64_t iter = 5; // Number of iterations in iterative procedure
- double detDiagW; // Determinant of diagonal matrix W
- double detDiagWk; // Determinant of diagonal matrix Wk
- double lambda = 0.0; // Volume
- double weightTotal = data->_weightTotal;
- double power = 1.0 / _pbDimension;
- double det = 0.0;
- double logDet = 0.0;
- double * W_k = new double[_pbDimension];
- double * Shape_k;
- double tmp;
- int64_t p;
-
-
- // Compute det[diag(W)]
- logDet = _W->determinant(minDeterminantDiagWValueError);
- detDiagW = powAndCheckIfNotNull(logDet,power);
-
- // Variance estimator for each of diagonal model
- switch(_modelType->_nameModel){
-
- //---------------------
- case (Gaussian_p_L_B) :
- case (Gaussian_pk_L_B) :
- //---------------------
- for (k=0; k<_nbCluster; k++){
- /*_tabLambda[k] = detDiagW / weightTotal;
- if ( _tabLambda[k]< minOverflow) {
- throw errorSigmaConditionNumber;
- }*/
- _tabSigma[k]->equalToMatrixDividedByDouble(_W, weightTotal); //_tabSigma[k] = _W/weightTtal
- }
- break;
-
-
- //---------------------
- case (Gaussian_p_L_Bk) :
- case (Gaussian_pk_L_Bk) :
- //----------------------
- for (k=0; k<_nbCluster; k++){
- det = _tabWk[k]->determinant(minDeterminantDiagWkValueError);
- detDiagWk = powAndCheckIfNotNull(det,power);
-
- _tabShape[k]->equalToMatrixDividedByDouble(_tabWk[k],detDiagWk); //_tabShape[k] = _tabW[k]/detWk
- lambda += detDiagWk;
- }
-
- for (k=0; k<_nbCluster; k++){
- _tabLambda[k] = lambda / weightTotal;
- if ( _tabLambda[k]< minOverflow)
- throw errorSigmaConditionNumber;
-
- _tabSigma[k]->equalToMatrixMultiplyByDouble(_tabShape[k],_tabLambda[k]); //tabSigma[k] = tabShape[k]*somme(lambda[k])
- }
- break;
-
-
- //---------------------
- case (Gaussian_p_Lk_Bk) :
- case (Gaussian_pk_Lk_Bk) :
- //----------------------
- for (k=0; k<_nbCluster; k++){
-
- logDet = _tabWk[k]->determinant( minDeterminantDiagWkValueError);
- detDiagWk = powAndCheckIfNotNull(logDet, power);
-
- _tabLambda[k] = detDiagWk / tabNk[k];
- if ( _tabLambda[k]< minOverflow)
- throw errorSigmaConditionNumber;
-
- _tabShape[k]->equalToMatrixDividedByDouble(_tabWk[k],detDiagWk); //_tabShape[k] = _tabW[k]/detWk
-
- _tabSigma[k]->equalToMatrixMultiplyByDouble(_tabShape[k],_tabLambda[k]); //tabSigma[k] = tabShape[k]*tabLambda[k]
-
- }
- break;
-
-
- //--------------------
- case (Gaussian_p_Lk_B) :
- case (Gaussian_pk_Lk_B) :
- //--------------------
- while (iter){
- /* Pb Overflow */
- for (k=0; k<_nbCluster; k++){
- if (_tabLambda[k] < minOverflow)
- throw errorSigmaConditionNumber;
- }
-
- /* Compute matrix B */
- (*B)=0.0;
- for (k=0; k<_nbCluster; k++){
- Bk->equalToMatrixDividedByDouble(_tabWk[k],_tabLambda[k]); //_tabShape[k] = _tabW[k]/_tabLambda[k]
- (*B) += Bk; // B->operator+=(Bk) : B=somme sur k des Bk
- }
- /* Compute det(B) */
- logDet = B->determinant(minDeterminantBValueError);
- detB = powAndCheckIfNotNull(logDet, power);
-
- /* Compute Shape[k] and Lambda[k] */
- for (k=0; k<_nbCluster; k++){
- // W_k = _tabWk[k]->getDiagonalValue();
- _tabWk[k]->putDiagonalValueInStore(W_k);
- _tabShape[k]->equalToMatrixDividedByDouble(B,detB);
- Shape_k = _tabShape[k]->getStore();
- tmp = 0.0;
- for(p=0; p<_pbDimension ; p++){
- tmp += W_k[p] / Shape_k[p];
- }
- tmp /= (_pbDimension*tabNk[k]);
-
- _tabLambda[k] = tmp;
- if ( _tabLambda[k] < minOverflow)
- throw errorSigmaConditionNumber;
- }
- iter--;
- }
-
- for (k=0; k<_nbCluster; k++){
- _tabSigma[k]->equalToMatrixMultiplyByDouble(_tabShape[k],_tabLambda[k]);
- }
-
- break;
-
-
- default :
- //------
- throw internalMixmodError;
- break;
- }
-
- updateTabInvSigmaAndDet() ;
- delete Bk;
- delete B;
- delete [] W_k;
-}
-
-
-
-
-
-//-------------------
-//getLogLikelihoodOne
-//-------------------
-double XEMGaussianDiagParameter::getLogLikelihoodOne() const{
-
- /* Compute log-likelihood for one cluster
- useful for NEC criterion */
-
- /* Initialization */
- int64_t nbSample = _model->getNbSample();
- int64_t i;
- XEMGaussianData * data = (XEMGaussianData*)(_model->getData());
- double logLikelihoodOne; // Log-likelihood for k=1
- double * Mean = new double[_pbDimension] ;
- double ** y = data->_yStore;
- double * yi;
- XEMDiagMatrix * Sigma = new XEMDiagMatrix(_pbDimension);
- XEMDiagMatrix * W = new XEMDiagMatrix(_pbDimension, 0.0);
- double norme, logDet, detDiagW ;
- double * weight = data->_weight;
- // Mean Estimator (empirical estimator)
- double totalWeight = data->_weightTotal;
- computeMeanOne(Mean,weight,y,nbSample,totalWeight);
- weight = data->_weight;
-
- /* Compute the Cluster Scattering Matrix W */
- int64_t p; // parcours
- double * xiMoinsMuk = data->getTmpTabOfSizePbDimension();
-
- for(i=0; iadd(xiMoinsMuk, weight[i]);
- }
-
- /* Compute determinant of diag(W) */
- logDet = W->detDiag(minDeterminantDiagWValueError); // virtual
- detDiagW = powAndCheckIfNotNull(logDet ,1.0/_pbDimension);
- Sigma->equalToMatrixDividedByDouble(W, totalWeight); // virtual
-
- // inverse of Sigma
- //XEMDiagMatrix * SigmaMoins1 = new XEMDiagMatrix(_pbDimension);
- XEMMatrix * SigmaMoins1 = NULL;
-// SigmaMoins1->inverse(Sigma);// virtual
- Sigma->inverse(SigmaMoins1); //construction de SigmaMoins1 dans la fonction inverse
-
- double detSigma = Sigma->determinant(minDeterminantSigmaValueError); // virtual
-
- // Compute the log-likelihood for one cluster (k=1)
- logLikelihoodOne = 0.0;
- for (i=0; inorme(xiMoinsMuk); // virtual
- logLikelihoodOne += norme * weight[i];
- }
-
- logLikelihoodOne += totalWeight * ( data->getPbDimensionLog2Pi() + log(detSigma)) ;
- logLikelihoodOne *= -0.5 ;
-
- delete W;
- delete Sigma;
- delete SigmaMoins1;
-
- delete[] Mean;
- return logLikelihoodOne;
-}
-
-
-
-//----------------
-//getFreeParameter
-//----------------
-int64_t XEMGaussianDiagParameter::getFreeParameter() const{
- int64_t nbParameter; // Number of parameters //
- int64_t k = _nbCluster; // Sample size //
- int64_t d = _pbDimension; // Sample dimension //
-
- int64_t alphaR = k*d; // alpha for for models with Restrainct proportions (Gaussian_p_...)
- int64_t alphaF = (k*d) + k-1; // alpha for models with Free proportions (Gaussian_pk_...)
-
- switch (_modelType->_nameModel) {
- case (Gaussian_p_L_B) :
- nbParameter = alphaR + d;
- break;
- case (Gaussian_p_Lk_B) :
- nbParameter = alphaR + d + k - 1;
- break;
- case (Gaussian_p_L_Bk) :
- nbParameter = alphaR + (k*d) - k + 1;
- break;
- case (Gaussian_p_Lk_Bk) :
- nbParameter = alphaR + (k*d);
- break;
- case (Gaussian_pk_L_B) :
- nbParameter = alphaF + d;
- break;
- case (Gaussian_pk_Lk_B) :
- nbParameter = alphaF + d + k - 1;
- break;
- case (Gaussian_pk_L_Bk) :
- nbParameter = alphaF + (k*d) - k + 1;
- break;
- case (Gaussian_pk_Lk_Bk) :
- nbParameter = alphaF + (k*d);
- break;
- default :
- throw internalMixmodError;
- break;
- }
- return nbParameter;
-}
-
diff --git a/lib/src/mixmod/MIXMOD/XEMGaussianDiagParameter.h b/lib/src/mixmod/MIXMOD/XEMGaussianDiagParameter.h
deleted file mode 100644
index 4f5091a..0000000
--- a/lib/src/mixmod/MIXMOD/XEMGaussianDiagParameter.h
+++ /dev/null
@@ -1,99 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMGaussianDiagParameter.h description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#ifndef XEMGaussianDiagParameter_H
-#define XEMGaussianDiagParameter_H
-
-#include "XEMGaussianEDDAParameter.h"
-#include "XEMMatrix.h"
-#include "XEMDiagMatrix.h"
-
-/**
- @brief Derived class of XEMGaussianParameter for Spherical Gaussian Model(s)
- @author F Langrognet & A Echenim
-*/
-
-class XEMGaussianDiagParameter : public XEMGaussianEDDAParameter{
-
- public:
-
- /// Default constructor
- XEMGaussianDiagParameter();
-
- /// Constructor
- // called by XEMModel
- XEMGaussianDiagParameter(XEMModel * iModel, XEMModelType * iModelType);
-
- /// Constructor (copy)
- XEMGaussianDiagParameter(const XEMGaussianDiagParameter * iParameter);
-
- /// Destructor
- virtual ~XEMGaussianDiagParameter();
-
- /// reset to default values
- virtual void reset();
-
- /** @brief Selector
- @return A copy of the model
- */
- XEMParameter * clone() const;
-
- /// initialisation USER
- void initUSER(XEMParameter * iParam);
-
- /// Compute table of sigmas of the samples of each cluster
- // NB : compute also lambda, shpae, orientation, wk, w
- void computeTabSigma();
-
-
- // SELECTORS
- // ------ / -------- //
-
- double * getTabLambda() const;
-
- //double ** getTabShape();
- XEMDiagMatrix** getTabShape() const;
-
- double getLogLikelihoodOne() const;
-
- int64_t getFreeParameter() const;
-
- protected :
- /// Table of volume of each cluster
- double * _tabLambda; /* Volume */
-
- //double ** _tabShape;
- XEMDiagMatrix** _tabShape;
-
-
-};
-
-inline double * XEMGaussianDiagParameter::getTabLambda() const{
- return _tabLambda;
-}
-
-inline XEMDiagMatrix** XEMGaussianDiagParameter::getTabShape() const{
- return _tabShape;
-}
-#endif
diff --git a/lib/src/mixmod/MIXMOD/XEMGaussianEDDAParameter.cpp b/lib/src/mixmod/MIXMOD/XEMGaussianEDDAParameter.cpp
deleted file mode 100644
index d09a41d..0000000
--- a/lib/src/mixmod/MIXMOD/XEMGaussianEDDAParameter.cpp
+++ /dev/null
@@ -1,583 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMGaussianEDDAParameter.cpp description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#include "XEMGaussianEDDAParameter.h"
-#include "XEMGaussianParameter.h"
-#include "XEMGaussianData.h"
-#include "XEMModel.h"
-
-
-
-#include "XEMGaussianGeneralParameter.h"
-#include "XEMUtil.h"
-#include "XEMGaussianData.h"
-#include "XEMRandom.h"
-
-#include "XEMMatrix.h"
-#include "XEMDiagMatrix.h"
-#include "XEMSphericalMatrix.h"
-#include "XEMSymmetricMatrix.h"
-#include "XEMGeneralMatrix.h"
-
-//-------------------- Constructors--------------
-
-XEMGaussianEDDAParameter::XEMGaussianEDDAParameter():XEMGaussianParameter(){
- throw wrongConstructorType;
-}
-
-
-//--------------
-// constructor
-//-------------
-XEMGaussianEDDAParameter::XEMGaussianEDDAParameter(XEMModel * iModel, XEMModelType * iModelType): XEMGaussianParameter(iModel, iModelType){
- _tabInvSqrtDetSigma = new double [_nbCluster];
- for (int64_t k=0; k<_nbCluster; k++){
- _tabInvSqrtDetSigma[k] = 1.0;
- }
- _tabInvSigma = new XEMMatrix* [_nbCluster];
- _tabSigma = new XEMMatrix* [_nbCluster];
-}
-
-
-
-
-//------------------------------------------
-// constructor called if USER_initialisation
-//-------------------------------------------
-XEMGaussianEDDAParameter::XEMGaussianEDDAParameter(int64_t iNbCluster, int64_t iPbDimension, XEMModelType * iModelType) : XEMGaussianParameter(iNbCluster, iPbDimension, iModelType){
- _tabInvSqrtDetSigma = new double [_nbCluster];
- for (int64_t k=0; k<_nbCluster; k++){
- _tabInvSqrtDetSigma[k] = 0.0;
- }
- _tabInvSigma = new XEMMatrix* [_nbCluster];
- _tabSigma = new XEMMatrix* [_nbCluster];
-}
-
-
-
-//-----------------
-// copy constructor
-//-----------------
-XEMGaussianEDDAParameter::XEMGaussianEDDAParameter(const XEMGaussianEDDAParameter * iParameter):XEMGaussianParameter(iParameter){
- int64_t k;
- _tabInvSqrtDetSigma = new double[_nbCluster];
- double * iTabInvSqrtDetSigma = iParameter->getTabInvSqrtDetSigma();
- for (k=0; k<_nbCluster; k++){
- _tabInvSqrtDetSigma[k] = iTabInvSqrtDetSigma[k];
- }
- _tabInvSigma = new XEMMatrix* [_nbCluster];
- _tabSigma = new XEMMatrix* [_nbCluster];
-}
-
-
-
-
-
-/**************/
-/* Destructor */
-/**************/
-XEMGaussianEDDAParameter::~XEMGaussianEDDAParameter(){
- if (_tabInvSqrtDetSigma){
- delete[] _tabInvSqrtDetSigma;
- _tabInvSqrtDetSigma = NULL;
- }
-
- if(_tabInvSigma){
- delete[] _tabInvSigma;
- _tabInvSigma = NULL;
- }
-
- if(_tabSigma){
- delete[] _tabSigma;
- _tabSigma = NULL;
- }
-
-}
-
-
-//------------------------
-// reset to default values
-//------------------------
-void XEMGaussianEDDAParameter::reset(){
- for (int64_t k=0; k<_nbCluster; k++){
- _tabInvSqrtDetSigma[k] = 0.0;
- *(_tabSigma[k]) = 1.0;
- *(_tabInvSigma[k]) = 1.0;
- }
- XEMGaussianParameter::reset();
-}
-
-
-
-
-void XEMGaussianEDDAParameter::getAllPdf(double** tabFik,double* tabProportion)const{
-
- XEMGaussianData * data = (XEMGaussianData*)(_model->getData());
- int64_t nbSample = _model->getNbSample();
- double * muk ; // to store mu_k
- double InvPiInvDetProp; // 1/(pi)^(d/2) * 1/det(\Sigma_k) ^ (.5) * p_k
- XEMMatrix * SigmakMoins1;
- double ** matrix = data->getYStore(); // to store x_i i=1,...,n
-
- double ** p_matrix; // pointeur : parcours de y
- double ** p_tabFik_i; // pointeur : parcours de tabFik
-
- double norme;
-
- // parcours
- int64_t i; // parcours des individus
- int64_t k; // cluster index
-
- double * xiMoinsMuk = data->getTmpTabOfSizePbDimension();
- double * xi;
- int64_t p;
-
- for(k=0 ; k<_nbCluster ; k++){
- InvPiInvDetProp = data->getInv2PiPow() * _tabInvSqrtDetSigma[k] * tabProportion[k];
- muk = _tabMean[k];
- SigmakMoins1 = _tabInvSigma[k];
-
- //computing
- p_matrix = matrix; // pointe sur la premiere ligne des donnees
- p_tabFik_i = tabFik;
- for(i=0 ; inorme(xiMoinsMuk); // virtual !!
- (*p_tabFik_i)[k] = InvPiInvDetProp * exp(-0.5 * norme);
-
- p_matrix++;
- p_tabFik_i++;
- } // end for i
- } // end for k
-
-}
-
-
-//-------
-// getPdf
-//-------
-double XEMGaussianEDDAParameter::getPdf(XEMSample * x, int64_t kCluster)const{
-
- XEMGaussianData * data = (XEMGaussianData*)(_model->getData());
-
- double normPdf = 0.0; /* Normal pdf */
- double invPi = data->getInv2PiPow() ; /* Inverse of (2*Pi)^(d/2) */
-
- /* Compute (x-m)*inv(S)*(x-m)' */
- double * ligne = ((XEMGaussianSample*)x)->getTabValue();
-
- XEMMatrix * Sigma_kMoins1 = _tabInvSigma[kCluster];
-
- double norme;
- double * muk = _tabMean[kCluster];
- double * xiMoinsMuk = data->getTmpTabOfSizePbDimension();
- int64_t p;
-
- for(p=0; p<_pbDimension; p++){
- xiMoinsMuk[p] = ligne[p] - muk[p];
- }
-
- // calcul de la norme (x-muk)'Sigma_k^{-1} (x-muk)
- norme = Sigma_kMoins1->norme(xiMoinsMuk); // virtual
-
- /* Compute normal density */
- normPdf = invPi * _tabInvSqrtDetSigma[kCluster] * exp( -0.5 * norme);
-
- return normPdf;
-}
-
-
-
-
-
-double XEMGaussianEDDAParameter::getPdf(int64_t iSample, int64_t kCluster)const{
-
- XEMGaussianData * data = (XEMGaussianData*)(_model->getData());
- double * x = (data->getYStore())[iSample];
-
- XEMMatrix * Sigma_kCluster_moins1 = _tabInvSigma[kCluster];
-
- double norme;
-
- double * xiMoinsMuk = data->getTmpTabOfSizePbDimension();
- double * muk = _tabMean[kCluster];
-
- int64_t p;
-
- for(p=0; p<_pbDimension; p++){
- xiMoinsMuk[p] = x[p] - muk[p];
- }
-
- // calcul de la norme (x-muk)'Sigma_k^{-1} (x-muk)
- norme = Sigma_kCluster_moins1->norme(xiMoinsMuk); // virtual
-
- double normPdf;
-
- // Compute normal density
- normPdf = data->getInv2PiPow() * _tabInvSqrtDetSigma[kCluster] * exp(-0.5*norme);
-
- return normPdf;
-}
-
-
-
-
-void XEMGaussianEDDAParameter::updateTabInvSigmaAndDet(){
- int64_t k;
- double detSigma;
- for (k=0; k<_nbCluster; k++){
- detSigma = _tabSigma[k]->determinant(minDeterminantSigmaValueError);
- _tabSigma[k]->inverse(_tabInvSigma[k]);
- _tabInvSqrtDetSigma[k] = 1.0 / sqrt(detSigma);
- }
-
-}
-
-
-//---------
-// MAP step
-//---------
-// update _tabWk and _W.
-// No need for this algo but if CV is the criterion, these values
-// must be updated
-/*void XEMGaussianEDDAParameter::MAPStep(){
- computeTabWkW();
- }*/
-
-
-
-//--------------------
-// init USER_PARTITION
-//--------------------
-void XEMGaussianEDDAParameter::initForInitUSER_PARTITION(int64_t & nbInitializedCluster, bool * tabNotInitializedCluster, XEMPartition * initPartition){
- // init : tabSigma, _tabInvSigma, _tabInvSqrtDetSigma,
- XEMDiagMatrix * matrixDataVar = new XEMDiagMatrix(_pbDimension,0.0);
- computeGlobalDiagDataVariance(matrixDataVar);
- // Vaiance matrix initialization to variance matrix of data (Symmetric, Diag or Spherical)
- for (int64_t k=0; k<_nbCluster; k++){
- (*_tabSigma[k]) = matrixDataVar;
- }
- updateTabInvSigmaAndDet();
- delete matrixDataVar;
-
- // _tabMean
- computeTabMeanInitUSER_PARTITION(nbInitializedCluster, tabNotInitializedCluster, initPartition);
-}
-
-
-
-
-
-void XEMGaussianEDDAParameter::initUSER(XEMParameter * iParam){
- /*
- updated : _tabMean, _tabWk, _tabSigma, _tabProportion, _tabShape, _tabOrientation, _tabInvSigma, _tabInvSqrtDetSigma
- */
-
- // recuperation
- // we got an XEMGaussianGeneralParameter
- // because of the implementation in class XEMStrategyType
- // in gaussian cases the init parameters are allways General
- XEMGaussianGeneralParameter * param = (XEMGaussianGeneralParameter *) iParam;
- double ** iTabMean = param->getTabMean();
- double * iTabProportion = param->getTabProportion();
- XEMMatrix ** iTabWk = param->getTabWk();
- XEMMatrix ** iTabSigma = param->getTabSigma();
- int64_t k;
- for (k=0; k<_nbCluster; k++){
- recopyTab(iTabMean[k], _tabMean[k], _pbDimension);
- (* _tabWk[k]) = iTabWk[k];
- (* _tabSigma[k]) = iTabSigma[k];
- if (hasFreeProportion(_modelType->_nameModel)) {
- _tabProportion[k] = iTabProportion[k];
- }
- else{
- _tabProportion[k] = 1.0 / _nbCluster;
- }
- }
-}
-
-
-
-void XEMGaussianEDDAParameter::input(ifstream & fi){
- // verify paramFile
- double tmp;
- int64_t cpt=0;
- fi.clear();
- fi.seekg(0, ios::beg);
- //while (!fi.eof()){
- while (fi>>tmp){ // retourne NULL si on est �la fin
- cpt++;
- }
- int64_t goodValue = _nbCluster * (1 + _pbDimension + _pbDimension*_pbDimension);
-
- if(cpt != goodValue){
- throw errorInitParameter;
- }
-
- fi.clear();
- fi.seekg(0, ios::beg);
-
- int64_t j,k;
- double * muk;
- for (k=0; k<_nbCluster; k++){
- muk = _tabMean[k];
- // Proportions //
- fi >> _tabProportion[k];
- // Center (mean) //
- for (j=0; j<_pbDimension; j++){
- fi >> muk[j];
- }
- // Variance matrix //
- _tabSigma[k]->input(fi) ; // virtual method
- } // end for k
-}
-
-
-
-
-
-void XEMGaussianEDDAParameter::updateForCV(XEMModel * originalModel, XEMCVBlock & CVBlock){
- XEMGaussianParameter::updateForCV(originalModel,CVBlock);
- computeTabSigma();
-}
-
-
-
-/*-------------------------------------------------------------------------------------------
- M step
- ------
-
- already updated :
- - _model (_tabFik, _tabTik, _tabCik, _tabNk) if Estep is done before
- - _model->_tabCik, _model->_tabNk if USER_PARTITION is done before (Disciminant analysis)
- In all cases, only _tabCik and _tabNk are needed
-
- updated in this method :
- - in XEMGaussianParameter::MStep :
- - _tabProportion
- - _tabMean
- - _tabWk
- - _W
- - here :
- - _tabSigma
- - _tabInvSigma
- - _tabInvSqrtDetSigma
- --------------------------------------------------------------------------------------------*/
-void XEMGaussianEDDAParameter::MStep(){
- XEMGaussianParameter::MStep();
- computeTabSigma();
-}
-
-
-
-
-
-
-//--------------------------------------------
-// initialize attributes before an InitRandom
-//--------------------------------------------
-void XEMGaussianEDDAParameter::initForInitRANDOM(){
- XEMDiagMatrix * matrixDataVar = new XEMDiagMatrix(_pbDimension,0.0);
- computeGlobalDiagDataVariance(matrixDataVar);
-
- // Vaiance matrix initialization to variance matrix of data (Symmetric, Diag or Spherical)
- // init : tabSigma, _tabInvSigma, _tabInvSqrtDetSigma,
- for (int64_t k=0; k<_nbCluster; k++){
- (*_tabSigma[k]) = matrixDataVar;
- }
- updateTabInvSigmaAndDet();
-
- delete matrixDataVar;
-}
-
-
-
-
-
-//-----------------------------------------
-// computeTik when underflow
-// -model->_tabSumF[i] pour ith sample = 0
-// i : 0 ->_nbSample-1
-//-----------------------------------------
-void XEMGaussianEDDAParameter::computeTikUnderflow(int64_t i, double ** tabTik){
-
- XEMGaussianData * data = (XEMGaussianData*)(_model->getData());
- int64_t * lnFk = new int64_t [_nbCluster];
- long double * lnFkPrim = new long double[_nbCluster];
- long double * fkPrim = new long double[_nbCluster];
- int64_t k,k0;
- long double lnFkMax, fkTPrim;
- long double detSigma;
-
- double * ligne = (data->getYStore())[i];
- long double norme;
- double ** p_tabMean;
- double * tabTik_i = tabTik[i];
-
- double * xiMoinsMuk = data->getTmpTabOfSizePbDimension();
- double * muk;
-
- p_tabMean = _tabMean;
- int64_t p;
- for (k=0; k<_nbCluster; k++){
- detSigma = (_tabSigma[k])->determinant(minDeterminantSigmaValueError); // virtual method
- muk = *p_tabMean;
- for(p=0; p<_pbDimension; p++){
- xiMoinsMuk[p] = ligne[p] - muk[p];
- }
- norme = (_tabInvSigma[k])->norme(xiMoinsMuk); // virtual method
-
- // compute lnFik[k]
-
- lnFk[k] = log(_tabProportion[k]) - data->getHalfPbDimensionLog2Pi() - 0.5*log(detSigma) - 0.5*norme;
-
- p_tabMean++;
- } // end for k
-
- lnFkMax = lnFk[0];
- for (k=1; k<_nbCluster; k++){
- if (lnFk[k] > lnFkMax){
- lnFkMax = lnFk[k];
- }
- }
-
-
- fkTPrim = 0.0;
- for (k=0; k<_nbCluster; k++){
- lnFkPrim[k] = lnFk[k] - lnFkMax;
- fkPrim[k] = exp(lnFkPrim[k]);
- fkTPrim += fkPrim[k];
- }
-
- // compute tabTik
- if (fkTPrim == 0){
- // reset tabTik
- initToZero(tabTik_i, _nbCluster);
- k0 = XEMGaussianParameter::computeClassAssigment(i);
- tabTik_i[k0] = 1.0;
- }
- else{
- for (k=0; k<_nbCluster; k++){
- tabTik_i[k] = fkPrim[k] / fkTPrim;
- }
- }
-
- delete [] lnFk;
- delete [] lnFkPrim;
- delete [] fkPrim;
-
-}
-
-
-
-
-void XEMGaussianEDDAParameter::edit(){
- int64_t k;
- for (k=0; k<_nbCluster; k++){
- cout<<"\tcomponent : " << k << endl;
- cout<<"\t\tproportion : " << _tabProportion[k]<edit(cout, "\t\t\t");
-
- cout<<"\t\tWk : "<edit(cout,"\t\t\t");
-
- cout<<"\t\tinvSigma : "<edit(cout,"\t\t\t");
-
- cout<<"\t\ttabInvSqrtDetSigma : "<<_tabInvSqrtDetSigma[k]<edit(cout, "\t\t");
-
-}
-
-
-
-
-
-
-//------
-// Edit
-//-------
-void XEMGaussianEDDAParameter::edit(ofstream & oFile, bool text){
- int64_t k;
- if (text){
- for (k=0; k<_nbCluster; k++){
- oFile<<"\t\t\tComponent "<edit(oFile, "\t\t\t\t");
-
- oFile<edit(oFile, "");
- oFile<getTabMean(), _tabMean, _nbCluster, _pbDimension);
- // _W->recopy(iParameter->getW());
- (*_W) = iParameter->getW();
-
- int64_t k;
- XEMMatrix ** iTabSigma = iParameter->getTabSigma();
- XEMMatrix ** iTabInvSigma = iParameter->getTabInvSigma();
- XEMMatrix ** iTabWk = iParameter->getTabWk();
- for(k=0; k<_nbCluster; k++){
- // _tabSigma[k]->recopy(iTabSigma[k]);
- (* _tabSigma[k]) = iTabSigma[k];
- // _tabInvSigma[k]->recopy(iTabInvSigma[k]);
- (*_tabInvSigma[k]) = iTabInvSigma[k];
- // _tabWk[k]->recopy(iTabWk[k]);
- (* _tabWk[k]) = iTabWk[k];
- }
- recopyTab(iParameter->getTabInvSqrtDetSigma(), _tabInvSqrtDetSigma, _nbCluster);
-
-}
-
diff --git a/lib/src/mixmod/MIXMOD/XEMGaussianEDDAParameter.h b/lib/src/mixmod/MIXMOD/XEMGaussianEDDAParameter.h
deleted file mode 100644
index d467896..0000000
--- a/lib/src/mixmod/MIXMOD/XEMGaussianEDDAParameter.h
+++ /dev/null
@@ -1,156 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMGaussianEDDAParameter.h description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#ifndef XEMGaussianEDDAParameter_H
-#define XEMGaussianEDDAParameter_H
-
-#include "XEMGaussianParameter.h"
-#include "XEMMatrix.h"
-#include "XEMGeneralMatrix.h"
-#include "XEMOldInput.h"
-
-class XEMGaussianEDDAParameter : public XEMGaussianParameter{
-
- public:
-
- /// Default constructor
- XEMGaussianEDDAParameter();
-
- /// Constructor
- // called by XEMModel
- XEMGaussianEDDAParameter(XEMModel * iModel, XEMModelType * iModelType);
-
-
- /// Constructor
- // called if USER initialisation
- XEMGaussianEDDAParameter(int64_t iNbCluster, int64_t iPbDimension, XEMModelType * iModelType);
-
-
- /// Constructor
- XEMGaussianEDDAParameter(const XEMGaussianEDDAParameter * iParameter);
-
- /// Destructor
- virtual ~XEMGaussianEDDAParameter();
-
- /// reset to default values
- virtual void reset();
-
-
- /** @brief Selector
- @return Table of inverse of sqrt of determinant of covariance matrix for each cluster
- */
- double * getTabInvSqrtDetSigma() const;
-
- /** @brief Selector
- @return Table of inverse of covariance matrix for each cluster
- */
- XEMMatrix ** getTabInvSigma() const;
-
-
- /** @brief Selector
- @return Table of covariance matrix for each cluster
- */
-
- XEMMatrix ** getTabSigma() const;
-
- /// Compute normal probability density function
- /// for iSample the sample and kCluster th cluster
- double getPdf(int64_t iSample, int64_t kCluster) const;
-
-
- /// compute normal probability density function
- /// for all i=1,..,n and k=1,..,K
- void getAllPdf(double ** tabFik,double * tabProportion) const;
-
-
- /// compute normal probability density function
- /// for the line x within the kCluster cluster
- double getPdf(XEMSample * x, int64_t kCluster) const;
-
- void updateTabInvSigmaAndDet() ;
-
- void computeTikUnderflow(int64_t i, double ** tabTik);
-
- void edit();
-
- void edit(ofstream & oFile, bool text=false);
-
- void recopy(XEMParameter * otherParameter);
-
- virtual int64_t getFreeParameter() const= 0;
- virtual void computeTabSigma() = 0;
-
- void updateForCV(XEMModel * originalModel, XEMCVBlock & CVBlock) ;
-
-
- virtual XEMParameter* clone() const = 0 ;
- void MStep();
- void MAPStep();
- virtual void input(ifstream & fi);
-
-
- //init
- //----
-
- /// User initialisation of the parameters of the model
- virtual void initUSER(XEMParameter * iParam);
-
-
- /// initialize attributes before an initRANDOM
- void initForInitRANDOM();
-
-
- /// initialize attributes for init USER_PARTITION
- /// outputs :
- /// - nbInitializedCluster
- /// - tabNotInitializedCluster (array of size _nbCluster)
- void initForInitUSER_PARTITION(int64_t & nbInitializedCluster, bool * tabNotInitializedCluster, XEMPartition * initPartition);
-
- protected :
- /// Table of inverse of covariance matrix of each cluster
- XEMMatrix ** _tabInvSigma;
-
- /// Table of covariance Matrix of each cluster
- XEMMatrix ** _tabSigma;
-
- /// 1/det(Sigma)
- double * _tabInvSqrtDetSigma;
-
-
-};
-
-
-inline double * XEMGaussianEDDAParameter::getTabInvSqrtDetSigma() const{
- return _tabInvSqrtDetSigma;
-}
-
-inline XEMMatrix ** XEMGaussianEDDAParameter::getTabInvSigma() const{
- return _tabInvSigma;
-}
-
-inline XEMMatrix ** XEMGaussianEDDAParameter::getTabSigma() const{
- return _tabSigma;
-}
-
-#endif
diff --git a/lib/src/mixmod/MIXMOD/XEMGaussianGeneralParameter.cpp b/lib/src/mixmod/MIXMOD/XEMGaussianGeneralParameter.cpp
deleted file mode 100644
index 3e44186..0000000
--- a/lib/src/mixmod/MIXMOD/XEMGaussianGeneralParameter.cpp
+++ /dev/null
@@ -1,981 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMGaussianGeneralParameter.cpp description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#include "XEMGaussianGeneralParameter.h"
-#include "XEMGaussianEDDAParameter.h"
-#include "XEMGaussianData.h"
-#include "XEMModel.h"
-#include "XEMDiagMatrix.h"
-#include "XEMSymmetricMatrix.h"
-#include "XEMGeneralMatrix.h"
-
-/****************/
-/* Constructors */
-/****************/
-
-//----------------------------------------------------
-//----------------------------------------------------
-XEMGaussianGeneralParameter::XEMGaussianGeneralParameter(){
- throw wrongConstructorType;
-}
-
-
-
-//-------------------------------------------------------------------------------------
-// constructor called by XEMModel
-//-------------------------------------------------------------------------------------
-XEMGaussianGeneralParameter::XEMGaussianGeneralParameter(XEMModel * iModel, XEMModelType * iModelType): XEMGaussianEDDAParameter(iModel, iModelType){
- int64_t k;
- _tabShape = new XEMDiagMatrix*[_nbCluster];
- _tabOrientation = new XEMGeneralMatrix*[_nbCluster];
- _tabLambda = new double [_nbCluster];
- _W = new XEMSymmetricMatrix(_pbDimension); //Id
-
- for (k=0; k<_nbCluster; k++){
- _tabShape[k] = new XEMDiagMatrix(_pbDimension); // Id
- _tabOrientation[k] = new XEMGeneralMatrix(_pbDimension);//Id
- _tabLambda[k] = 1.0;
- _tabSigma[k] = new XEMSymmetricMatrix(_pbDimension); //Id
- _tabInvSigma[k] = new XEMSymmetricMatrix(_pbDimension); //Id
- _tabWk[k] = new XEMSymmetricMatrix(_pbDimension); //Id
- }
-
- __storeDim = _pbDimension * (_pbDimension + 1) / 2;
-}
-
-
-
-
-//-------------------------------------------------------------------------------------
-//constructeur avec une initialisation USER
-//-------------------------------------------------------------------------------------
-XEMGaussianGeneralParameter::XEMGaussianGeneralParameter(int64_t iNbCluster, int64_t iPbDimension, XEMModelType * iModelType, string & iFileName) : XEMGaussianEDDAParameter(iNbCluster, iPbDimension, iModelType){
- int64_t k;
- __storeDim = _pbDimension * (_pbDimension + 1) / 2;
- _tabShape = new XEMDiagMatrix*[_nbCluster];
- _tabOrientation = new XEMGeneralMatrix*[_nbCluster];
- _tabLambda = new double [_nbCluster];
-
- for (k=0; k<_nbCluster; k++){
- _tabShape[k] = new XEMDiagMatrix(_pbDimension); //Id
- _tabOrientation[k] = new XEMGeneralMatrix(_pbDimension); //Id
- _tabLambda[k] = 1.0;
-
- // _tabSigma, _tabInvSigma, _tabWk will be initialized in XEMGaussianEDDAparameter
- _tabInvSigma[k] = new XEMSymmetricMatrix(_pbDimension); //Id
- _tabSigma[k] = new XEMSymmetricMatrix(_pbDimension); // Id
- _tabWk[k] = new XEMSymmetricMatrix(_pbDimension); //Id
- *_tabWk[k] = 1.0;
- }
- _W = new XEMSymmetricMatrix(_pbDimension); //Id
-
- // read parameters in file iFileName//
- if (iFileName.compare("") != 0){
- ifstream paramFile(iFileName.c_str(), ios::in);
- if (! paramFile.is_open()){
- throw wrongParamFileName;
- }
- input(paramFile);
- paramFile.close();
- }
-
- updateTabInvSigmaAndDet(); // method of XEMGaussianParameter
-
-}
-
-
-
-//-------------------------------------------------------------------------------------
-//constructeur par copie
-//-------------------------------------------------------------------------------------
-XEMGaussianGeneralParameter::XEMGaussianGeneralParameter(const XEMGaussianGeneralParameter * iParameter):XEMGaussianEDDAParameter(iParameter){
- int64_t k;
- __storeDim = _pbDimension * (_pbDimension + 1) / 2;
- _tabShape = new XEMDiagMatrix*[_nbCluster];
- _tabOrientation = new XEMGeneralMatrix*[_nbCluster];
- _tabLambda = new double[_nbCluster];
- XEMDiagMatrix ** iTabShape = iParameter->getTabShape();
- XEMGeneralMatrix ** iTabOrientation = iParameter->getTabOrientation();
- double * iTabLambda = iParameter->getTabLambda();
- XEMMatrix ** iTabSigma = iParameter->getTabSigma();
- XEMMatrix ** iTabInvSigma = iParameter->getTabInvSigma();
- XEMMatrix ** iTabWk = iParameter->getTabWk();
- _W = new XEMSymmetricMatrix((XEMSymmetricMatrix *)(iParameter->getW()));
-
- for (k=0; k<_nbCluster; k++){
- //_tabInvSigma[k] = NULL;
- _tabShape[k] = new XEMDiagMatrix(iTabShape[k]); // copy constructor
- _tabOrientation[k] = new XEMGeneralMatrix(iTabOrientation[k]); // copy constructor
- _tabLambda[k] = iTabLambda[k];
-
- _tabSigma[k] = new XEMSymmetricMatrix(_pbDimension);
- (* _tabSigma[k]) = iTabSigma[k];
- _tabWk[k] = new XEMSymmetricMatrix(_pbDimension);
- (* _tabWk[k]) = iTabWk[k];
- _tabInvSigma[k] = new XEMSymmetricMatrix(_pbDimension);
- (* _tabInvSigma[k]) = iTabInvSigma[k];
- }
-}
-
-
-
-/**************/
-/* Destructor */
-/**************/
-XEMGaussianGeneralParameter::~XEMGaussianGeneralParameter(){
- int64_t k;
-
- if (_tabShape){
- for (k=0;k<_nbCluster;k++){
- delete _tabShape[k];
- _tabShape[k] = NULL;
- }
- delete[] _tabShape;
- _tabShape = NULL;
- }
-
- if (_tabOrientation){
- for (k=0;k<_nbCluster;k++){
- delete _tabOrientation[k];
- _tabOrientation[k] = NULL;
- }
- delete[] _tabOrientation;
- _tabOrientation = NULL;
- }
-
- if (_tabLambda){
- delete[] _tabLambda;
- _tabLambda = NULL;
- }
-
- if (_tabInvSigma){
- for(k=0; k<_nbCluster; k++){
- delete _tabInvSigma[k];
- _tabInvSigma[k] = NULL;
- }
- }
-
- if (_tabSigma){
- //cout<<"destructeur gaussgene"<getTabShape();
- XEMGeneralMatrix ** iTabOrientation = param->getTabOrientation();
- double * iTabLambda = param->getTabLambda();
- int64_t k;
- for (k=0; k<_nbCluster; k++){
- (*_tabShape[k]) = (iTabShape[k]); // copy constructor
- (* _tabOrientation[k]) = iTabOrientation[k]; // copy constructor
- _tabLambda[k] = iTabLambda[k];
- }
-
-}
-
-
-
-/***********************/
-/* computeTabSigma_L_C */
-/***********************/
-void XEMGaussianGeneralParameter::computeTabSigma_L_C(){
- int64_t k;
- double totalWeight = ((XEMGaussianData*)(_model->getData()))->_weightTotal;
- for(k=0; k<_nbCluster; k++){
- _tabSigma[k]->equalToMatrixDividedByDouble(_W, totalWeight); // :: Sigma_k = W / totalWeight
- }
-}
-
-
-
-/*************************/
-/* computeTabSigma_Lkk_Ck */
-/*************************/
-void XEMGaussianGeneralParameter::computeTabSigma_Lk_Ck(){
- int64_t k;
- double * tabNk = _model->getTabNk();
- for (k=0; k<_nbCluster; k++){
- _tabSigma[k]->equalToMatrixDividedByDouble(_tabWk[k], tabNk[k]); // :: Sigma_k = W_k / n_k
- }
-}
-
-
-
-/************************/
-/* computeTabSigma_L_Ck */
-/************************/
-void XEMGaussianGeneralParameter::computeTabSigma_L_Ck(){
- double lambda = 0.0;
- int64_t k;
- double logDet;
- XEMGaussianData * data = (XEMGaussianData*)(_model->getData());
- double totalWeight = data->_weightTotal;
- double detWk_k,tmp;
- double * detWk = new double[_nbCluster];
-
- try {
- for (k=0; k<_nbCluster; k++){
- logDet = _tabWk[k]->determinant(minDeterminantWkValueError );
- detWk_k = powAndCheckIfNotNull(logDet ,1.0/_pbDimension);
- lambda += detWk_k;
- detWk[k] = detWk_k;
- }
-
- lambda /= totalWeight; // lambda = somme sur k detWk^(1/pbDimension) / totalWeight
- if (lambda < minOverflow){
- throw errorSigmaConditionNumber;
- }
-
- for (k=0; k<_nbCluster; k++){
- _tabLambda[k] = lambda;
- tmp = detWk[k] / lambda;
- _tabSigma[k]->equalToMatrixDividedByDouble(_tabWk[k], tmp); // :: Sigma_k = Wk / tmp
- }
-
- delete [] detWk;
- detWk = NULL;
- }
- catch (XEMErrorType errorType) {
- delete [] detWk;
- detWk = NULL;
- throw errorType;
- }
-}
-
-
-
-/*****************************/
-/* computeTabSigma_L_Dk_A_Dk */
-/*****************************/
-void XEMGaussianGeneralParameter::computeTabSigma_L_Dk_A_Dk(){
-
- //Ma modification
- int64_t k;
- XEMGaussianData * data = (XEMGaussianData*)(_model->getData());
- XEMDiagMatrix* Shape = new XEMDiagMatrix(_pbDimension, 0.0); //0.0
- double logDet;
- double detShape;
- double totalWeight = data->_weightTotal;
- double * tabShape_k_store;
- // SVD Decomposition of Cluster Scattering Matrix (Symmetric) Wk: Wk=U*S*U'
-
-
-
- try{
-
- for (k=0; k<_nbCluster; k++){
- _tabWk[k]->computeSVD(_tabShape[k], _tabOrientation[k]);
- tabShape_k_store = _tabShape[k]->getStore();
- (*Shape) += _tabShape[k];
- }
-
- // Compute determinant of Shape matrix
- logDet = Shape->determinant(minDeterminantShapeValueError);
- detShape = powAndCheckIfNotNull(logDet , 1.0/_pbDimension);
-
- for (k=0; k<_nbCluster; k++){
- _tabLambda[k] = detShape / totalWeight; //determinant de la somme sur k de _tabShape[k] / totalWeight
- if (_tabLambda[k]equalToMatrixDividedByDouble(Shape,detShape); // _tabShape[k] = somme sur k de _tabShape / determinant de la somme sur k de _tabShape[k]
- _tabSigma[k]->compute_as__multi_O_S_O(_tabLambda[k], _tabOrientation[k], _tabShape[k]) ;
- }
-
- delete Shape;
- }
- catch(XEMErrorType errorType){
- delete Shape;
- throw errorType;
- }
-}
-
-
-
-/********************************/
-/* computeTabSigma_Lk_Dk_A_Dk() */
-/********************************/
-void XEMGaussianGeneralParameter::computeTabSigma_Lk_Dk_A_Dk(){
- int64_t k;
- int64_t iter = 5;
-
- XEMDiagMatrix* Shape = new XEMDiagMatrix(_pbDimension); //Id
- XEMDiagMatrix* S = new XEMDiagMatrix(_pbDimension);//Id
- double * tabNk = _model->getTabNk();
-
- XEMDiagMatrix ** tabS = new XEMDiagMatrix*[_nbCluster];
- XEMGeneralMatrix ** tabU = new XEMGeneralMatrix*[_nbCluster];
- for(k=0; k<_nbCluster ; k++){
- tabS[k] = new XEMDiagMatrix(_pbDimension); //Id
- tabU[k] = new XEMGeneralMatrix(_pbDimension); //Id
- }
-
-
- double logDet;
- double detShape;
-
- try{
- // SVD Decomposition of Cluster Scattering Matrix (Symmetric) Wk: Wk=U*S*U'
- for (k=0; k<_nbCluster; k++){
- _tabWk[k] ->computeSVD(tabS[k],tabU[k]);
- }
- while (iter){
-
- (*Shape)=0.0;
- for (k=0; k<_nbCluster; k++){
- (S)->equalToMatrixDividedByDouble(tabS[k],_tabLambda[k]);
- (*Shape) += S;
- }
-
-
- //Compute determinant of Shape matrix
- logDet = Shape->determinant(minDeterminantShapeValueError);
- detShape = powAndCheckIfNotNull(logDet,1.0/_pbDimension);
-
- for (k=0; k<_nbCluster; k++){
- _tabShape[k]->equalToMatrixDividedByDouble(Shape,detShape); //A=_tabShape = somme sur k 1/lambdak*Shape
- _tabLambda[k] = _tabWk[k]->trace_this_O_Sm1_O(tabU[k], _tabShape[k]);
- _tabLambda[k] /= (_pbDimension*tabNk[k]);
-
- if (_tabLambda[k]trace_this_O_Sm1_O(_tabOrientation[k], _tabShape[k]); //trace(Wk Dk A^(-1) Dk') / pbDimension*tabNk
- _tabLambda[k] /= (_pbDimension*tabNk[k]);
-
- if (_tabLambda[k]compute_as__multi_O_S_O(_tabLambda[k], _tabOrientation[k], _tabShape[k]);
- }
-
-
- for(k=0; k<_nbCluster; k++){
- delete tabS[k];
- tabS[k] = NULL;
- delete tabU[k];
- tabU[k] = NULL;
- }
- delete [] tabU;
- delete [] tabS;
- delete Shape;
- delete S;
-
- }
-
- catch(XEMErrorType errorType){
- for(k=0; k<_nbCluster; k++){
- delete tabS[k];
- tabS[k] = NULL;
- delete tabU[k];
- tabU[k] = NULL;
- }
- delete [] tabU;
- delete [] tabS;
- delete Shape;
- delete S;
- throw errorType;
- }
-}
-
-
-
-
-
-
-/************************/
-/* computeTabSigma_Lk_C */
-/************************/
-void XEMGaussianGeneralParameter::computeTabSigma_Lk_C(){
-
- int64_t k;
- double * tabNk = _model->getTabNk();
- XEMSymmetricMatrix * C = new XEMSymmetricMatrix(_pbDimension); //Id
- XEMMatrix * R = new XEMSymmetricMatrix(_pbDimension);
- XEMMatrix * Inv;
- Inv = new XEMSymmetricMatrix(_pbDimension);
- double logDet, detR;
- int64_t iter = 5;
-
- try {
-
- while (iter){
- // Inv = NULL;
- (*R)=0.0;
- for (k=0; k<_nbCluster; k++){
- R->compute_product_Lk_Wk(_tabWk[k], _tabLambda[k]); // on calcule Wk/Lk
- }
-
- logDet = R->determinant(minDeterminantRValueError);
- detR = powAndCheckIfNotNull(logDet,1.0/_pbDimension);
- C->equalToMatrixDividedByDouble(R,detR); // C = somme sur k de 1/lambdak * Wk / det (somme sur k de 1/lambdak * Wk )^(1/d)
-
- C->inverse(Inv);
- for (k=0; k<_nbCluster; k++){
- _tabLambda[k] = _tabWk[k]->compute_trace_W_C(Inv); //lambdak = trace(Wk C^(-1)) / pbDimension * tabNk
- _tabLambda[k] /= (_pbDimension*tabNk[k]);
-
- if (_tabLambda[k]getData());
- int64_t k;
- int64_t iter = 5;
- double diff, FOld, F;
- double sumTraceM = 0.0;
- double lambda, logDet, DetDiagQtmp;
- XEMDiagMatrix * Shape_0 = new XEMDiagMatrix(_pbDimension);//Id
-
- /* Flury algorithm */
- /* SVD Decomposition of Cluster Scattering Matrix (Symmetric) Wk: Wk=U*S*U' */
- /* Orientation matrix initialisation: _tabOrientation */
- try{
- (*Shape_0) = (_tabShape[0]);
- _tabWk[0]-> computeSVD(_tabShape[0], _tabOrientation[0]);
- (*_tabShape[0]) = ( Shape_0);
-
- /* Compute Lambda: _lambda */
- /* Mtmp = D*Ak^-1*D'*Wk */
-
-
- for (k=0; k<_nbCluster; k++){
- sumTraceM += _tabWk[k]->trace_this_O_Sm1_O(_tabOrientation[0], _tabShape[k]);
- }
- lambda = sumTraceM/(_pbDimension*data->_weightTotal); // lambda = somme sur k de trace(D Ak^(-1) D' Wk) / weightTotal
-
- diff = 10.0;
- FOld = 0.0;
- F = 0.0;
- /* Iterative procedure 1 */
- while ((iter) && (diff > defaultFluryEpsilon)){
- /* Compute Shape matrix: _tabShape */
- /* Mtmp = D*Wk^-1*D' */
-
- F = 0.0;
-
- for (k=0; k<_nbCluster; k++){
- F += _tabWk[k]->trace_this_O_Sm1_O(_tabOrientation[0], _tabShape[k]);
-
- _tabWk[k]->computeShape_as__diag_Ot_this_O(_tabShape[k], _tabOrientation[0]);
- logDet = _tabShape[k]->determinant(minDeterminantDiagQtmpValueError);
- DetDiagQtmp = powAndCheckIfNotNull(logDet ,1.0/_pbDimension);
- (*_tabShape[k]) /= DetDiagQtmp; // tabShpaek = Ak = diag(D' Wk D) / det(diag(D' Wk D))^(1/d)
- }
-
- FOld = F;
- F = flury(F);
- diff = fabs(F-FOld);
- iter--;
- }
-
- if (lambda < minOverflow){
- throw errorSigmaConditionNumber;
- }
-
- for (k=0; k<_nbCluster; k++){
- _tabLambda[k] = lambda;
- _tabSigma[k]->compute_as__multi_O_S_O(lambda, _tabOrientation[0], _tabShape[k]);
- }
- delete Shape_0;
-
- }
- catch(XEMErrorType errorType){
- delete Shape_0;
- Shape_0 = NULL;
- throw errorType;
- }
-
-}
-
-
-
-
-/*****************************/
-/* computeTabSigma_Lk_D_Ak_D */
-/*****************************/
-void XEMGaussianGeneralParameter::computeTabSigma_Lk_D_Ak_D(){
-
- double * tabNk = _model->getTabNk();
- int64_t k, iter = 5;
-
- /* Flury algorithm */
- /* SVD Decomposition of Cluster Scattering Matrix (Symmetric) Wk: Wk=U*S*U' */
- //((XEMSymmetricMatrix *) _tabWk[0]) -> computeSVD(_tabShape[0], _tabOrientation[0]);
-
- _tabWk[0] -> computeSVD(_tabShape[0], _tabOrientation[0]);
-
- double diff = 10.0;
- double F = 0.0;
- double FOld, logDet;
-
- /* Iterative procedure 1 */
- while ((iter) && (diff > defaultFluryEpsilon)){
- /* Compute Shape matrix: _tabShape[k] */
- /* Qtmp = D*Wk^-1*D' */
- for (k=0; k<_nbCluster; k++){
- _tabWk[k]->computeShape_as__diag_Ot_this_O(_tabShape[k], _tabOrientation[0], tabNk[k]); //
- logDet = (_tabShape[k])->determinant(minDeterminantShapeValueError);
- }
- FOld = F;
- F = flury(F);
- diff = fabs(F-FOld);
- iter--;
- }
-
- for (k=0; k<_nbCluster; k++){
-
- (*_tabOrientation[k]) = _tabOrientation[0];
- //((XEMSymmetricMatrix *)_tabSigma[k])->compute_as__multi_O_S_O(1.0, _tabOrientation[k], _tabShape[k]);
- _tabSigma[k]->compute_as__multi_O_S_O(1.0, _tabOrientation[k], _tabShape[k]);
- }
-
-}
-
-
-
-
-/*******************/
-/* computeTabSigma */
-/*******************/
-//------------------------------------------------/
-/* Variance estimator for each of general model */
-//------------------------------------------------/
-void XEMGaussianGeneralParameter::computeTabSigma(){
- switch (_modelType->_nameModel) {
-
- case (Gaussian_p_L_C) :
- case (Gaussian_pk_L_C) :
- computeTabSigma_L_C();
- break;
- case (Gaussian_p_Lk_Ck) :
- case (Gaussian_pk_Lk_Ck) :
- computeTabSigma_Lk_Ck();
- break;
- case (Gaussian_p_L_Ck) :
- case (Gaussian_pk_L_Ck) :
- computeTabSigma_L_Ck();
- break;
- case (Gaussian_p_L_Dk_A_Dk) :
- case (Gaussian_pk_L_Dk_A_Dk) :
- computeTabSigma_L_Dk_A_Dk();
- break;
- case (Gaussian_p_Lk_Dk_A_Dk) :
- case (Gaussian_pk_Lk_Dk_A_Dk) :
- computeTabSigma_Lk_Dk_A_Dk();
- break;
- case (Gaussian_p_Lk_C) :
- case (Gaussian_pk_Lk_C) :
- computeTabSigma_Lk_C();
- break;
- case (Gaussian_p_L_D_Ak_D) :
- case (Gaussian_pk_L_D_Ak_D) :
- computeTabSigma_L_D_Ak_D();
- break;
- case (Gaussian_p_Lk_D_Ak_D) :
- case (Gaussian_pk_Lk_D_Ak_D) :
- computeTabSigma_Lk_D_Ak_D();
- break;
- default :
- throw internalMixmodError;
- break;
- }
-
- updateTabInvSigmaAndDet();
-}
-
-
-
-//----------------
-// Flury Algorithm
-//----------------
-double XEMGaussianGeneralParameter::flury(double F){
-
- // coefficients de la matrice a� diagonaliser
- // [ a b ]
- // M = [ ]
- // [ b c ]
- double a, b, c;
- // plus petite valeur propre
- double eigenValueMoins;
- // vecteur propre associe a� eigenValueMoins
- double eigenVectorMoins_1;
- double eigenVectorMoins_2;
- int64_t k,il,im, iter=0;
- double diff = 10;
- double FOld, tmp, tmp2, termesHorsDiag,termesDiag;
- double * Wk_store;
- int64_t p,q,r;
- double * tabShape_k_store;
- double * Ori = _tabOrientation[0]->getStore();
- XEMGaussianData * data = (XEMGaussianData*)(_model->getData());
- double * D_im = data->getTmpTabOfSizePbDimension();
- double * D_il = new double[_pbDimension];
- int64_t il_p, im_p;
- double iSim_iSil_k;
-
- while ((iter < maxFluryIter) && ( diff > defaultFluryEpsilon)){
-
- /* Compute orientation matrix: _tabOrientation */
- for (il=0; il<_pbDimension; il++){
- for (im=il+1; im<_pbDimension; im++){
-
- /* Extraction of colums dl and dm */
- il_p = il;
- im_p = im;
- for(p=0;p<_pbDimension;p++){
- D_il[p] = Ori[il_p]; // remplir la colonne 1 :
- D_im[p] = Ori[im_p]; // remplir la colonne 2 :
-
- il_p += _pbDimension;
- im_p += _pbDimension;
- }
-
- /* Compute matrix R where q1 is the second eigen vector */
- //RR_store = RR.Store();
- //initToZero(RR_store,4);
-
- a = 0.0;
- b = 0.0;
- c = 0.0;
-
- for (k=0; k<_nbCluster; k++){
-
- //Wk_store = ((XEMSymmetricMatrix *)_tabWk[k])->getStore();
- Wk_store = _tabWk[k]->getSymmetricStore();
- tabShape_k_store = _tabShape[k]->getStore();
- iSim_iSil_k = 1.0 / tabShape_k_store[il] - 1.0 / tabShape_k_store[im] ;
-
- // remplir RR(1,1)
- termesHorsDiag = termesDiag = 0.0;
- for(p=0,r=0; p<_pbDimension ; p++,r++){
- tmp = D_il[p] ; // Dlm(p+1,1);
- for(q=0; qtrace_this_O_Sm1_O(_tabOrientation[0], _tabShape[k]);
- F += _tabWk[k]->trace_this_O_Sm1_O(_tabOrientation[0], _tabShape[k]);
- }
-
- diff = fabs(F-FOld);
-
- iter++;
- }
-
- delete[] D_il;
- return F;
-}
-
-
-
-//-------------------
-//getLogLikelihoodOne
-//-------------------
-double XEMGaussianGeneralParameter::getLogLikelihoodOne() const{
- /* Compute log-likelihood for one cluster
- useful for NEC criterion */
-
- /* Initialization */
- int64_t nbSample = _model->getNbSample();
- int64_t i;
- XEMGaussianData * data = (XEMGaussianData*)(_model->getData());
- double logLikelihoodOne; // Log-likelihood for k=1
- double * Mean = new double[_pbDimension] ;
- double ** y = data->_yStore;
- double * yi;
- XEMSymmetricMatrix * Sigma = new XEMSymmetricMatrix(_pbDimension); //Id
- XEMSymmetricMatrix * W = new XEMSymmetricMatrix(_pbDimension, 0.0); // 0.0
- double norme, logDet, detDiagW ;
- double * weight = data->_weight;
-
- // Mean Estimator (empirical estimator)
- double totalWeight = data->_weightTotal;
- computeMeanOne(Mean,weight,y,nbSample,totalWeight);
- weight = data->_weight;
-
- /* Compute the Cluster Scattering Matrix W */
-
- int64_t p; // parcours
- double * xiMoinsMuk = data->getTmpTabOfSizePbDimension();
-
- for(i=0; iadd(xiMoinsMuk, weight[i]);
- }
- /* Compute determinant of diag(W) */
- logDet = W->detDiag(minDeterminantDiagWValueError); // virtual
- detDiagW = powAndCheckIfNotNull(logDet ,1.0/_pbDimension);
-
- Sigma->equalToMatrixDividedByDouble(W, totalWeight); // virtual
-
- // inverse of Sigma
-
- //XEMSymmetricMatrix * SigmaMoins1 = new XEMSymmetricMatrix(_pbDimension);
- XEMMatrix* SigmaMoins1 = NULL;
- //SigmaMoins1->inverse(Sigma);// virtual
- Sigma->inverse(SigmaMoins1);
-
- double detSigma = Sigma->determinant(minDeterminantSigmaValueError); // virtual
-
-
- // Compute the log-likelihood for one cluster (k=1)
- logLikelihoodOne = 0.0;
- for (i=0; inorme(xiMoinsMuk); // virtual
- logLikelihoodOne += norme * weight[i];
- }
-
- logLikelihoodOne += totalWeight * ( data->getPbDimensionLog2Pi() + log(detSigma)) ;
- logLikelihoodOne *= -0.5 ;
-
- delete W;
- delete Sigma;
- delete SigmaMoins1;
-
- delete[] Mean;
-
- return logLikelihoodOne;
-}
-
-
-//----------------
-//getFreeParameter
-//----------------
-int64_t XEMGaussianGeneralParameter::getFreeParameter() const{
- int64_t nbParameter; // Number of parameters //
- int64_t k = _nbCluster; // Sample size //
- int64_t d = _pbDimension; // Sample dimension //
-
- int64_t alphaR = k*d; // alpha for for models with Restrainct proportions (Gaussian_p_...)
- int64_t alphaF = (k*d) + k-1; // alpha for models with Free proportions (Gaussian_pk_...)
- int64_t beta = d*(d+1)/2;
-
- switch (_modelType->_nameModel) {
- case (Gaussian_p_L_C):
- nbParameter = alphaR + beta;
- break;
- case (Gaussian_p_Lk_C):
- nbParameter = alphaR + beta + (k-1);
- break;
- case (Gaussian_p_L_D_Ak_D):
- nbParameter = alphaR + beta + (k-1)*(d-1);
- break;
- case (Gaussian_p_Lk_D_Ak_D):
- nbParameter = alphaR+beta + (k-1)*d;
- break;
- case (Gaussian_p_L_Dk_A_Dk):
- nbParameter = alphaR + (k*beta) - (k-1)*d;
- break;
- case (Gaussian_p_Lk_Dk_A_Dk):
- nbParameter = alphaR + (k*beta) - (k-1)*(d-1);
- break;
- case (Gaussian_p_L_Ck):
- nbParameter = alphaR + (k*beta) - (k-1);
- break;
- case (Gaussian_p_Lk_Ck):
- nbParameter = alphaR + (k*beta);
- break;
- case (Gaussian_pk_L_C):
- nbParameter = alphaF + beta;
- break;
- case (Gaussian_pk_Lk_C):
- nbParameter = alphaF + beta + (k-1);
- break;
- case (Gaussian_pk_L_D_Ak_D):
- nbParameter = alphaF + beta + (k-1)*(d-1);
- break;
- case (Gaussian_pk_Lk_D_Ak_D):
- nbParameter = alphaF + beta + (k-1)*d;
- break;
- case (Gaussian_pk_L_Dk_A_Dk):
- nbParameter = alphaF + (k*beta) - (k-1)*d;
- break;
- case (Gaussian_pk_Lk_Dk_A_Dk):
- nbParameter = alphaF + (k*beta) - (k-1)*(d-1);
- break;
- case (Gaussian_pk_L_Ck):
- nbParameter = alphaF + (k*beta) - (k-1);
- break;
- case (Gaussian_pk_Lk_Ck):
- nbParameter = alphaF + (k*beta);
- break;
- default :
- throw internalMixmodError;
- break;
- }
- return nbParameter;
-}
-
-
diff --git a/lib/src/mixmod/MIXMOD/XEMGaussianGeneralParameter.h b/lib/src/mixmod/MIXMOD/XEMGaussianGeneralParameter.h
deleted file mode 100644
index 7c0a43d..0000000
--- a/lib/src/mixmod/MIXMOD/XEMGaussianGeneralParameter.h
+++ /dev/null
@@ -1,138 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMGaussianGeneralParameter.h description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#ifndef XEMGaussianGeneralParameter_H
-#define XEMGaussianGeneralParameter_H
-
-#include "XEMGaussianEDDAParameter.h"
-#include "XEMMatrix.h"
-#include "XEMGeneralMatrix.h"
-#include "XEMSymmetricMatrix.h"
-
-/**
- @brief Derived class of XEMGaussianParameter for Spherical Gaussian Model(s)
- @author F Langrognet & A Echenim
-*/
-
-class XEMGaussianGeneralParameter : public XEMGaussianEDDAParameter{
-
- public:
-
- /// Default constructor
- XEMGaussianGeneralParameter();
-
- /// Constructor
- // called by XEMModel
- XEMGaussianGeneralParameter(XEMModel * iModel, XEMModelType * iModelType);
-
- // Constructor
- // called by XEMStrategyType if initialization is USER
- XEMGaussianGeneralParameter(int64_t iNbCluster, int64_t iPbDimension, XEMModelType * iModelType, string & iFileName);
-
- /// Constructor (copy)
- XEMGaussianGeneralParameter(const XEMGaussianGeneralParameter * iParameter);
-
-
- /// Destructor
- virtual ~XEMGaussianGeneralParameter();
-
- /// reset to default values
- virtual void reset();
-
- /** @brief Selector
- @return A copy of the model
- */
- XEMParameter * clone() const;
-
- void initUSER(XEMParameter * iParam);
-
- /// Compute table of sigmas of the samples of each cluster
- // NB : compute also lambda, shape, orientation, wk, w
- void computeTabSigma();
-
- /// Flury Algorithm
- /// return the value of Flury function
- double flury(double F);
-
-
- // SELECTORS
- // ------ / -------- //
- double * getTabLambda() const;
-
- /** @brief Selector
- @return Table of shape matrix for each cluster
- */
- XEMDiagMatrix ** getTabShape() const;
-
- /** @brief Selector
- @return Table of orientation matrix for each cluster
- */
- XEMGeneralMatrix ** getTabOrientation() const;
-
-
- double getLogLikelihoodOne() const;
-
- protected :
- /// Table of volume of each cluster
- double * _tabLambda; /* Volume */
-
-
- /// Table of shape matrix of each cluster
- XEMDiagMatrix ** _tabShape; /* Shape */
-
- // Table of orientation matrix of each cluster
- XEMGeneralMatrix ** _tabOrientation; /* Orientation */
-
- int64_t __storeDim;
-
- // model dependant methods for computing _tabSigma
- void computeTabSigma_L_C();
- void computeTabSigma_Lk_Ck();
- void computeTabSigma_L_Ck();
- void computeTabSigma_L_Dk_A_Dk();
- void computeTabSigma_Lk_Dk_A_Dk();
- void computeTabSigma_Lk_C();
- void computeTabSigma_L_D_Ak_D();
- void computeTabSigma_Lk_D_Ak_D();
-
- //void recopySymmetricMatrixInMatrix(SymmetricMatrix & sym, Matrix& mat, double facteur);
-
- int64_t getFreeParameter() const;
-
-};
-
-
-inline XEMDiagMatrix ** XEMGaussianGeneralParameter::getTabShape() const{
- return _tabShape;
-}
-
-inline XEMGeneralMatrix ** XEMGaussianGeneralParameter::getTabOrientation() const{
- return _tabOrientation;
-}
-
-inline double * XEMGaussianGeneralParameter::getTabLambda() const{
- return _tabLambda;
-}
-
-#endif
diff --git a/lib/src/mixmod/MIXMOD/XEMGaussianHDDAParameter.cpp b/lib/src/mixmod/MIXMOD/XEMGaussianHDDAParameter.cpp
deleted file mode 100644
index 8ff7a12..0000000
--- a/lib/src/mixmod/MIXMOD/XEMGaussianHDDAParameter.cpp
+++ /dev/null
@@ -1,1677 +0,0 @@
-/***************************************************************************
- SRC/MIXMOD/XEMGaussianHDDAParameter.cpp description
- copyright : (C) MIXMOD Team - 2001-2011
- email : contact@mixmod.org
-***************************************************************************/
-
-/***************************************************************************
- This file is part of MIXMOD
-
- MIXMOD 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.
-
- MIXMOD 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 GNU General Public License
- along with MIXMOD. If not, see .
-
- All informations available on : http://www.mixmod.org
-***************************************************************************/
-#include "XEMGaussianHDDAParameter.h"
-#include "XEMGaussianParameter.h"
-#include "XEMGaussianData.h"
-#include "XEMModel.h"
-
-
-
-#include "XEMGaussianGeneralParameter.h"
-#include "XEMUtil.h"
-#include "XEMGaussianData.h"
-#include "XEMRandom.h"
-
-#include "XEMMatrix.h"
-#include "XEMDiagMatrix.h"
-#include "XEMSphericalMatrix.h"
-#include "XEMGeneralMatrix.h"
-
-/****************/
-/* Constructors */
-/****************/
-
-//----------------------------------------------------
-//----------------------------------------------------
-XEMGaussianHDDAParameter::XEMGaussianHDDAParameter():XEMGaussianParameter(){
- throw wrongConstructorType;
-}
-
-
-
-//-------------------------------------------------------------------------------------
-// constructor called by XEMModel
-//-------------------------------------------------------------------------------------
-XEMGaussianHDDAParameter::XEMGaussianHDDAParameter( XEMModel * iModel, XEMModelType * iModelType): XEMGaussianParameter(iModel, iModelType){
- int64_t k;
-
- _tabAkj = new double*[_nbCluster];
- _tabBk = new double[_nbCluster];
- _tabShape = new XEMDiagMatrix*[_nbCluster];
- _tabQk = new XEMGeneralMatrix*[_nbCluster];
- _W = new XEMSymmetricMatrix(_pbDimension); //Id
- _tabDk = new int64_t [_nbCluster];
- _tabGammak = NULL;
- _Gammak = NULL;
-
-
- for (k=0; k<_nbCluster; k++){
- _tabShape[k] = new XEMDiagMatrix(_pbDimension); //Id
- _tabQk[k] = new XEMGeneralMatrix(_pbDimension); //Id
- _tabWk[k] = new XEMSymmetricMatrix(_pbDimension); // Id
- _tabDk[k] = 0;
- }
- __storeDim = _pbDimension * (_pbDimension + 1) / 2;
-
-
- if ((iModelType->_tabSubDimensionFree !=NULL) &&
- isFreeSubDimension(iModelType->_nameModel)){
- for (k=0;k<_nbCluster;k++){
- _tabDk[k] = iModelType->_tabSubDimensionFree[k];
- }
- }
- else if ((iModelType->_subDimensionEqual !=0) &&
- !isFreeSubDimension(iModelType->_nameModel)){
- for (k=0;k<_nbCluster;k++){
- _tabDk[k] = iModelType->_subDimensionEqual;
- }
- }
-
- for (k=0;k<_nbCluster;k++){
- _tabAkj[k] = new double[_tabDk[k]];
- for (int64_t j=0;j<_tabDk[k];j++){
- _tabAkj[k][j] = 1.0;
- }
- _tabBk[k] = 1.0; // /(_pbDimension - _tabDk[k]);
- }
-
-}
-
-
-
-//constructeur avec une initialisation USER
-//-------------------------------------------------------------------------------------
-//-------------------------------------------------------------------------------------
-XEMGaussianHDDAParameter::XEMGaussianHDDAParameter(int64_t iNbCluster, int64_t iPbDimension, XEMModelType * iModelType, string & iFileName) : XEMGaussianParameter(iNbCluster, iPbDimension, iModelType){
- int64_t k;
- _tabAkj = new double*[_nbCluster];
- _tabBk = new double[_nbCluster];
- _tabDk = new int64_t [_nbCluster];
- _tabGammak = NULL;
- _Gammak = NULL;
- __storeDim = _pbDimension * (_pbDimension + 1) / 2;
-
- _tabShape = new XEMDiagMatrix*[_nbCluster];
- _tabQk = new XEMGeneralMatrix*[_nbCluster];
-
-
- for (k=0; k<_nbCluster; k++){
- _tabShape[k] = new XEMDiagMatrix(_pbDimension); //Id
- _tabQk[k] = new XEMGeneralMatrix(_pbDimension); //Id
- _tabWk[k] = new XEMSymmetricMatrix(_pbDimension); //Id
- _tabAkj[k] = NULL;
- }
- _W = new XEMSymmetricMatrix(_pbDimension); //Id
-
- // read parameters in file iFileName//
- if (iFileName.compare("") != 0){
- ifstream paramFile(iFileName.c_str(), ios::in);
- if (! paramFile.is_open()){
- throw wrongParamFileName;
- }
- input(paramFile);
- paramFile.close();
- }
-
-}
-
-
-//constructeur par copie
-//-------------------------------------------------------------------------------------
-//-------------------------------------------------------------------------------------
-XEMGaussianHDDAParameter::XEMGaussianHDDAParameter(const XEMGaussianHDDAParameter * iParameter):XEMGaussianParameter(iParameter){
- int64_t k;
- __storeDim = _pbDimension * (_pbDimension + 1) / 2;
- int64_t * iTabD = iParameter->getTabD();
- double ** iTabA = iParameter->getTabA();
- double *iTabB = iParameter->getTabB();
-
- _tabShape = new XEMDiagMatrix*[_nbCluster];
- _tabQk = new XEMGeneralMatrix*[_nbCluster];
- _tabDk = new int64_t [_nbCluster];
- _tabAkj = new double*[_nbCluster];
- _tabBk = new double[_nbCluster];
-
- XEMDiagMatrix ** iTabShape = iParameter->getTabShape();
- XEMGeneralMatrix ** iTabQ = iParameter->getTabQ();
- XEMMatrix ** iTabWk = iParameter->getTabWk();
-
- _tabGammak = NULL;
- _Gammak = NULL;
- _W = new XEMSymmetricMatrix(_pbDimension);//Id
- (* _W) = iParameter->getW();
-
- recopyTab(iTabD,_tabDk,_nbCluster);
- recopyTab(iTabB,_tabBk,_nbCluster);
-
- for (k=0; k<_nbCluster; k++){
- _tabAkj[k] = new double[_tabDk[k]];
- recopyTab(iTabA[k], _tabAkj[k],_tabDk[k]);
- _tabShape[k] = new XEMDiagMatrix(iTabShape[k]);
- _tabQk[k] = new XEMGeneralMatrix(iTabQ[k]); // copy constructor
- _tabWk[k] = new XEMSymmetricMatrix(_pbDimension); //Id
- (* _tabWk[k]) = iTabWk[k];
- }
-
-}
-
-
-
-
-
-
-/**************/
-/* Destructor */
-/**************/
-XEMGaussianHDDAParameter::~XEMGaussianHDDAParameter(){
- int64_t k;
- if (_tabShape){
- for (k=0;k<_nbCluster;k++){
- delete _tabShape[k];
- _tabShape[k] = NULL;
- }
- delete[] _tabShape;
- _tabShape = NULL;
- }
-
- if (_tabQk){
- for (k=0;k<_nbCluster;k++){
- delete _tabQk[k];
- _tabQk[k] = NULL;
- }
- delete[] _tabQk;
- _tabQk = NULL;
- }
-
- if (_tabAkj){
- for (k=0;k<_nbCluster;k++){
- delete[] _tabAkj[k];
- _tabAkj[k] = NULL;
- }
- delete[] _tabAkj;
- _tabAkj = NULL;
- }
-
- if(_tabBk){
- delete[] _tabBk;
- _tabBk = NULL;
- }
-
- if (_tabDk){
- delete[] _tabDk;
- _tabDk = NULL;
- }
-
- if (_Gammak){
- for (k=0;k<_nbCluster;k++){
- delete[] _Gammak[k];
- _Gammak[k] = NULL;
- }
- delete[] _Gammak;
- _Gammak = NULL;
- }
-
- if(_tabGammak){
- for(k=0; k<_nbCluster; k++){
- delete _tabGammak[k];
- }
- delete[] _tabGammak;
- _tabGammak = NULL;
- }
-}
-
-
-
-
-//------------------------
-// reset to default values
-//------------------------
-void XEMGaussianHDDAParameter::reset(){
- throw internalMixmodError;
- // faire ensuite : XEMGaussianParameter::reset();
-}
-
-
-
-
-/*********/
-/* clone */
-/*********/
-XEMParameter * XEMGaussianHDDAParameter::clone() const{
- XEMGaussianHDDAParameter * newParam = new XEMGaussianHDDAParameter(this);
- return(newParam);
-}
-
-
-
-
-/*********/
-/* MStep */
-/********/
-void XEMGaussianHDDAParameter::MStep(){
- XEMGaussianParameter::MStep();
-
- switch(_modelType->_nameModel){
- case (Gaussian_HD_pk_AkjBkQkDk) :
- case (Gaussian_HD_p_AkjBkQkDk) :
- case (Gaussian_HD_pk_AkjBkQkD) :
- case (Gaussian_HD_p_AkjBkQkD) :
- computeAkjBkQk();
- break;
-
- case (Gaussian_HD_pk_AkBkQkDk) :
- case (Gaussian_HD_p_AkBkQkDk) :
- case (Gaussian_HD_pk_AkBkQkD) :
- case (Gaussian_HD_p_AkBkQkD) :
- computeAkBkQk();
- break;
-
- case (Gaussian_HD_p_AkjBQkD):
- case (Gaussian_HD_pk_AkjBQkD):
- computeAkjBQk();
- break;
-
- case (Gaussian_HD_p_AkBQkD):
- case (Gaussian_HD_pk_AkBQkD):
- computeAkBQk();
- break;
-
- case (Gaussian_HD_p_AjBkQkD):
- case (Gaussian_HD_pk_AjBkQkD):
- computeAjBkQk();
- break;
-
- case (Gaussian_HD_p_AjBQkD):
- case (Gaussian_HD_pk_AjBQkD):
- computeAjBQk();
- break;
- default : throw internalMixmodError;
- }
-
-}
-
-
-
-//**********/
-/* getPdf */
-/**********/
-// returns the density of a sample using the cost function property Cost = -2LL
-double XEMGaussianHDDAParameter::getPdf(int64_t iSample, int64_t kCluster)const{
-
- XEMGaussianData * data = (XEMGaussianData*)(_model->getData());
- double * ligne = (data->getYStore())[iSample];
- double normPdf = 0, K = 0;
-
- XEMParameter * parameter = _model->getParameter();
- XEMGaussianParameter* gparameter = (XEMGaussianParameter*)parameter;
- double ** tabMean = gparameter->getTabMean();
- double* tabProportion = gparameter->getTabProportion();
- double * xiMoinsMuk = new double[_pbDimension];
- double* tabShapek_store = new double[_pbDimension];
-
-
- //----------------calcul le produit Q*t(Q)---------------
-
- XEMSymmetricMatrix* Pk = new XEMSymmetricMatrix(_pbDimension); //Id
- Pk->compute_as_M_tM(_tabQk[kCluster], _tabDk[kCluster]);
-
- //---------------calcul du produit A = Qi_L^-1_t(Qi)------------------
-
- XEMSymmetricMatrix* A = new XEMSymmetricMatrix(_pbDimension); //id
-
- double sum_lambda = 0.0;
- for (int64_t j=0;j<_tabDk[kCluster];j++){
- tabShapek_store[j] = 1/_tabAkj[kCluster][j];
- sum_lambda += log(_tabAkj[kCluster][j]);
- }
- for (int64_t j=_tabDk[kCluster]; j<_pbDimension;j++){
- tabShapek_store[j] = 0.0;
- }
- A->compute_as_O_S_O(_tabQk[kCluster], tabShapek_store);
-
- double constante = sum_lambda + (_pbDimension-_tabDk[kCluster])*log(_tabBk[kCluster])-2*log(tabProportion[kCluster])+_pbDimension*log(2*XEMPI);
-
-
-
- //-----------soustraction des moyennes
- for (int64_t j=0;j<_pbDimension;j++){
- xiMoinsMuk[j] = ligne[j] - tabMean[kCluster][j];
- }
-
- //-----------produit Qt(Q)(x-muk)----------------------
-
- XEMSymmetricMatrix * Pi = new XEMSymmetricMatrix(_pbDimension); //Id
- Pi->compute_as_M_V(Pk,xiMoinsMuk);
- double * storePi = Pi->getStore();
-
- //----------------norme(muk-Pi(x))_A = t(muk-Pi(x))_A_(muk-Pi(x))-----------
- double normeA = A->norme(xiMoinsMuk);
-
- //----------------norme(x-Pi(x))------------------------------
- double norme = 0.0;
- for (int64_t i=0;i<_pbDimension;i++){
- storePi[i] += tabMean[kCluster][i];
- norme += (ligne[i] - storePi[i])*(ligne[i] - storePi[i]);
- }
-
- //----------------calcul de normPdf---------
- K = normeA + 1.0/_tabBk[kCluster]*norme + constante;
- normPdf = exp(-1/2*K);
-
- delete Pk;
- delete A;
- delete Pi;
- delete[] xiMoinsMuk;
- xiMoinsMuk = NULL;
- delete[] tabShapek_store;
- tabShapek_store = NULL;
-
- return normPdf;
-}
-
-
-
-
-
-/*************/
-/* getAllPdf */
-/*************/
-// returns the density of all the data (tabFik)
-void XEMGaussianHDDAParameter::getAllPdf(double ** tabFik,double * tabProportion)const{
- //cout<<"XEMGaussianHDDAParameter::getAllPdf"<getNbSample();
-
- for (int64_t i=0; igetTabValue();
- double normPdf = 0, K = 0;
- XEMParameter * parameter = _model->getParameter();
- XEMGaussianParameter* gparameter = (XEMGaussianParameter*)parameter;
- double ** tabMean = gparameter->getTabMean();
- double* tabProportion = gparameter->getTabProportion();
- double * xiMoinsMuk = new double[_pbDimension];
- double* tabShapek_store = new double[_pbDimension];
-
-
- //----------------calcul le produit Q*t(Q)---------------
-
- XEMSymmetricMatrix* Pk = new XEMSymmetricMatrix(_pbDimension); //Id
- Pk->compute_as_M_tM(_tabQk[kCluster],_tabDk[kCluster]);
-
- //---------------calcul du produit A = Qi_L^-1_t(Qi)------------------
-
- XEMSymmetricMatrix* A = new XEMSymmetricMatrix(_pbDimension); //Id
-
-
- double sum_lambda = 0.0;
- for (int64_t j=0;j<_tabDk[kCluster];j++){
- tabShapek_store[j] = 1/_tabAkj[kCluster][j];
- sum_lambda += log(_tabAkj[kCluster][j]);
- }
- for (int64_t j=_tabDk[kCluster];j<_pbDimension;j++){
- tabShapek_store[j] = 0.0;
- }
-
- A->compute_as_O_S_O(_tabQk[kCluster], tabShapek_store);
-
- double constante = sum_lambda + (_pbDimension-_tabDk[kCluster])*log(_tabBk[kCluster])-2*log(tabProportion[kCluster])+_pbDimension*log(2*XEMPI);
-
- //-----------soustraction des moyennes
- for (int64_t j=0;j<_pbDimension;j++){
- xiMoinsMuk[j] = ligne[j] - tabMean[kCluster][j];
- }
-
- //-----------produit Qt(Q)(x-muk)----------------------
-
- XEMSymmetricMatrix * Pi = new XEMSymmetricMatrix(_pbDimension); //Id
- Pi->compute_as_M_V(Pk,xiMoinsMuk);
- double * storePi = Pi->getStore();
-
- //----------------norme(muk-Pi(x))_A = t(muk-Pi(x))_A_(muk-Pi(x))-----------
- double normeA = A->norme(xiMoinsMuk);
-
- //----------------norme(x-Pi(x))------------------------------
- double norme = 0.0;
- for (int64_t i=0;i<_pbDimension;i++){
- storePi[i] += tabMean[kCluster][i];
- norme += (ligne[i] - storePi[i])*(ligne[i] - storePi[i]);
- }
-
- //----------------calcul de normPdf---------
- K = normeA + 1.0/_tabBk[kCluster]*norme + constante;
- normPdf = exp(-1/2*K);
-
- delete Pk;
- delete A;
- delete Pi;
- delete[] xiMoinsMuk;
- xiMoinsMuk = NULL;
- delete[] tabShapek_store;
- tabShapek_store = NULL;
-
- return normPdf;
-}
-
-
-
-
-/*********/
-/* input */
-/*********/
-// reads the input data file for HDDA models
-void XEMGaussianHDDAParameter::input(ifstream & fi){
- int64_t j,k;
- for (k=0; k<_nbCluster; k++){
- // Proportions //
- fi >> _tabProportion[k];
- // Center (mean) //
- for (j=0; j<_pbDimension; j++){
- fi >> _tabMean[k][j];
- }
-
- // Sub Dimension //
- fi >> _tabDk[k];
- if (_tabAkj[k]){
- //cout<<_tabAkj[k][0]<> _tabAkj[k][j];
- }
-
- // Parameters Bk //
- fi >> _tabBk[k];
- // Orientation matrix //
- _tabQk[k]->input(fi,_tabDk[k]) ; // virtual method
- } // end for k
-
-}
-
-
-
-
-/*****************/
-/* computeTabWkW */
-/*****************/
-//computes W_k and W but also _tabGammak and Gammak
-void XEMGaussianHDDAParameter:: computeTabWkW(){
- double* tabNk = _model->getTabNk();
- double ** tabCik = _model->getTabCik();
- int64_t nbSample = _model->getNbSample();
- XEMGaussianData * data = (XEMGaussianData*)(_model->getData());
- double * weight = data->_weight;
- int64_t i;
- double ** matrix = data->getYStore(); // to store x_i i=1,...,n
- int64_t j, k, l, dimStoreG;
-
- k = 0;
- while(k < _nbCluster){
- if(tabNk[k]<_pbDimension){
- _tabGammak = new XEMSymmetricMatrix*[_nbCluster];
- k = _nbCluster;
- }
- else{
- k++;
- }
- }
-
-
- XEMGaussianParameter::computeTabWkW();
-
- for (k=0;k<_nbCluster;k++){
- if ((tabNk[k]<_pbDimension) && (_tabDk[k] < (tabNk[k]+1))){
- l = 0;
- double test = floor(tabNk[k]);
- if (test==tabNk[k]){
- _Gammak = new double*[_nbCluster];
- int64_t nk = (int64_t ) tabNk[k];
- _tabGammak[k] = new XEMSymmetricMatrix(nk); //Id
- dimStoreG = nk * _pbDimension;
- _Gammak[k] = new double[dimStoreG];
- for (i=0;icompute_M_tM(_Gammak[k],dimStoreG);
- }
- else{
- throw tabNkNotInteger;
- } // end test
- }
- } // fin for k
-}
-
-
-
-
-
-//-------------------------------------------
-// initialize attributes before an InitRandom
-//-------------------------------------------
-void XEMGaussianHDDAParameter::initForInitRANDOM(){
- throw internalMixmodError;
-}
-
-
-
-
-//---------------------------
-// initForInitUSER_PARTITION
-//--------------------------
-void XEMGaussianHDDAParameter::initForInitUSER_PARTITION(int64_t & nbInitializedCluster, bool * tabNotInitializedCluster, XEMPartition * initPartition){
- computeTabMeanInitUSER_PARTITION(nbInitializedCluster, tabNotInitializedCluster, initPartition);
- // initialization of _tabAkj, ;... :
-
- XEMDiagMatrix * matrixDataVar = new XEMDiagMatrix(_pbDimension,0.0);
- computeGlobalDiagDataVariance(matrixDataVar);
- matrixDataVar->sortDiagMatrix();
-
- double * store = matrixDataVar->getStore();
- double sum_lambda = 0.0;
-
- for (int64_t k=0; k<_nbCluster; k++){
- *_tabQk[k] = 1.0;
- }
-
- for (int64_t j=0;j<_tabDk[0];j++){
- _tabAkj[0][j] = store[j];
- sum_lambda += _tabAkj[0][j] ;
- }
-
- double trace = matrixDataVar->computeTrace();
-
- _tabBk[0] = 1.0/(_pbDimension-_tabDk[0])*(trace-sum_lambda);
-
-
- for (int64_t k=1; k<_nbCluster; k++){
- for (int64_t j=0;j<_tabDk[k];j++){
- _tabAkj[k][j] = store[j];
- }
- _tabBk[k] = _tabBk[0];
- }
-
- if (nbInitializedCluster != _nbCluster){
- /* ca ne devrait pas arriver puisque on a v�rifi� dans XEMOldInput que :
- - la partition doit �tre 'compl�te' (y compris que toutes les classes sont repr�sent�es) quand on a l'algo M
- - que l'on autorise que M ou MAP pour les mod�les HD
- - que MAP => USER (et donc pas 'USER_PARTITION' avec MAP)
- */
- throw internalMixmodError;
- }
- delete matrixDataVar;
-}
-
-
-
-
-/**********************/
-/* initUSER_PARTITION */
-/**********************/
-/*void XEMGaussianHDDAParameter:: initUSER_PARTITION(){
- computeTabMeanInitUSER_PARTITION();
-
-
- _model->computeNk();
- computeTabWkW();
- if (_tabDk[0]==0){
- computeTabDk();
- }
-
- switch(_modelType){
- case (Gaussian_HD_pk_AkjBkQkDk) :
- case (Gaussian_HD_p_AkjBkQkDk) :
- case (Gaussian_HD_pk_AkjBkQkD) :
- case (Gaussian_HD_p_AkjBkQkD) :
- computeAkjBkQk();
- break;
-
- case (Gaussian_HD_pk_AkBkQkDk) :
- case (Gaussian_HD_p_AkBkQkDk) :
- case (Gaussian_HD_pk_AkBkQkD) :
- case (Gaussian_HD_p_AkBkQkD) :
- computeAkBkQk();
- break;
-
- case (Gaussian_HD_p_AkjBQkD):
- case (Gaussian_HD_pk_AkjBQkD):
- computeAkjBQk();
- break;
-
- case (Gaussian_HD_p_AkBQkD):
- case (Gaussian_HD_pk_AkBQkD):
- computeAkBQk();
- break;
-
- case (Gaussian_HD_p_AjBkQkD):
- case (Gaussian_HD_pk_AjBkQkD):
- computeAjBkQk();
- break;
-
- case (Gaussian_HD_p_AjBQkD):
- case (Gaussian_HD_pk_AjBQkD):
- computeAjBQk();
- break;
- }
- }*/
-
-
-
-
-
-/************/
-/* initUSER */
-/***********/
-void XEMGaussianHDDAParameter::initUSER(XEMParameter * iParam){
- XEMGaussianHDDAParameter * param = (XEMGaussianHDDAParameter *) iParam;
-
- double ** iTabMean = param->getTabMean();
- double * iTabProportion = param->getTabProportion();
- XEMMatrix ** iTabWk = param->getTabWk();
- int64_t * iTabDk = param->getTabD();
- double * iTabBk = param->getTabB();
- double ** iTabAkj = param->getTabA();
- XEMGeneralMatrix** iTabQk = param->getTabQ();
-
-
-
-
- int64_t k;
- recopyTab(iTabBk, _tabBk,_nbCluster);
- for (k=0; k<_nbCluster; k++){
- if (_tabDk[k] != iTabDk[k]) {
- throw differentSubDimensionsWithMAP;
- }
- else {
- recopyTab(iTabMean[k], _tabMean[k], _pbDimension);
- recopyTab(iTabAkj[k],_tabAkj[k],_tabDk[k]);
- // recopy _tabWk
- (*_tabWk[k]) = iTabWk[k];
-
- if (!hasFreeProportion(_modelType->_nameModel)){
- _tabProportion[k] = 1.0 / _nbCluster;
- }
- else {
- _tabProportion[k] = iTabProportion[k];
- }
- (*_tabQk[k]) = iTabQk[k];
- }
-
- }
-}
-
-
-
-
-/***********/
-/* recopy */
-/***********/
-void XEMGaussianHDDAParameter::recopy(XEMParameter * otherParameter){
-
- XEMGaussianParameter * iParameter = (XEMGaussianParameter *)otherParameter;
-
- recopyTab(iParameter->getTabMean(), _tabMean, _nbCluster, _pbDimension);
- //_W->recopy(iParameter->getW());
- (* _W) = iParameter->getW();
-
- int64_t k;
- XEMMatrix ** iTabWk = iParameter->getTabWk();
- for(k=0; k<_nbCluster; k++){
- //_tabWk[k]->recopy(iTabWk[k]);
- (* _tabWk[k]) = iTabWk[k];
- }
-}
-
-
-//-------------------
-//getLogLikelihoodOne
-//-------------------
-double XEMGaussianHDDAParameter::getLogLikelihoodOne()const{
- cout<<"XEMGaussianHDDAParameter::getLogLikelihoodOne"<getNbSample();
- int64_t i;
- XEMGaussianData * data = (XEMGaussianData*)(_model->getData());
- int64_t j, n, l;
- double logLikelihoodOne; // Log-likelihood for k=1
- double * Mean = new double[_pbDimension] ;
- double ** y = data->_yStore;
- double * yi;
- // XEMGeneralMatrix * Sigma = new XEMGeneralMatrix(_pbDimension);
- XEMSymmetricMatrix * W = new XEMSymmetricMatrix(_pbDimension);
- *W = 0.0;
- double norme, logDet, detDiagW ;
- double * weight = data->_weight;
-
- // Mean Estimator (empirical estimator)
- double totalWeight = data->_weightTotal;
- computeMeanOne(Mean,weight,y,nbSample,totalWeight);
- weight = data->_weight;
- // Compute the Cluster Scattering Matrix W
- int64_tp; // parcours
- double * xiMoinsMuk = data->getTmpTabOfSizePbDimension();
- double tmp;
- for(i=0; iadd(xiMoinsMuk, weight[i]);
- }
-
- //Compute determinant of diag(W)
- // logDet = W->detDiag(minDeterminantDiagWValueError); // virtual
- detDiagW = powAndCheckIfNotNull(logDet ,1.0/_pbDimension);
-
-
- Sigma->divise(W, totalWeight); // virtual
-
-
- // inverse of Sigma
- XEMGeneralMatrix * SigmaMoins1 = new XEMGeneralMatrix(_pbDimension);
- SigmaMoins1->inverse(Sigma);// virtual
-
- double detSigma = Sigma->determinant(minDeterminantSigmaValueError); // virtual
-
- // Compute the log-likelihood for one cluster (k=1)
- logLikelihoodOne = 0.0;
- for (i=0; inorme(xiMoinsMuk); // virtual
- logLikelihoodOne += norme * weight[i];
- }
-
- logLikelihoodOne += totalWeight * ( data->getPbDimensionLog2Pi() + log(detSigma)) ;
- logLikelihoodOne *= -0.5 ;
-
-
- delete W;
- delete Sigma;
- delete SigmaMoins1;
-
- delete[] Mean;
-
- return logLikelihoodOne;*/
- return 0.0;
-}
-
-
-
-//Calcule la log-vraismeblance avec la fonction de coût
-/*************************/
-/* computeLoglikelihoodK */
-/*************************/
-double* XEMGaussianHDDAParameter::computeLoglikelihoodK(double**K){
- int64_t i,k;
- int64_t nbSample = _model->getNbSample();
- int64_t ** tabZikKnown = _model->getTabZikKnown();
-
- double* L = new double[_nbCluster];
- for (k=0;k<_nbCluster;k++){
- L[k] = 0.0;
- }
-
- for (i=0;i_nameModel) {
- case (Gaussian_HD_pk_AkjBkQkDk):
- for (int64_t cls=0;clsedit(cout, "\t\t\t"," ",_tabDk[k]);
-
- cout<<"\t\tWk : "<edit(cout,"\t\t\t");
- }
-
- cout<<"\tW : "<edit(cout, "\t\t");
-
-}
-
-
-
-
-/********/
-/* edit */
-/********/
-void XEMGaussianHDDAParameter:: edit(ofstream & oFile, bool text){
-
- int64_t k;
-
- if (text){
- for (k=0; k<_nbCluster; k++){
- oFile<<"\t\t\tComponent "<edit(oFile, "\t\t\t\t\t"," ",_tabDk[k]);
- oFile<edit(oFile, ""," ",_tabDk[k]);
- oFile<getTabNk();
- int64_tk;
- // double* tabShapek_store = new double[_pbDimension];
- for (k=0;k<_nbCluster;k++){
- if (tabNk[k]<_pbDimension){
- int64_tnk = (int) tabNk[k];
- XEMGeneralMatrix * tabQk = new XEMGeneralMatrix(nk);
- tabQk->resetToZero();
- _tabQk[k]->resetToZero();
- //W_k = (XEMSymmetricMatrix* )(_tabGammak[k]);
- //W_k = _tabGammak[k];
- //W_k -> computeSVD(_tabShape[k], tabQk);
- _tabGammak[k] -> computeSVD(_tabShape[k], tabQk);
- _tabQk[k]->multiply(_Gammak[k], nk, tabQk);
- delete tabQk;
- delete _tabGammak[k];
- _tabGammak[k] = NULL;
- }
- else{
- //W_k = (XEMSymmetricMatrix* )(_tabWk[k]);
- //W_k = _tabWk[k];
- //W_k -> computeSVD(_tabShape[k], _tabQk[k]);
- _tabWk[k] -> computeSVD(_tabShape[k], _tabQk[k]);
- }
- double * tabShapek_store = _tabShape[k]->getStore();
-
-
-
- double max_shape = (tabShapek_store[0] - tabShapek_store[1])/tabNk[k];
- int64_ti = 2;
- double seuil = 0.1; // � d�cider
- while (i<_pbDimension){
- double diff = (tabShapek_store[i-1] - tabShapek_store[i])/tabNk[k];
-
- int64_tres = 0;
- if( diff < seuil*max_shape){
- int64_tr = i+1;
- while (r < _pbDimension){
- double diff2 = (tabShapek_store[r-1] - tabShapek_store[r])/tabNk[k];
- if (diff2getTabNk();
- int64_tnbSample = _model->getNbSample();
- int64_tk;
-
- for (k=0;k<_nbCluster;k++){
- _tabDk[k] = 1;
- }
-
- double nbFreeParameter;// = getFreeParameter();
-
- switch (_modelType){
- case Gaussian_HD_pk_AkjBkQkDk:
- case Gaussian_HD_p_AkjBkQkDk:
- case Gaussian_HD_pk_AkjBkQkD:
- case Gaussian_HD_p_AkjBkQkD:
- computeAkjBkQk();
- break;
- case Gaussian_HD_pk_AkBkQkDk:
- case Gaussian_HD_p_AkBkQkDk:
- case Gaussian_HD_pk_AkBkQkD:
- case Gaussian_HD_p_AkBkQkD:
- computeAkBkQk();
- break;
- case Gaussian_HD_pk_AkjBQkD:
- case Gaussian_HD_p_AkjBQkD:
- computeAkjBQk();
- break;
- case Gaussian_HD_pk_AkBQkD:
- case Gaussian_HD_p_AkBQkD:
- computeAkBQk();
- break;
- case Gaussian_HD_pk_AjBQkD:
- case Gaussian_HD_p_AjBQkD:
- computeAjBQk();
- break;
- case Gaussian_HD_pk_AjBkQkD:
- case Gaussian_HD_p_AjBkQkD:
- computeAjBkQk();
- break;
- } // end switch
-
- double** Cost1 = computeCost(_tabQk);
- double* Lk_min = computeLoglikelihoodK(Cost1);
- double * BIC_min = new double[_nbCluster];
- double* BIC = new double[_nbCluster];
- int64_t* tabD_min = new int[_nbCluster];
- double* Lk;
- double** Cost;
-
- for (int64_tk=0;k<_nbCluster;k++){
- nbFreeParameter = (_pbDimension+1.0-1.0/_nbCluster) + _tabDk[k]*(_pbDimension-(_tabDk[k]+1.0)/2.0)+_tabDk[k]+2.0;
- BIC_min[k] = (-2*Lk_min[k] + nbFreeParameter * log(tabNk[k])) / (tabNk[k]);
- // cout<<"d : "<<1<<" BIC_min["<getParameter();
- XEMGaussianParameter* gparameter = (XEMGaussianParameter*)parameter;
- double ** tabMean = gparameter->getTabMean();
- double* tabProportion = gparameter->getTabProportion();
- int64_t nbSample = _model->getNbSample();
- XEMGaussianData * data = (XEMGaussianData *)(_model->getData());
- double ** matrix = data->getYStore();
- double * xiMoinsMuk = new double[_pbDimension];
- XEMSymmetricMatrix* Pk = new XEMSymmetricMatrix(_pbDimension); //Id
- XEMSymmetricMatrix* A = new XEMSymmetricMatrix(_pbDimension); //Id
- XEMSymmetricMatrix * Pi = new XEMSymmetricMatrix(_pbDimension); //id
-
- for (int64_t classe=0; classe<_nbCluster; classe++){
- double* tabShapek_store = new double[_pbDimension];
- K[classe] = new double[nbSample];
- //----------------calcul le produit Q*t(Q)---------------
-
- Pk->compute_as_M_tM(tabQ[classe],_tabDk[classe]);
- //---------------calcul du produit A = Qi_L^-1_t(Qi)------------------
- double sum_lambda = 0.0;
- for (j=0;j<_tabDk[classe];j++){
- tabShapek_store[j] = 1.0/_tabAkj[classe][j];
- sum_lambda += log(_tabAkj[classe][j]);
- }
- for (j=_tabDk[classe];j<_pbDimension;j++){
- tabShapek_store[j] = 0.0;
- }
- A->compute_as_O_S_O(tabQ[classe],tabShapek_store);
- double constante = sum_lambda + (_pbDimension-_tabDk[classe])*log(_tabBk[classe])-2*log(tabProportion[classe])+_pbDimension*log(2*XEMPI);
-
- for (int64_t sample=0; sample