current version :)

cbaselib/Example.h

@@ -1,7 +1,7 @@
 #ifndef EXAMPLEINCLUDE
 #define EXAMPLEINCLUDE
 
-/** 
+/**
 * @file Example.h
 * @brief data caching of several feature images and many more
 * @author Erik Rodner
@@ -15,121 +15,121 @@ namespace OBJREC {
 /** Cached example with window information */
 class Example
 {
-    public:
-	/** weight of example */
-	double weight;
-
-	/** complete image example */
-	CachedExample *ce;
-
-	/** x position of window */
-	long long int x;
-	/** y position of window */
-	long long int y;
-
-	/** width of window */
-	int width;
-	/** height of window */
-	int height;
-	
-	//! if some examples are related, they have the same position
-	int position;
-	
-	//! scale of the feature
-	double scale;
-	
-	/** evil workaround: store simple vector example */
-	NICE::Vector *vec;
-	
-	/** evil workaround: store sparse vector example */
-	NICE::SparseVector *svec;
-
-	/** simple constructor, initialize everything with 0
-	 */
-	Example();
-
-	/** constructor using a window covering the whole image
-	    @param ce associated image data
-	*/
-	Example ( CachedExample *_ce );
-
-	/** constructor 
-	    @param ce associated image data
-	    @param x x position of window
-	    @param y y position of window
-	    @param weight weight of example
-	*/
-	Example ( CachedExample *ce, int x, int y, double weight = 1.0 );
-	
-	/** constructor 
-	    @param ce associated image data
-	    @param x x position of window
-	    @param y y position of window
-	    @param weight weight of example
-	*/
-	Example ( CachedExample *ce, int x, int y, int width, int height, double weight = 1.0 );
-
-	/** evil workaround: simple constructors for std::vector examples 
-	    @param vec simple feature vector
-	    @param weight optional weight
-	*/
-	Example ( NICE::Vector *vec, double weight = 1.0 );
-
-	/**
-	 * copy constructor
-	 * @param ex input Example
-	 */
-	Example ( const Example &ex);
-	
-	/**
-	 * copies the values of ex in this example
-	 * @param ex input Example
-	 */
-	void copy ( const Example &ex);
-	
-	/** destructor */
-	~Example ();
-
-	/** delete all data associated with this example
-	    (sparse vector, vector, cached example, etc.) */
-	void clean ();
-	
-	/**
-	 * load from file (warning: no cachedexamples supported yet
-	 * @param is file
-	 * @param format 
-	 */
-	void restore (std::istream & is, int format = 0);
-	
-	
-	/**
-	 * write to file (warning: no cachedexamples supported yet
-	 * @param is file
-	 * @param format 
-	 */
-	void store (std::ostream & os, int format = 0) const;
+  public:
+    /** weight of example */
+    double weight;
+
+    /** complete image example */
+    CachedExample *ce;
+
+    /** x position of window */
+    long long int x;
+    /** y position of window */
+    long long int y;
+
+    /** width of window */
+    int width;
+    /** height of window */
+    int height;
+
+    //! if some examples are related, they have the same position
+    int position;
+
+    //! scale of the feature
+    double scale;
+
+    /** evil workaround: store simple vector example */
+    NICE::Vector *vec;
+
+    /** evil workaround: store sparse vector example */
+    NICE::SparseVector *svec;
+
+    /** simple constructor, initialize everything with 0
+     */
+    Example();
+
+    /** constructor using a window covering the whole image
+        @param ce associated image data
+    */
+    Example ( CachedExample *_ce );
+
+    /** constructor
+        @param ce associated image data
+        @param x x position of window
+        @param y y position of window
+        @param weight weight of example
+    */
+    Example ( CachedExample *ce, int x, int y, double weight = 1.0 );
+
+    /** constructor
+        @param ce associated image data
+        @param x x position of window
+        @param y y position of window
+        @param weight weight of example
+    */
+    Example ( CachedExample *ce, int x, int y, int width, int height, double weight = 1.0 );
+
+    /** evil workaround: simple constructors for std::vector examples
+        @param vec simple feature vector
+        @param weight optional weight
+    */
+    Example ( NICE::Vector *vec, double weight = 1.0 );
+
+    /**
+     * copy constructor
+     * @param ex input Example
+     */
+    Example ( const Example &ex );
+
+    /**
+     * copies the values of ex in this example
+     * @param ex input Example
+     */
+    void copy ( const Example &ex );
+
+    /** destructor */
+    ~Example ();
+
+    /** delete all data associated with this example
+        (sparse vector, vector, cached example, etc.) */
+    void clean ();
+
+    /**
+     * load from file (warning: no cachedexamples supported yet
+     * @param is file
+     * @param format
+     */
+    void restore ( std::istream & is, int format = 0 );
+
+
+    /**
+     * write to file (warning: no cachedexamples supported yet
+     * @param is file
+     * @param format
+     */
+    void store ( std::ostream & os, int format = 0 ) const;
 };
 
-/** labeled pair of Example 
+/** labeled pair of Example
     @see Example */
 class Examples : public std::vector< std::pair<int, Example> >
 {
-    public:
-		std::string filename;
+  public:
+    std::string filename;
 
-    inline int getMaxClassNo () const 
+    inline int getMaxClassNo () const
     {
-		int maxClassno = -1;
-		for ( const_iterator i = begin(); i != end(); i++ )
-			if ( i->first > maxClassno )
-			maxClassno = i->first;
-		
-		return maxClassno;
+      int maxClassno = -1;
+      for ( const_iterator i = begin(); i != end(); i++ )
+        if ( i->first > maxClassno )
+          maxClassno = i->first;
+
+      return maxClassno;
     }
 
-	/** delete all data associated with all examples
-	    (sparse vector, vector, cached example, etc.) */
-	void clean ();
+    /** delete all data associated with all examples
+        (sparse vector, vector, cached example, etc.) */
+    void clean ();
 };
 
 

classifier/fpclassifier/logisticregression/FPCSMLR.cpp

@@ -9,153 +9,153 @@ using namespace NICE;
 
 FPCSMLR::FPCSMLR ()
 {
-	inpic = false;
+  inpic = false;
 }
 
