atio.hpp

Go to the documentation of this file.

#ifndef ATIO_HPP
#define ATIO_HPP

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

// Typedef declaration
typedef Eigen::SparseMatrix<double, Eigen::RowMajor> SpMatCSR;
typedef Eigen::SparseMatrix<double, Eigen::ColMajor> SpMatCSC;
typedef Eigen::Triplet<double> Tri;


// Declarations
void Eigen2Pajek(FILE *, const SpMatCSR *, const char *);
void Eigen2Compressed(FILE *, const SpMatCSC *, const char *);
void Eigen2Compressed(FILE *, const SpMatCSR *, const char *);
SpMatCSC * pajek2Eigen(FILE *, const char *);
SpMatCSR * Compressed2Eigen(FILE *, const char *);
gsl_matrix * Compressed2EdgeList(FILE *, const char *);
void Compressed2EdgeList(FILE *, FILE *, const char *);
SpMatCSR * CSC2CSR(const SpMatCSC *T);
SpMatCSC * CSR2CSC(const SpMatCSR *T);
std::vector<Tri> Eigen2Triplet(const SpMatCSC *);
std::vector<Tri> Eigen2Triplet(const SpMatCSR *);
void fprintfEigen(FILE *, const SpMatCSR *, const char *);
size_t lineCount(FILE *);

// Definitions
void
Eigen2Pajek(FILE *fp, const SpMatCSR *P, const char *dataFormat="%lf"){
  int N = P->rows();
  int E = P->nonZeros();

  // Write vertices
  fprintf(fp, "*Vertices %d\n", N);
  for (int k = 0; k < N; k++)
    fprintf(fp, "%d \"%d\"\n", k, k);

  // Write Edges
  fprintf(fp, "Edges %d\n", E);
  for (int i = 0; i < P->rows(); i++){
    for (SpMatCSR::InnerIterator it(*P, i); it; ++it){
      fprintf(fp, "%d %d ", i, it.col());
      fprintf(fp, dataFormat, it.value());
      fprintf(fp, "\n");
    }
  }

  return;
}

SpMatCSC *
pajek2Eigen(FILE *fp, const char *dataFormat="%lf"){
  char label[20];
  int N, E;
  int i, j, i0;
  double x;

  std::vector<Tri> tripletList;

  // Read vertices
  fscanf(fp, "%s %d", label, &N);

  // Define sparse matrix
  SpMatCSC *T = new SpMatCSC(N, N);

  // Read first (assume monotonous)
  fscanf(fp, "%d %s", &i0, label);
  for (int k = 1; k < N; k++){
    fscanf(fp, "%d %s", &i, label);
  }

  // Read Edges
  fscanf(fp, "%s %d", label, &E);

  // Reserve triplet capacity
  tripletList.reserve(E);

  for (int k = 0; k < E; k++){
    fscanf(fp, "%d %d", &i, &j);
    fscanf(fp, dataFormat, &x);
    tripletList.push_back(Tri(i - i0, j - i0, x));
  }
  
  // Assign matrix elements
  T->setFromTriplets(tripletList.begin(), tripletList.end());

  return T;
}

void
Eigen2Compressed(FILE *fp, const SpMatCSC *P, const char *dataFormat="%lf"){
  char sparseType[] = "CSC";

  // Write type, inner size, outer size and number of non-zeros
  fprintf(fp, "%s %d %d %d\n", sparseType, P->innerSize(), P->outerSize(), P->nonZeros());

  // Write values
  for (int nz = 0; nz < P->nonZeros(); nz++){
    fprintf(fp, dataFormat, (P->valuePtr())[nz]);
    fprintf(fp, " ");
  }
  fprintf(fp, "\n");

  // Write row indices
  for (int nz = 0; nz < P->nonZeros(); nz++)
    fprintf(fp, "%d ", (P->innerIndexPtr())[nz]);
  fprintf(fp, "\n");

  // Write first element of column pointer
  for (int outer = 0; outer < P->outerSize()+1; outer++)
    fprintf(fp, "%d ", (P->outerIndexPtr())[outer]);
  fprintf(fp, "\n");
  
  return;
}

