tmatrix.h

Go to the documentation of this file.

 
#ifndef _T_MATRIX_H
#define _T_MATRIX_H
 
#ifdef HAS_GSL
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#endif
 
#include <string.h>
#include <sstream>
#include <string>
#include <assert.h>
#include <vector>
#include "output.h"
 
using namespace std;
 
class TMatrix {
private:
  
  unsigned int _rows, _cols, _length;
  
  double* _val;
  
public:
        
  TMatrix  () : _rows(0), _cols(0), _length(0), _val(0) { }
 
  TMatrix  (const TMatrix& mat) : _rows(0), _cols(0), _length(0), _val(0)
  {
    copy(mat);
  } 
 
  TMatrix ( unsigned int rows, unsigned int cols ) : _rows(0), _cols(0), _length(0), _val(0)
  {
    _length = rows*cols;
    _val = new double [_length];
    _rows = rows; _cols = cols;
  }
  
  ~TMatrix () {if(_val != NULL) delete [] _val;}
  
  void copy (const TMatrix& mat)
  {
    _rows = mat._rows;
    _cols = mat._cols;
    _length = _rows*_cols;
    if(_val) delete [] _val;
    _val = new double [_length];
    memcpy(_val,mat._val,_length*sizeof(double));
  }
 
  void copy_recycle (const TMatrix& mat)
  {
    if(_rows == mat._rows && _cols == mat._cols)
      copy(mat);
    else if(_val == NULL || _length == 0)
      error("TMatrix::copy_recycle::matrix must be allocated before call\n");
    else {
      for(unsigned int i = 0; i < _rows; ++i)
         for(unsigned int j = 0; j < _cols; ++j)
           set(i, j, mat._val[ (i%mat._rows)*mat._cols + (j%mat._cols) ]); //mat.get( i%nrow, j%ncol )); _val[i*_cols + j]
    }
  }
  void set (unsigned int i, unsigned int j, double val) {
    if( i*j < _length)
      _val[i*_cols + j] = val;
    else
      error("TMatrix::set overflow, i*j > length!\n");
  }
 
  void set_row (unsigned int i, double val) {
    if( i < _rows)
      for(unsigned int j = 0, stride = i*_cols; j < _cols; ++j)
        _val[stride + j] = val;
    else
      error("TMatrix::set_row overflow, i > num rows!\n");
  }
 
  void set_row (unsigned int i, vector<double> vec) {
 
    if(vec.size() < _cols)
      fatal("TMatrix::set_row copy from vector too small to fit\n");
 
    if( i < _rows)
      for(unsigned int j = 0, stride = i*_cols; j < _cols; ++j)
        _val[stride + j] = vec[j];
    else
      error("TMatrix::set_row overflow, i > num rows!\n");
  }
 
  void set_col (unsigned int i, double val) {
    if( i < _cols)
      for(unsigned int j = 0; j < _rows; ++j)
        _val[j*_cols + i] = val;
    else
      error("TMatrix::set_col overflow, i > num cols!\n");
  }
 
  void set_col (unsigned int i, vector<double> vec) {
 
    if(vec.size() < _rows)
      fatal("TMatrix::set_col copy from vector too small to fit\n");
 
    if( i < _cols)
      for(unsigned int j = 0; j < _rows; ++j)
        _val[j*_cols + i] = vec[j];
    else
      error("TMatrix::set_col overflow, i > num cols!\n");
  }
 
  void assign (double val)
  {
    for(unsigned int i = 0; i < _length; ++i) _val[i] = val;
  }
 
  void reset (unsigned int rows, unsigned int cols) {
    _length = rows * cols;
    if(_length == 0) error("TMatrix::attempt to reset a matrix with size = 0!\n");
    if(_val != NULL) delete [] _val;
    _val = new double [_length];
    memset(_val, 0, _length*sizeof(double));
    _rows = rows; _cols = cols;
  }
 
  void reset (unsigned int rows, unsigned int cols, double value) {
    reset(rows, cols);
    assign(value);
  }
 
  void reset (unsigned int rows, unsigned int cols, const double* array) {
    reset(rows, cols);
    memcpy(_val, array, _length * sizeof(double));
  }  
 
  void reset ( ) 
  {
    _rows = 0; 
    _cols = 0; 
    _length = 0;
    if(_val != NULL) delete [] _val; 
    _val = NULL;
  }
  