-FPCSMLR::FPCSMLR( const Config *_conf, string section ) : conf(_conf)
+FPCSMLR::FPCSMLR ( const Config *_conf, string section ) : conf ( _conf )
 {
-	confsection = section;
-	inpic = conf->gB(section, "inpic", false );
+  confsection = section;
+  inpic = conf->gB ( section, "inpic", false );
 }
 
 FPCSMLR::~FPCSMLR()
 {
-	//clean up
+  //clean up
 }
 
 ClassificationResult FPCSMLR::classify ( Example & pce )
 {
-	FullVector overall_distribution(maxClassNo+1);
-	overall_distribution[maxClassNo] = 0.0;
-
-	double maxp = -numeric_limits<double>::max();
-	int classno = 0;
-	
-	double sum  = 0.0;
-		
-	for(int i = 0; i < maxClassNo; i++)
-	{
-		overall_distribution[i] = classifiers[i].classify(pce);
-		
-		sum += overall_distribution[i];
-		
-		if(maxp < overall_distribution[i])
-		{
-			classno = i;
-			maxp = overall_distribution[i];
-		}
-	}
-	for(int i = 0; i < maxClassNo; i++)
-	{
-		overall_distribution[i] /= sum;
-	}
-
-	return ClassificationResult ( classno, overall_distribution );
+  FullVector overall_distribution ( maxClassNo + 1 );
+  overall_distribution[maxClassNo] = 0.0;
+
+  double maxp = -numeric_limits<double>::max();
+  int classno = 0;
+
+  double sum  = 0.0;
+
+  for ( int i = 0; i < maxClassNo; i++ )
+  {
+    overall_distribution[i] = classifiers[i].classify ( pce );
+
+    sum += overall_distribution[i];
+
+    if ( maxp < overall_distribution[i] )
+    {
+      classno = i;
+      maxp = overall_distribution[i];
+    }
+  }
+  for ( int i = 0; i < maxClassNo; i++ )
+  {
+    overall_distribution[i] /= sum;
+  }
+
+  return ClassificationResult ( classno, overall_distribution );
 }
 
 void FPCSMLR::train ( FeaturePool & _fp, Examples & examples )
 {
 
-	cout << "start train" << endl;
-	fp = FeaturePool(_fp);
-    
-	// Anzahl von Merkmalen
-	int fanz = examples.size();
-
-	maxClassNo = -1;
-	for(int i = 0; i < fanz; i++)
-	{
-		maxClassNo = std::max(maxClassNo, examples[i].first);
-	}
-	maxClassNo++;
-	
-	assert(fanz >= maxClassNo);
-	
-	classifiers.resize(maxClassNo);
-	for(int i = 0; i < maxClassNo; i++)
-	{
-		cout << "classifier no " << i << " training starts" << endl;
-		classifiers[i] = SLR(conf, confsection);
-
-		if(inpic)
-		{
-	
-			vector<bool> classinpic;
-			
-			for(int j = 0; j < (int)examples.size(); j++)
-			{
-				if(examples[j].first == i)
-				{
-					if(examples[j].second.position < (int)classinpic.size())
-						classinpic[examples[j].second.position] = true;
-					else if(examples[j].second.position == (int)classinpic.size())
-						classinpic.push_back(true);
-					else
-					{
-						while(examples[j].second.position > (int)classinpic.size())
-						{
-							classinpic.push_back(false);
-						}
-						classinpic.push_back(true);
-					}
-				}
-			}
-			
-			Examples ex2;
-
-			for(int j = 0; j < (int)examples.size(); j++)
-			{
-				if(examples[j].second.position >= (int)classinpic.size())
-						continue;
-				if(classinpic[examples[j].second.position])
-				{
-					Example e;
-					e.svec = examples[j].second.svec;
-					e.vec = examples[j].second.vec;
-					ex2.push_back ( pair<int, Example> ( examples[j].first, e ) );
-				}
-			}
-			cout << "examples for class " << i << ": " << ex2.size() << endl;
-			
-			if(ex2.size() <= 2)
-				continue;
-			
-			classifiers[i].train(_fp, ex2, i);
-			
-			for(int j = 0; j < (int)ex2.size(); j++)
-			{
-				ex2[j].second.svec = NULL;
-			}
-		}
-		else
-		{
-			classifiers[i].train(_fp, examples, i);
-		}
-	}
-
-	cout << "end train" << endl;
+  cout << "start train" << endl;
+  fp = FeaturePool ( _fp );
+
+  // Anzahl von Merkmalen
+  int fanz = examples.size();
+
+  maxClassNo = -1;
+  for ( int i = 0; i < fanz; i++ )
+  {
+    maxClassNo = std::max ( maxClassNo, examples[i].first );
+  }
+  maxClassNo++;
+
+  assert ( fanz >= maxClassNo );
+
+  classifiers.resize ( maxClassNo );
+  for ( int i = 0; i < maxClassNo; i++ )
+  {
+    cout << "classifier no " << i << " training starts" << endl;
+    classifiers[i] = SLR ( conf, confsection );
+
+    if ( inpic )
+    {
+
+      vector<bool> classinpic;
+
+      for ( int j = 0; j < ( int ) examples.size(); j++ )
+      {
+        if ( examples[j].first == i )
+        {
+          if ( examples[j].second.position < ( int ) classinpic.size() )
+            classinpic[examples[j].second.position] = true;
+          else if ( examples[j].second.position == ( int ) classinpic.size() )
+            classinpic.push_back ( true );
+          else
+          {
+            while ( examples[j].second.position > ( int ) classinpic.size() )
+            {
+              classinpic.push_back ( false );
+            }
+            classinpic.push_back ( true );
+          }
+        }
+      }
+
+      Examples ex2;
+
+      for ( int j = 0; j < ( int ) examples.size(); j++ )
+      {
+        if ( examples[j].second.position >= ( int ) classinpic.size() )
+          continue;
+        if ( classinpic[examples[j].second.position] )
+        {
+          Example e;
+          e.svec = examples[j].second.svec;
+          e.vec = examples[j].second.vec;
+          ex2.push_back ( pair<int, Example> ( examples[j].first, e ) );
+        }
+      }
+      cout << "examples for class " << i << ": " << ex2.size() << endl;
+
+      if ( ex2.size() <= 2 )
+        continue;
+
+      classifiers[i].train ( _fp, ex2, i );
+
+      for ( int j = 0; j < ( int ) ex2.size(); j++ )
+      {
+        ex2[j].second.svec = NULL;
+      }
+    }
+    else
+    {
+      classifiers[i].train ( _fp, examples, i );
+    }
+  }
+
+  cout << "end train" << endl;
 }
 
 
-void FPCSMLR::restore (istream & is, int format)
+void FPCSMLR::restore ( istream & is, int format )
 {
-	is >> maxClassNo;
-	classifiers.resize(maxClassNo);
-	for(int i = 0; i < maxClassNo; i++)
-	{
-		classifiers[i].restore(is, format);
-	}
+  is >> maxClassNo;
+  classifiers.resize ( maxClassNo );
+  for ( int i = 0; i < maxClassNo; i++ )
+  {
+    classifiers[i].restore ( is, format );
+  }
 }
 
-void FPCSMLR::store (ostream & os, int format) const
+void FPCSMLR::store ( ostream & os, int format ) const
 {
 
-	if(format!=-9999) os << maxClassNo;
+  if ( format != -9999 ) os << maxClassNo;
 
-	for(int i = 0; i < maxClassNo; i++)
-	{
-		classifiers[i].store(os, format);
-	}
+  for ( int i = 0; i < maxClassNo; i++ )
+  {
+    classifiers[i].store ( os, format );
+  }
 }
 
 void FPCSMLR::clear ()
@@ -165,15 +165,15 @@ void FPCSMLR::clear ()
 
 FeaturePoolClassifier *FPCSMLR::clone () const
 {
-	//TODO: wenn alle Variablen der Klasse bekannt sind, dann übergebe diese mit
-	FPCSMLR *o = new FPCSMLR ( conf, confsection );
+  //TODO: wenn alle Variablen der Klasse bekannt sind, dann übergebe diese mit
+  FPCSMLR *o = new FPCSMLR ( conf, confsection );
 
-	o->maxClassNo = maxClassNo;
+  o->maxClassNo = maxClassNo;
 
-	return o;
+  return o;
 }
 
 void FPCSMLR::setComplexity ( int size )
 {
-	cerr << "FPCSMLR: no complexity to set" << endl;
+  cerr << "FPCSMLR: no complexity to set" << endl;
 }

classifier/fpclassifier/logisticregression/SLR.cpp

@@ -18,115 +18,115 @@ SLR::SLR ()
 {
 }
 
