// file: $isip_ifc/class/algo/Rps/rps_06.cc
// version: $Id: rps_06.cc 10484 2006-03-14 23:50:33Z srinivas $
//

// isip include files
//
#include "Rps.h"

// method: computeTimeDelayEmbedding
//
// arguments:
//  MatrixFloat& rps: (output) covariance coefficients
//  const VectorFloat& x: (input) input data
//
// return: a bool8 value indicating status
//
// this method computes RPS matrix from a scalar time-series using
// time-delay embedding
//
bool8 Rps::computeTimeDelayEmbedding(MatrixFloat& rps_a,
				       const VectorFloat& x_a) {

  int32 L = x_a.length();
  int32 delay =  (int32)Integral::round(delay_d * sample_freq_d);
  int32 N =  L - ((int32)embed_dimension_d - 1) * delay;
  int32 i, j;

  if (embed_dimension_d < (int32)1) {
    return Error::handle(name(), L"computeTimeDelayEmbedding", ERR_EMBEDDIM,
			 __FILE__, __LINE__);
  }

  if (N < (int32)1) {
    return Error::handle(name(), L"computeTimeDelayEmbedding", ERR_INP,
			 __FILE__, __LINE__);
  }

  rps_a.setDimensions(N, embed_dimension_d + (int32)1, false, Integral::FULL);
  
  // time-delay embedding
  //
  for (i = 0; i < N; i++) {
    for (j = 0; j < embed_dimension_d; j++) {
      rps_a.setValue(i, j, x_a(i + j * delay));
    }
    
    // Append a last column of ones for flagging as "no discontinuity"
    //
    rps_a.setValue(i, embed_dimension_d, (float32)1.0);
  }
  // Flag last column as the end of this RPS matrix
  //
  rps_a.setValue(N - 1, embed_dimension_d, (float32)0.0);
  
  return true;
}

// method: computeTimeDelayEmbedding
//
// arguments:
//  MatrixFloat& rps: (output) covariance coefficients
//  const MatrixFloat& x: (input) input data
//
// return: a bool8 value indicating status
//
// this method computes RPS matrix from a vector time-series using
// time-delay embedding
//
bool8 Rps::computeTimeDelayEmbedding(MatrixFloat& rps_a,
				       const MatrixFloat& x_a) {
  
  int32 L = x_a.getNumRows();
  int32 M = x_a.getNumColumns();
  int32 delay = (int32)Integral::round(delay_d * sample_freq_d);
  int32 N = L - ((int32)embed_dimension_d - 1) * delay;
  int32 i, j, k;
  int32 num_rows_concat = (int32)(embed_dimension_d / M);
  int32 num_last_row_elements = (int32)(embed_dimension_d % M);
  
  if (embed_dimension_d < (int32)1) {
    return Error::handle(name(), L"computeTimeDelayEmbedding", ERR_EMBEDDIM,
			 __FILE__, __LINE__);
  }

  if (N < 1) {
    return Error::handle(name(), L"computeTimeDelayEmbedding", ERR_INP,
			 __FILE__, __LINE__);
  }

  if (embed_dimension_d < M ) {
    return Error::handle(name(), L"computeTimeDelayEmbedding", ERR_EMBEDDIM,
			 __FILE__, __LINE__);
  }    
  rps_a.setDimensions(N, embed_dimension_d, false, Integral::FULL);
  
  for (i = 0; i < N; i++) {
    for (j = 0; j < num_rows_concat; j++) {
      for (k = 0; k < M; k++) {
	rps_a.setValue(i, j * M + k, x_a.getValue(i + j * delay, k));
      }
    }
    for (j = 0; j < num_last_row_elements; j++) {
      rps_a.setValue(i, num_rows_concat * M + j, x_a.getValue(i + num_rows_concat * delay, j));
    }
  }
  
  rps_a.setDimensions(N, embed_dimension_d + (int32)1, true, Integral::FULL);
  
  // Append a last column of ones for flagging as "no discontinuity"
  //
  for (i = 0; i < N; i++) {
    rps_a.setValue(i, embed_dimension_d, (float32)1.0);
  }
  // Flag last column as the end of this RPS matrix
  //
  rps_a.setValue(N - 1, embed_dimension_d, (float32)0.0);
  
  return true;
}

