transferOperator.hpp

Go to the documentation of this file.

#ifndef TRANSFEROPERATOR_HPP
#define TRANSFEROPERATOR_HPP

#include <iostream>
#include <vector>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_uint.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_matrix_uint.h>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <ATSuite/atmatrix.hpp>
#include <ATSuite/atmarkov.hpp>
#include <ATSuite/atio.hpp>
#if defined (WITH_OMP) && WITH_OMP == 1
#include <omp.h>
#endif

/* 
 * Typedef declaration
 */

typedef Eigen::Triplet<double> triplet;
typedef std::vector<triplet> tripletVector;
typedef Eigen::Triplet<size_t> tripletUInt;
typedef std::vector<tripletUInt> tripletUIntVector;
typedef Eigen::SparseMatrix<double, Eigen::ColMajor> SpMatCSC;
typedef Eigen::SparseMatrix<double, Eigen::RowMajor> SpMatCSR;


/*
 * Class declarations
 */

class Grid {
  void allocate(gsl_vector_uint *);
  void getRectGrid(gsl_vector_uint *, const gsl_vector *, const gsl_vector *);

public:
  size_t dim;
  size_t N;
  gsl_vector_uint *nx;
  std::vector<gsl_vector *> *gridBounds;

  Grid(){}
  Grid(gsl_vector_uint *nx_){ allocate(nx_); }
  Grid(gsl_vector_uint *, const gsl_vector *, const gsl_vector *);
  Grid(size_t, size_t, double, double);
  ~Grid();
  
  int printGrid(const char *, const char *, bool);
};

class transferOperator {
  void allocate(size_t);
  void buildFromMembership(const gsl_matrix_uint *);
    
public:
  size_t N;              
  SpMatCSR *P;           
  SpMatCSR *Q;           
  gsl_vector *initDist;  
  gsl_vector *finalDist; 

  transferOperator(){}
  transferOperator(size_t gridSize){ allocate(gridSize); }
  transferOperator(const gsl_matrix_uint *, size_t);
  transferOperator(const gsl_vector_uint *, size_t); 
  transferOperator(const gsl_matrix *, const gsl_matrix *, const Grid *);
  transferOperator(const gsl_matrix *, const Grid *, size_t);
  ~transferOperator();

  // Output methods
  int printForwardTransition(const char *, const char *);
  int printBackwardTransition(const char *, const char *);
  int printInitDist(const char *, const char *);
  int printFinalDist(const char *, const char *);
};


/*
 *  Functions declarations
 */

gsl_matrix_uint *getGridMemMatrix(const gsl_matrix *, const gsl_matrix *, const Grid *);
gsl_matrix_uint *getGridMemMatrix(const gsl_matrix *, const Grid *, const size_t);
gsl_vector_uint *getGridMemVector(const gsl_matrix *, const Grid *);
gsl_matrix_uint *memVector2memMatrix(const gsl_vector_uint *, size_t);
gsl_matrix_uint * memVectorList2memMatrix(const std::vector<gsl_vector_uint *> *, size_t);
int getBoxMembership(const gsl_vector *, const Grid *);
tripletUIntVector *getTransitionCountTriplet(const gsl_matrix_uint *, size_t);


/*
 * Constructors and destructors definitions
 */

transferOperator::transferOperator(const gsl_matrix_uint *gridMem, size_t gridSize)
{
  // Allocate
  allocate(gridSize);

  // Get transition matrices and distributions from grid membership matrix
  buildFromMembership(gridMem);

  return;
}

transferOperator::transferOperator(const gsl_matrix *initStates,
                   const gsl_matrix *finalStates,
                   const Grid *grid)
{
  gsl_matrix_uint *gridMem;

  // Allocate
  N = grid->N;
  allocate(N);

  // Get grid membership matrix
  gridMem = getGridMemMatrix(initStates, finalStates, grid);

  // Get transition matrices and distributions from grid membership matrix
  buildFromMembership(gridMem);

  // Free
  gsl_matrix_uint_free(gridMem);
  
  return;
}

transferOperator::transferOperator(const gsl_matrix *states, const Grid *grid,
                   const size_t tauStep)
{
  gsl_matrix_uint *gridMem;

  // Allocate
  N = grid->N;
  allocate(N);

  // Get grid membership matrix from a single long trajectory
  gridMem = getGridMemMatrix(states, grid, tauStep);

  // Get transition matrices and distributions from grid membership matrix
  buildFromMembership(gridMem);

  // Free
  gsl_matrix_uint_free(gridMem);
  
  return;
}