-SLR::SLR( const Config *_conf, string section ) : conf(_conf)
+SLR::SLR ( const Config *_conf, string section ) : conf ( _conf )
 {
-	maxiter = conf->gI(section, "maxiter", 20000 );
-	resamp_decay = conf->gD(section, "resamp_decay", 0.5 );
-	convergence_tol = conf->gD(section, "convergence_tol", 1e-7 );
-	min_resamp =  conf->gD(section, "min_resamp", 0.001 );
-	lambda = conf->gD(section, "lambda", 0.1 );
-	samplesperclass = conf->gD(section, "samplesperclass", 200.0);
-	weight.setDim(0);
-	fdim = 0;
+  maxiter = conf->gI ( section, "maxiter", 20000 );
+  resamp_decay = conf->gD ( section, "resamp_decay", 0.5 );
+  convergence_tol = conf->gD ( section, "convergence_tol", 1e-7 );
+  min_resamp =  conf->gD ( section, "min_resamp", 0.001 );
+  lambda = conf->gD ( section, "lambda", 0.1 );
+  samplesperclass = conf->gD ( section, "samplesperclass", 200.0 );
+  weight.setDim ( 0 );
+  fdim = 0;
 }
 
 SLR::~SLR()
 {
-	//clean up
+  //clean up
 }
 
 double SLR::classify ( Example & pce )
 {
-	double result;
-	SparseVector *svec;
+  double result;
+  SparseVector *svec;
 
-	if(weight.getDim() == 0)
-		return 0.0;
+  if ( weight.getDim() == 0 )
+    return 0.0;
 
-	bool newvec = false;
-	
-	if(pce.svec != NULL)
-	{
-		svec = pce.svec;
-	}
-	else
-	{
-		Vector x;
+  bool newvec = false;
 
-		x = *(pce.vec);
+  if ( pce.svec != NULL )
+  {
+    svec = pce.svec;
+  }
+  else
+  {
+    Vector x;
+
+    x = * ( pce.vec );
 
 #ifdef featnorm
-		for(int m = 0; m < (int)x.size(); m++)
-		{
-			x[m] = (x[m]-minval[m]) /(maxval[m] - minval[m]);
-		}
+    for ( int m = 0; m < ( int ) x.size(); m++ )
+    {
+      x[m] = ( x[m] - minval[m] ) / ( maxval[m] - minval[m] );
+    }
 #endif
-		svec = new SparseVector(x);
+    svec = new SparseVector ( x );
+
+    svec->setDim ( x.size() );
 
-		svec->setDim(x.size());
+    newvec = true;
 
-		newvec = true;
+  }
 
-	}
-		
-	if(weight.size() == 0)
-		result = 0.0;
-	else
-	{
-		result = 1.0/(1.0+exp(-svec->innerProduct(weight)));
-	}
+  if ( weight.size() == 0 )
+    result = 0.0;
+  else
+  {
+    result = 1.0 / ( 1.0 + exp ( -svec->innerProduct ( weight ) ) );
+  }
 
-	if(newvec)
-		delete svec;
+  if ( newvec )
+    delete svec;
 
-	return result;
+  return result;
 }
 
 void SLR::train ( FeaturePool & _fp, Examples & examples, int classno )
 {
 
-	cout << "start train" << endl;
-	fp = FeaturePool(_fp);
-    
-	// Anzahl von Merkmalen
-	int fanz = examples.size();
-	assert(fanz >= 2);
-
-	// Merkmalsdimension bestimmen
-	Vector x;
-	fp.calcFeatureVector(examples[0].second, x);
-	fdim = x.size();
-	assert(fdim > 0);
-
-	stepwise_regression(examples, classno);	
-	cout << "end train" << endl;
+  cout << "start train" << endl;
+  fp = FeaturePool ( _fp );
+
+  // Anzahl von Merkmalen
+  int fanz = examples.size();
+  assert ( fanz >= 2 );
+
+  // Merkmalsdimension bestimmen
+  Vector x;
+  fp.calcFeatureVector ( examples[0].second, x );
+  fdim = x.size();
+  assert ( fdim > 0 );
+
+  stepwise_regression ( examples, classno );
+  cout << "end train" << endl;
 }
 
-void SLR::restore (istream & is, int format)
+void SLR::restore ( istream & is, int format )
 {
-	weight.restore(is, format);
-	fdim = (int) weight.getDim();
-	maxval.resize(fdim);
-	minval.resize(fdim);
-	for(int i = 0; i < fdim; i++)
-	{
-		is >> maxval[i];
-		is >> minval[i];
-	}
+  weight.restore ( is, format );
+  fdim = ( int ) weight.getDim();
+  maxval.resize ( fdim );
+  minval.resize ( fdim );
+  for ( int i = 0; i < fdim; i++ )
+  {
+    is >> maxval[i];
+    is >> minval[i];
+  }
 }
 
-void SLR::store (ostream & os, int format) const
+void SLR::store ( ostream & os, int format ) const
 
 {
-	if(format!=-9999){
-		weight.store(os, format);
-		for(int i = 0; i < fdim; i++)
-		{
-			os << maxval[i] << " " << minval[i] << endl;
-		}
-	}else{
-		weight.store(os, format);
-	}
+  if ( format != -9999 ) {
+    weight.store ( os, format );
+    for ( int i = 0; i < fdim; i++ )
+    {
+      os << maxval[i] << " " << minval[i] << endl;
+    }
+  } else {
+    weight.store ( os, format );
+  }
 }
 
 void SLR::clear ()
@@ -134,385 +134,385 @@ void SLR::clear ()
 //TODO: einbauen
 }
 
