matrix.h

这是一款用C++编写的实现gmm算法的程序

/*
  Copyright (c) 2008 Eric Sauser, 
  LASA Lab, EPFL, CH-1015 Lausanne, Switzerland, 
  http://lasa.epfl.ch
*/

#ifndef MATRIX_H
#define MATRIX_H

#ifdef WIN32
#include "windows.h"
#endif

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>

#include "Macros.h"
#include "Vector.h"

#ifdef  USE_T_EXTENSIONS
template<unsigned int ROW> class TMatrix;
#endif
    
class Matrix
{
  friend class Vector;
#ifdef  USE_T_EXTENSIONS
  template<unsigned int ROW> friend class TMatrix;
#endif
  
protected: 
  static int bInverseOk;
  
  unsigned int  row;
  unsigned int  column;
  float        *_;

public:

  inline Matrix() {
    row    = 0;
    column = 0;
    _      = NULL;
  }
  
  inline virtual ~Matrix(){
    Release(); 
  }

  inline Matrix(const Matrix &matrix)
  {
    row    = 0;
    column = 0;
    _      = NULL;
    Resize(matrix.row,matrix.column,false);
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        _[j*column+i] = matrix._[j*column+i];
  }

  inline Matrix(unsigned int rowSize, unsigned int colSize, bool clear = true)
  {
    row    = 0;
    column = 0;
    _      = NULL;
    Resize(rowSize,colSize,false);
    if(clear)
      Zero();
  }
  
  inline Matrix(const float _[], unsigned int rowSize, unsigned int colSize)
  {
    row       = 0;
    column    = 0;
    this->_   = NULL;
    Resize(rowSize,colSize,false);
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        this->_[j*column+i] = _[j*column+i];
  }

#ifdef  USE_T_EXTENSIONS
  template<unsigned int ROW> inline Matrix(const TMatrix<ROW> &matrix)
  {
    row    = 0;
    column = 0;
    _      = NULL;
    Resize(ROW,ROW,false);
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        _[j*column+i] = matrix._[j*column+i];
  }
#endif
   
  inline Matrix& Zero()
  {
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        _[j*column+i] = 0.0f;
    return *this;
  }
    
  inline unsigned int RowSize() const{
    return row;
  }
  inline unsigned int ColumnSize() const{
    return column;
  } 
  inline float *Array() const{
    return _;
  }

  inline float& operator() (const unsigned int row, const unsigned int col)
  {
    if((row<this->row)&&(col<this->column))
      return _[row*column+col];
    return Vector::undef; 
  }

  inline Vector GetRow(const unsigned int row) const
  {
    Vector result(column,false);    
    return GetRow(row,result);     
  }

  inline Vector& GetRow(const unsigned int row, Vector& result) const
  {
    result.Resize(column,false);
    for (unsigned int i = 0; i < column; i++)
      result._[i] = _[row*column+i];
    return result;     
  }

  inline Vector GetColumn(const unsigned int col) const
  {
    Vector result(row,false);    
    return GetColumn(col,result);     
  }

  inline Vector& GetColumn(const unsigned int col, Vector& result) const
  {
    result.Resize(row,false);
    if(col<column){
      for (unsigned int j = 0; j < row; j++)
        result._[j] = _[j*column+col];
    }else{
      result.Zero();
    }
    return result;     
  }

  inline Matrix GetColumnSpace(const unsigned int col, const unsigned int len) const  
  {
    if(len>0){
      Matrix result(row,len,false);    
      return GetColumnSpace(col,len,result);
    }else
      return Matrix();     
  }

  inline Matrix GetRowSpace(const unsigned int row, const unsigned int len) const
  {
    if(len>0){
      Matrix result(len,column,false);    
      return GetRowSpace(row,len,result);
    }else
      return Matrix();     
  }