transferOperator::~transferOperator()
{
  delete P;
  delete Q;
  gsl_vector_free(initDist);
  gsl_vector_free(finalDist);
}

Grid::Grid(gsl_vector_uint *nx_, const gsl_vector *xmin, const gsl_vector *xmax)
{
  // Allocate and build uniform rectangular grid
  getRectGrid(nx_, xmin, xmax);
}

Grid::Grid(size_t dim_, size_t inx, double dxmin, double dxmax)
{
  // Convert to uniform vectors to call getRectGrid.
  gsl_vector_uint *nx_ = gsl_vector_uint_alloc(dim_);
  gsl_vector *xmin_ = gsl_vector_alloc(dim_);
  gsl_vector *xmax_ = gsl_vector_alloc(dim_);
  gsl_vector_uint_set_all(nx_, inx);
  gsl_vector_set_all(xmin_, dxmin);
  gsl_vector_set_all(xmax_, dxmax);

  // Allocate and build uniform rectangular grid
  getRectGrid(nx_, xmin_, xmax_);

  // Free
  gsl_vector_uint_free(nx_);
  gsl_vector_free(xmin_);
  gsl_vector_free(xmax_);
}

Grid::~Grid()
{
  gsl_vector_uint_free(nx);
  for (size_t d = 0; d < dim; d++)
    gsl_vector_free((*gridBounds)[d]);
  delete gridBounds;
}


/*
 * Methods definitions
 */

void
transferOperator::allocate(size_t gridSize)
{
  N = gridSize;
  P = new SpMatCSR(N, N);
  Q = new SpMatCSR(N, N);
  initDist = gsl_vector_alloc(N);
  finalDist = gsl_vector_alloc(N);
  
  return;
}

void
transferOperator::buildFromMembership(const gsl_matrix_uint *gridMem)
{
  // Get transition count triplets
  tripletUIntVector *T = getTransitionCountTriplet(gridMem, N);
  
  // Convert to CSR matrix
  P->setFromTriplets(T->begin(), T->end());
  
  // Get initial and final distribution
  getRowSum(P, initDist);
  getColSum(P, finalDist);
  
  // Get forward and backward transition matrices
  *Q = SpMatCSR(P->transpose());
  toLeftStochastic(P);
  toLeftStochastic(Q);
  normalizeVector(initDist);
  normalizeVector(finalDist);

  // Free
  delete T;

  return;
}

int
transferOperator::printForwardTransition(const char *path, const char *dataFormat="%lf")
{
  FILE *fp;

  // Open file
  if ((fp = fopen(path, "w")) == NULL){
    fprintf(stderr, "Can't open %s for printing the forward transition matrix!\n", path);
    return(EXIT_FAILURE);
  }

  // Print
  Eigen2Compressed(fp, P, dataFormat);

  // Close
  fclose(fp);

  return 0;
}

int
transferOperator::printBackwardTransition(const char *path, const char *dataFormat="%lf")
{
  FILE *fp;

  // Open file
  if ((fp = fopen(path, "w")) == NULL){
    fprintf(stderr, "Can't open %s for printing the backward transition matrix!\n", path);
    return(EXIT_FAILURE);
  }

  // Print
  Eigen2Compressed(fp, Q, dataFormat);

  // Close
  fclose(fp);

  return 0;
}

int
transferOperator::printInitDist(const char *path, const char *dataFormat="%lf")
{
  FILE *fp;

  // Open file
  if ((fp = fopen(path, "w")) == NULL){
    fprintf(stderr, "Can't open %s for printing the initial distribution!\n", path);
    return(EXIT_FAILURE);
  }

  // Print
  gsl_vector_fprintf(fp, initDist, dataFormat);

  // Close
  fclose(fp);

  return 0;
}

int
transferOperator::printFinalDist(const char *path, const char *dataFormat="%lf")
{
  FILE *fp;

  // Open file
  if ((fp = fopen(path, "w")) == NULL){
    fprintf(stderr, "Can't open %s for printing the final distribution!\n", path);
    return(EXIT_FAILURE);
  }

  // Print
  gsl_vector_fprintf(fp, finalDist, dataFormat);

  // Close
  fclose(fp);

  return 0;
}

