|
@@ -17,760 +17,760 @@ using namespace std;
|
|
|
|
|
|
GMM::GMM()
|
|
GMM::GMM()
|
|
{
|
|
{
|
|
- gaussians = 3;
|
|
|
|
- maxiter = 200;
|
|
|
|
- featsperclass = 0;
|
|
|
|
- tau = 10.0;
|
|
|
|
- srand ( time (NULL) );
|
|
|
|
- comp = false;
|
|
|
|
|
|
+ gaussians = 3;
|
|
|
|
+ maxiter = 200;
|
|
|
|
+ featsperclass = 0;
|
|
|
|
+ tau = 10.0;
|
|
|
|
+ srand ( time ( NULL ) );
|
|
|
|
+ comp = false;
|
|
}
|
|
}
|
|
|
|
|
|
-GMM::GMM(int _gaussians):gaussians(_gaussians)
|
|
|
|
|
|
+GMM::GMM ( int _gaussians ) : gaussians ( _gaussians )
|
|
{
|
|
{
|
|
- maxiter = 200;
|
|
|
|
- featsperclass = 0;
|
|
|
|
- tau = 0.0;
|
|
|
|
- srand ( time (NULL) );
|
|
|
|
- pyramid = false;
|
|
|
|
- comp = false;
|
|
|
|
|
|
+ maxiter = 200;
|
|
|
|
+ featsperclass = 0;
|
|
|
|
+ tau = 0.0;
|
|
|
|
+ srand ( time ( NULL ) );
|
|
|
|
+ pyramid = false;
|
|
|
|
+ comp = false;
|
|
}
|
|
}
|
|
|
|
|
|
-GMM::GMM(const Config *conf, int _gaussians):gaussians(_gaussians)
|
|
|
|
|
|
+GMM::GMM ( const Config *conf, int _gaussians ) : gaussians ( _gaussians )
|
|
{
|
|
{
|
|
- this->conf = conf;
|
|
|
|
-
|
|
|
|
- if(gaussians < 2)
|
|
|
|
- gaussians = conf->gI("GMM", "gaussians", 2 );
|
|
|
|
-
|
|
|
|
- maxiter = conf->gI("GMM", "maxiter", 200 );
|
|
|
|
-
|
|
|
|
- featsperclass = conf->gI("GMM", "featsperclass", 0 );
|
|
|
|
-
|
|
|
|
- tau = conf->gD("GMM", "tau", 100.0);
|
|
|
|
-
|
|
|
|
- pyramid = conf->gB("GMM", "pyramid", false);
|
|
|
|
-
|
|
|
|
- comp = false;
|
|
|
|
-
|
|
|
|
- srand ( time (NULL) );
|
|
|
|
|
|
+ this->conf = conf;
|
|
|
|
+
|
|
|
|
+ if ( gaussians < 2 )
|
|
|
|
+ gaussians = conf->gI ( "GMM", "gaussians", 2 );
|
|
|
|
+
|
|
|
|
+ maxiter = conf->gI ( "GMM", "maxiter", 200 );
|
|
|
|
+
|
|
|
|
+ featsperclass = conf->gI ( "GMM", "featsperclass", 0 );
|
|
|
|
+
|
|
|
|
+ tau = conf->gD ( "GMM", "tau", 100.0 );
|
|
|
|
+
|
|
|
|
+ pyramid = conf->gB ( "GMM", "pyramid", false );
|
|
|
|
+
|
|
|
|
+ comp = false;
|
|
|
|
+
|
|
|
|
+ srand ( time ( NULL ) );
|
|
}
|
|
}
|
|
|
|
|
|
-void GMM::computeMixture(Examples examples)
|
|
|
|
|
|
+void GMM::computeMixture ( Examples examples )
|
|
{
|
|
{
|
|
- int fsize = (int)examples.size();
|
|
|
|
- assert(fsize >= gaussians);
|
|
|
|
-
|
|
|
|
- dim = examples[0].second.vec->size();
|
|
|
|
-
|
|
|
|
- int samples = fsize;
|
|
|
|
- if(featsperclass > 0)
|
|
|
|
- {
|
|
|
|
- samples = featsperclass*gaussians;
|
|
|
|
- samples = std::min(samples, fsize);
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- // Copy data in Matrix
|
|
|
|
- VVector dataset;
|
|
|
|
-
|
|
|
|
- cout << "reduced training data for GMM from " << fsize << " features to " << samples << " features.";
|
|
|
|
- for(int i = 0; i < samples; i++)
|
|
|
|
- {
|
|
|
|
- int k = rand() % fsize;
|
|
|
|
- NICE::Vector *vec = examples[k].second.vec;
|
|
|
|
- dataset.push_back(*vec);
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- computeMixture(dataset);
|
|
|
|
|
|
+ int fsize = ( int ) examples.size();
|
|
|
|
+ assert ( fsize >= gaussians );
|
|
|
|
+
|
|
|
|
+ dim = examples[0].second.vec->size();
|
|
|
|
+
|
|
|
|
+ int samples = fsize;
|
|
|
|
+ if ( featsperclass > 0 )
|
|
|
|
+ {
|
|
|
|
+ samples = featsperclass * gaussians;
|
|
|
|
+ samples = std::min ( samples, fsize );
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ // Copy data in Matrix
|
|
|
|
+ VVector dataset;
|
|
|
|
+
|
|
|
|
+ cout << "reduced training data for GMM from " << fsize << " features to " << samples << " features.";
|
|
|
|
+ for ( int i = 0; i < samples; i++ )
|
|
|
|
+ {
|
|
|
|
+ int k = rand() % fsize;
|
|
|
|
+ NICE::Vector *vec = examples[k].second.vec;
|
|
|
|
+ dataset.push_back ( *vec );
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ computeMixture ( dataset );
|
|
}
|
|
}
|
|
|
|
|
|
-void GMM::computeMixture(const VVector &DataSet)
|
|
|
|
|
|
+void GMM::computeMixture ( const VVector &DataSet )
|
|
{
|
|
{
|
|
- // Learn the GMM model
|
|
|
|
- assert(DataSet.size() >= (uint)gaussians);
|
|
|
|
- //initEMkMeans(DataSet); // initialize the model
|
|
|
|
- srand ( time (NULL) );
|
|
|
|
-
|
|
|
|
- bool poweroftwo = false;
|
|
|
|
- int power = 1;
|
|
|
|
- while(power <= gaussians)
|
|
|
|
- {
|
|
|
|
- if(gaussians == power)
|
|
|
|
- poweroftwo = true;
|
|
|
|
- power *= 2;
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- if(poweroftwo && pyramid)
|
|
|
|
- {
|
|
|
|
- initEM(DataSet); // initialize the model
|
|
|
|
- int g = 1;
|
|
|
|
-
|
|
|
|
- while(g <= gaussians)
|
|
|
|
- {
|
|
|
|
- cout << "g = " << g << endl;
|
|
|
|
- doEM(DataSet, g);
|
|
|
|
-
|
|
|
|
- if(2*g <= gaussians)
|
|
|
|
- {
|
|
|
|
- for(int i = g; i < g*2; i++)
|
|
|
|
- {
|
|
|
|
- mu[i] = mu[i-g];
|
|
|
|
- sparse_sigma[i] = sparse_sigma[i-g];
|
|
|
|
- sparse_inv_sigma[i] = sparse_inv_sigma[i-g];
|
|
|
|
- log_det_sigma[i] = log_det_sigma[i-g];
|
|
|
|
- priors[i] = priors[i-g];
|
|
|
|
-
|
|
|
|
- double interval = 1.0;
|
|
|
|
- for(int k = 0; k < (int)mu[i].size(); k++)
|
|
|
|
- {
|
|
|
|
- interval = mu[i][k];
|
|
|
|
- interval = std::max(interval, 0.1);
|
|
|
|
- double r = (interval*((double)rand()/(double)RAND_MAX))-interval/2.0;
|
|
|
|
- mu[i][k] += r;
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
- g *= 2;
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
- else
|
|
|
|
- {
|
|
|
|
- initEMkMeans(DataSet); // initialize the model
|
|
|
|
- doEM(DataSet, gaussians);
|
|
|
|
- }
|
|
|
|
- // performs EM
|
|
|
|
|
|
+ // Learn the GMM model
|
|
|
|
+ assert ( DataSet.size() >= ( uint ) gaussians );
|
|
|
|
+ //initEMkMeans(DataSet); // initialize the model
|
|
|
|
+ srand ( time ( NULL ) );
|
|
|
|
+
|
|
|
|
+ bool poweroftwo = false;
|
|
|
|
+ int power = 1;
|
|
|
|
+ while ( power <= gaussians )
|
|
|
|
+ {
|
|
|
|
+ if ( gaussians == power )
|
|
|
|
+ poweroftwo = true;
|
|
|
|
+ power *= 2;
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ if ( poweroftwo && pyramid )
|
|
|
|
+ {
|
|
|
|
+ initEM ( DataSet ); // initialize the model
|
|
|
|
+ int g = 1;
|
|
|
|
+
|
|
|
|
+ while ( g <= gaussians )
|
|
|
|
+ {
|
|
|
|
+ cout << "g = " << g << endl;
|
|
|
|
+ doEM ( DataSet, g );
|
|
|
|
+
|
|
|
|
+ if ( 2*g <= gaussians )
|
|
|
|
+ {
|
|
|
|
+ for ( int i = g; i < g*2; i++ )
|
|
|
|
+ {
|
|
|
|
+ mu[i] = mu[i-g];
|
|
|
|
+ sparse_sigma[i] = sparse_sigma[i-g];
|
|
|
|
+ sparse_inv_sigma[i] = sparse_inv_sigma[i-g];
|
|
|
|
+ log_det_sigma[i] = log_det_sigma[i-g];
|
|
|
|
+ priors[i] = priors[i-g];
|
|
|
|
+
|
|
|
|
+ double interval = 1.0;
|
|
|
|
+ for ( int k = 0; k < ( int ) mu[i].size(); k++ )
|
|
|
|
+ {
|
|
|
|
+ interval = mu[i][k];
|
|
|
|
+ interval = std::max ( interval, 0.1 );
|
|
|
|
+ double r = ( interval * ( ( double ) rand() / ( double ) RAND_MAX ) ) - interval / 2.0;
|
|
|
|
+ mu[i][k] += r;
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ g *= 2;
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ else
|
|
|
|
+ {
|
|
|
|
+ initEMkMeans ( DataSet ); // initialize the model
|
|
|
|
+ doEM ( DataSet, gaussians );
|
|
|
|
+ }
|
|
|
|
+ // performs EM
|
|
}
|
|
}
|
|
|
|
|
|
-inline double diagDeterminant(const NICE::Vector &sparse_mat)
|
|
|
|
|
|
+inline double diagDeterminant ( const NICE::Vector &sparse_mat )
|
|
{
|
|
{
|
|
- double det = 1.0;
|
|
|
|
- for(int i = 0; i < (int)sparse_mat.size(); i++)
|
|
|
|
- {
|
|
|
|
- det *= sparse_mat[i];
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- return det;
|
|
|
|
|
|
+ double det = 1.0;
|
|
|
|
+ for ( int i = 0; i < ( int ) sparse_mat.size(); i++ )
|
|
|
|
+ {
|
|
|
|
+ det *= sparse_mat[i];
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ return det;
|
|
}
|
|
}
|
|
|
|
|
|
-inline double logdiagDeterminant(const NICE::Vector &sparse_mat)
|
|
|
|
|
|
+inline double logdiagDeterminant ( const NICE::Vector &sparse_mat )
|
|
{
|
|
{
|
|
- double det = 0.0;
|
|
|
|
- for(int i = 0; i < (int)sparse_mat.size(); i++)
|
|
|
|
- {
|
|
|
|
- det += log(sparse_mat[i]);
|
|
|
|
- }
|
|
|
|
|
|
+ double det = 0.0;
|
|
|
|
+ for ( int i = 0; i < ( int ) sparse_mat.size(); i++ )
|
|
|
|
+ {
|
|
|
|
+ det += log ( sparse_mat[i] );
|
|
|
|
+ }
|
|
|
|
|
|
- return det;
|
|
|
|
|
|
+ return det;
|
|
}
|
|
}
|
|
|
|
|
|
-inline NICE::Vector diagInverse(const NICE::Vector &sparse_mat)
|
|
|
|
|
|
+inline NICE::Vector diagInverse ( const NICE::Vector &sparse_mat )
|
|
{
|
|
{
|
|
- NICE::Vector inv(sparse_mat.size());
|
|
|
|
- for(int i = 0; i < (int)sparse_mat.size(); i++)
|
|
|
|
- {
|
|
|
|
- inv[i] = 1.0/sparse_mat[i];
|
|
|
|
- }
|
|
|
|
- return inv;
|
|
|
|
|
|
+ NICE::Vector inv ( sparse_mat.size() );
|
|
|
|
+ for ( int i = 0; i < ( int ) sparse_mat.size(); i++ )
|
|
|
|
+ {
|
|
|
|
+ inv[i] = 1.0 / sparse_mat[i];
|
|
|
|
+ }
|
|
|
|
+ return inv;
|
|
}
|
|
}
|
|
|
|
|
|
-void GMM::initEMkMeans(const VVector &DataSet)
|
|
|
|
|
|
+void GMM::initEMkMeans ( const VVector &DataSet )
|
|
{
|
|
{
|
|
- /*init GMM with k-Means*/
|
|
|
|
- KMeans k( gaussians );
|
|
|
|
- VVector means;
|
|
|
|
- vector<double> weights;
|
|
|
|
- vector<int> assignment;
|
|
|
|
- k.cluster ( DataSet, means, weights, assignment );
|
|
|
|
-
|
|
|
|
- int nData = DataSet.size();
|
|
|
|
- this->dim = DataSet[0].size();
|
|
|
|
- cdimval = dim*log(2*3.14159);
|
|
|
|
- vector<int> pop(gaussians, 0);
|
|
|
|
- priors.resize(gaussians);
|
|
|
|
- mu = VVector(gaussians, dim);
|
|
|
|
- log_det_sigma.clear();
|
|
|
|
- vector<int> allk;
|
|
|
|
-
|
|
|
|
- NICE::Vector globmean(dim);
|
|
|
|
- globmean.set(0.0);
|
|
|
|
-
|
|
|
|
- for(int n = 0; n < (int)DataSet.size(); n++) /* getting the max value for time */
|
|
|
|
- {
|
|
|
|
- globmean = globmean + DataSet[n];
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- globmean *= (1.0/nData);
|
|
|
|
-
|
|
|
|
- NICE::Vector sparse_globsigma;
|
|
|
|
-
|
|
|
|
- sparse_globsigma.resize(dim);
|
|
|
|
- sparse_globsigma.set(0.0);
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < (int)DataSet.size(); i++) // Covariances updates
|
|
|
|
- {
|
|
|
|
- // nur diagonal Elemente berechnen
|
|
|
|
- for(int d = 0; d < dim; d++)
|
|
|
|
- {
|
|
|
|
- double diff = (DataSet[i][d] - globmean[d]);
|
|
|
|
- sparse_globsigma[d] += diff*diff;
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- sparse_globsigma *= (1.0/DataSet.size());
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < gaussians; i++)
|
|
|
|
- {
|
|
|
|
- NICE::Vector _inv_sigma = diagInverse(sparse_globsigma);
|
|
|
|
-
|
|
|
|
- sparse_sigma.push_back(sparse_globsigma);
|
|
|
|
- sparse_inv_sigma.push_back(_inv_sigma);
|
|
|
|
- log_det_sigma.push_back(logdiagDeterminant(sparse_globsigma));
|
|
|
|
- mu[i] = means[i];
|
|
|
|
- //priors[i]=1.0/(double)gaussians; // set equi-probables states
|
|
|
|
- priors[i]=weights[i]; // set kMeans weights
|
|
|
|
- }
|
|
|
|
|
|
+ /*init GMM with k-Means*/
|
|
|
|
+ KMeans k ( gaussians );
|
|
|
|
+ VVector means;
|
|
|
|
+ vector<double> weights;
|
|
|
|
+ vector<int> assignment;
|
|
|
|
+ k.cluster ( DataSet, means, weights, assignment );
|
|
|
|
+
|
|
|
|
+ int nData = DataSet.size();
|
|
|
|
+ this->dim = DataSet[0].size();
|
|
|
|
+ cdimval = dim * log ( 2 * 3.14159 );
|
|
|
|
+ vector<int> pop ( gaussians, 0 );
|
|
|
|
+ priors.resize ( gaussians );
|
|
|
|
+ mu = VVector ( gaussians, dim );
|
|
|
|
+ log_det_sigma.clear();
|
|
|
|
+ vector<int> allk;
|
|
|
|
+
|
|
|
|
+ NICE::Vector globmean ( dim );
|
|
|
|
+ globmean.set ( 0.0 );
|
|
|
|
+
|
|
|
|
+ for ( int n = 0; n < ( int ) DataSet.size(); n++ ) /* getting the max value for time */
|
|
|
|
+ {
|
|
|
|
+ globmean = globmean + DataSet[n];
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ globmean *= ( 1.0 / nData );
|
|
|
|
+
|
|
|
|
+ NICE::Vector sparse_globsigma;
|
|
|
|
+
|
|
|
|
+ sparse_globsigma.resize ( dim );
|
|
|
|
+ sparse_globsigma.set ( 0.0 );
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < ( int ) DataSet.size(); i++ ) // Covariances updates
|
|
|
|
+ {
|
|
|
|
+ // nur diagonal Elemente berechnen
|
|
|
|
+ for ( int d = 0; d < dim; d++ )
|
|
|
|
+ {
|
|
|
|
+ double diff = ( DataSet[i][d] - globmean[d] );
|
|
|
|
+ sparse_globsigma[d] += diff * diff;
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ sparse_globsigma *= ( 1.0 / DataSet.size() );
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < gaussians; i++ )
|
|
|
|
+ {
|
|
|
|
+ NICE::Vector _inv_sigma = diagInverse ( sparse_globsigma );
|
|
|
|
+
|
|
|
|
+ sparse_sigma.push_back ( sparse_globsigma );
|
|
|
|
+ sparse_inv_sigma.push_back ( _inv_sigma );
|
|
|
|
+ log_det_sigma.push_back ( logdiagDeterminant ( sparse_globsigma ) );
|
|
|
|
+ mu[i] = means[i];
|
|
|
|
+ //priors[i]=1.0/(double)gaussians; // set equi-probables states
|
|
|
|
+ priors[i] = weights[i]; // set kMeans weights
|
|
|
|
+ }
|
|
}
|
|
}
|
|
|
|
|
|
-void GMM::initEM(const VVector &DataSet)
|
|
|
|
|
|
+void GMM::initEM ( const VVector &DataSet )
|
|
{
|
|
{
|
|
|
|
|
|
- /* init the GaussianMixture by using randomized meanvectors */
|
|
|
|
- int nData = DataSet.size();
|
|
|
|
- this->dim = DataSet[0].size();
|
|
|
|
- cdimval = dim*log(2*3.14159);
|
|
|
|
- vector<int> pop(gaussians, 0);
|
|
|
|
- priors.resize(gaussians);
|
|
|
|
- mu = VVector(gaussians, dim);
|
|
|
|
- log_det_sigma.clear();
|
|
|
|
- vector<int> allk;
|
|
|
|
-
|
|
|
|
- NICE::Vector globmean(dim);
|
|
|
|
- globmean.set(0.0);
|
|
|
|
-
|
|
|
|
- for(int n = 0; n < (int)DataSet.size(); n++) /* getting the max value for time */
|
|
|
|
- {
|
|
|
|
- globmean = globmean + DataSet[n];
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- globmean *= (1.0/nData);
|
|
|
|
-
|
|
|
|
- NICE::Vector sparse_globsigma;
|
|
|
|
-
|
|
|
|
- sparse_globsigma.resize(dim);
|
|
|
|
- sparse_globsigma.set(0.0);
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < (int)DataSet.size(); i++) // Covariances updates
|
|
|
|
- {
|
|
|
|
- // nur diagonal Elemente berechnen
|
|
|
|
- for(int d = 0; d < dim; d++)
|
|
|
|
- {
|
|
|
|
- double diff = (DataSet[i][d] - globmean[d]);
|
|
|
|
- sparse_globsigma[d] += diff*diff;
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- sparse_globsigma *= (1.0/DataSet.size());
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < gaussians; i++)
|
|
|
|
- {
|
|
|
|
- priors[i]=1.0/(double)gaussians; // set equi-probables states
|
|
|
|
- NICE::Vector _inv_sigma = diagInverse(sparse_globsigma);
|
|
|
|
-
|
|
|
|
- sparse_sigma.push_back(sparse_globsigma);
|
|
|
|
- sparse_inv_sigma.push_back(_inv_sigma);
|
|
|
|
- log_det_sigma.push_back(logdiagDeterminant(sparse_globsigma));
|
|
|
|
- bool newv = false;
|
|
|
|
- int k;
|
|
|
|
- while(!newv)
|
|
|
|
- {
|
|
|
|
- newv = true;
|
|
|
|
- k = rand() % nData;
|
|
|
|
-
|
|
|
|
- for(int nk = 0; nk < (int)allk.size(); nk++)
|
|
|
|
- if(allk[nk] == k)
|
|
|
|
- {
|
|
|
|
- newv = false;
|
|
|
|
- }
|
|
|
|
- if(newv)
|
|
|
|
- allk.push_back(k);
|
|
|
|
- }
|
|
|
|
- mu[i] = DataSet[k];
|
|
|
|
- }
|
|
|
|
|
|
+ /* init the GaussianMixture by using randomized meanvectors */
|
|
|
|
+ int nData = DataSet.size();
|
|
|
|
+ this->dim = DataSet[0].size();
|
|
|
|
+ cdimval = dim * log ( 2 * 3.14159 );
|
|
|
|
+ vector<int> pop ( gaussians, 0 );
|
|
|
|
+ priors.resize ( gaussians );
|
|
|
|
+ mu = VVector ( gaussians, dim );
|
|
|
|
+ log_det_sigma.clear();
|
|
|
|
+ vector<int> allk;
|
|
|
|
+
|
|
|
|
+ NICE::Vector globmean ( dim );
|
|
|
|
+ globmean.set ( 0.0 );
|
|
|
|
+
|
|
|
|
+ for ( int n = 0; n < ( int ) DataSet.size(); n++ ) /* getting the max value for time */
|
|
|
|
+ {
|
|
|
|
+ globmean = globmean + DataSet[n];
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ globmean *= ( 1.0 / nData );
|
|
|
|
+
|
|
|
|
+ NICE::Vector sparse_globsigma;
|
|
|
|
+
|
|
|
|
+ sparse_globsigma.resize ( dim );
|
|
|
|
+ sparse_globsigma.set ( 0.0 );
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < ( int ) DataSet.size(); i++ ) // Covariances updates
|
|
|
|
+ {
|
|
|
|
+ // nur diagonal Elemente berechnen
|
|
|
|
+ for ( int d = 0; d < dim; d++ )
|
|
|
|
+ {
|
|
|
|
+ double diff = ( DataSet[i][d] - globmean[d] );
|
|
|
|
+ sparse_globsigma[d] += diff * diff;
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ sparse_globsigma *= ( 1.0 / DataSet.size() );
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < gaussians; i++ )
|
|
|
|
+ {
|
|
|
|
+ priors[i] = 1.0 / ( double ) gaussians; // set equi-probables states
|
|
|
|
+ NICE::Vector _inv_sigma = diagInverse ( sparse_globsigma );
|
|
|
|
+
|
|
|
|
+ sparse_sigma.push_back ( sparse_globsigma );
|
|
|
|
+ sparse_inv_sigma.push_back ( _inv_sigma );
|
|
|
|
+ log_det_sigma.push_back ( logdiagDeterminant ( sparse_globsigma ) );
|
|
|
|
+ bool newv = false;
|
|
|
|
+ int k;
|
|
|
|
+ while ( !newv )
|
|
|
|
+ {
|
|
|
|
+ newv = true;
|
|
|
|
+ k = rand() % nData;
|
|
|
|
+
|
|
|
|
+ for ( int nk = 0; nk < ( int ) allk.size(); nk++ )
|
|
|
|
+ if ( allk[nk] == k )
|
|
|
|
+ {
|
|
|
|
+ newv = false;
|
|
|
|
+ }
|
|
|
|
+ if ( newv )
|
|
|
|
+ allk.push_back ( k );
|
|
|
|
+ }
|
|
|
|
+ mu[i] = DataSet[k];
|
|
|
|
+ }
|
|
}
|
|
}
|
|
|
|
|
|
-inline void sumRow(NICE::Matrix mat, NICE::Vector &res)
|
|
|
|
|
|
+inline void sumRow ( NICE::Matrix mat, NICE::Vector &res )
|
|
{
|
|
{
|
|
- int row = mat.rows();
|
|
|
|
- int column = mat.cols();
|
|
|
|
- res =NICE::Vector(column);
|
|
|
|
- res.set(1e-5f);
|
|
|
|
|
|
+ int row = mat.rows();
|
|
|
|
+ int column = mat.cols();
|
|
|
|
+ res = NICE::Vector ( column );
|
|
|
|
+ res.set ( 1e-5f );
|
|
//#pragma omp parallel for
|
|
//#pragma omp parallel for
|
|
- for(int i = 0; i < column; i++){
|
|
|
|
- for(int j = 0; j < row; j++){
|
|
|
|
- res[i] += mat(j, i);
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
|
|
+ for ( int i = 0; i < column; i++ ) {
|
|
|
|
+ for ( int j = 0; j < row; j++ ) {
|
|
|
|
+ res[i] += mat ( j, i );
|
|
|
|
+ }
|
|
|
|
+ }
|
|
}
|
|
}
|
|
|
|
|
|
-double GMM::logpdfState(const NICE::Vector &Vin,int state)
|
|
|
|
|
|
+double GMM::logpdfState ( const NICE::Vector &Vin, int state )
|
|
{
|
|
{
|
|
- /* get the probability density for a given state and a given vector */
|
|
|
|
- double p;
|
|
|
|
- NICE::Vector dif;
|
|
|
|
-
|
|
|
|
- dif = Vin;
|
|
|
|
- dif -= mu[state];
|
|
|
|
-
|
|
|
|
- p = 0.0;
|
|
|
|
- for(int i = 0; i < (int)dif.size(); i++)
|
|
|
|
- {
|
|
|
|
- p += dif[i] * dif[i] * sparse_inv_sigma[state][i];
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- p = -0.5*(p + cdimval + log_det_sigma[state]);
|
|
|
|
- return p;
|
|
|
|
|
|
+ /* get the probability density for a given state and a given vector */
|
|
|
|
+ double p;
|
|
|
|
+ NICE::Vector dif;
|
|
|
|
+
|
|
|
|
+ dif = Vin;
|
|
|
|
+ dif -= mu[state];
|
|
|
|
+
|
|
|
|
+ p = 0.0;
|
|
|
|
+ for ( int i = 0; i < ( int ) dif.size(); i++ )
|
|
|
|
+ {
|
|
|
|
+ p += dif[i] * dif[i] * sparse_inv_sigma[state][i];
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ p = -0.5 * ( p + cdimval + log_det_sigma[state] );
|
|
|
|
+ return p;
|
|
}
|
|
}
|
|
|
|
|
|
-int GMM::doEM(const VVector &DataSet, int nbgaussians)
|
|
|
|
|
|
+int GMM::doEM ( const VVector &DataSet, int nbgaussians )
|
|
{
|
|
{
|
|
/* perform Expectation/Maximization on the given Dataset :
|
|
/* perform Expectation/Maximization on the given Dataset :
|
|
- Matrix DataSet(nSamples,Dimensions).
|
|
|
|
- The GaussianMixture Object must be initialised before
|
|
|
|
- (see initEM or initEMkMeans methods ) */
|
|
|
|
-
|
|
|
|
|
|
+ Matrix DataSet(nSamples,Dimensions).
|
|
|
|
+ The GaussianMixture Object must be initialised before
|
|
|
|
+ (see initEM or initEMkMeans methods ) */
|
|
|
|
+
|
|
#ifdef DEBUG
|
|
#ifdef DEBUG
|
|
- cerr << "GMM::start EM" << endl;
|
|
|
|
-#endif
|
|
|
|
-
|
|
|
|
- int nData = DataSet.size();
|
|
|
|
- int iter = 0;
|
|
|
|
- double log_lik;
|
|
|
|
- double log_lik_threshold = 1e-6f;
|
|
|
|
- double log_lik_old = -1e10f;
|
|
|
|
-
|
|
|
|
- NICE::Matrix unity(dim,dim,0.0);
|
|
|
|
- for(int k = 0; k < dim; k++)
|
|
|
|
- unity(k, k)=1.0;
|
|
|
|
-
|
|
|
|
- //EM loop
|
|
|
|
- while(true)
|
|
|
|
- {
|
|
|
|
|
|
+ cerr << "GMM::start EM" << endl;
|
|
|
|
+#endif
|
|
|
|
+
|
|
|
|
+ int nData = DataSet.size();
|
|
|
|
+ int iter = 0;
|
|
|
|
+ double log_lik;
|
|
|
|
+ double log_lik_threshold = 1e-6f;
|
|
|
|
+ double log_lik_old = -1e10f;
|
|
|
|
+
|
|
|
|
+ NICE::Matrix unity ( dim, dim, 0.0 );
|
|
|
|
+ for ( int k = 0; k < dim; k++ )
|
|
|
|
+ unity ( k, k ) = 1.0;
|
|
|
|
+
|
|
|
|
+ //EM loop
|
|
|
|
+ while ( true )
|
|
|
|
+ {
|
|
#ifdef DEBUGGMM
|
|
#ifdef DEBUGGMM
|
|
- cerr << "GMM::EM: iteration: " << iter << endl;
|
|
|
|
|
|
+ cerr << "GMM::EM: iteration: " << iter << endl;
|
|
#endif
|
|
#endif
|
|
-
|
|
|
|
- vector<double> sum_p;
|
|
|
|
- sum_p.resize(nData);
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < nData; i++)
|
|
|
|
- {
|
|
|
|
- sum_p[i] = 0.0;
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- NICE::Matrix pxi(nData,gaussians);
|
|
|
|
- pxi.set(0.0);
|
|
|
|
- NICE::Matrix pix(nData,gaussians);
|
|
|
|
- pix.set(0.0);
|
|
|
|
- NICE::Vector E;
|
|
|
|
-
|
|
|
|
- iter++;
|
|
|
|
- if (iter>maxiter){
|
|
|
|
- cerr << "EM stops here. Max number of iterations (" << maxiter << ") has been reached." << endl;
|
|
|
|
- return iter;
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- double sum_log = 0.0;
|
|
|
|
-
|
|
|
|
- // computing expectation
|
|
|
|
- double maxp = -numeric_limits<double>::max();
|
|
|
|
- vector<double> maxptmp(nData,-numeric_limits<double>::max());
|
|
|
|
|
|
+
|
|
|
|
+ vector<double> sum_p;
|
|
|
|
+ sum_p.resize ( nData );
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < nData; i++ )
|
|
|
|
+ {
|
|
|
|
+ sum_p[i] = 0.0;
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ NICE::Matrix pxi ( nData, gaussians );
|
|
|
|
+ pxi.set ( 0.0 );
|
|
|
|
+ NICE::Matrix pix ( nData, gaussians );
|
|
|
|
+ pix.set ( 0.0 );
|
|
|
|
+ NICE::Vector E;
|
|
|
|
+
|
|
|
|
+ iter++;
|
|
|
|
+ if ( iter > maxiter ) {
|
|
|
|
+ cerr << "EM stops here. Max number of iterations (" << maxiter << ") has been reached." << endl;
|
|
|
|
+ return iter;
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ double sum_log = 0.0;
|
|
|
|
+
|
|
|
|
+ // computing expectation
|
|
|
|
+ double maxp = -numeric_limits<double>::max();
|
|
|
|
+ vector<double> maxptmp ( nData, -numeric_limits<double>::max() );
|
|
#pragma omp parallel for
|
|
#pragma omp parallel for
|
|
- for(int i = 0; i < nData; i++)
|
|
|
|
- {
|
|
|
|
- for(int j = 0; j < nbgaussians; j++)
|
|
|
|
- {
|
|
|
|
- double p = logpdfState(DataSet[i],j); // log(P(x|i))
|
|
|
|
- maxptmp[i] = std::max(maxptmp[i], p);
|
|
|
|
- pxi(i, j) = p;
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < nData; i++)
|
|
|
|
- {
|
|
|
|
- maxp = std::max(maxp, maxptmp[i]);
|
|
|
|
- }
|
|
|
|
|
|
+ for ( int i = 0; i < nData; i++ )
|
|
|
|
+ {
|
|
|
|
+ for ( int j = 0; j < nbgaussians; j++ )
|
|
|
|
+ {
|
|
|
|
+ double p = logpdfState ( DataSet[i], j ); // log(P(x|i))
|
|
|
|
+ maxptmp[i] = std::max ( maxptmp[i], p );
|
|
|
|
+ pxi ( i, j ) = p;
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < nData; i++ )
|
|
|
|
+ {
|
|
|
|
+ maxp = std::max ( maxp, maxptmp[i] );
|
|
|
|
+ }
|
|
|
|
|
|
#pragma omp parallel for
|
|
#pragma omp parallel for
|
|
- for(int i = 0; i < nData; i++)
|
|
|
|
- {
|
|
|
|
- sum_p[i] = 0.0;
|
|
|
|
- for(int j = 0; j < nbgaussians; j++)
|
|
|
|
- {
|
|
|
|
- double p = pxi(i, j) - maxp; // log(P(x|i))
|
|
|
|
- p = exp(p); // P(x|i)
|
|
|
|
-
|
|
|
|
- if(p < 1e-40)
|
|
|
|
- p = 1e-40;
|
|
|
|
-
|
|
|
|
- pxi(i, j) = p;
|
|
|
|
- sum_p[i] += p*priors[j];
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < nData; i++)
|
|
|
|
- {
|
|
|
|
- sum_log += log(sum_p[i]);
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
|
|
+ for ( int i = 0; i < nData; i++ )
|
|
|
|
+ {
|
|
|
|
+ sum_p[i] = 0.0;
|
|
|
|
+ for ( int j = 0; j < nbgaussians; j++ )
|
|
|
|
+ {
|
|
|
|
+ double p = pxi ( i, j ) - maxp; // log(P(x|i))
|
|
|
|
+ p = exp ( p ); // P(x|i)
|
|
|
|
+
|
|
|
|
+ if ( p < 1e-40 )
|
|
|
|
+ p = 1e-40;
|
|
|
|
+
|
|
|
|
+ pxi ( i, j ) = p;
|
|
|
|
+ sum_p[i] += p * priors[j];
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < nData; i++ )
|
|
|
|
+ {
|
|
|
|
+ sum_log += log ( sum_p[i] );
|
|
|
|
+ }
|
|
|
|
+
|
|
#pragma omp parallel for
|
|
#pragma omp parallel for
|
|
- for(int j = 0; j < nbgaussians; j++)
|
|
|
|
- {
|
|
|
|
- for(int i = 0; i < nData; i++)
|
|
|
|
- {
|
|
|
|
- pix(i, j) = (pxi(i, j)*priors[j])/sum_p[i]; // then P(i|x)
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- // here we compute the log likehood
|
|
|
|
- log_lik = sum_log/nData;
|
|
|
|
-#ifdef DEBUGGMM
|
|
|
|
- cout << "diff: " << fabs((log_lik/log_lik_old)-1) << " thresh: " << log_lik_threshold << " sum: " << sum_log << endl;
|
|
|
|
- //cout << "logold: " << log_lik_old << " lognew: " << log_lik << endl;
|
|
|
|
|
|
+ for ( int j = 0; j < nbgaussians; j++ )
|
|
|
|
+ {
|
|
|
|
+ for ( int i = 0; i < nData; i++ )
|
|
|
|
+ {
|
|
|
|
+ pix ( i, j ) = ( pxi ( i, j ) * priors[j] ) / sum_p[i]; // then P(i|x)
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ // here we compute the log likehood
|
|
|
|
+ log_lik = sum_log / nData;
|
|
|
|
+#ifdef DEBUGGMM
|
|
|
|
+ cout << "diff: " << fabs ( ( log_lik / log_lik_old ) - 1 ) << " thresh: " << log_lik_threshold << " sum: " << sum_log << endl;
|
|
|
|
+ //cout << "logold: " << log_lik_old << " lognew: " << log_lik << endl;
|
|
#endif
|
|
#endif
|
|
-
|
|
|
|
- if(fabs((log_lik/log_lik_old)-1) < log_lik_threshold )
|
|
|
|
- {
|
|
|
|
- //if log likehood hasn't move enough, the algorithm has converged, exiting the loop
|
|
|
|
- return iter;
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- log_lik_old = log_lik;
|
|
|
|
-
|
|
|
|
- // Update Step
|
|
|
|
- sumRow(pix,E);
|
|
|
|
-
|
|
|
|
|
|
+
|
|
|
|
+ if ( fabs ( ( log_lik / log_lik_old ) - 1 ) < log_lik_threshold )
|
|
|
|
+ {
|
|
|
|
+ //if log likehood hasn't move enough, the algorithm has converged, exiting the loop
|
|
|
|
+ return iter;
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ log_lik_old = log_lik;
|
|
|
|
+
|
|
|
|
+ // Update Step
|
|
|
|
+ sumRow ( pix, E );
|
|
|
|
+
|
|
#pragma omp parallel for
|
|
#pragma omp parallel for
|
|
- for(int j = 0; j < nbgaussians; j++)
|
|
|
|
- {
|
|
|
|
- priors[j] = (E[j]+tau)/(nData+tau*nbgaussians); // new priors
|
|
|
|
-
|
|
|
|
- NICE::Vector tmu(dim);
|
|
|
|
- tmu.set(0.0);
|
|
|
|
-
|
|
|
|
- NICE::Vector sparse_tmsigma(dim);
|
|
|
|
- sparse_tmsigma.set(0.0);
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < nData; i++) // Means update loop
|
|
|
|
- {
|
|
|
|
- tmu = tmu + (DataSet[i]*pix(i,j));
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- NICE::Vector tmu2 = mu[j];
|
|
|
|
- mu[j] = tmu + tau*tmu2;
|
|
|
|
- mu[j] = mu[j]*(1.0/(E[j]+tau));
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < nData; i++) // Covariances updates
|
|
|
|
- {
|
|
|
|
- // nur diagonal Elemente berechnen
|
|
|
|
- for(int d = 0; d < dim; d++)
|
|
|
|
- {
|
|
|
|
- sparse_tmsigma[d] += DataSet[i][d]*DataSet[i][d]*pix(i, j);
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- NICE::Vector sparse_tmsigma2(dim);
|
|
|
|
- for(int d = 0; d < dim; d++)
|
|
|
|
- {
|
|
|
|
- sparse_tmsigma[d] += 1e-5f;
|
|
|
|
- sparse_tmsigma2[d] = sparse_sigma[j][d]+tmu2[d]*tmu2[d];
|
|
|
|
- sparse_sigma[j][d] = (sparse_tmsigma[d]+tau*sparse_tmsigma2[d])/(E[j]+tau)-(mu[j][d]*mu[j][d]);
|
|
|
|
- }
|
|
|
|
- sparse_inv_sigma[j] = diagInverse(sparse_sigma[j]);
|
|
|
|
-
|
|
|
|
- log_det_sigma[j] = logdiagDeterminant(sparse_sigma[j]);
|
|
|
|
- }
|
|
|
|
- if(comp)
|
|
|
|
- {
|
|
|
|
- double dist = compare();
|
|
|
|
- cout << "dist for iteration " << iter << endl;
|
|
|
|
- cout << "d: " << dist << endl;
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
|
|
+ for ( int j = 0; j < nbgaussians; j++ )
|
|
|
|
+ {
|
|
|
|
+ priors[j] = ( E[j] + tau ) / ( nData + tau * nbgaussians ); // new priors
|
|
|
|
+
|
|
|
|
+ NICE::Vector tmu ( dim );
|
|
|
|
+ tmu.set ( 0.0 );
|
|
|
|
+
|
|
|
|
+ NICE::Vector sparse_tmsigma ( dim );
|
|
|
|
+ sparse_tmsigma.set ( 0.0 );
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < nData; i++ ) // Means update loop
|
|
|
|
+ {
|
|
|
|
+ tmu = tmu + ( DataSet[i] * pix ( i, j ) );
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ NICE::Vector tmu2 = mu[j];
|
|
|
|
+ mu[j] = tmu + tau * tmu2;
|
|
|
|
+ mu[j] = mu[j] * ( 1.0 / ( E[j] + tau ) );
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < nData; i++ ) // Covariances updates
|
|
|
|
+ {
|
|
|
|
+ // nur diagonal Elemente berechnen
|
|
|
|
+ for ( int d = 0; d < dim; d++ )
|
|
|
|
+ {
|
|
|
|
+ sparse_tmsigma[d] += DataSet[i][d] * DataSet[i][d] * pix ( i, j );
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ NICE::Vector sparse_tmsigma2 ( dim );
|
|
|
|
+ for ( int d = 0; d < dim; d++ )
|
|
|
|
+ {
|
|
|
|
+ sparse_tmsigma[d] += 1e-5f;
|
|
|
|
+ sparse_tmsigma2[d] = sparse_sigma[j][d] + tmu2[d] * tmu2[d];
|
|
|
|
+ sparse_sigma[j][d] = ( sparse_tmsigma[d] + tau * sparse_tmsigma2[d] ) / ( E[j] + tau ) - ( mu[j][d] * mu[j][d] );
|
|
|
|
+ }
|
|
|
|
+ sparse_inv_sigma[j] = diagInverse ( sparse_sigma[j] );
|
|
|
|
+
|
|
|
|
+ log_det_sigma[j] = logdiagDeterminant ( sparse_sigma[j] );
|
|
|
|
+ }
|
|
|
|
+ if ( comp )
|
|
|
|
+ {
|
|
|
|
+ double dist = compare();
|
|
|
|
+ cout << "dist for iteration " << iter << endl;
|
|
|
|
+ cout << "d: " << dist << endl;
|
|
|
|
+ }
|
|
|
|
+ }
|
|
#ifdef DEBUG
|
|
#ifdef DEBUG
|
|
- cerr << "GMM::finished EM after " << iter << " iterations" << endl;
|
|
|
|
|
|
+ cerr << "GMM::finished EM after " << iter << " iterations" << endl;
|
|
#endif
|
|
#endif
|
|
- return iter;
|
|
|
|
|
|
+ return iter;
|
|
}
|
|
}
|
|
|
|
|
|
-int GMM::getBestClass(const NICE::Vector &v, double *bprob)
|
|
|
|
|
|
+int GMM::getBestClass ( const NICE::Vector &v, double *bprob )
|
|
{
|
|
{
|
|
- int bestclass = 0;
|
|
|
|
- double maxprob = logpdfState(v,0)+log(priors[0]); // log(P(x|i))
|
|
|
|
-
|
|
|
|
- for(int i = 1; i < gaussians; i++)
|
|
|
|
- {
|
|
|
|
- double prob = logpdfState(v,i)+log(priors[i]); // log(P(x|i))
|
|
|
|
-
|
|
|
|
- if(prob > maxprob)
|
|
|
|
- {
|
|
|
|
- maxprob = prob;
|
|
|
|
- bestclass = i;
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- if(bprob != NULL)
|
|
|
|
- *bprob = maxprob;
|
|
|
|
-
|
|
|
|
- return bestclass;
|
|
|
|
|
|
+ int bestclass = 0;
|
|
|
|
+ double maxprob = logpdfState ( v, 0 ) + log ( priors[0] ); // log(P(x|i))
|
|
|
|
+
|
|
|
|
+ for ( int i = 1; i < gaussians; i++ )
|
|
|
|
+ {
|
|
|
|
+ double prob = logpdfState ( v, i ) + log ( priors[i] ); // log(P(x|i))
|
|
|
|
+
|
|
|
|
+ if ( prob > maxprob )
|
|
|
|
+ {
|
|
|
|
+ maxprob = prob;
|
|
|
|
+ bestclass = i;
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ if ( bprob != NULL )
|
|
|
|
+ *bprob = maxprob;
|
|
|
|
+
|
|
|
|
+ return bestclass;
|
|
}
|
|
}
|
|
-
|
|
|
|
-void GMM::getProbs(const NICE::Vector &vin, SparseVector &outprobs)
|
|
|
|
|
|
+
|
|
|
|
+void GMM::getProbs ( const NICE::Vector &vin, SparseVector &outprobs )
|
|
{
|
|
{
|
|
- outprobs.clear();
|
|
|
|
- outprobs.setDim(gaussians);
|
|
|
|
- Vector probs;
|
|
|
|
- getProbs(vin, probs);
|
|
|
|
- for(int i = 0; i < gaussians; i++)
|
|
|
|
- {
|
|
|
|
- if(probs[i] > 1e-7f )
|
|
|
|
- outprobs[i] = probs[i];
|
|
|
|
- }
|
|
|
|
|
|
+ outprobs.clear();
|
|
|
|
+ outprobs.setDim ( gaussians );
|
|
|
|
+ Vector probs;
|
|
|
|
+ getProbs ( vin, probs );
|
|
|
|
+ for ( int i = 0; i < gaussians; i++ )
|
|
|
|
+ {
|
|
|
|
+ if ( probs[i] > 1e-7f )
|
|
|
|
+ outprobs[i] = probs[i];
|
|
|
|
+ }
|
|
}
|
|
}
|
|
-
|
|
|
|
-void GMM::getProbs(const NICE::Vector &vin, Vector &outprobs)
|
|
|
|
|
|
+
|
|
|
|
+void GMM::getProbs ( const NICE::Vector &vin, Vector &outprobs )
|
|
{
|
|
{
|
|
- Vector probs;
|
|
|
|
- probs.resize(gaussians);
|
|
|
|
- probs.set(0.0);
|
|
|
|
-
|
|
|
|
- double probsum = 0.0;
|
|
|
|
- double maxp = -numeric_limits<double>::max();
|
|
|
|
- for(int i = 0; i < gaussians; i++)
|
|
|
|
- {
|
|
|
|
- probs[i] = logpdfState(vin, i) + log(priors[i]); // log(P(x|i))
|
|
|
|
- maxp = std::max(maxp, probs[i]);
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < gaussians; i++)
|
|
|
|
- {
|
|
|
|
- probs[i] -= maxp;
|
|
|
|
- probs[i] = exp(probs[i]);
|
|
|
|
- probsum += probs[i];
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- // normalize probs
|
|
|
|
|
|
+ Vector probs;
|
|
|
|
+ probs.resize ( gaussians );
|
|
|
|
+ probs.set ( 0.0 );
|
|
|
|
+
|
|
|
|
+ double probsum = 0.0;
|
|
|
|
+ double maxp = -numeric_limits<double>::max();
|
|
|
|
+ for ( int i = 0; i < gaussians; i++ )
|
|
|
|
+ {
|
|
|
|
+ probs[i] = logpdfState ( vin, i ) + log ( priors[i] ); // log(P(x|i))
|
|
|
|
+ maxp = std::max ( maxp, probs[i] );
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < gaussians; i++ )
|
|
|
|
+ {
|
|
|
|
+ probs[i] -= maxp;
|
|
|
|
+ probs[i] = exp ( probs[i] );
|
|
|
|
+ probsum += probs[i];
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ // normalize probs
|
|
#pragma omp parallel for
|
|
#pragma omp parallel for
|
|
- for(int i = 0; i < gaussians; i++)
|
|
|
|
- {
|
|
|
|
- probs[i] /= probsum;
|
|
|
|
- }
|
|
|
|
- outprobs = probs;
|
|
|
|
|
|
+ for ( int i = 0; i < gaussians; i++ )
|
|
|
|
+ {
|
|
|
|
+ probs[i] /= probsum;
|
|
|
|
+ }
|
|
|
|
+ outprobs = probs;
|
|
}
|
|
}
|
|
|
|
|
|
-void GMM::getFisher(const NICE::Vector &vin, SparseVector &outprobs)
|
|
|
|
|
|
+void GMM::getFisher ( const NICE::Vector &vin, SparseVector &outprobs )
|
|
{
|
|
{
|
|
- int size = gaussians*2*dim;
|
|
|
|
- outprobs.clear();
|
|
|
|
- outprobs.setDim(size);
|
|
|
|
-
|
|
|
|
- int counter = 0;
|
|
|
|
-
|
|
|
|
- NICE::Vector classprobs;
|
|
|
|
- classprobs.resize(gaussians);
|
|
|
|
- classprobs.set(0.0);
|
|
|
|
-
|
|
|
|
- double maxp = -numeric_limits<double>::max();
|
|
|
|
- for(int i = 0; i < gaussians; i++)
|
|
|
|
- {
|
|
|
|
- classprobs[i] = logpdfState(vin, i) + log(priors[i]); // log(P(x|i))
|
|
|
|
-
|
|
|
|
- maxp = std::max(maxp, classprobs[i]);
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < gaussians; i++)
|
|
|
|
- {
|
|
|
|
- double p = classprobs[i] - maxp;
|
|
|
|
-
|
|
|
|
- p = exp(p);
|
|
|
|
-
|
|
|
|
- for(int d = 0; d < dim; d++)
|
|
|
|
- {
|
|
|
|
- double diff = vin[d]-mu[i][d];
|
|
|
|
-
|
|
|
|
- double sigma2 = sparse_sigma[i][d]*sparse_sigma[i][d];
|
|
|
|
- double sigma3 = sigma2*sparse_sigma[i][d];
|
|
|
|
- double normmu = sqrt(priors[i]/sigma2);
|
|
|
|
- double normsig = sqrt((2.0*priors[i])/sigma2);
|
|
|
|
-
|
|
|
|
- double gradmu = (p * (diff/sigma2))/normmu;
|
|
|
|
- double gradsig = (p * ((diff*diff)/sigma3 - 1.0/sparse_sigma[i][d])) /normsig;
|
|
|
|
- if(fabs(gradmu) > 1e-7f)
|
|
|
|
- outprobs[counter] = gradmu;
|
|
|
|
- counter++;
|
|
|
|
-
|
|
|
|
- if(fabs(gradsig) > 1e-7f)
|
|
|
|
- outprobs[counter] = gradsig;
|
|
|
|
- counter++;
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
|
|
+ int size = gaussians * 2 * dim;
|
|
|
|
+ outprobs.clear();
|
|
|
|
+ outprobs.setDim ( size );
|
|
|
|
+
|
|
|
|
+ int counter = 0;
|
|
|
|
+
|
|
|
|
+ NICE::Vector classprobs;
|
|
|
|
+ classprobs.resize ( gaussians );
|
|
|
|
+ classprobs.set ( 0.0 );
|
|
|
|
+
|
|
|
|
+ double maxp = -numeric_limits<double>::max();
|
|
|
|
+ for ( int i = 0; i < gaussians; i++ )
|
|
|
|
+ {
|
|
|
|
+ classprobs[i] = logpdfState ( vin, i ) + log ( priors[i] ); // log(P(x|i))
|
|
|
|
+
|
|
|
|
+ maxp = std::max ( maxp, classprobs[i] );
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < gaussians; i++ )
|
|
|
|
+ {
|
|
|
|
+ double p = classprobs[i] - maxp;
|
|
|
|
+
|
|
|
|
+ p = exp ( p );
|
|
|
|
+
|
|
|
|
+ for ( int d = 0; d < dim; d++ )
|
|
|
|
+ {
|
|
|
|
+ double diff = vin[d] - mu[i][d];
|
|
|
|
+
|
|
|
|
+ double sigma2 = sparse_sigma[i][d] * sparse_sigma[i][d];
|
|
|
|
+ double sigma3 = sigma2 * sparse_sigma[i][d];
|
|
|
|
+ double normmu = sqrt ( priors[i] / sigma2 );
|
|
|
|
+ double normsig = sqrt ( ( 2.0 * priors[i] ) / sigma2 );
|
|
|
|
+
|
|
|
|
+ double gradmu = ( p * ( diff / sigma2 ) ) / normmu;
|
|
|
|
+ double gradsig = ( p * ( ( diff * diff ) / sigma3 - 1.0 / sparse_sigma[i][d] ) ) / normsig;
|
|
|
|
+ if ( fabs ( gradmu ) > 1e-7f )
|
|
|
|
+ outprobs[counter] = gradmu;
|
|
|
|
+ counter++;
|
|
|
|
+
|
|
|
|
+ if ( fabs ( gradsig ) > 1e-7f )
|
|
|
|
+ outprobs[counter] = gradsig;
|
|
|
|
+ counter++;
|
|
|
|
+ }
|
|
|
|
+ }
|
|
}
|
|
}
|
|
|
|
|
|
void GMM::cluster ( const VVector & features, VVector & prototypes, vector<double> & weights, vector<int> & assignment )
|
|
void GMM::cluster ( const VVector & features, VVector & prototypes, vector<double> & weights, vector<int> & assignment )
|
|
{
|
|
{
|
|
- computeMixture(features);
|
|
|
|
-
|
|
|
|
- prototypes.clear();
|
|
|
|
- weights.clear();
|
|
|
|
- assignment.clear();
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < (int)features.size(); i++)
|
|
|
|
- {
|
|
|
|
- int c = getBestClass(features[i]);
|
|
|
|
- assignment.push_back(c);
|
|
|
|
-
|
|
|
|
- weights.push_back(priors[c]);
|
|
|
|
- }
|
|
|
|
- for(int c = 0; c < gaussians; c++)
|
|
|
|
- prototypes.push_back(mu[c]);
|
|
|
|
-
|
|
|
|
- cout << "tau: " << tau << endl;
|
|
|
|
|
|
+ computeMixture ( features );
|
|
|
|
+
|
|
|
|
+ prototypes.clear();
|
|
|
|
+ weights.clear();
|
|
|
|
+ assignment.clear();
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < ( int ) features.size(); i++ )
|
|
|
|
+ {
|
|
|
|
+ int c = getBestClass ( features[i] );
|
|
|
|
+ assignment.push_back ( c );
|
|
|
|
+
|
|
|
|
+ weights.push_back ( priors[c] );
|
|
|
|
+ }
|
|
|
|
+ for ( int c = 0; c < gaussians; c++ )
|
|
|
|
+ prototypes.push_back ( mu[c] );
|
|
|
|
+
|
|
|
|
+ cout << "tau: " << tau << endl;
|
|
}
|
|
}
|
|
|
|
|
|
-void GMM::saveData(const std::string filename)
|
|
|
|
|
|
+void GMM::saveData ( const std::string filename )
|
|
{
|
|
{
|
|
- ofstream fout( filename.c_str() );
|
|
|
|
-
|
|
|
|
- fout << gaussians << " " << dim << endl;
|
|
|
|
-
|
|
|
|
- mu >> fout;
|
|
|
|
- fout << endl;
|
|
|
|
-
|
|
|
|
- for(int n = 0; n < gaussians; n++)
|
|
|
|
- {
|
|
|
|
- fout << sparse_sigma[n] << endl;
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- for(int n = 0; n < gaussians; n++)
|
|
|
|
- {
|
|
|
|
- fout << priors[n] << " ";
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- fout.close();
|
|
|
|
|
|
+ ofstream fout ( filename.c_str() );
|
|
|
|
+
|
|
|
|
+ fout << gaussians << " " << dim << endl;
|
|
|
|
+
|
|
|
|
+ mu >> fout;
|
|
|
|
+ fout << endl;
|
|
|
|
+
|
|
|
|
+ for ( int n = 0; n < gaussians; n++ )
|
|
|
|
+ {
|
|
|
|
+ fout << sparse_sigma[n] << endl;
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ for ( int n = 0; n < gaussians; n++ )
|
|
|
|
+ {
|
|
|
|
+ fout << priors[n] << " ";
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ fout.close();
|
|
}
|
|
}
|
|
|
|
|
|
-bool GMM::loadData(const std::string filename)
|
|
|
|
|
|
+bool GMM::loadData ( const std::string filename )
|
|
{
|
|
{
|
|
- cerr << "read GMM Data" << endl;
|
|
|
|
- ifstream fin( filename.c_str() );
|
|
|
|
-
|
|
|
|
- if(fin.fail())
|
|
|
|
- {
|
|
|
|
- cerr << "GMM: Datei " << filename << " nicht gefunden!" << endl;
|
|
|
|
- return false;
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- fin >> gaussians;
|
|
|
|
- fin >> dim;
|
|
|
|
- cdimval = log(pow(2*3.14159,dim));
|
|
|
|
- mu.clear();
|
|
|
|
-
|
|
|
|
- for(int i = 0; i < gaussians; i++)
|
|
|
|
- {
|
|
|
|
- NICE::Vector tmp;
|
|
|
|
- fin >> tmp;
|
|
|
|
- mu.push_back(tmp);
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- log_det_sigma.clear();
|
|
|
|
-
|
|
|
|
- for(int n = 0; n < gaussians; n++)
|
|
|
|
- {
|
|
|
|
- NICE::Matrix _sigma;
|
|
|
|
- NICE::Vector _sparse_sigma;
|
|
|
|
-
|
|
|
|
- _sparse_sigma =NICE::Vector(dim);
|
|
|
|
- fin >> _sparse_sigma;
|
|
|
|
- sparse_sigma.push_back(_sparse_sigma);
|
|
|
|
-
|
|
|
|
- sparse_inv_sigma.push_back(diagInverse(sparse_sigma[n]));
|
|
|
|
- log_det_sigma.push_back(logdiagDeterminant(sparse_sigma[n]));
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- for(int n = 0; n < gaussians; n++)
|
|
|
|
- {
|
|
|
|
- double tmpd;
|
|
|
|
- fin >> tmpd;
|
|
|
|
- priors.push_back(tmpd);
|
|
|
|
- }
|
|
|
|
- fin.close();
|
|
|
|
-
|
|
|
|
- cerr << "reading GMM Data finished" << endl;
|
|
|
|
- return true;
|
|
|
|
|
|
+ cerr << "read GMM Data" << endl;
|
|
|
|
+ ifstream fin ( filename.c_str() );
|
|
|
|
+
|
|
|
|
+ if ( fin.fail() )
|
|
|
|
+ {
|
|
|
|
+ cerr << "GMM: Datei " << filename << " nicht gefunden!" << endl;
|
|
|
|
+ return false;
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ fin >> gaussians;
|
|
|
|
+ fin >> dim;
|
|
|
|
+ cdimval = log ( pow ( 2 * 3.14159, dim ) );
|
|
|
|
+ mu.clear();
|
|
|
|
+
|
|
|
|
+ for ( int i = 0; i < gaussians; i++ )
|
|
|
|
+ {
|
|
|
|
+ NICE::Vector tmp;
|
|
|
|
+ fin >> tmp;
|
|
|
|
+ mu.push_back ( tmp );
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ log_det_sigma.clear();
|
|
|
|
+
|
|
|
|
+ for ( int n = 0; n < gaussians; n++ )
|
|
|
|
+ {
|
|
|
|
+ NICE::Matrix _sigma;
|
|
|
|
+ NICE::Vector _sparse_sigma;
|
|
|
|
+
|
|
|
|
+ _sparse_sigma = NICE::Vector ( dim );
|
|
|
|
+ fin >> _sparse_sigma;
|
|
|
|
+ sparse_sigma.push_back ( _sparse_sigma );
|
|
|
|
+
|
|
|
|
+ sparse_inv_sigma.push_back ( diagInverse ( sparse_sigma[n] ) );
|
|
|
|
+ log_det_sigma.push_back ( logdiagDeterminant ( sparse_sigma[n] ) );
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ for ( int n = 0; n < gaussians; n++ )
|
|
|
|
+ {
|
|
|
|
+ double tmpd;
|
|
|
|
+ fin >> tmpd;
|
|
|
|
+ priors.push_back ( tmpd );
|
|
|
|
+ }
|
|
|
|
+ fin.close();
|
|
|
|
+
|
|
|
|
+ cerr << "reading GMM Data finished" << endl;
|
|
|
|
+ return true;
|
|
}
|
|
}
|
|
|
|
|
|
-void GMM::getParams(VVector &mean, VVector &sSigma, vector<double> &p)
|
|
|
|
|
|
+void GMM::getParams ( VVector &mean, VVector &sSigma, vector<double> &p )
|
|
{
|
|
{
|
|
- mean = mu;
|
|
|
|
- sSigma.resize(gaussians);
|
|
|
|
- p.clear();
|
|
|
|
- for(int i = 0; i < gaussians; i++)
|
|
|
|
- {
|
|
|
|
- sSigma[i] = sparse_sigma[i];
|
|
|
|
- p.push_back(priors[i]);
|
|
|
|
- }
|
|
|
|
|
|
+ mean = mu;
|
|
|
|
+ sSigma.resize ( gaussians );
|
|
|
|
+ p.clear();
|
|
|
|
+ for ( int i = 0; i < gaussians; i++ )
|
|
|
|
+ {
|
|
|
|
+ sSigma[i] = sparse_sigma[i];
|
|
|
|
+ p.push_back ( priors[i] );
|
|
|
|
+ }
|
|
}
|
|
}
|
|
|
|
|
|
-void GMM::setCompareGM(VVector mean, VVector sSigma, vector<double> p)
|
|
|
|
|
|
+void GMM::setCompareGM ( VVector mean, VVector sSigma, vector<double> p )
|
|
{
|
|
{
|
|
- mu2 = mean;
|
|
|
|
- sparse_sigma2 = sSigma;
|
|
|
|
- priors2 = vector<double>(p);
|
|
|
|
|
|
+ mu2 = mean;
|
|
|
|
+ sparse_sigma2 = sSigma;
|
|
|
|
+ priors2 = vector<double> ( p );
|
|
}
|
|
}
|
|
-
|
|
|
|
-double GMM::kPPK(NICE::Vector sigma1, NICE::Vector sigma2, NICE::Vector mu1, NICE::Vector mu2, double p)
|
|
|
|
|
|
+
|
|
|
|
+double GMM::kPPK ( NICE::Vector sigma1, NICE::Vector sigma2, NICE::Vector mu1, NICE::Vector mu2, double p )
|
|
{
|
|
{
|
|
- double d = mu1.size();
|
|
|
|
-
|
|
|
|
- double det1 = 1.0;
|
|
|
|
- double det2 = 1.0;
|
|
|
|
- double det3 = 1.0;
|
|
|
|
- double eval = 0.0;
|
|
|
|
- for(int i = 0; i < d; i++)
|
|
|
|
- {
|
|
|
|
- det1 *= sigma1[i];
|
|
|
|
- det2 *= sigma2[i];
|
|
|
|
- double sigma = 1.0/((1.0/sigma1[i]+1.0/sigma2[i])*p);
|
|
|
|
- det3 *= sigma;
|
|
|
|
- double mu = p*(mu1[i]*sigma1[i]+mu2[i]*sigma2[i]);
|
|
|
|
- eval += 0.5*mu*sigma*mu - (p*mu1[i]*mu1[i])/(2.0*sigma1[i]) - (p*mu2[i]*mu2[i]) / (2.0*sigma2[i]);
|
|
|
|
- }
|
|
|
|
-
|
|
|
|
- return (pow(2.0*M_PI, ((1.0-2.0*p)*d)/2.0) * pow(det1, -p/2.0) * pow(det2, -p/2.0) * pow(det3, 0.5) * exp(eval));
|
|
|
|
|
|
+ double d = mu1.size();
|
|
|
|
+
|
|
|
|
+ double det1 = 1.0;
|
|
|
|
+ double det2 = 1.0;
|
|
|
|
+ double det3 = 1.0;
|
|
|
|
+ double eval = 0.0;
|
|
|
|
+ for ( int i = 0; i < d; i++ )
|
|
|
|
+ {
|
|
|
|
+ det1 *= sigma1[i];
|
|
|
|
+ det2 *= sigma2[i];
|
|
|
|
+ double sigma = 1.0 / ( ( 1.0 / sigma1[i] + 1.0 / sigma2[i] ) * p );
|
|
|
|
+ det3 *= sigma;
|
|
|
|
+ double mu = p * ( mu1[i] * sigma1[i] + mu2[i] * sigma2[i] );
|
|
|
|
+ eval += 0.5 * mu * sigma * mu - ( p * mu1[i] * mu1[i] ) / ( 2.0 * sigma1[i] ) - ( p * mu2[i] * mu2[i] ) / ( 2.0 * sigma2[i] );
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ return ( pow ( 2.0*M_PI, ( ( 1.0 - 2.0*p ) *d ) / 2.0 ) * pow ( det1, -p / 2.0 ) * pow ( det2, -p / 2.0 ) * pow ( det3, 0.5 ) * exp ( eval ) );
|
|
}
|
|
}
|
|
|
|
|
|
double GMM::compare()
|
|
double GMM::compare()
|
|
{
|
|
{
|
|
- double distkij = 0.0;
|
|
|
|
- double distkjj = 0.0;
|
|
|
|
- double distkii = 0.0;
|
|
|
|
- for(int i = 0; i < gaussians; i++)
|
|
|
|
- {
|
|
|
|
- for(int j = 0; j < gaussians; j++)
|
|
|
|
- {
|
|
|
|
- double kij = kPPK(sparse_sigma[i],sparse_sigma2[j], mu[i], mu2[j], 0.5);
|
|
|
|
- double kii = kPPK(sparse_sigma[i],sparse_sigma[j], mu[i], mu[j], 0.5);
|
|
|
|
- double kjj = kPPK(sparse_sigma2[i],sparse_sigma2[j], mu2[i], mu2[j], 0.5);
|
|
|
|
- kij = priors[i]*priors2[j]*kij;
|
|
|
|
- kii = priors[i]*priors[j]*kii;
|
|
|
|
- kjj = priors2[i]*priors2[j]*kjj;
|
|
|
|
- distkij += kij;
|
|
|
|
- distkii += kii;
|
|
|
|
- distkjj += kjj;
|
|
|
|
- }
|
|
|
|
- }
|
|
|
|
- return (distkij / (sqrt(distkii)*sqrt(distkjj)));
|
|
|
|
|
|
+ double distkij = 0.0;
|
|
|
|
+ double distkjj = 0.0;
|
|
|
|
+ double distkii = 0.0;
|
|
|
|
+ for ( int i = 0; i < gaussians; i++ )
|
|
|
|
+ {
|
|
|
|
+ for ( int j = 0; j < gaussians; j++ )
|
|
|
|
+ {
|
|
|
|
+ double kij = kPPK ( sparse_sigma[i], sparse_sigma2[j], mu[i], mu2[j], 0.5 );
|
|
|
|
+ double kii = kPPK ( sparse_sigma[i], sparse_sigma[j], mu[i], mu[j], 0.5 );
|
|
|
|
+ double kjj = kPPK ( sparse_sigma2[i], sparse_sigma2[j], mu2[i], mu2[j], 0.5 );
|
|
|
|
+ kij = priors[i] * priors2[j] * kij;
|
|
|
|
+ kii = priors[i] * priors[j] * kii;
|
|
|
|
+ kjj = priors2[i] * priors2[j] * kjj;
|
|
|
|
+ distkij += kij;
|
|
|
|
+ distkii += kii;
|
|
|
|
+ distkjj += kjj;
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ return ( distkij / ( sqrt ( distkii ) *sqrt ( distkjj ) ) );
|
|
}
|
|
}
|
|
|
|
|
|
-void GMM::comparing(bool c)
|
|
|
|
|
|
+void GMM::comparing ( bool c )
|
|
{
|
|
{
|
|
- comp = c;
|
|
|
|
|
|
+ comp = c;
|
|
}
|
|
}
|