diff --git a/matrix.hpp b/matrix.hpp
deleted file mode 100644
index 89b77d85..00000000
--- a/matrix.hpp
+++ /dev/null
@@ -1,560 +0,0 @@
-/***************************************************************************
- *  $Id: matrix.hpp,v 1.11 2006/05/13 10:05:53 nojhan Exp $
- *  Copyright : Free Software Foundation
- *  Author : Johann Dréo <nojhan@gmail.com>
- ****************************************************************************/
-
-/*
- *  This program is free software; you can redistribute it and/or modify
- *  it under the terms of the GNU Lesser General Public License as published by
- *  the Free Software Foundation; either version 2 of the License, or
- *  (at your option) any later version.
- *
- *  This program is distributed in the hope that it will be useful,
- *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *  GNU Lesser General Public License for more details.
- *
- *  You should have received a copy of the GNU Lesser General Public License
- *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
- */
- 
-#ifndef MATRIX
-#define MATRIX
-
-#include <vector>
-#include <sstream>
-
-#include "Exception_oMetah.hpp"
-
-using namespace std;
-
-namespace ometah {
-
-//! Test if a vector is comprised in bounds
-template<class T>
-bool isInBounds( T aVector, T mins, T maxs)
-{
-  unsigned int i;
-  for(i=0; i<aVector.size(); i++ ){
-    // too low
-    if( aVector[i] < mins[i] ){
-      return false;
-    // too high
-    }else if( aVector[i] > maxs[i] ){
-      return false;
-    }
-  }
-  return true;
-}
-
-
-//! Force a vector to be in bounds
-template<class T>
-T forceBounds( T aVector, T mins, T maxs)
-{
-  T CastedVector=aVector;
-    
-  unsigned int i;
-  for(i=0; i<aVector.size(); i++ ){
-    
-    if( aVector[i] < mins[i] ){
-      CastedVector[i]=mins[i];
-            
-    }else if( aVector[i] > maxs[i] ){
-      CastedVector[i]=maxs[i];
-    }
-  }
-  return CastedVector;
-}
-
-//! Create a 2D matrix filled with values
-/* 
-    if we want a vector<vector<double> > :
-    T stand for double
-    V stand for vector<vector<double> >
-*/
-template <class T, class U>
-U matrixFilled( unsigned int dimL, unsigned int dimC, T fillValue )
-{
-  unsigned int i;
-
-  // make the vector<double> possible at this step
-  typename U::value_type vec(dimC, fillValue);
-
-  U mat;
-  for(i=0; i<dimL; i++){
-    mat.push_back(vec);
-  }
-
-  return mat;
-}
-
-
-template <class T>
-vector<vector< T > > matrixFilled( unsigned int dimL, unsigned int dimC, T fillValue )
-{
-  unsigned int i;
-
-  // make the vector<double> possible at this step
-  vector< T > vec(dimC, fillValue);
-
-  vector<vector< T > >  mat;
-  for(i=0; i<dimL; i++){
-    mat.push_back(vec);
-  }
-
-  return mat;
-}
-
-
-//! Multipliate two 2D matrix
-template<class T>
-T multiply( T matA, T matB)
-{
-
-  T newMat;
-
-  unsigned int Al=matA.size();
-  unsigned int Ac=matA[0].size();
-  unsigned int Bl=matB.size();
-  unsigned int Bc=matB[0].size();
-
-  newMat=matrixFilled( Al,Bc,0.0);
-
-  if(Ac!=Bl) {
-    throw Exception_Size_Match("Cannot multiply matrices, sizes does not match", EXCEPTION_INFOS );
-  }
-
-  for( unsigned int i=0; i<Al; i++ ) {
-    for( unsigned int j=0; j<Bc; j++ ) {
-      for( unsigned int k=0; k<Ac ;k++ ) {
-         newMat[i][j] += matA[i][k]*matB[k][j];
-      }
-    }
-  }
-  return newMat;
-}
-
-//! Multiply each term of a vector by a scalar
-template<class T, class U>
-U multiply(U aVector, T aNb)
-{
-  U res;
-
-  res.reserve( aVector.size() );
-  
-  unsigned int i;
-  for(i=0; i<aVector.size(); i++){
-    double x=aVector[i]*aNb;
-    res.push_back(x);
-  }
-  return res;
-}
-
-
-//! Cholesky factorization
-template<class T>
-T cholesky( T A)
-{
-
-  // FIXME : vérifier que A est symétrique définie positive
-  
-  T B;
-  unsigned int Al=A.size();
-  unsigned int Ac=A[0].size();
-  B = matrixFilled(Al, Ac, 0.0);
-
-  unsigned int i,j,k;
-  
-  // first column
-  i=0;
-
-  // diagonal
-  j=0;
-  B[0][0]=sqrt(A[0][0]);
-  
-  // end of the column
-  for(j=1;j<Ac;j++) {
-    B[j][0] = A[0][j] / B[0][0];
-  }
-
-  // end of the matrix
-  for(i=1;i<Al;i++){ // each column
-
-    // diagonal
-    double sum=0.0;
-    for(k=0; k<i; k++) {
-      sum += B[i][k]*B[i][k];
-    }
-    
-    // Check for math error
-    if( (A[i][i]-sum) <= 0 ) {
-        ostringstream msg;
-        msg << "Error: Cannot compute the Cholesky decomposition, matrix may not be positive definite (A[";
-        msg << i << "][" << i << "]-sum(B[i][k]^2) = " << A[i][i]-sum << ").";
-        throw Exception_Math(msg.str(), EXCEPTION_INFOS );
-    }
-    
-    B[i][i] = sqrt( A[i][i] - sum );
-
-    for(j=i+1;j<Al;j++){ // rows
-
-      // one element
-      sum = 0.0;
-      for(k=0; k<i; k++) {
-          sum += B[j][k]*B[i][k];
-      }
-      B[j][i] = (A[j][i] - sum) / B[i][i];
-      
-    }
-  }
-  return B;
-}
-
-//! Transposition of a matrix
-template<class T>
-T transpose( T &mat)
-{
-  unsigned int iSize=mat.size();
-  unsigned int jSize=mat[0].size();
-
-  if ( iSize == 0 || jSize == 0 ) { 
-    ostringstream msg;
-    msg << "ErrorSize: matrix not defined "
-        << "(iSize:" << iSize << ", jSize:" << jSize << ")";
-    throw Exception_Size( msg.str(), EXCEPTION_INFOS );
-  }
-
-  typename T::value_type  aVector;
-  T newMat;
-
-  unsigned int i, j;
-
-  for (j=0; j<jSize; j++) {
-    for(i=0; i<iSize; i++) {
-      
-      if ( mat[i].size() != jSize ) { 
-        ostringstream msg;
-        msg << "ErrorSize: matrix not defined "
-            << "(iSize:" << iSize << ", jSize:" << jSize << ", matrix[" << i << "].size:" << mat[i].size() << ")";
-        throw Exception_Size(msg.str(), EXCEPTION_INFOS );
-      }
-
-      aVector.push_back(mat[i][j]);
-    }//j
-
-    newMat.push_back(aVector);
-    aVector.clear();
-  }//i
-
-  return newMat;
-}
-
-//! Calculate the mean vector of a matrix
-template<class T>
-vector<T> mean( vector<vector<T> > mat)
-{
-  vector<T> moyDim;
-  moyDim.reserve(mat.size());
-
-  unsigned int i,a; 
-  a=mat.size();
-
-  for(i=0;i<a;i++) {
-    moyDim.push_back( mean(mat[i]) ); 
-  }
-  return moyDim;
-}
-
-//! Calculate the mean of a vector
-template<class T>
-T mean( vector<T> aVector, unsigned int begin=0, unsigned int during=0)
-{
-  if (during==0) {
-    during = aVector.size() - begin; // if no end : take all
-  }
-
-  T aSum, aMean;
-
-  aSum = sum(aVector, begin, during); // Sum
-  aMean = aSum / (during - begin); // Mean
-
-  return aMean;
-}
-
-//! Calculate a variance-covariance matrix from a list of vector
-/*!
-    For a population of p points on n dimensions :
-        if onRow==true, the matrix should have p rows and n columns.
-        if onRow==false, the matrix should have n rows and p columns.
-*/
-template<class U>
-U varianceCovariance( U pop, bool onRow = true)
-{
-/*
-  // vector of means
-  typename U::value_type  vecMeanCentered;
-  if(onRow) {
-    vecMeanCentered = mean( transpose(pop) ); // p rows and n columns => means of p
-  } else {
-    vecMeanCentered = mean( pop ); // n rows and p columns => means of n
-  }
-  
-  // centered population 
-  // same size as the initial matrix
-  U popMeanCentered = matrixFilled(pop.size(),pop[0].size(), 0.0);
-
-  // centering
-  // rows
-  for(unsigned int i=0;i<pop.size();i++) {
-    // columns
-    for(unsigned int j=0;j<pop[i].size();j++) {
-      popMeanCentered[i][j] = (pop[i][j] - vecMeanCentered[j]);
-    }
-  }
-*/
-  // no centering
-  U popMeanCentered = pop;
-
-  // transposition of the centered matrix
-  U popMeanCenteredT;
-  popMeanCenteredT = transpose(popMeanCentered);
-
-  // final variance/covariance matrix
-  U popVar;
-  if(onRow) {
-      popVar = multiply( popMeanCenteredT, popMeanCentered ); // if p rows and n columns => covariance of p
-  } else {
-      popVar = multiply( popMeanCentered, popMeanCenteredT ); // if n rows and p columns => covariance of n
-  }
-
-  // multiplication by 1/n :
-  for(unsigned int i=0;i<popVar.size();i++) {
-    for(unsigned int j=0;j<popVar[i].size();j++) {
-      popVar[i][j]=popVar[i][j]/(pop.size());
-    }
-  }
-
-  return popVar;
-}
-
-//! Calculate the sum of a vector
-template<class T>
-T sum(vector<T> aVector, unsigned int begin=0, unsigned int during=0)
-{
-  if ( begin > aVector.size() || during > aVector.size() ) {
-    ostringstream msg;
-    msg << "ErrorSize: parameters are out of vector bounds "
-        << "(begin:" << begin << ", during:" << during
-        << ", size:" << aVector.size() << ")";
-    throw Exception_Size_Index( msg.str(), EXCEPTION_INFOS );
-  }
-
-  if (during==0) {
-    during = aVector.size() - begin;
-  }
-
-  T aSum=0;
-
-  for (unsigned int j=begin; j<during; j++) {
-    aSum = aSum + aVector[j]; // sum
-  }//for (j)
-
-  return aSum;
-
-}
-
-//! Calculate the standard deviation of a vector
-template<class T>
-T stdev(vector<T> aVector, unsigned int begin=0, unsigned int during=0)
-{
-  if ( begin > aVector.size() || during > aVector.size() ) {
-    ostringstream msg;
-    msg << "ErrorSize: parameters are out of vector bounds "
-        << "(begin:" << begin << ", during:" << during
-        << ", size:" << aVector.size() << ")";
-    throw Exception_Size_Index( msg.str(), EXCEPTION_INFOS );
-  }
-
-  if (during==0) {
-    during = aVector.size() - begin;
-  }
-
-  vector<T> deviation;
-  T aMean, aDev, aStd;
-
-  aMean = mean(aVector, begin, during); // mean
-
-  for (unsigned int j=begin; j<during; j++) {
-    aDev = aMean - aVector[j];
-    deviation.push_back(aDev*aDev);
-  }//for (j)
-
-  aStd = sqrt( mean(deviation, begin, during) );
-
-  return aStd;
-
-}
-
-//! Find the minimum value of a vector
-template<class T>
-typename T::value_type min(T aVector, unsigned int begin=0, unsigned int during=0)
-{
-  if ( begin > aVector.size() || during > aVector.size() ) {
-    ostringstream msg;
-    msg << "ErrorSize: parameters are out of vector bounds "
-        << "(begin:" << begin << ", during:" << during
-        << ", size:" << aVector.size() << ")";
-    throw Exception_Size_Index( msg.str(), EXCEPTION_INFOS );
-  }
-
-  if (during==0) {
-    during = aVector.size() - begin;
-  }
-
-  typename T::value_type aMin = aVector[begin];
-
-  for (unsigned int i=begin+1; i<during; i++) {
-    if ( aVector[i] < aMin ) {
-      aMin = aVector[i];
-    }
-  }
-  
-  return aMin;
-}
-
-//! Find the minimums values of a matrix, for each row
-template<class T>
-vector<T> mins(vector<vector< T > > aMatrix) 
-{
-  vector<T> mins;
-
-  for( unsigned int i=0; i < aMatrix.size(); i++ ) {
-    mins.push_back( min(aMatrix[i]) );
-  }
-  
-  return mins;
-}
-
-//! Find the maximums values of a matrix, for each row
-template<class T>
-vector<T> maxs(vector<vector< T > > aMatrix) 
-{
-  vector<T> maxs;
-
-  for( unsigned int i=0; i < aMatrix.size(); i++ ) {
-    maxs.push_back( max(aMatrix[i]) );
-  }
-  
-  return maxs;
-}
-
-//! Find the maximum value of a vector
-template<class T>
-typename T::value_type max(T aVector, unsigned int begin=0, unsigned int during=0)
-{
-  if ( begin > aVector.size() || during > aVector.size() ) {
-    ostringstream msg;
-    msg << "ErrorSize: parameters are out of vector bounds "
-        << "(begin:" << begin << ", during:" << during
-        << ", size:" << aVector.size() << ")";
-    throw Exception_Size_Index( msg.str(), EXCEPTION_INFOS );
-  }
-
-  if (during==0) {
-    during = aVector.size() - begin;
-  }
-
-  typename T::value_type aMax = aVector[begin];
-
-  for (unsigned int i=begin+1; i<during; i++) {
-    if ( aVector[i] > aMax ) {
-      aMax = aVector[i];
-    }
-  }
-  
-  return aMax;
-}
-
-//! Substraction of two vectors, terms by terms
-template<class T>
-T substraction(T from, T that)
-{
-  T res;
-  
-  res.reserve(from.size());
-
-  for(unsigned int i=0; i<from.size(); i++){
-    res.push_back( from[i]-that[i] );
-  }
-  return res;
-}
-
-//! Addition of two vectors, terms by terms
-template<class T>
-T addition(T from, T that)
-{
-  T res;
-
-  res.reserve( from.size() );
-
-  for(unsigned int i=0; i<from.size(); i++){
-    res.push_back( from[i]+that[i] );
-  }
-  return res;
-}
-
-//! Calculate the absolute values of a vector
-template<class T>
-T absolute(T aVector)
-{
-  for(unsigned int i=0; i<aVector.size(); i++){
-    aVector[i] = abs(aVector[i]);
-  }
-    
-  return aVector;
-}
-
-
-template<class T>
-vector<T> gravityCenter( vector<vector<T> > points, vector<T> weights )
-{
-
-  // if we have only one weight, we use it for all items
-  if ( weights.size() == 1 ) {
-    for ( unsigned int i=1; i < points.size(); i++ ) {
-      weights.push_back( weights[0] );
-    }
-  }
-
-  // if sizes does not match : error
-  if ( points.size() != weights.size() ) {
-    ostringstream msg;
-    msg << "ErrorSize: "
-        << "points size (" << points.size() << ")" 
-        << " does not match weights size (" << weights.size() << ")";
-    throw Exception_Size_Match( msg.str(), EXCEPTION_INFOS );
-  }
-
-  T weightsSum = sum(weights);
-
-  vector<vector< T > > pointsT = transpose( points );
-
-  vector<T> gravity;
-
-  for ( unsigned int i=0; i < pointsT.size(); i++ ) { // dimensions
-    T g = 0;
-    for ( unsigned int j=0; j < pointsT[i].size(); j++ ) { // points
-      g += ( pointsT[i][j] * weights[j] ) / weightsSum;
-    }
-    gravity.push_back( g );
-  }
-  
-  return gravity;
-}
-
-} // ometah
-
-#endif // MATRIX