void
Grid::allocate(gsl_vector_uint *nx_)
{
  dim = nx_->size;
  
  nx = gsl_vector_uint_alloc(dim);
  gsl_vector_uint_memcpy(nx, nx_);
  
  N = 1;
  gridBounds = new std::vector<gsl_vector *>(dim);
  for (size_t d = 0; d < dim; d++){
    N *= gsl_vector_uint_get(nx, d);
    (*gridBounds)[d] = gsl_vector_alloc(gsl_vector_uint_get(nx, d) + 1);
  }

  return;
}

void
Grid::getRectGrid(gsl_vector_uint *nx_, const gsl_vector *xmin, const gsl_vector *xmax)
{
  double delta;
  
  // Allocate
  allocate(nx_);

  // Build uniform grid bounds
  for (size_t d = 0; d < dim; d++) {
    // Get spatial step
    delta = (gsl_vector_get(xmax, d) - gsl_vector_get(xmin, d))
      / gsl_vector_uint_get(nx, d);
    // Set grid bounds
    gsl_vector_set((*gridBounds)[d], 0, gsl_vector_get(xmin, d));
    for (size_t i = 1; i < gsl_vector_uint_get(nx, d) + 1; i++)
      gsl_vector_set((*gridBounds)[d], i,
             gsl_vector_get((*gridBounds)[d], i-1) + delta);
  }

  return;
}

int
Grid::printGrid(const char *path, const char *dataFormat="%lf", bool verbose=false)
{
  gsl_vector *bounds;
  FILE *fp;

  // Open file
  if ((fp = fopen(path, "w")) == NULL){
    fprintf(stderr, "Can't open %s for printing the grid!\n", path);
    return(EXIT_FAILURE);
  }

  if (verbose)
    std::cout << "Domain grid (min, max, n):" << std::endl;

  // Print grid
  for (size_t d = 0; d < dim; d++) {
    bounds = (*gridBounds)[d];
    if (verbose) {
      std::cout << "dim " << d+1 << ": ("
        << gsl_vector_get(bounds, 0) << ", "
        << gsl_vector_get(bounds, bounds->size - 1) << ", "
        << (bounds->size - 1) << ")" << std::endl;
    }
    
    for (size_t i = 0; i < bounds->size; i++){
      fprintf(fp, dataFormat, gsl_vector_get((*gridBounds)[d], i));
      fprintf(fp, " ");
    }
    fprintf(fp, "\n");
  }

  // Close
  fclose(fp);

  return 0;
}


/*
 * Function definitions
 */

gsl_matrix_uint *
getGridMemMatrix(const gsl_matrix *initStates, const gsl_matrix *finalStates,
         const Grid *grid)
{
  const size_t nTraj = initStates->size1;
  const size_t dim = initStates->size2;
  gsl_matrix_uint *gridMem;

  // Allocate
  gridMem = gsl_matrix_uint_alloc(grid->N, 2);
  
  // Assign a pair of source and destination boxes to each trajectory
#pragma omp parallel
  {
    gsl_vector *X = gsl_vector_alloc(dim);
    
#pragma omp for
    for (size_t traj = 0; traj < nTraj; traj++) {
      // Find initial box
      gsl_matrix_get_row(X, initStates, traj);
      gsl_matrix_uint_set(gridMem, traj, 0, getBoxMembership(X, grid));
      
      // Find final box
      gsl_matrix_get_row(X, finalStates, traj);
      gsl_matrix_uint_set(gridMem, traj, 1, getBoxMembership(X, grid));
    }
    gsl_vector_free(X);
  }
  
  return gridMem;
}

gsl_vector_uint *
getGridMemVector(const gsl_matrix *states, const Grid *grid)
{
  const size_t nStates = states->size1;
  const size_t dim = states->size2;
  gsl_vector_uint *gridMem = gsl_vector_uint_alloc(nStates);

  // Assign a pair of source and destination boxes to each trajectory
#pragma omp parallel
  {
    gsl_vector *X = gsl_vector_alloc(dim);
    
#pragma omp for
    for (size_t traj = 0; traj < nStates; traj++) {
      // Find initial box
      gsl_matrix_get_row(X, states, traj);
#pragma omp critical
      {
    gsl_vector_uint_set(gridMem, traj, getBoxMembership(X, grid));
      }
    }
    gsl_vector_free(X);
  }
  
  return gridMem;
}