// method: computeSVDEmbedding
//
// arguments:
//  MatrixFloat& rps: (output) covariance coefficients
//  const VectorFloat& x: (input) input data
//
// return: a bool8 value indicating status
//
// this method computes RPS matrix from a scalar time-series using
// SVD embedding
//
bool8 Rps::computeSVDEmbedding(MatrixFloat& rps_a,
				 const VectorFloat& x_a) {
  
  int32 L = x_a.length();
  
  // Setting the delay to the specified value, default = 1 sample
  //
  int32 delay = (int32)Integral::round(delay_d * sample_freq_d);
  int32 N = L - ((int32)svd_window_size_d - 1) * delay;
  int32 i, j;
  MatrixFloat u, v, w;
  MatrixFloat rps1;
  
  if (svd_window_size_d < (int32)1) {
    return Error::handle(name(), L"computeSVDEmbedding", ERR_SVDDIM,
			 __FILE__, __LINE__);
  }

  if (N < 1) {
    return Error::handle(name(), L"computeSVDEmbedding", ERR_INP,
			 __FILE__, __LINE__);
  }

  rps1.setDimensions(N, svd_window_size_d, false, Integral::FULL);
  
  // first-step is similar to time-delay embedding but
  // on a higher dimension, svd_window_size_d
  //
  for (i = 0; i < N; i++) {
    for (j = 0; j < svd_window_size_d; j++) {
      rps1.setValue(i, j, x_a(i + j * delay));
    }
  }
  
  // use SVD to get a matrix of dimension embed_dimension_d
  //
  rps1.decompositionSVD(u, w, v);

  // If embedding dimension is not specified estimate it by
  // comparing SVD threshold value with the singular values
  //
  if (embed_dimension_d < (int32)1) {
    int32 min_mn = ((svd_window_size_d < N)? (int32)svd_window_size_d:N);

    embed_dimension_d = 0;
    for (i = 0; i < min_mn; i++) {
      if (w.getValue(min_mn, min_mn) < svd_threshold_d) {
	i--;
	break;
      }
    }
    embed_dimension_d = i + 1;

    // Error if embedding dimension is estimated as zero
    //
    if (embed_dimension_d == (int32)0) {
      return Error::handle(name(), L"computeSVDEmbedding", ERR_EMBEDDIM,
			   __FILE__, __LINE__);
    }
  }

  v.setDimensions(svd_window_size_d, embed_dimension_d, true);
  rps_a.mult(rps1, v);
  
  rps_a.setDimensions(N, embed_dimension_d + (int32)1, true, Integral::FULL);
  
  // Append a last column of ones for flagging as "no discontinuity"
  //
  for (i = 0; i < N; i++) {
    rps_a.setValue(i, embed_dimension_d, (float32)1.0);
  }
  // Flag last column as the end of this RPS matrix
  //
  rps_a.setValue(N - 1, embed_dimension_d, (float32)0.0);
  
  return true;
}

// method: computeSVDEmbedding
//
// arguments:
//  MatrixFloat& rps: (output) covariance coefficients
//  const MatrixFloat& x: (input) input data
//
// return: a bool8 value indicating status
//
// this method computes RPS matrix from a vector time-series using
// SVD embedding
//
bool8 Rps::computeSVDEmbedding(MatrixFloat& rps_a,
				 const MatrixFloat& x_a) {
  
  int32 L = x_a.getNumRows();
  int32 M = x_a.getNumColumns();
  // Setting the delay to the specified value, default = 1 sample
  //
  int32 delay = (int32)Integral::round(delay_d * sample_freq_d);
  int32 N = L - ((int32)svd_window_size_d - 1) * delay;
  int32 i, j, k;
  int32 num_rows_concat = (int32)(svd_window_size_d / M);
  int32 num_last_row_elements = (int32)(svd_window_size_d % M);
  MatrixFloat u, v, w;
  MatrixFloat rps1;
  
  if (svd_window_size_d < (int32)1) {
    return Error::handle(name(), L"computeSVDEmbedding", ERR_SVDDIM,
			 __FILE__, __LINE__);
  }

  if (N < 1) {
    return Error::handle(name(), L"computeSVDEmbedding", ERR_INP,
			 __FILE__, __LINE__);
  }

  if (svd_window_size_d < M ) {
    return Error::handle(name(), L"computeSVDEmbedding", ERR_SVDDIM,
			 __FILE__, __LINE__);
  }    
  rps1.setDimensions(N, svd_window_size_d, false, Integral::FULL);
  
  // first-step is similar to time-delay embedding but
  // on a higher dimension, svd_window_size_d
  //
  for (i = 0; i < N; i++) {
    for (j = 0; j < num_rows_concat; j++) {
      for (k = 0; k < M; k++) {
	rps1.setValue(i, j * M + k, x_a.getValue(i + j * delay, k));
      }
    }
    for (j = 0; j < num_last_row_elements; j++) {
      rps1.setValue(i, num_rows_concat * M + j, x_a.getValue(i + num_rows_concat * delay, j));
    }
  }
  
  // use SVD to get a matrix of dimension embed_dimension_d
  //
  rps1.decompositionSVD(u, w, v);

  // If embedding dimension is not specified estimate it by
  // comparing SVD threshold value with the singular values
  //
  if (embed_dimension_d < (int32)1) {
    int32 min_mn = ((svd_window_size_d < N)? (int32)svd_window_size_d : N);

    embed_dimension_d = 0;
    for (i = 0; i < min_mn; i++) {
      if (w.getValue(min_mn, min_mn) < svd_threshold_d) {
	i--;
	break;
      }
    }
    embed_dimension_d = i + 1;

    // Error if embedding dimension is estimated as zero
    //
    if (embed_dimension_d == (int32)0) {
      return Error::handle(name(), L"computeSVDEmbedding", ERR_EMBEDDIM,
			   __FILE__, __LINE__);
    }
  }

  v.setDimensions(svd_window_size_d, embed_dimension_d, true);
  rps_a.mult(rps1, v);
  
  rps_a.setDimensions(N, embed_dimension_d + (int32)1, true, Integral::FULL);
  
  // Append a last column of ones for flagging as "no discontinuity"
  //
  for (i = 0; i < N; i++) {
    rps_a.setValue(i, embed_dimension_d, (float32)1.0);
  }
  // Flag last column as the end of this RPS matrix
  //
  rps_a.setValue(N - 1, embed_dimension_d, (float32)0.0);
  
  return true;
}