atmatrix.hpp

Go to the documentation of this file.

#ifndef ATMATRIX_HPP
#define ATMATRIX_HPP

#include <vector>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_int.h>
#include <gsl/gsl_matrix.h>
#include <Eigen/Dense>
#include <Eigen/Sparse>

// Typedef definitions
typedef Eigen::SparseMatrix<double, Eigen::RowMajor> SpMatCSR;
typedef Eigen::SparseMatrix<int, Eigen::RowMajor> SpMatIntCSR;
typedef Eigen::SparseMatrix<double, Eigen::ColMajor> SpMatCSC;
typedef Eigen::SparseMatrix<bool, Eigen::ColMajor> SpMatCSCBool;

// Declarations
gsl_vector* getRowSum(SpMatCSR *);
void getRowSum(SpMatCSR *, gsl_vector *);
gsl_vector_int* getRowSum(SpMatIntCSR *);
gsl_vector* getColSum(SpMatCSR *);
void getColSum(SpMatCSR *, gsl_vector *);
gsl_vector* getColSum(SpMatCSC *);
double getSum(SpMatCSR *);
double sumVectorElements(gsl_vector *);
void normalizeVector(gsl_vector *);
void normalizeRows(SpMatCSR *, gsl_vector *);
SpMatCSCBool * cwiseGT(SpMatCSC *, double);
SpMatCSCBool * cwiseLT(SpMatCSC *, double);
bool any(SpMatCSCBool *);
double max(SpMatCSC *);
double min(SpMatCSC *);
std::vector<int> * argmax(SpMatCSC *);
std::vector<int> * argmin(SpMatCSC *);


// Definitions
gsl_vector*
getRowSum(SpMatCSR *T)
{
  int N = T->rows();
  gsl_vector *rowSum = gsl_vector_calloc(N);

  // Calculate row sums
  for (int i = 0; i < N; i++)
    for (SpMatCSR::InnerIterator it(*T, i); it; ++it)
      gsl_vector_set(rowSum, i, gsl_vector_get(rowSum, i)+it.value()); 
  
  return rowSum;
}

void
getRowSum(SpMatCSR *T, gsl_vector *rowSum)
{
  int N = T->rows();

  // Calculate row sums
  gsl_vector_set_all(rowSum, 0.);
  for (int i = 0; i < N; i++)
    for (SpMatCSR::InnerIterator it(*T, i); it; ++it)
      gsl_vector_set(rowSum, i, gsl_vector_get(rowSum, i)+it.value()); 
  
  return;
}

gsl_vector_int *
getRowSum(SpMatIntCSR *T)
{
  int N = T->rows();
  gsl_vector_int *rowSum = gsl_vector_int_calloc(N);
  // Calculate row sums
  for (int i = 0; i < N; i++)
    for (SpMatIntCSR::InnerIterator it(*T, i); it; ++it)
      gsl_vector_int_set(rowSum, i, gsl_vector_int_get(rowSum, i)+it.value()); 
  
  return rowSum;
}

gsl_vector *
getColSum(SpMatCSR *T)
{
  int N = T->rows();
  gsl_vector *colSum = gsl_vector_calloc(N);

  // Calculate col sums
  for (int irow = 0; irow < N; irow++)
    for (SpMatCSR::InnerIterator it(*T, irow); it; ++it)
      gsl_vector_set(colSum, it.col(),
             gsl_vector_get(colSum, it.col()) + it.value()); 
  
  return colSum;
}

void
getColSum(SpMatCSR *T, gsl_vector *colSum)
{
  int N = T->rows();

  // Calculate col sums
  gsl_vector_set_all(colSum, 0.);
  for (int irow = 0; irow < N; irow++)
    for (SpMatCSR::InnerIterator it(*T, irow); it; ++it)
      gsl_vector_set(colSum, it.col(),
             gsl_vector_get(colSum, it.col()) + it.value()); 
  
  return;
}

gsl_vector *
getColSum(SpMatCSC *T)
{
  int N = T->rows();
  gsl_vector *colSum = gsl_vector_calloc(N);
  // Calculate col sums
  for (int icol = 0; icol < N; icol++)
    for (SpMatCSC::InnerIterator it(*T, icol); it; ++it)
      gsl_vector_set(colSum, icol, gsl_vector_get(colSum, icol) + it.value()); 
  
  return colSum;
}