  inline Matrix& GetColumnSpace(const unsigned int col, const unsigned int len, Matrix &result) const
  {
    if(len>0){
      const unsigned int end  = col+len-1;    
      const unsigned int size = len; 
      result.Resize(row,size,false);
      
      if(col<column){
        const unsigned int k = (end+1<=column?end+1:column);  
        
        for (unsigned int i = col; i < k; i++)
          for (unsigned int j = 0; j < row; j++)
            result._[j*size+(i-col)] = _[j*column+i];
        for (unsigned int i = k; i < end+1; i++)
          for (unsigned int j = 0; j < row; j++)
           result._[j*size+(i-col)] = 0.0f;            
      }else{
        result.Zero();
      }
    }else{
      result.Resize(0,0,false);
    }
    return result;     
  }

  inline Matrix& GetRowSpace(const unsigned int row, const unsigned int len, Matrix &result) const
  {      
    if(len>0){
      const unsigned int end  = row+len-1;
      const unsigned int size = len; 
      result.Resize(size,column,false);
      
      if(row<this->row){
        const unsigned int k = (end+1<=this->row?end+1:this->row);  
        
        for (unsigned int j = 0; j < column; j++)
          for (unsigned int i = row; i < k; i++)
            result._[(i-row)*column+j] = _[i*column+j];
        for (unsigned int j = 0; j < column; j++)
          for (unsigned int i = k; i < end+1; i++)
           result._[(i-row)*column+j] = 0.0f;            
      }else{
        result.Zero();
      }
    }else{
      result.Resize(0,0,false);
    }
    return result;     
  }

  inline Matrix& SetRow(const Vector &vector, const unsigned int row)
  {
    if(row<this->row){    
      const unsigned int ki = (column<=vector.row?column:vector.row);
      for (unsigned int i = 0; i < ki; i++)
        _[row*column+i] = vector._[i]; 
    }
    return *this;     
  }

  inline Matrix& SetColumn(const Vector &vector, const unsigned int col)
  {
    if(col<this->column){    
      const unsigned int kj = (row<=vector.row?row:vector.row);
      for (unsigned int j = 0; j < kj; j++)
        _[j*column+col] = vector._[j];
    }
    return *this;
  }

  inline Matrix& SetColumnSpace(const Matrix &matrix, const unsigned int col)
  {
    if(col<this->column){    
      const unsigned int kj = (row<=matrix.row?row:matrix.row);
      const unsigned int ki = (col+matrix.column<=this->column?col+matrix.column:this->column);
      for (unsigned int j = 0; j < kj; j++)
        for (unsigned int i = col; i < ki; i++)
          _[j*column+i] = matrix._[j*matrix.column+(i-col)];
    }
    return *this;
  }

  inline Matrix& SetRowSpace(const Matrix &matrix, const unsigned int row)
  {
    if(row<this->row){
      const unsigned int ki = (column<=matrix.column?column:matrix.column);
      const unsigned int kj = (row+matrix.row<=this->row?row+matrix.row:this->row);
      for (unsigned int j = row; j < kj; j++)
        for (unsigned int i = 0; i < ki; i++)
          _[j*column+i] = matrix._[(j-row)*matrix.column+i]; 
    }
    return *this;     
  }

  inline Matrix GetRowSpace(const Vector &ids) const
  {
    Matrix result(ids.Size(),column);
    return GetRowSpace(ids,result);
  }

  inline Matrix GetColumnSpace(const Vector &ids) const
  {
    Matrix result(row,ids.Size());
    return GetColumnSpace(ids,result);
  }

  inline Matrix GetMatrixSpace(const Vector &rowIds,const Vector &colIds) const
  {
    Matrix result(rowIds.Size(),colIds.Size());
    return GetMatrixSpace(rowIds,colIds,result);
  }