gsl_matrix_uint *
memVector2memMatrix(gsl_vector_uint *gridMemVect, size_t tauStep)
{
  const size_t nStates = gridMemVect->size;
  gsl_matrix_uint *gridMem = gsl_matrix_uint_alloc(nStates - tauStep, 2);

  // Get membership matrix from vector
  for (size_t traj = 0; traj < (nStates - tauStep); traj++) {
    gsl_matrix_uint_set(gridMem, traj, 0,
               gsl_vector_uint_get(gridMemVect, traj));
    gsl_matrix_uint_set(gridMem, traj, 1,
               gsl_vector_uint_get(gridMemVect, traj + tauStep));
  }

  return gridMem;
}

gsl_matrix_uint *
getGridMemMatrix(const gsl_matrix *states, const Grid *grid, const size_t tauStep)
{
  // Get membership vector
  gsl_vector_uint *gridMemVect = getGridMemVector(states, grid);

  // Get membership matrix from vector
  gsl_matrix_uint *gridMem = memVector2memMatrix(gridMemVect, tauStep);

  // Free
  gsl_vector_uint_free(gridMemVect);
  
  return gridMem;
}

gsl_matrix_uint *
memVectorList2memMatrix(const std::vector<gsl_vector_uint *> *memList, size_t tauStep)
{
  size_t nStatesTot = 0;
  size_t count;
  const size_t listSize = memList->size();
  gsl_matrix_uint *gridMem, *gridMemMatrixL;

  // Get total number of states and allocate grid membership matrix
  for (size_t l = 0; l < listSize; l++)
    nStatesTot += (memList->at(l))->size;
  gridMem = gsl_matrix_uint_alloc(nStatesTot - tauStep * listSize, 2);
  
  // Get membership matrix from list of membership vectors
  count = 0;
  for (size_t l = 0; l < listSize; l++) {
    gridMemMatrixL = memVector2memMatrix(memList->at(l), tauStep);
    for (size_t t = 0; t < gridMemMatrixL->size1; t++) {
      gsl_matrix_uint_set(gridMem, count, 0,
              gsl_matrix_uint_get(gridMemMatrixL, t, 0));
      gsl_matrix_uint_set(gridMem, count, 1,
              gsl_matrix_uint_get(gridMemMatrixL, t, 1));
      count++;
    }
    gsl_matrix_uint_free(gridMemMatrixL);
  }
  
  return gridMem;
}

int
getBoxMembership(const gsl_vector *state, const Grid *grid)
{
  const size_t dim = state->size;
  size_t inBox, nBoxDir;
  size_t foundBox;
  size_t subbp, subbn, ids;
  gsl_vector *bounds;

  // Get box
  foundBox = grid->N;
  for (size_t box = 0; box < grid->N; box++){
    inBox = 0;
    subbp = box;
    for (size_t d = 0; d < dim; d++){
      bounds = grid->gridBounds->at(d);
      nBoxDir = bounds->size - 1;
      subbn = (size_t) (subbp / nBoxDir);
      ids = subbp - subbn * nBoxDir;
      inBox += (size_t) ((gsl_vector_get(state, d)
              >= gsl_vector_get(bounds, ids))
             & (gsl_vector_get(state, d)
                < gsl_vector_get(bounds, ids+1)));
      subbp = subbn;
    }
    if (inBox == dim){
      foundBox = box;
      break;
    }
  }
  
  return foundBox;
}

tripletUIntVector *
getTransitionCountTriplet(const gsl_matrix_uint *gridMem, size_t N)
{
  const size_t nTraj = gridMem->size1;
  size_t box0, boxf;
  size_t nOut = 0;
  tripletUIntVector *T = new tripletUIntVector;
  T->reserve(nTraj);

  for (size_t traj = 0; traj < nTraj; traj++) {
    box0 = gsl_matrix_uint_get(gridMem, traj, 0);
    boxf = gsl_matrix_uint_get(gridMem, traj, 1);
    
    // Add transition triplet
    if ((box0 < N) && (boxf < N))
      T->push_back(tripletUInt(box0, boxf, 1));
    else
      nOut++;
  }
  std::cout <<  nOut * 100. / nTraj
        << "% of the trajectories ended up out of the domain." << std::endl;

  return T;
}

#endif