void
Eigen2Compressed(FILE *fp, const SpMatCSR *P, const char *dataFormat="%lf"){
  char sparseType[] = "CSR";

  // Write type, inner size, outer size and number of non-zeros
  fprintf(fp, "%s %d %d %d\n", sparseType, P->innerSize(), P->outerSize(), P->nonZeros());

  // Write values
  for (int nz = 0; nz < P->nonZeros(); nz++){
    fprintf(fp, dataFormat, (P->valuePtr())[nz]);
    fprintf(fp, " ");
  }
  fprintf(fp, "\n");

  // Write column indices
  for (int nz = 0; nz < P->nonZeros(); nz++)
    fprintf(fp, "%d ", (P->innerIndexPtr())[nz]);
  fprintf(fp, "\n");

  // Write first element of row pointer
  for (int outer = 0; outer < P->outerSize()+1; outer++)
    fprintf(fp, "%d ", (P->outerIndexPtr())[outer]);
  fprintf(fp, "\n");
  
  return;
}

SpMatCSR *
Compressed2Eigen(FILE *fp, const char *dataFormat="%lf")
{
  int innerSize, outerSize, nnz;
  char sparseType[4];
  double *nzval;
  int *innerIndexPtr, *outerIndexPtr;
  SpMatCSR *T;
  std::vector<Tri> tripletList;

  // Read type, inner size, outer size and number of non-zeros and allocate
  fscanf(fp, "%s %d %d %d", sparseType, &outerSize, &innerSize, &nnz);
  nzval = new double [nnz];
  innerIndexPtr = new int [nnz];
  outerIndexPtr = new int[outerSize+1];
  T = new SpMatCSR(outerSize, innerSize);
  Eigen::VectorXf innerNNZ(outerSize);

  // Read values
  for (int nz = 0; nz < nnz; nz++)
    fscanf(fp, dataFormat, &nzval[nz]);

  // Read inner indices (column)
  for (int nz = 0; nz < nnz; nz++)
    fscanf(fp, "%d ", &innerIndexPtr[nz]);

  // Read first element of column pointer
  fscanf(fp, "%d ", &outerIndexPtr[0]);
  for (int outer = 1; outer < outerSize+1; outer++){
    fscanf(fp, "%d ", &outerIndexPtr[outer]);
    innerNNZ(outer-1) = outerIndexPtr[outer] - outerIndexPtr[outer-1];
  }
  T->reserve(innerNNZ);

  // Insert elements
  for (int outer = 0; outer < outerSize; outer++)
     for (int nzInner = outerIndexPtr[outer]; nzInner < outerIndexPtr[outer+1]; nzInner++)
       T->insertBackUncompressed(outer, innerIndexPtr[nzInner]) =  nzval[nzInner];
  //T.sumupDuplicates();

  delete nzval;
  delete innerIndexPtr;
  delete outerIndexPtr;
  
  return T;
}

gsl_matrix *
Compressed2EdgeList(FILE *fp, const char *dataFormat="%lf")
{
  int innerSize, outerSize, nnz;
  char sparseType[4];
  double *nzval;
  int *idx, *ptr;
  gsl_matrix *EdgeList;
  int iOuter, iInner;
  
  // Read type, inner size, outer size and number of non-zeros and allocate
  fscanf(fp, "%s %d %d %d", sparseType, &innerSize, &outerSize, &nnz);
  nzval = new double [nnz];
  idx = new int [nnz];
  ptr = new int[outerSize+1];
  EdgeList = gsl_matrix_alloc(nnz, 3);
  if (strcmp(sparseType, "CSR") == 0)
    iOuter = 0;
  else if (strcmp(sparseType, "CSR") == 0)
    iOuter = 1;
  else{
    std::cerr << "Invalid sparse matrix type." << std::endl;
    exit(EXIT_FAILURE);
  }
  iInner = (iOuter + 1)%2;

  // Read values
  for (int nz = 0; nz < nnz; nz++)
    fscanf(fp, dataFormat, &nzval[nz]);

  // Read row indices
  for (int nz = 0; nz < nnz; nz++)
    fscanf(fp, "%d ", &idx[nz]);

  // Read first element of column pointer
  for (int outer = 0; outer < outerSize+1; outer++)
    fscanf(fp, "%d ", &ptr[outer]);

  int nz = 0;
  for (int outer = 0; outer < outerSize; outer++){
    for (int inner = ptr[outer]; inner < ptr[outer+1]; inner++){
    gsl_matrix_set(EdgeList, nz, iOuter, outer);
    gsl_matrix_set(EdgeList, nz, iInner, idx[inner]);
    gsl_matrix_set(EdgeList, nz, 2, nzval[inner]);
    nz++;
    }
  }

  delete(nzval);
  delete(idx);
  delete(ptr);
  return EdgeList;
}