double
getSum(SpMatCSR *T)
{
  int N = T->rows();
  double sum = 0.;
  // Calculate col sums
  for (int irow = 0; irow < N; irow++)
    for (SpMatCSR::InnerIterator it(*T, irow); it; ++it)
      sum += it.value();
  
  return sum;
}


double
sumVectorElements(gsl_vector *v)
{
  double sum = 0.;
  
  for (size_t i = 0; i < v->size; i++)
    sum += gsl_vector_get(v, i);
  
  return sum;
}

void
normalizeVector(gsl_vector *v)
{
  double sum;

  sum = sumVectorElements(v);
  gsl_vector_scale(v, 1. / sum);
  
  return;
}

void
normalizeRows(SpMatCSR *T, gsl_vector *rowSum)
{
  double rowSumi;
  for (int i = 0; i < T->rows(); i++){
    rowSumi = gsl_vector_get(rowSum, i);
    for (SpMatCSR::InnerIterator it(*T, i); it; ++it)
      if (rowSumi > 0.)
    it.valueRef() /= rowSumi;
  }
  return ;
}

SpMatCSCBool *
cwiseGT(SpMatCSC *T, double ref)
{
  int j;
  SpMatCSCBool *cmpT = new SpMatCSCBool(T->rows(), T->cols());
  for (j = 0; j < T->cols(); j++)
    for (SpMatCSC::InnerIterator it(*T, j); it; ++it)
      cmpT->insert(it.row(), j) = it.value() > ref;

  return cmpT;
}
  
SpMatCSCBool *
cwiseLT(SpMatCSC *T, double ref)
{
  int j;
  SpMatCSCBool *cmpT = new SpMatCSCBool(T->rows(), T->cols());
  for (j = 0; j < T->cols(); j++)
    for (SpMatCSC::InnerIterator it(*T, j); it; ++it)
      cmpT->insert(it.row(), j) = it.value() < ref;

  return cmpT;
}

bool
any(SpMatCSCBool *T)
{
  int j;
  for (j = 0; j < T->cols(); j++)
    for (SpMatCSCBool::InnerIterator it(*T, j); it; ++it)
      if (it.value())
    return true;

  return false;
}

double
max(SpMatCSC *T)
{
  int j;
  SpMatCSC::InnerIterator it(*T, 0);
  double maxValue = it.value();
  for (j = 0; j < T->cols(); j++)
    for (SpMatCSC::InnerIterator it(*T, j); it; ++it)
      if (it.value() > maxValue)
    maxValue = it.value();

  return maxValue;
}


double
min(SpMatCSC *T)
{
  int j;
  SpMatCSC::InnerIterator it(*T, 0);
  double minValue = it.value();
  for (j = 0; j < T->cols(); j++)
    for (SpMatCSC::InnerIterator it(*T, j); it; ++it)
      if (it.value() < minValue)
    minValue = it.value();

  return minValue;
}


std::vector<int> *
argmax(SpMatCSC *T)
{
  int j;
  std::vector<int> *argmax = new std::vector<int>(2);
  SpMatCSC::InnerIterator it(*T, 0);
  double maxValue = it.value();
  for (j = 0; j < T->cols(); j++){
    for (SpMatCSC::InnerIterator it(*T, j); it; ++it){
      if (it.value() > maxValue){
    argmax->at(0) = it.row();
    argmax->at(1) = it.col();
    maxValue = it.value();
      }
    }
  }

  return argmax;
}

std::vector<int> *
argmin(SpMatCSC *T)
{
  int j;
  std::vector<int> *argmin = new std::vector<int>(2);
  SpMatCSC::InnerIterator it(*T, 0);
  double minValue = it.value();
  for (j = 0; j < T->cols(); j++){
    for (SpMatCSC::InnerIterator it(*T, j); it; ++it){
      if (it.value() < minValue){
    argmin->at(0) = it.row();
    argmin->at(1) = it.col();
    minValue = it.value();
      }
    }
  }

  return argmin;
}


#endif