  double get (unsigned int i, unsigned int j) const {
    if( !((i+1)*(j+1) > _length) )
      return _val[i*_cols + j];
    else
      fatal("TMatrix::get overflow!\n"); 
    return 0;
  }
  double* get () const {return _val;}
  double* getValArray () const {return _val;}
  unsigned int get_dims (unsigned int* dims) const {
    if(dims != NULL) { dims[0] = _rows; dims[1] = _cols; }
    return _length;
  }
  unsigned int getNbRows ( ) const {return _rows;}
  unsigned int nrows () const {return _rows;}
  unsigned int getNbCols ( ) const {return _cols;}
  unsigned int ncols () const {return _cols;}
  unsigned int length    ( ) const {return _length;}
  void getColumnView (unsigned int col, unsigned int n, double* array)
  {
    if(col > _cols-1) {
      error("TMatrix::getColumnView: not that many columns in matrix\n");
      return;
    }
    if(n != _rows) {
      error("TMatrix::getColumnView: array size not equal to number of rows in matrix\n");
      return;
    }
    for(unsigned int i = 0; i < _rows; ++i)
      array[i] = _val[i*_cols + col];
  }
  void getRowView (unsigned int row, unsigned int n, double* array)
  {
    if(row > _rows-1) {
      error("TMatrix::getRowView: not that many rows in matrix\n");
      return;
    }
    if(n != _cols) {
      error("TMatrix::getRowView: array size not equal to number of columns in matrix\n");
      return;
    }
    for(unsigned int i = 0, stride = row*_cols; i < _cols; ++i)
      array[i] = _val[stride + i];
  }
  void plus (unsigned int i, unsigned int j, double value)
  {
    if( i*j < _length)
      _val[i*_cols + j] += value;
    else
      error("TMatrix::plus overflow!\n");
  }
  void matrix_increment (double value)
  {
     for(unsigned int i = 0; i < _rows; ++i){
       for(unsigned int j = 0; j < _cols; ++j){
         plus(i,j,value);
       }
     }
  }
  void sweep_plus (double value)
  {
     for(unsigned int i = 0; i < _rows; ++i){
       for(unsigned int j = 0; j < _cols; ++j){
       plus(i,j,value);
       }
     }
  }
  void minus (unsigned int i, unsigned int j, double value)
  {
    if( i*j < _length)
      _val[i*_cols + j] -= value;
    else
      error("TMatrix::minus overflow!\n");
  }
  void sweep_minus (double value)
  {
    for(unsigned int i = 0; i < _rows; ++i){
      for(unsigned int j = 0; j < _cols; ++j){
        minus(i,j,value);
      }
    }
  }
  void multi (unsigned int i, unsigned int j, double value)
  {
    if( i*j < _length)
      _val[i*_cols + j] *= value;
    else
      error("TMatrix::multi overflow!\n");
  }
  void sweep_multiply (double value)
  {
    for(unsigned int i = 0; i < _rows; ++i){
      for(unsigned int j = 0; j < _cols; ++j){
        multi(i,j,value);
      }
    }
  }
  void divide (unsigned int i, unsigned int j, double value)
  {
    if( i*j < _length)
      _val[i*_cols + j] /= value;
    else
      error("TMatrix::divide overflow!\n");
  }
  void sweep_divide (double value)
  {
    for(unsigned int i = 0; i < _rows; ++i){
      for(unsigned int j = 0; j < _cols; ++j){
        divide(i,j,value);
      }
    }
  }
  void transpose ()
  {
    TMatrix tmp(_cols, _rows);
    
    for(unsigned int i = 0; i < _rows; i++)
      for(unsigned int j = 0; j < _cols; j++)
        tmp.set(j, i, get(i, j));
    
    reset(_cols, _rows, tmp.get());
  }
  double colSum (unsigned int col)
  {
    assert(col < _cols); //has to be [0, _cols-1]
    double sum = 0;
    for (unsigned int i = 0; i < _rows; ++i) {
      sum += _val[i*_cols + col];
    }
    return sum;
  }
  double rowSum (unsigned int row)
  {
    assert(row < _rows); // has to be [0, _rows-1]
    double sum = 0;
    for (unsigned int i = 0; i < _cols; ++i) {
      sum += _val[row*_cols + i];
    }
    return sum;
  }
  
  void show_up()
  {
    message("TMatrix dimensions: \nrows = %i, columns = %i, length = %i\n",_rows,_cols, _length);
    for(unsigned int i = 0; i < _rows; i++) {
      for(unsigned int j = 0; j < _cols; j++)
        message("%.3f ",_val[i*_cols + j]);
      message("\n");
    }
  }
  string to_string ()
  {
    ostringstream OUT;
    
    OUT << "{";
    for(unsigned int i = 0; i < _rows; i++) {
      OUT<<"{";
      for(unsigned int j = 0; j < _cols-1; j++) {
        OUT<<_val[i*_cols + j]<<",";
      }
      OUT<<_val[i*_cols + _cols-1]<<"}";
    }
    OUT << "}";
    
    return OUT.str();
  }
  
#ifdef HAS_GSL
  void get_gsl_matrix(gsl_matrix* M)
  {
    if(M->size1 != _rows)
      fatal("TMatrix::get_gsl_matrix row size of input matrix doesn't match!\n");
    if(M->size2 != _cols)
      fatal("TMatrix::get_gsl_matrix col size of input matrix doesn't match!\n");
    for(unsigned int i = 0; i < _rows; ++i)
      for(unsigned int j = 0; j < _cols; ++j)
        gsl_matrix_set(M,i,j,_val[i*_cols + j]);
  }
  void set_from_gsl_matrix (gsl_matrix* M)
  {
    if(_val != 0) delete [] _val;
    _rows = M->size1;
    _cols = M->size2;
    _length = _rows*_cols;
    _val = new double [_length];
    memcpy(_val, M->data, _length*sizeof(double));
  }
  void inverse ()
  {
    if(_rows != _cols) {
      error("TMatrix::inverse matrix must be square to invert!\n");
      return;
    }
    
    gsl_matrix *S = gsl_matrix_alloc(_rows,_cols);
    get_gsl_matrix(S);
    gsl_matrix *LU = gsl_matrix_alloc(_rows,_cols);
    gsl_permutation *P = gsl_permutation_alloc(_rows);
    int snum;
    
    gsl_matrix_memcpy(LU,S);
    gsl_linalg_LU_decomp(LU,P,&snum);
    gsl_linalg_LU_invert(LU,P,S);
    
    set_from_gsl_matrix(S);
    
    gsl_matrix_free(S);
    gsl_matrix_free(LU);
    gsl_permutation_free(P);
  }
#endif
};
 
#endif