- useless comments and files removed
This commit is contained in:
Caner Candan
parent
91fa58ccaf
commit
a70a083735
1 changed files with 0 additions and 560 deletions
560
matrix.hpp
560
matrix.hpp
|
|
@ -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
|
||||
Reference in a new issue