-int SLR::stepwise_regression(Examples &x, int classno)
+int SLR::stepwise_regression ( Examples &x, int classno )
 {
-	// initialize randomization
-	srand (time(NULL));
-	//srand (1);
-	
-	cout << "start regression" << endl;
-	
-	// get the number of features
-	int fnum = x.size();
-	
-	// create Datamatrix
-	GMSparseVectorMatrix X;
-	
-	GMSparseVectorMatrix Y;
-	
-	vector<int> count(2, 0);
-
-	maxval.resize(fdim);
-	minval.resize(fdim);
-
-	if(x[0].second.svec == NULL) //input normal vectors
-	{
-		Vector feat;
-
-		for(int i = 0; i < fdim; i++)
-		{
-			maxval[i] = numeric_limits<double>::min();
-			minval[i] = numeric_limits<double>::max();
-		}
-		
-		for(int i = 0; i < fnum; i++)
-		{
-			int pos = 0;
-			if(x[i].first != classno)
-				pos = 1;
-			
-			++count[pos];
-
-			fp.calcFeatureVector(x[i].second, feat);
+  // initialize randomization
+  srand ( time ( NULL ) );
+  //srand (1);
+
+  cout << "start regression" << endl;
+
+  // get the number of features
+  int fnum = x.size();
+
+  // create Datamatrix
+  GMSparseVectorMatrix X;
+
+  GMSparseVectorMatrix Y;
+
+  vector<int> count ( 2, 0 );
+
+  maxval.resize ( fdim );
+  minval.resize ( fdim );
+
+  if ( x[0].second.svec == NULL ) //input normal vectors
+  {
+    Vector feat;
+
+    for ( int i = 0; i < fdim; i++ )
+    {
+      maxval[i] = numeric_limits<double>::min();
+      minval[i] = numeric_limits<double>::max();
+    }
+
+    for ( int i = 0; i < fnum; i++ )
+    {
+      int pos = 0;
+      if ( x[i].first != classno )
+        pos = 1;
+
+      ++count[pos];
+
+      fp.calcFeatureVector ( x[i].second, feat );
 #ifdef featnorm
-			for(int m = 0; m < (int)feat.size(); m++)
-			{
-				minval[m] = std::min(minval[m], feat[m]);
-				maxval[m] = std::max(maxval[m], feat[m]);
-			}
+      for ( int m = 0; m < ( int ) feat.size(); m++ )
+      {
+        minval[m] = std::min ( minval[m], feat[m] );
+        maxval[m] = std::max ( maxval[m], feat[m] );
+      }
 #endif
-			fdim = feat.size();
-		}
-	}
-	else //input Sparse Vectors
-	{
-		for(int i = 0; i < fnum; i++)
-		{
-			int pos = 0;
-			if(x[i].first != classno)
-				pos = 1;
-			
-			++count[pos];
-		}
-		fdim = x[0].second.svec->getDim();
-	}
-	
-	if(count[0] == 0 || count[1] == 0)
-	{
-		cerr << "not enought samples " << count[0] << " : " << count[1] << endl;
-		weight.setDim(0);
-		return -1;
-	}
-
-	double samples = std::min(count[0], count[1]);
-	samples = std::min(samples, samplesperclass);
-	//samples = std::max(samples, 200);
-	//samples = samplesperclass;
-	
-	vector<double> rands(2,0.0);
-	vector<double> allsizes(2,0.0);
-	for(int i = 0; i < 2; i++)
-	{
-		rands[i] = 1.0;								//use all samples (default if samples>count[i])
-		if(samples>0 && samples<1) rands[i] = samples;				//use relative number of samples wrt. class size
-		if(samples>1 && samples<=count[i]) rands[i] = samples/(double)count[i]; //use (approximately) fixed absolute number of samples for each class
-		allsizes[i]=count[i];
-		count[i] = 0;
-	}
-
-	for(int i = 0; i < fnum; i++)
-	{
-		int pos = 0;
-		if(x[i].first != classno)
-			pos = 1;
-			
-			
-		double r = (double)rand()/(double)RAND_MAX;	
-			
-		if( r > rands[pos])
-			continue;
-		++count[pos];
-
-		if(x[0].second.svec == NULL) 
-		{
-			Vector feat;
-			fp.calcFeatureVector(x[i].second, feat);
+      fdim = feat.size();
+    }
+  }
+  else //input Sparse Vectors
+  {
+    for ( int i = 0; i < fnum; i++ )
+    {
+      int pos = 0;
+      if ( x[i].first != classno )
+        pos = 1;
+
+      ++count[pos];
+    }
+    fdim = x[0].second.svec->getDim();
+  }
+
+  if ( count[0] == 0 || count[1] == 0 )
+  {
+    cerr << "not enought samples " << count[0] << " : " << count[1] << endl;
+    weight.setDim ( 0 );
+    return -1;
+  }
+
+  double samples = std::min ( count[0], count[1] );
+  samples = std::min ( samples, samplesperclass );
+  //samples = std::max(samples, 200);
+  //samples = samplesperclass;
+
+  vector<double> rands ( 2, 0.0 );
+  vector<double> allsizes ( 2, 0.0 );
+  for ( int i = 0; i < 2; i++ )
+  {
+    rands[i] = 1.0;        //use all samples (default if samples>count[i])
+    if ( samples > 0 && samples < 1 ) rands[i] = samples;    //use relative number of samples wrt. class size
+    if ( samples > 1 && samples <= count[i] ) rands[i] = samples / ( double ) count[i]; //use (approximately) fixed absolute number of samples for each class
+    allsizes[i] = count[i];
+    count[i] = 0;
+  }
+
+  for ( int i = 0; i < fnum; i++ )
+  {
+    int pos = 0;
+    if ( x[i].first != classno )
+      pos = 1;
+
+
+    double r = ( double ) rand() / ( double ) RAND_MAX;
+
+    if ( r > rands[pos] )
+      continue;
+    ++count[pos];
+
+    if ( x[0].second.svec == NULL )
+    {
+      Vector feat;
+      fp.calcFeatureVector ( x[i].second, feat );
 #ifdef featnorm
-			for(int m = 0; m < (int)feat.size(); m++)
-			{
-				feat[m] = (feat[m]-minval[m]) /(maxval[m] - minval[m]);
-			}
+      for ( int m = 0; m < ( int ) feat.size(); m++ )
+      {
+        feat[m] = ( feat[m] - minval[m] ) / ( maxval[m] - minval[m] );
+      }
 #endif
 
-			X.addRow(feat);
-		}
-		else
-		{
-			X.addRow(x[i].second.svec);
-		}
-		SparseVector *v = new SparseVector(2);
-		(*v)[pos] = 1.0;
-		Y.addRow(v);
-	}
-	Y.setDel();
-
-	if(x[0].second.svec == NULL)
-		X.setDel();
-		
-	for(int i = 0; i < 2; i++)
-	{
-		cerr << "Examples for class " << i << ": " << count[i] << " out of " << allsizes[i] << " with p = " << rands[i] << endl;
-	}	
-	
+      X.addRow ( feat );
+    }
+    else
+    {
+      X.addRow ( x[i].second.svec );
+    }
+    SparseVector *v = new SparseVector ( 2 );
+    ( *v ) [pos] = 1.0;
+    Y.addRow ( v );
+  }
+  Y.setDel();
+
+  if ( x[0].second.svec == NULL )
+    X.setDel();
+
+  for ( int i = 0; i < 2; i++ )
+  {
+    cerr << "Examples for class " << i << ": " << count[i] << " out of " << allsizes[i] << " with p = " << rands[i] << endl;
+  }
+
 #undef NORMALIZATION
 #ifdef NORMALIZATION
-	GMSparseVectorMatrix Xred;
-	Xred.resize(X.rows(), X.cols());
-
-	for(int r = 0; r < (int)Xred.rows(); r++)
-	{
-		for(int c = 0; c < (int)Xred.cols(); c++)
-		{
-			double tmp = X[r].get(c);
-			
-			if(Y[r].get(0) == 1)
-				tmp *= count[0]/fnum;
-			else
-				tmp *= count[1]/fnum;
-			
-			if(fabs(tmp) > 10e-7)
-				Xred[r][c] = tmp;
-		}
-	}
+  GMSparseVectorMatrix Xred;
+  Xred.resize ( X.rows(), X.cols() );
+
+  for ( int r = 0; r < ( int ) Xred.rows(); r++ )
+  {
+    for ( int c = 0; c < ( int ) Xred.cols(); c++ )
+    {
+      double tmp = X[r].get ( c );
+
+      if ( Y[r].get ( 0 ) == 1 )
+        tmp *= count[0] / fnum;
+      else
+        tmp *= count[1] / fnum;
+
+      if ( fabs ( tmp ) > 10e-7 )
+        Xred[r][c] = tmp;
+    }
+  }
 #endif
 
-	fnum = X.rows();
-	
-	GMSparseVectorMatrix xY;
+  fnum = X.rows();
+
+  GMSparseVectorMatrix xY;
 #ifdef NORMALIZATION
-	Xred.mult(Y, xY, true);
+  Xred.mult ( Y, xY, true );
 #else
-	X.mult(Y, xY, true);
+  X.mult ( Y, xY, true );
 #endif
 
-	weight.setDim(fdim);
-
-	// for faster Computing Xw = X*w	
-	GMSparseVectorMatrix Xw;
-	X.mult(weight, Xw, false, true);
-	SparseVector S(fnum);
-
-	for(int r = 0; r < fnum; r++)
-	{
-		S[r] = 0.0;
-		for(int c = 0; c < 2; c++)
-		{
-			S[r] += exp(Xw[r].get(c));
-		}
-	}
-
-	// for faster computing ac[i] = (maxClassNo-1)/(2*maxClassNo) * Sum_j (x_i (j))^2
-	
-	Vector ac(fdim);
-	Vector lm_2_ac(fdim);
-	
-	for(int f = 0; f < fdim; f++)
-	{
-		ac[f] = 0.0;
-		for(int a = 0; a < fnum; a++)
-		{
-			ac[f] += X[a].get(f)*X[a].get(f);
-		}
-		ac[f] *= 0.25;
-		lm_2_ac[f] = (lambda/2.0)/ac[f];
-	}
-	
-	// initialize the iterative optimization
-	double incr = numeric_limits<double>::max();
-	long non_zero = 0;
-	long wasted_basis = 0;
-	long needed_basis = 0;
-	
-	// prob of resample each weight
-	vector<double> p_resamp;
-	p_resamp.resize(fdim);
-
-	// loop over cycles
-	long cycle = 0;
-	
-	for (cycle = 0; cycle < maxiter; cycle++)
-	{
+  weight.setDim ( fdim );
+
+  // for faster Computing Xw = X*w
+  GMSparseVectorMatrix Xw;
+  X.mult ( weight, Xw, false, true );
+  SparseVector S ( fnum );
+
+  for ( int r = 0; r < fnum; r++ )
+  {
+    S[r] = 0.0;
+    for ( int c = 0; c < 2; c++ )
+    {
+      S[r] += exp ( Xw[r].get ( c ) );
+    }
+  }
+
+  // for faster computing ac[i] = (maxClassNo-1)/(2*maxClassNo) * Sum_j (x_i (j))^2
+
+  Vector ac ( fdim );
+  Vector lm_2_ac ( fdim );
+
+  for ( int f = 0; f < fdim; f++ )
+  {
+    ac[f] = 0.0;
+    for ( int a = 0; a < fnum; a++ )
+    {
+      ac[f] += X[a].get ( f ) * X[a].get ( f );
+    }
+    ac[f] *= 0.25;
+    lm_2_ac[f] = ( lambda / 2.0 ) / ac[f];
+  }
+
+  // initialize the iterative optimization
+  double incr = numeric_limits<double>::max();
+  long non_zero = 0;
+  long wasted_basis = 0;
+  long needed_basis = 0;
+
+  // prob of resample each weight
+  vector<double> p_resamp;
+  p_resamp.resize ( fdim );
+
+  // loop over cycles
+  long cycle = 0;
+
+  for ( cycle = 0; cycle < maxiter; cycle++ )
+  {
 #ifdef SLRDEBUG
-		cerr << "iteration: " << cycle << " of " << maxiter << endl;
+    cerr << "iteration: " << cycle << " of " << maxiter << endl;
 #endif
-    	// zero out the diffs for assessing change
-		double sum2_w_diff = 0.0;
-		double sum2_w_old = 0.0;
-		wasted_basis = 0;
-		if (cycle==1)
-			needed_basis = 0;		
-		
-    	// update each weight
+    // zero out the diffs for assessing change
+    double sum2_w_diff = 0.0;
+    double sum2_w_old = 0.0;
+    wasted_basis = 0;
+    if ( cycle == 1 )
+      needed_basis = 0;
+
+    // update each weight
 //#pragma omp parallel for
-		for (int basis = 0; basis < fdim; basis++) // über alle Dimensionen
-		{
-			int changed = 0;
-			// get the starting weight
-			double w_old = weight.get(basis);
-
-			// set the p_resamp if it's the first cycle
-			if (cycle == 0)
-			{
-				p_resamp[basis] = 1.0;
-			}
-
-			// see if we're gonna update
-			double rval = (double)rand()/(double)RAND_MAX;			
-			
-			if ((w_old != 0) || (rval < p_resamp[basis]))
-			{
-	  			// calc the probability
-				double XdotP = 0.0;
-				for (int i = 0; i < fnum; i++)
-				{
+    for ( int basis = 0; basis < fdim; basis++ ) // über alle Dimensionen
+    {
+      int changed = 0;
+      // get the starting weight
+      double w_old = weight.get ( basis );
+
+      // set the p_resamp if it's the first cycle
+      if ( cycle == 0 )
+      {
+        p_resamp[basis] = 1.0;
+      }
+
+      // see if we're gonna update
+      double rval = ( double ) rand() / ( double ) RAND_MAX;
+
+      if ( ( w_old != 0 ) || ( rval < p_resamp[basis] ) )
+      {
+        // calc the probability
+        double XdotP = 0.0;
+        for ( int i = 0; i < fnum; i++ )
+        {
 #ifdef NORMALIZATION
-					double e = Xred[i].get(basis) * exp(Xw[i].get(0))/S[i];
+          double e = Xred[i].get ( basis ) * exp ( Xw[i].get ( 0 ) ) / S[i];
 #else
-					double e = X[i].get(basis) * exp(Xw[i].get(0))/S[i];
+          double e = X[i].get ( basis ) * exp ( Xw[i].get ( 0 ) ) / S[i];
 #endif
-					if(finite(e))
-						XdotP += e;
+          if ( finite ( e ) )
+            XdotP += e;
 #ifdef SLRDEBUG
-					else
-						throw "numerical problems";
+          else
+            throw "numerical problems";
 #endif
-				}
-					
-				// get the gradient
-				double grad = xY[basis].get(0) - XdotP;
-	  			// set the new weight
-				double w_new = w_old + grad/ac[basis];
-
-				// test that we're within bounds
-				if (w_new > lm_2_ac[basis])
-				{
-	    			// more towards bounds, but keep it
-					w_new -= lm_2_ac[basis];
-					changed = 1;
-
-	    			// umark from being zero if necessary
-					if (w_old == 0.0)
-					{
-						non_zero += 1;
-
-	      				// reset the p_resample
-						p_resamp[basis] = 1.0;
-
-	      				// we needed the basis
-						needed_basis += 1;
-					}
-				}
-				else if (w_new < -lm_2_ac[basis])
-				{
-	    			// more towards bounds, but keep it
-					w_new += lm_2_ac[basis];
-					changed = 1;
-
-	    			// umark from being zero if necessary
-					if (w_old == 0.0)
-					{
-						non_zero += 1;
-
-	      				// reset the p_resample
-						p_resamp[basis] = 1.0;
-
-	      				// we needed the basis
-						needed_basis += 1;
-					}
-				}
-				else
-				{
-	    			// gonna zero it out
-					w_new = 0.0;
-
-	    			// decrease the p_resamp
-					p_resamp[basis] -= (p_resamp[basis] - min_resamp) * resamp_decay;
-
-	    			// set the number of non-zero
-					if (w_old == 0.0)
-					{
-	     				// we didn't change
-						changed = 0;
-
-	      				// and wasted a basis
-						wasted_basis += 1;
-					}
-					else
-					{
-	      				// we changed
-						changed = 1;
-
-	      				// must update num non_zero
-						non_zero -= 1;
-					}
-				}
-	  			// process changes if necessary
-				if (changed == 1)
-				{
-	   				// update the expected values
-					double w_diff = w_new - w_old;
-					for (int i = 0; i < fnum; i++)
-					{
-						double E_old = exp(Xw[i].get(0));
-						double val = X[i].get(basis)*w_diff;
-						if(Xw[i].get(0) == 0.0)
-						{
-							if(fabs(val) > 10e-7)
-								Xw[i][0] = val;
-						}
-						else
-							Xw[i][0] += X[i].get(basis)*w_diff;
-						double E_new_m = exp(Xw[i].get(0));			
-						
-						S[i] += E_new_m - E_old;
-					}
-
-	    			// update the weight
-					if( w_new == 0.0)
-					{
-						if(weight.get(basis) != 0.0)
-						{
-							weight.erase(basis);
-						}
-					}
-					else
-						weight[basis] = w_new;
-	    			// keep track of the sqrt sum squared diffs
+        }
+
+        // get the gradient
+        double grad = xY[basis].get ( 0 ) - XdotP;
+        // set the new weight
+        double w_new = w_old + grad / ac[basis];
+
+        // test that we're within bounds
+        if ( w_new > lm_2_ac[basis] )
+        {
+          // more towards bounds, but keep it
+          w_new -= lm_2_ac[basis];
+          changed = 1;
+
+          // umark from being zero if necessary
+          if ( w_old == 0.0 )
+          {
+            non_zero += 1;
+
+            // reset the p_resample
+            p_resamp[basis] = 1.0;
+
+            // we needed the basis
+            needed_basis += 1;
+          }
+        }
+        else if ( w_new < -lm_2_ac[basis] )
+        {
+          // more towards bounds, but keep it
+          w_new += lm_2_ac[basis];
+          changed = 1;
+
+          // umark from being zero if necessary
+          if ( w_old == 0.0 )
+          {
+            non_zero += 1;
+
+            // reset the p_resample
+            p_resamp[basis] = 1.0;
+
+            // we needed the basis
+            needed_basis += 1;
+          }
+        }
+        else
+        {
+          // gonna zero it out
+          w_new = 0.0;
+
+          // decrease the p_resamp
+          p_resamp[basis] -= ( p_resamp[basis] - min_resamp ) * resamp_decay;
+
+          // set the number of non-zero
+          if ( w_old == 0.0 )
+          {
+            // we didn't change
+            changed = 0;
+
+            // and wasted a basis
+            wasted_basis += 1;
+          }
+          else
+          {
+            // we changed
+            changed = 1;
+
+            // must update num non_zero
+            non_zero -= 1;
+          }
+        }
+        // process changes if necessary
+        if ( changed == 1 )
+        {
+          // update the expected values
+          double w_diff = w_new - w_old;
+          for ( int i = 0; i < fnum; i++ )
+          {
+            double E_old = exp ( Xw[i].get ( 0 ) );
+            double val = X[i].get ( basis ) * w_diff;
+            if ( Xw[i].get ( 0 ) == 0.0 )
+            {
+              if ( fabs ( val ) > 10e-7 )
+                Xw[i][0] = val;
+            }
+            else
+              Xw[i][0] += X[i].get ( basis ) * w_diff;
+            double E_new_m = exp ( Xw[i].get ( 0 ) );
+
+            S[i] += E_new_m - E_old;
+          }
+
+          // update the weight
+          if ( w_new == 0.0 )
+          {
+            if ( weight.get ( basis ) != 0.0 )
+            {
+              weight.erase ( basis );
+            }
+          }
+          else
+            weight[basis] = w_new;
+          // keep track of the sqrt sum squared diffs
 //#pragma omp critical
-					sum2_w_diff += w_diff*w_diff;
-				}
-	  			// no matter what we keep track of the old
+          sum2_w_diff += w_diff * w_diff;
+        }
+        // no matter what we keep track of the old
 //#pragma omp critical
-				sum2_w_old += w_old*w_old;
-			}
-		}
+        sum2_w_old += w_old * w_old;
+      }
+    }
 
-    	// finished a cycle, assess convergence
-		incr = sqrt(sum2_w_diff) / (sqrt(sum2_w_old)+numeric_limits<double>::epsilon());
+    // finished a cycle, assess convergence
+    incr = sqrt ( sum2_w_diff ) / ( sqrt ( sum2_w_old ) + numeric_limits<double>::epsilon() );
 #ifdef SLRDEBUG
-		cout << "sum2_w_diff: " << sum2_w_diff << " sum2_w_old " << sum2_w_old << endl;
-		cout << "convcrit: " << incr << " tol: " << convergence_tol << endl;
-		//cout << "sum2_w_wold = " << sum2_w_old << " sum2_w_diff = " << sum2_w_diff << endl;
+    cout << "sum2_w_diff: " << sum2_w_diff << " sum2_w_old " << sum2_w_old << endl;
+    cout << "convcrit: " << incr << " tol: " << convergence_tol << endl;
+    //cout << "sum2_w_wold = " << sum2_w_old << " sum2_w_diff = " << sum2_w_diff << endl;
 #endif
-		if (incr < convergence_tol)
-		{
-      		// we converged!!!
-			break;
-		}
-	}
-
-  	// finished updating weights
-  	// assess convergence
-	cerr << "end regression after " << cycle << " steps"  << endl;
-	return cycle;
+    if ( incr < convergence_tol )
+    {
+      // we converged!!!
+      break;
+    }
+  }
+
+  // finished updating weights
+  // assess convergence
+  cerr << "end regression after " << cycle << " steps"  << endl;
+  return cycle;
 }