void
Compressed2EdgeList(FILE *src, FILE *dst, const char *dataFormat="%lf")
{
  gsl_matrix *EdgeList = Compressed2EdgeList(src);
  size_t nnz = EdgeList->size1;

  for (size_t nz = 0; nz < nnz; nz++){
    fprintf(dst, "%d\t%d\t", (int) gsl_matrix_get(EdgeList, nz, 0),
        (int) gsl_matrix_get(EdgeList, nz, 1));
    fprintf(dst, dataFormat, gsl_matrix_get(EdgeList, nz, 2));
    fprintf(dst, "\n");
  }

  gsl_matrix_free(EdgeList);
  return;
}

SpMatCSR *
CSC2CSR(const SpMatCSC *T){
  int N = T->rows();

  // Define sparse matrix
  SpMatCSR *TCSR = new SpMatCSR(N, N);

  // Get triplet list
  std::vector<Tri> tripletList = Eigen2Triplet(T);

  // Assign matrix elements
  TCSR->setFromTriplets(tripletList.begin(), tripletList.end());

  return TCSR;
}

SpMatCSC *
CSR2CSC(const SpMatCSR *T){
  int N = T->rows();

  // Define sparse matrix
  SpMatCSC *TCSC = new SpMatCSC(N, N);

  // Get triplet list
  std::vector<Tri> tripletList = Eigen2Triplet(T);

  // Assign matrix elements
  TCSC->setFromTriplets(tripletList.begin(), tripletList.end());

  return TCSC;
}

std::vector<Tri>
Eigen2Triplet(const SpMatCSC *T)
{
  // Works for column major only
  int nnz = T->nonZeros();
  std::vector<Tri> tripletList;
  tripletList.reserve(nnz);

  for (int beta = 0; beta < T->outerSize(); ++beta)
    for (SpMatCSC::InnerIterator it(*T, beta); it; ++it)
      tripletList.push_back(Tri(it.row(), it.col(), it.value()));

  return tripletList;
}

std::vector<Tri>
Eigen2Triplet(const SpMatCSR *T)
{
  // Works for column major only
  int nnz = T->nonZeros();
  std::vector<Tri> tripletList;
  tripletList.reserve(nnz);

  for (int beta = 0; beta < T->outerSize(); ++beta)
    for (SpMatCSR::InnerIterator it(*T, beta); it; ++it)
      tripletList.push_back(Tri(it.row(), it.col(), it.value()));

  return tripletList;
}

void
fprintfEigen(FILE *fp, const SpMatCSR *T, const char *dataFormat)
{
  int count;
  for (int irow = 0; irow < T->outerSize(); ++irow){
    count = 0;
    for (SpMatCSR::InnerIterator it(*T, irow); it; ++it){
      while (count < it.col()){
    fprintf(fp, "---\t");
    count++;
      }
      fprintf(fp, dataFormat, it.value());
      fprintf(fp, "\t");
      count++;
    }
    while (count < T->innerSize()){
      fprintf(fp, "---\t");
      count++;
    }
    fprintf(fp, "\n");
  }
  return;
}

size_t
lineCount(FILE *fp)
{
  size_t count = 0;
  int ch;

  // Count lines
  do {
    ch = fgetc(fp);
    if (ch == '\n')
      count++;
  } while (ch != EOF);

  return count;
}

#endif