  inline Matrix& GetColumnSpace(const Vector &ids, Matrix &result) const
  {
    const unsigned int k=ids.Size();
    result.Resize(row,k);
    for(unsigned int i=0;i<k;i++){
      const unsigned int g = (unsigned int)(fabs(ROUND(ids._[i])));
      if(g<column){
        for(unsigned int j=0;j<row;j++)
          result._[j*k+i] = _[j*column+g];
      }else{
        for(unsigned int j=0;j<row;j++)
          result._[j*k+i] = 0.0f;        
      }
    }
    return result;     
  }

  inline Matrix& GetRowSpace(const Vector &ids, Matrix &result) const
  {
    const unsigned int k=ids.Size();
    result.Resize(k,column);
    for(unsigned int i=0;i<k;i++){
      const unsigned int g = (unsigned int)(fabs(ROUND(ids._[i])));
      if(g<row){
        for(unsigned int j=0;j<column;j++)
          result._[i*column+j] = _[g*column+j];
      }else{
        for(unsigned int j=0;j<column;j++)
          result._[i*column+j] = 0.0f;
      }
    }
    return result;     
  }

  inline Matrix& GetMatrixSpace(const Vector &rowIds,const Vector &colIds, Matrix &result) const
  {
    const unsigned int k1=rowIds.Size();
    const unsigned int k2=colIds.Size();
    result.Resize(k1,k2);
    for(unsigned int i=0;i<k1;i++){
      const unsigned int g1 = (unsigned int)(fabs(ROUND(rowIds._[i])));
      if(g1<row){
        for(unsigned int j=0;j<k2;j++){      
          const unsigned int g2 = (unsigned int)(fabs(ROUND(colIds._[j])));
          if(g2<column){
            result._[i*k2+j] = _[g1*column+g2];            
          }else{
            result._[i*k2+j] = 0.0f;
          }
        }
      }else{
        for(unsigned int j=0;j<k2;j++)
          result._[i*k2+j] = 0.0f;
      }
    }
    return result;     
  }

  inline Matrix operator - () const
  {
    Matrix result(row,column,false);
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        result._[j*column+i] = -_[j*column+i];
    return result;
  }
  
  inline virtual Matrix& operator = (const Matrix &matrix)
  {
    Resize(matrix.row,matrix.column,false);
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        _[j*column+i] = matrix._[j*column+i];
    return *this;    
  }

  inline Matrix& operator += (const Matrix &matrix)
  {
    const unsigned int kj = (row<=matrix.row?row:matrix.row);
    const unsigned int ki = (column<=matrix.column?column:matrix.column);
    for (unsigned int j = 0; j < kj; j++)
      for (unsigned int i = 0; i < ki; i++)
        _[j*column+i] += matrix._[j*column+i];
    return *this;
  }

  inline Matrix& operator -= (const Matrix &matrix)
  {
    const unsigned int kj = (row<=matrix.row?row:matrix.row);
    const unsigned int ki = (column<=matrix.column?column:matrix.column);
    for (unsigned int j = 0; j < kj; j++)
      for (unsigned int i = 0; i < ki; i++)
        _[j*column+i] -= matrix._[j*column+i];
    return *this;
  }

  inline Matrix& operator ^= (const Matrix &matrix)
  {
    const unsigned int kj = (row<=matrix.row?row:matrix.row);
    const unsigned int ki = (column<=matrix.column?column:matrix.column);
    for (unsigned int j = 0; j < kj; j++)
      for (unsigned int i = 0; i < ki; i++)
        _[j*column+i] *= matrix._[j*column+i];
    return *this;
  }

  inline Matrix& operator /= (const Matrix &matrix)
  {
    const unsigned int kj = (row<=matrix.row?row:matrix.row);
    const unsigned int ki = (column<=matrix.column?column:matrix.column);
    for (unsigned int j = 0; j < kj; j++)
      for (unsigned int i = 0; i < ki; i++)
        _[j*column+i] /= matrix._[j*column+i];
    return *this;
  }

  inline Matrix& operator += (float scalar)
  {
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        _[j*column+i] += scalar;
    return *this;
  }

  inline Matrix& operator -= (float scalar)
  {
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)