classifier/fpclassifier/logisticregression/SLR.h

@@ -1,4 +1,4 @@
-/** 
+/**
  * @file SLR.h
  * @brief implementation of Sparse Logistic Regression (SMLR) Classififier (one vs. all)
  * @author Björn Fröhlich
@@ -20,98 +20,98 @@ namespace OBJREC {
 
 class SLR : public NICE::Persistent
 {
-  //protected:
-	public:
-	//! the configuration file
-	const NICE::Config *conf;
-
-	//! section in the configfile
-	std::string confsection;
-	
-	//! weight vectors
-	SparseVector weight;
-	
-	//! the featurepool
-	FeaturePool fp;
-	
-	//! maximum number of iterations
-	int maxiter;
-	
-	//! Decay rate in the probability of resampling a zero weight.1.0 will immediately decrease to the min_resamp from 1.0, 0.0 will never decrease from 1.0.
-	double resamp_decay;
-	
-	//! convergence criterium for stepwise regression
-	double convergence_tol;
-	
-	//! Minimum resampling probability for zeroed weights"
-	double min_resamp;
-	
-	//! Feature Dimension
-	int fdim;
-	
-	//! The penalty term lambda. Larger values will give rise to	more sparsification
-	double lambda;
-	
-	//! how many samples per class should be used
-	double samplesperclass;
-	
-	//! for normalization
-	std::vector<double> minval, maxval;
-	
-    public:
-  
-	
-	/**
-	 * standard constructor
-	 * @param conf configfile
-	 * @param section section name in configfile for classifier
-	 */
-	SLR( const NICE::Config *conf, std::string section="SMLR");
-      
-	
-	/**
-	 * simple constructor -> does nothing
-	 */
-	SLR ();
-
-	/**
-	 * simple destructor
-	 */
-	~SLR();
-
-	/**
-	 * main classification function
-	 * @param pce input feature
-	 * @return a classification result
-	 */
-	double classify ( Example & pce );
-
-	/**
-	 * start training
-	 * @param fp a featurepool (how to handle which features...)
-	 * @param classno which class should be learned
-	 * @param examples input features
-	 */
-	void train ( FeaturePool & _fp, Examples & examples, int classno = 1 );
-
-	/**
-	 * clone this object
-	 * @return a copy of this object
-	 */
-	FeaturePoolClassifier *clone () const;
-
-	/**
-	 * perform the stepwise regression
-	 * @param x input features
-	 * @param classno which class should be learned
-	 * @return number of iterations
-	 */
-	int stepwise_regression(Examples &x, int classno = 1);
-	
-	/** IO functions */
-	void restore (std::istream & is, int format = 0);
-	void store (std::ostream & os, int format = 0) const;
-	void clear ();
+    //protected:
+  public:
+    //! the configuration file
+    const NICE::Config *conf;
+
+    //! section in the configfile
+    std::string confsection;
+
+    //! weight vectors
+    SparseVector weight;
+
+    //! the featurepool
+    FeaturePool fp;
+
+    //! maximum number of iterations
+    int maxiter;
+
+    //! Decay rate in the probability of resampling a zero weight.1.0 will immediately decrease to the min_resamp from 1.0, 0.0 will never decrease from 1.0.
+    double resamp_decay;
+
+    //! convergence criterium for stepwise regression
+    double convergence_tol;
+
+    //! Minimum resampling probability for zeroed weights"
+    double min_resamp;
+
+    //! Feature Dimension
+    int fdim;
+
+    //! The penalty term lambda. Larger values will give rise to more sparsification
+    double lambda;
+
+    //! how many samples per class should be used
+    double samplesperclass;
+
+    //! for normalization
+    std::vector<double> minval, maxval;
+
+  public:
+
+
+    /**
+     * standard constructor
+     * @param conf configfile
+     * @param section section name in configfile for classifier
+     */
+    SLR ( const NICE::Config *conf, std::string section = "SMLR" );
+
+
+    /**
+     * simple constructor -> does nothing
+     */
+    SLR ();
+
+    /**
+     * simple destructor
+     */
+    ~SLR();
+
+    /**
+     * main classification function
+     * @param pce input feature
+     * @return a classification result
+     */
+    double classify ( Example & pce );
+
+    /**
+     * start training
+     * @param fp a featurepool (how to handle which features...)
+     * @param classno which class should be learned
+     * @param examples input features
+     */
+    void train ( FeaturePool & _fp, Examples & examples, int classno = 1 );
+
+    /**
+     * clone this object
+     * @return a copy of this object
+     */
+    FeaturePoolClassifier *clone () const;
+
+    /**
+     * perform the stepwise regression
+     * @param x input features
+     * @param classno which class should be learned
+     * @return number of iterations
+     */
+    int stepwise_regression ( Examples &x, int classno = 1 );
+
+    /** IO functions */
+    void restore ( std::istream & is, int format = 0 );
+    void store ( std::ostream & os, int format = 0 ) const;
+    void clear ();
 };
 
 } // namespace

math/cluster/GMM.cpp

@@ -17,760 +17,760 @@ using namespace std;
 
 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
-	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 :
-	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
-	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
-		cerr << "GMM::EM: iteration: " << iter << endl;
+    cerr << "GMM::EM: iteration: " << iter << endl;
 #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
-		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
-		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
-		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
-		
-		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
-		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
-	cerr << "GMM::finished EM after " << iter << " iterations" << endl;
+  cerr << "GMM::finished EM after " << iter << " iterations" << endl;
 #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
-	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 )
 {
-	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 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;
 }

math/cluster/GMM.h

@@ -1,12 +1,12 @@
-/** 
+/**
  * @file GMM.h
- * @brief Gaussian Mixture Model based on   
+ * @brief Gaussian Mixture Model based on
   article{Calinon07SMC,
-  title="On Learning, Representing and Generalizing a Task in a Humanoid 
+  title="On Learning, Representing and Generalizing a Task in a Humanoid
   Robot",
   author="S. Calinon and F. Guenter and A. Billard",
-  journal="IEEE Transactions on Systems, Man and Cybernetics, Part B. 
-  Special issue on robot learning by observation, demonstration and 
+  journal="IEEE Transactions on Systems, Man and Cybernetics, Part B.
+  Special issue on robot learning by observation, demonstration and
   imitation",
   year="2007",
   volume="37",
@@ -22,7 +22,7 @@
 
 #include "core/vector/VectorT.h"
 #include "core/vector/MatrixT.h"
- 
+
 #include "vislearning/cbaselib/MultiDataset.h"
 #include "vislearning/cbaselib/LocalizationResult.h"
 #include "vislearning/cbaselib/CachedExample.h"
@@ -35,206 +35,210 @@ namespace OBJREC {
 
 class GMM : public ClusterAlgorithm
 {
-	protected:
-		
-		//! number of gaussians
-		int gaussians;
-		
-		//! dimension of each feature
-		int dim;
-
-		//! mean vectors
-		NICE::VVector mu;
-	
-		//! sparse sigma vectors
-		NICE::VVector sparse_sigma;
-			
-		//! sparse inverse sigma vectors (if usediag)
-		NICE::VVector sparse_inv_sigma;
-		
-		//! save det_sigma for fast computing
-		std::vector<double> log_det_sigma;
-				
-		//! the configfile
-		const NICE::Config *conf;
-		
-		//! the weight for each gaussian
-		std::vector<double> priors;
-		
-		//! parameters for other GM to compare
-		NICE::VVector mu2;
-		NICE::VVector sparse_sigma2;
-		std::vector<double> priors2;
-		bool comp;
-		
-		//! maximum number of iterations for EM
-		int maxiter;
-		
-		//! how many features to use, use 0 for alle input features
-		int featsperclass;
-				
-		//! for faster computing cdimval = dim*2*PI
-		double cdimval;
-		
-		//! parameter for the map estimation
-		double tau;
-		
-		//! whether to use pyramid initialisation or not
-		bool pyramid;
-	public:
-		/**
-	 	 * simplest constructor
-		 */
-		GMM();
-		
-		/**
-		 * simple constructor
-		 * @param _no_classes 
-		 */
-		GMM(int _no_classes);
-		
-		/**
-		 * standard constructor
-		 * @param conf a Configfile
-		 * @param _no_classes number of gaussian
-		 */
-		GMM(const NICE::Config *conf, int _no_classes = -1);
-		
-		/**
-		 * standard destructor
-		 */
-		~GMM(){std::cerr << "dadada" << std::endl;};
-		
-		/**
-		 * computes the mixture
-		 * @param examples the input features
-		 */
-		void computeMixture(Examples examples);
-		
-		/**
-		 * computes the mixture
-		 * @param DataSet the input features
-		 */
-		void computeMixture(const NICE::VVector &DataSet);
-		
-		/**
-		 * returns the probabilities for each gaussian in a sparse vector
-		 * @param vin input vector
-		 * @param probs BoV output vector
-		 */
-		void getProbs(const NICE::Vector &vin, SparseVector &probs);
-		
-		/**
-		 * returns the probabilities for each gaussian
-		 * @param vin input vector
-		 * @param probs BoV output vector
-		 */
-		void getProbs(const NICE::Vector &vin, NICE::Vector &probs);
-		
-		/**
-		 * returns the fisher score for the gmm
-		 * @param vin input vector
-		 * @param probs Fisher score output vector
-		 */
-		void getFisher(const NICE::Vector &vin, SparseVector &probs);
-		
-		/**
-		 *   init the GaussianMixture by selecting randomized mean vectors and using the coovariance of all features
-		 * @param DataSet input Matrix
-		 */
-		void initEM(const NICE::VVector &DataSet);
-				
-		/**
-		 *   alternative for initEM: init the GaussianMixture with a K-Means clustering
-		 * @param DataSet input Matrix
-		 */
-		void initEMkMeans(const NICE::VVector &DataSet);
-		
-		/**
-		 * performs Expectation Maximization on the Dataset, 	in order to obtain a nState GMM Dataset is a Matrix(no_classes,nDimensions)
-		 * @param DataSet input Matrix
-		 * @param gaussians number gaussians to use
-		 * @return number of iterations
-		 */
-		int doEM(const NICE::VVector &DataSet, int nbgaussians);
-			
-		/**
-		 * Compute log probabilty of vector v for the given state. *
-		 * @param Vin 
-		 * @param state 
-		 * @return 
-				 */
-		double logpdfState(const NICE::Vector &Vin,int state);
-		
-		/**
-		 * determine the best mixture for the input feature
-		 * @param v input feature
-		 * @param bprob probability of the best mixture
-		 * @return numer of the best mixture
-		 */
-		int getBestClass(const NICE::Vector &v, double *bprob = NULL);
-		
-		/**
-		 * Cluster a given Set of features and return the labels for each feature
-		 * @param features input features
-		 * @param prototypes mean of the best gaussian
-		 * @param weights weight of the best gaussian
-		 * @param assignment number of the best gaussian
-		 */
-		void cluster ( const NICE::VVector & features, NICE::VVector & prototypes, std::vector<double> & weights, std::vector<int> & assignment );
-		
-		/**
-		 * save GMM data
-		 * @param filename filename
-		 */
-		void saveData(const std::string filename);
-		
-        /**
-		 * load GMM data
-		 * @param filename filename
-		 * @return true if everything works fine
-	 	 */
-		bool loadData(const std::string filename);
-		
-		/**
-		 * return the parameter of the mixture
-		 * @param mu 
-		 * @param sSigma 
-		 * @param p 
-		 */
-		void getParams(NICE::VVector &mean, NICE::VVector &sSigma, std::vector<double> &p);
-		
-		/**
-		 * Set the parameters of an other mixture for comparing with this one
-		 * @param mean mean vectors
-		 * @param sSigma diagonal covariance Matrixs
-		 * @param p weights
-		 */
-		void setCompareGM(NICE::VVector mean, NICE::VVector sSigma, std::vector<double> p);
-		
-		/**
-		 * probability product kernel
-		 * @param sigma1 
-		 * @param sigma2 
-		 * @param mu1 
-		 * @param mu2 
-		 * @param p 
-		 * @return 
-		 */
-		double kPPK(NICE::Vector sigma1, NICE::Vector sigma2, NICE::Vector mu1, NICE::Vector mu2, double p);
-		
-		/**
-		 * starts a comparison between this Mixture and a other one seted bei "setComparGM"
-		 */
-		double compare();
-		
-		/**
-		 * whether to compare or not
-		 * @param c 
-		 */
-		void comparing(bool c = true);
-		
-		int getSize(){return gaussians;}
+  protected:
+
+    //! number of gaussians
+    int gaussians;
+
+    //! dimension of each feature
+    int dim;
+
+    //! mean vectors
+    NICE::VVector mu;
+
+    //! sparse sigma vectors
+    NICE::VVector sparse_sigma;
+
+    //! sparse inverse sigma vectors (if usediag)
+    NICE::VVector sparse_inv_sigma;
+
+    //! save det_sigma for fast computing
+    std::vector<double> log_det_sigma;
+
+    //! the configfile
+    const NICE::Config *conf;
+
+    //! the weight for each gaussian
+    std::vector<double> priors;
+
+    //! parameters for other GM to compare
+    NICE::VVector mu2;
+    NICE::VVector sparse_sigma2;
+    std::vector<double> priors2;
+    bool comp;
+
+    //! maximum number of iterations for EM
+    int maxiter;
+
+    //! how many features to use, use 0 for alle input features
+    int featsperclass;
+
+    //! for faster computing cdimval = dim*2*PI
+    double cdimval;
+
+    //! parameter for the map estimation
+    double tau;
+
+    //! whether to use pyramid initialisation or not
+    bool pyramid;
+  public:
+    /**
+      * simplest constructor
+     */
+    GMM();
+
+    /**
+     * simple constructor
+     * @param _no_classes
+     */
+    GMM ( int _no_classes );
+
+    /**
+     * standard constructor
+     * @param conf a Configfile
+     * @param _no_classes number of gaussian
+     */
+    GMM ( const NICE::Config *conf, int _no_classes = -1 );
+
+    /**
+     * standard destructor
+     */
+    ~GMM() {
+      std::cerr << "dadada" << std::endl;
+    };
+
+    /**
+     * computes the mixture
+     * @param examples the input features
+     */
+    void computeMixture ( Examples examples );
+
+    /**
+     * computes the mixture
+     * @param DataSet the input features
+     */
+    void computeMixture ( const NICE::VVector &DataSet );
+
+    /**
+     * returns the probabilities for each gaussian in a sparse vector
+     * @param vin input vector
+     * @param probs BoV output vector
+     */
+    void getProbs ( const NICE::Vector &vin, SparseVector &probs );
+
+    /**
+     * returns the probabilities for each gaussian
+     * @param vin input vector
+     * @param probs BoV output vector
+     */
+    void getProbs ( const NICE::Vector &vin, NICE::Vector &probs );
+
+    /**
+     * returns the fisher score for the gmm
+     * @param vin input vector
+     * @param probs Fisher score output vector
+     */
+    void getFisher ( const NICE::Vector &vin, SparseVector &probs );
+
+    /**
+     *   init the GaussianMixture by selecting randomized mean vectors and using the coovariance of all features
+     * @param DataSet input Matrix
+     */
+    void initEM ( const NICE::VVector &DataSet );
+
+    /**
+     *   alternative for initEM: init the GaussianMixture with a K-Means clustering
+     * @param DataSet input Matrix
+     */
+    void initEMkMeans ( const NICE::VVector &DataSet );
+
+    /**
+     * performs Expectation Maximization on the Dataset,  in order to obtain a nState GMM Dataset is a Matrix(no_classes,nDimensions)
+     * @param DataSet input Matrix
+     * @param gaussians number gaussians to use
+     * @return number of iterations
+     */
+    int doEM ( const NICE::VVector &DataSet, int nbgaussians );
+
+    /**
+     * Compute log probabilty of vector v for the given state. *
+     * @param Vin
+     * @param state
+     * @return
+       */
+    double logpdfState ( const NICE::Vector &Vin, int state );
+
+    /**
+     * determine the best mixture for the input feature
+     * @param v input feature
+     * @param bprob probability of the best mixture
+     * @return numer of the best mixture
+     */
+    int getBestClass ( const NICE::Vector &v, double *bprob = NULL );
+
+    /**
+     * Cluster a given Set of features and return the labels for each feature
+     * @param features input features
+     * @param prototypes mean of the best gaussian
+     * @param weights weight of the best gaussian
+     * @param assignment number of the best gaussian
+     */
+    void cluster ( const NICE::VVector & features, NICE::VVector & prototypes, std::vector<double> & weights, std::vector<int> & assignment );
+
+    /**
+     * save GMM data
+     * @param filename filename
+     */
+    void saveData ( const std::string filename );
+
+    /**
+    * load GMM data
+    * @param filename filename
+    * @return true if everything works fine
+    */
+    bool loadData ( const std::string filename );
+
+    /**
+     * return the parameter of the mixture
+     * @param mu
+     * @param sSigma
+     * @param p
+     */
+    void getParams ( NICE::VVector &mean, NICE::VVector &sSigma, std::vector<double> &p );
+
+    /**
+     * Set the parameters of an other mixture for comparing with this one
+     * @param mean mean vectors
+     * @param sSigma diagonal covariance Matrixs
+     * @param p weights
+     */
+    void setCompareGM ( NICE::VVector mean, NICE::VVector sSigma, std::vector<double> p );
+
+    /**
+     * probability product kernel
+     * @param sigma1
+     * @param sigma2
+     * @param mu1
+     * @param mu2
+     * @param p
+     * @return
+     */
+    double kPPK ( NICE::Vector sigma1, NICE::Vector sigma2, NICE::Vector mu1, NICE::Vector mu2, double p );
+
+    /**
+     * starts a comparison between this Mixture and a other one seted bei "setComparGM"
+     */
+    double compare();
+
+    /**
+     * whether to compare or not
+     * @param c
+     */
+    void comparing ( bool c = true );
+
+    int getSize() {
+      return gaussians;
+    }
 };