text:
/* Routines for manipulating evolutionary rate matrices.
 *
 * There is no specific object for this module.  Rate matrix
 * operations use square nxn ESL_DMATRIX data objects. (The rmx
 * module essentially subclasses the dmx module.)
 *
 * An instantaneous rate matrix is usually denoted by Q.  A
 * conditional probability matrix (for a specific t) is usually
 * denoted by P. An exchangeability matrix is denoted by E.
 * A stationary residue probability vector is denoted by pi.
 *
 * Two important relations among these:
 *
 *    Q in terms of E and pi:
 *      $Q_{ij} = E_{ij} \pi_j$ for $i \neq j$;
 *      $Q_{ii} = -\sum_{j \neq i} Q_{ij}$
 *
 *    P in terms of Q and t:
 *      $P = e^{tQ}$
 *
 * Contents:
 *   1. Setting standard rate matrix models.
 *   2. Debugging routines for validating or dumping rate matrices.
 *   3. Other routines in the exposed ratematrix API.
 *   4. Benchmark driver.
 *   5. Regression test driver.
 *   6. Unit tests.
 *   7. Test driver.
 *   8. Example.
 *   9. Copyright/license information.
 *
 *
 *
 * See also:
 *   paml   - i/o of rate matrices from/to data files in PAML format
 *
 */
#include "esl_config.h"

#include <math.h>

#include "easel.h"
#include "esl_composition.h"
#include "esl_dmatrix.h"
#include "esl_vectorops.h"
#include "esl_ratematrix.h"

/*****************************************************************
 * 1. Setting standard rate matrix models.
 *****************************************************************/

/* Function:  esl_rmx_SetWAG()
 * Incept:    SRE, Thu Mar  8 18:00:00 2007 [Janelia]
 *
 * Purpose:   Sets a $20 \times 20$ rate matrix <Q> to WAG parameters.
 *            The caller allocated <Q>.
 *
 *            If <pi> is non-<NULL>, it provides a vector of 20 amino
 *            acid stationary probabilities in Easel alphabetic order,
 *            A..Y, and the WAG stationary probabilities are set to
 *            these desired $\pi_i$. If <pi> is <NULL>, the default
 *            WAG stationary probabilities are used.
 *
 *            The WAG parameters are a maximum likelihood
 *            parameterization obtained by Whelan and Goldman
 *            \citep{WhelanGoldman01}.
 *
 * Note:      The data table was reformatted from wag.dat by the UTILITY1
 *            executable in the paml module. The wag.dat file was obtained from
 *            \url{http://www.ebi.ac.uk/goldman/WAG/wag.dat}. A copy
 *            is in formats/wag.dat.
 *
 * Args:      Q   - a 20x20 rate matrix to set, allocated by caller.
 *            pi  - desired stationary probabilities A..Y, or
 *                  NULL to use WAG defaults.
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslEINVAL> if <Q> isn't a 20x20 general matrix; and
 *            the state of <Q> is undefined.
 */
int
esl_rmx_SetWAG(ESL_DMATRIX *Q, double *pi)
{
  static double wagE[190] = {
    1.027040, 0.738998, 0.030295, 1.582850, 0.021352, 6.174160, 0.210494, 0.398020, 0.046730, 0.081134,
    1.416720, 0.306674, 0.865584, 0.567717, 0.049931, 0.316954, 0.248972, 0.930676, 0.570025, 0.679371,
    0.249410, 0.193335, 0.170135, 0.039437, 0.127395, 1.059470, 0.030450, 0.138190, 0.906265, 0.074034,
    0.479855, 2.584430, 0.088836, 0.373558, 0.890432, 0.323832, 0.397915, 0.384287, 0.084805, 0.154263,
    2.115170, 0.061304, 0.499462, 3.170970, 0.257555, 0.893496, 0.390482, 0.103754, 0.315124, 1.190630,
    0.174100, 0.404141, 4.257460, 0.934276, 4.854020, 0.509848, 0.265256, 5.429420, 0.947198, 0.096162,
    1.125560, 3.956290, 0.554236, 3.012010, 0.131528, 0.198221, 1.438550, 0.109404, 0.423984, 0.682355,
    0.161444, 0.243570, 0.696198, 0.099929, 0.556896, 0.415844, 0.171329, 0.195081, 0.908598, 0.098818,
    0.616783, 5.469470, 0.099921, 0.330052, 4.294110, 0.113917, 3.894900, 0.869489, 1.545260, 1.543640,
    0.933372, 0.551571, 0.528191, 0.147304, 0.439157, 0.102711, 0.584665, 2.137150, 0.186979, 5.351420,
    0.497671, 0.683162, 0.635346, 0.679489, 3.035500, 3.370790, 1.407660, 1.071760, 0.704939, 0.545931,
    1.341820, 0.740169, 0.319440, 0.967130, 0.344739, 0.493905, 3.974230, 1.613280, 1.028870, 1.224190,
    2.121110, 0.512984, 0.374866, 0.822765, 0.171903, 0.225833, 0.473307, 1.458160, 1.386980, 0.326622,
    1.516120, 2.030060, 0.795384, 0.857928, 0.554413, 4.378020, 2.006010, 1.002140, 0.152335, 0.588731,
    0.649892, 0.187247, 0.118358, 7.821300, 0.305434, 1.800340, 2.058450, 0.196246, 0.314887, 0.301281,
    0.251849, 0.232739, 1.388230, 0.113133, 0.717070, 0.129767, 0.156557, 1.529640, 0.336983, 0.262569,
    0.212483, 0.137505, 0.665309, 0.515706, 0.071917, 0.139405, 0.215737, 1.163920, 0.523742, 0.110864,
    0.365369, 0.240735, 0.543833, 0.325711, 0.196303, 6.454280, 0.103604, 3.873440, 0.420170, 0.133264,
    0.398618, 0.428437, 1.086000, 0.216046, 0.227710, 0.381533, 0.786993, 0.291148, 0.314730, 2.485390};
  static double wagpi[20];
  int i,j,z;

  if (Q->m != 20 || Q->n != 20 || Q->type != eslGENERAL)
    ESL_EXCEPTION(eslEINVAL, "Q must be a 20x20 general matrix");
  esl_composition_WAG(wagpi);

  /* 1. Transfer the wag E lower triagonal matrix directly into Q. */
  z = 0;
  for (i = 0; i < 20; i++)
    {
      Q->mx[i][i] = 0.; /* code below depends on this zero initialization */
      for (j = 0; j < i; j++) {
	Q->mx[i][j] = wagE[z++];
	Q->mx[j][i] = Q->mx[i][j];
      }
    }

  /* 2. Set offdiagonals Q_ij = E_ij * pi_j */
  for (i = 0; i < 20; i++)
    for (j = 0; j < 20; j++)
      if (pi != NULL) Q->mx[i][j] *= pi[j];
      else            Q->mx[i][j] *= wagpi[j];

  /* 3. Set diagonal Q_ii to -\sum_{i \neq j} Q_ij */
  for (i = 0; i < 20; i++)
    Q->mx[i][i] = -1. * esl_vec_DSum(Q->mx[i], 20);

  /* 4. Renormalize matrix to units of 1 substitution/site. */
  if (pi != NULL) esl_rmx_ScaleTo(Q, pi,    1.0);
  else            esl_rmx_ScaleTo(Q, wagpi, 1.0);

  return eslOK;
}

/* Function:  esl_rmx_SetJukesCantor()
 * Incept:    SRE, Thu Mar 15 13:04:56 2007 [Janelia]
 *
 * Purpose:   Sets a 4x4 rate matrix to a Jukes-Cantor model,
 *            scaled to units of 1t = 1.0 substitutions/site.
 *
 * Note:     eigenvalues of Q are 0, -4\alpha, -4\alpha, -4\alpha
 */
int
esl_rmx_SetJukesCantor(ESL_DMATRIX *Q)
{
  int    i,j;
  double pi[4] = { 0.25, 0.25, 0.25, 0.25 };

  if (Q->m != 4 || Q->n != 4 || Q->type != eslGENERAL)
    ESL_EXCEPTION(eslEINVAL, "Q must be a 4x4 general matrix");

  for (i = 0; i < 4; i++) {
    for (j = 0; j < 4; j++)
      {
	if (i != j) Q->mx[i][j] = 1.0;
	else        Q->mx[i][j] = 0.0;
      }
    Q->mx[i][i] =  -1. * esl_vec_DSum(Q->mx[i], 4);
  }
  esl_rmx_ScaleTo(Q, pi, 1.0);
  return eslOK;
}

/* Function:  esl_rmx_SetKimura()
 * Incept:    SRE, Thu Mar 15 13:08:08 2007 [Janelia]
 *
 * Purpose:   Sets a 4x4 rate matrix to a Kimura 2-parameter
 *            model, given transition and transversion
 *            relative rates <alpha> and <beta>, respectively,
 *            scaled to units of 1t = 1.0 substitutions/site.
 *
 * Note:     eigenvalues of Q are 0, -4\alpha, -2(\alpha+\beta), -2(\alpha+\beta)
 */
int
esl_rmx_SetKimura(ESL_DMATRIX *Q, double alpha, double beta)
{
  int i,j;
  double pi[4] = { 0.25, 0.25, 0.25, 0.25 };

  if (Q->m != 4 || Q->n != 4 || Q->type != eslGENERAL)
    ESL_EXCEPTION(eslEINVAL, "Q must be a 4x4 general matrix");

  for (i = 0; i < 4; i++) {
    for (j = 0; j < 4; j++)
      {
	if (i != j) Q->mx[i][j] = ((i+j)%2)? beta : alpha; /* even=0=transition;odd=1=transversion */
	else        Q->mx[i][j] = 0.0;
      }
    Q->mx[i][i] =  -1. * esl_vec_DSum(Q->mx[i], 4);
  }
  esl_rmx_ScaleTo(Q, pi, 1.0);
  return eslOK;
}

/* Function:  esl_rmx_SetF81()
 * Incept:    SRE, Thu Mar 15 13:33:30 2007 [Janelia]
 *
 * Purpose:   Sets a 4x4 rate matrix to the F81 model (aka
 *            equal-input model) given stationary base
 *            compositions <pi>,
 *            scaled to units of 1t = 1.0 substitutions/site.
 */
int
esl_rmx_SetF81(ESL_DMATRIX *Q, double *pi)
{
  int i,j;

  if (Q->m != 4 || Q->n != 4 || Q->type != eslGENERAL)
    ESL_EXCEPTION(eslEINVAL, "Q must be a 4x4 general matrix");

  for (i = 0; i < 4; i++) {
    for (j = 0; j < 4; j++)
      {
	if (i != j) Q->mx[i][j] = pi[j];
	else        Q->mx[i][j] = 0.0;
      }
    Q->mx[i][i] =  -1. * esl_vec_DSum(Q->mx[i], 4);
  }
  esl_rmx_ScaleTo(Q, pi, 1.0);
  return eslOK;
}

/* Function:  esl_rmx_SetHKY()
 * Incept:    SRE, Thu Aug 12 08:26:39 2004 [St. Louis]
 *
 * Purpose:   Given stationary base composition <pi> for ACGT, and
 *            transition and transversion relative rates <alpha> and
 *            <beta> respectively, sets the matrix <Q> to be the
 *            corresponding HKY (Hasegawa/Kishino/Yano) DNA rate
 *            matrix, scaled in units of 1t= 1.0 substitutions/site
 *            \citep{Hasegawa85}.
 *
 * Args:      pi     - stationary base composition A..T
 *            alpha  - relative transition rate
 *            beta   - relative transversion rate
 *
 *
 * Returns:   <eslOK>
 *
 * Xref:
 */
int
esl_rmx_SetHKY( ESL_DMATRIX *Q, double *pi, double alpha, double beta)
{
  int i,j;

  if (Q->m != 4 || Q->n != 4 || Q->type != eslGENERAL)
    ESL_EXCEPTION(eslEINVAL, "Q must be a 4x4 general matrix");

  for (i = 0; i < 4; i++) {
    for (j = 0; j < 4; j++)
      {
	if (i != j)  Q->mx[i][j] = ((i+j)%2)? pi[j]*beta : pi[j]*alpha; /* even=0=transition;odd=1=transversion */
	else         Q->mx[i][j] = 0.;
      }
    Q->mx[i][i] =  -1. * esl_vec_DSum(Q->mx[i], 4);
  }
  esl_rmx_ScaleTo(Q, pi, 1.0);
  return eslOK;
}

/*****************************************************************
 * 2. Debugging routines for validating or dumping rate matrices.
 *****************************************************************/

/* Function:  esl_rmx_ValidateP()
 * Incept:    SRE, Sun Mar 11 10:30:50 2007 [Janelia]
 *
 * Purpose:   Validates a conditional probability matrix <P>, whose
 *            elements $P_{ij}$ represent conditional probabilities
 *            $P(j \mid i)$; for example in a first-order Markov
 *            chain, or a continuous-time Markov transition process
 *            where <P> is for a particular $t$.
 *
 *            Rows must sum to one, and each element $P_{ij}$ is a
 *            probability $0 \leq P_{ij} \leq 1$.
 *
 *            <tol> specifies the floating-point tolerance to which
 *            the row sums must equal one: <fabs(sum-1.0) <= tol>.
 *
 *            <errbuf> is an optional error message buffer. The caller
 *            may pass <NULL> or a pointer to a buffer of at least
 *            <eslERRBUFSIZE> characters.
 *
 * Args:      P      - matrix to validate
 *            tol    - floating-point tolerance (0.00001, for example)
 *            errbuf - OPTIONAL: ptr to an error buffer of at least
 *                     <eslERRBUFSIZE> characters.
 *
 * Returns:   <eslOK> on successful validation.
 *            <eslFAIL> on failure, and if a non-<NULL> <errbuf> was
 *            provided by the caller, a message describing
 *            the reason for the failure is put there.
 *
 * Throws:    (no abnormal error conditions)
 */
int
esl_rmx_ValidateP(ESL_DMATRIX *P, double tol, char *errbuf)
{
  int    i,j;
  double sum;

  if (P->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "P must be type eslGENERAL to be validated");

  for (i = 0; i < P->n; i++)
    {
      sum = esl_vec_DSum(P->mx[i], P->m);
      if (fabs(sum-1.0) > tol) ESL_FAIL(eslFAIL, errbuf, "row %d does not sum to 1.0", i);

      for (j = 0; j < P->m; j++)
	if (P->mx[i][j] < 0.0 || P->mx[i][j] > 1.0)
	  ESL_FAIL(eslFAIL, errbuf, "element %d,%d is not a probability (%f)", i,j,P->mx[i][j]);
    }
  return eslOK;
}

/* Function:  esl_rmx_ValidateQ()
 * Incept:    SRE, Sun Mar 11 10:30:50 2007 [Janelia]
 *
 * Purpose:   Validates an instantaneous rate matrix <Q> for a
 *            continuous-time Markov process, whose elements $q_{ij}$
 *            represent instantaneous transition rates $i \rightarrow
 *            j$.
 *
 *            Rows satisfy the condition that
 *            $q_{ii} = -\sum_{i \neq j} q_{ij}$, and also
 *            that $q_{ij} \geq 0$ for all $j \neq i$.
 *
 *            <tol> specifies the floating-point tolerance to which
 *            that condition must hold: <fabs(sum-q_ii) <= tol>.
 *
 *            <errbuf> is an optional error message buffer. The caller
 *            may pass <NULL> or a pointer to a buffer of at least
 *            <eslERRBUFSIZE> characters.
 *
 * Args:      Q      - rate matrix to validate
 *            tol    - floating-point tolerance (0.00001, for example)
 *            errbuf - OPTIONAL: ptr to an error buffer of at least
 *                     <eslERRBUFSIZE> characters.
 *
 * Returns:   <eslOK> on successful validation.
 *            <eslFAIL> on failure, and if a non-<NULL> <errbuf> was
 *            provided by the caller, a message describing
 *            the reason for the failure is put there.
 *
 * Throws:    (no abnormal error conditions)
 */
int
esl_rmx_ValidateQ(ESL_DMATRIX *Q, double tol, char *errbuf)
{
  int    i,j;
  double qi;

  if (Q->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "Q must be type eslGENERAL to be validated");
  if (Q->n    != Q->m)       ESL_EXCEPTION(eslEINVAL, "a rate matrix Q must be square");

  for (i = 0; i < Q->n; i++)
    {
      qi = 0.;
      for (j = 0; j < Q->m; j++)
	{
	  if (i != j) {
	    if (Q->mx[i][j] < 0.)       ESL_FAIL(eslFAIL, errbuf, "offdiag elem %d,%d < 0",i,j);
	    qi += Q->mx[i][j];
	  } else {
	    if (Q->mx[i][j] > 0.)       ESL_FAIL(eslFAIL, errbuf, "diag elem %d,%d < 0", i,j);
	  }
	}
      if (fabs(qi + Q->mx[i][i]) > tol) ESL_FAIL(eslFAIL, errbuf, "row %d does not sum to 0.0", i);
    }
  return eslOK;
}

/*****************************************************************
 * 3. Other routines in the exposed ratematrix API.
 *****************************************************************/

/* Function:  esl_rmx_ScaleTo()
 * Incept:    SRE, Tue Jul 13 16:05:16 2004 [St. Louis]
 *
 * Purpose:   Rescales rate matrix <Q> so that expected substitution
 *            rate per dt is <unit>.
 *
 *            Expected substitution rate is:
 *               $\sum_i \sum_j pi_i Q_ij  \forall i \neq j$
 *
 *            <unit> typically taken to be 1.0, so time units are substitutions/site.
 *            An exception is PAM, where <unit> = 0.01 for 1 PAM unit.
 *
 * Args:      Q     - rate matrix to normalize
 *            pi    - stationary residue frequencies
 *            unit  - expected subsitution rate per dt
 *                    (1.0 = substitutions/site; 0.01 = PAMs)
 *
 * Returns:   <eslOK> on success, and matrix Q is rescaled.
 *
 * Xref:      STL8/p56.
 */
int
esl_rmx_ScaleTo(ESL_DMATRIX *Q, double *pi, double unit)
{
  int     i,j;
  double  sum = 0.;

  if (Q->m != Q->n || Q->type != eslGENERAL)
    ESL_EXCEPTION(eslEINVAL, "Q must be a square general matrix");

  for (i = 0; i < Q->m; i++)
    for (j = 0; j < Q->n; j++)
      if (i != j) sum += pi[i] * Q->mx[i][j];

  for (i = 0; i < Q->m; i++)
    for (j = 0; j < Q->n; j++)
      Q->mx[i][j] *= (unit / sum);

  return eslOK;
}

/* Function:  esl_rmx_E2Q()
 * Incept:    SRE, Tue Jul 13 15:52:41 2004 [St. Louis]
 *
 * Purpose:   Given a lower triangular matrix ($j<i$) of
 *            residue exchangeabilities <E>, and a stationary residue
 *            frequency vector <pi>; assuming $E_{ij} = E_{ji}$;
 *            calculates a rate matrix <Q> as
 *
 *            $Q_{ij} = E_{ij} * \pi_j$
 *
 *            The resulting <Q> is not normalized to any particular
 *            number of substitutions/site/time unit. See
 *            <esl_rmx_ScaleTo()> for that.
 *
 * Args:      E     - symmetric residue "exchangeabilities";
 *                    only lower triangular entries are used.
 *            pi    - residue frequencies at stationarity.
 *            Q     - RETURN: rate matrix, square (NxN).
 *                    Caller allocates the memory for this.
 *
 * Returns:   <eslOK> on success; Q is calculated and filled in.
 *
 * Xref:      STL8/p56.
 */
int
esl_rmx_E2Q(ESL_DMATRIX *E, double *pi, ESL_DMATRIX *Q)
{
  int          i,j;

  if (E->n != Q->n) ESL_EXCEPTION(eslEINVAL, "E and Q sizes differ");

  /* Scale all off-diagonals to pi[j] * E[i][j].
   */
  for (i = 0; i < E->n; i++)
    for (j = 0; j < i; j++)	/* only look at lower triangle of E. */
      {
	Q->mx[i][j] = pi[j] * E->mx[i][j];
	Q->mx[j][i] = pi[i] * E->mx[i][j];
      }

  /* Set diagonal to  -\sum of all j != i.
   */
  for (i = 0; i < Q->n; i++)
    {
      Q->mx[i][i] = 0.;		/* makes the vector sum work for j != i */
      Q->mx[i][i] = -1. * esl_vec_DSum(Q->mx[i], Q->n);
    }
  return eslOK;
}

/* Function:  esl_rmx_RelativeEntropy()
 * Incept:    SRE, Fri Mar 23 09:18:26 2007 [Janelia]
 *
 * Purpose:   Given a conditional substitution probability matrix <P>,
 *            with stationary probabilities <pi>, calculate its
 *            relative entropy $H$:
 *
 *               $H_t = \sum_{ij} P(j \mid i,t) \pi_i \log_2 \frac{P(j \mid i,t)} {\pi_j}$
 *
 *            This assumes that the stationary probabilities are the
 *            same as the background (null model) probabilities.
 *
 * Returns:   the relative entropy, $H$, in bits
 */
double
esl_rmx_RelativeEntropy(ESL_DMATRIX *P, double *pi)
{
  double H = 0.;
  int    i,j;

  for (i = 0; i < P->m; i++)
    for (j = 0; j < P->n; j++)
      H += P->mx[i][j] * pi[i] * log(P->mx[i][j] / pi[j]);
  return H / eslCONST_LOG2;
}

/* Function:  esl_rmx_ExpectedScore()
 * Incept:    SRE, Fri Mar 23 09:32:05 2007 [Janelia]
 *
 * Purpose:   Given a conditional substitution probability matrix <P>
 *            with stationary probabilities <pi>, calculate its
 *            expected score:
 *
 *               $ = \sum_{ij} \pi_j \pi_i \log_2 \frac{P(j \mid i,t)} {\pi_j}$
 *
 *            This assumes that the stationary probabilities are the
 *            same as the background (null model) probabilities.
 *
 * Returns:   the expected score, in bits
 */
double
esl_rmx_ExpectedScore(ESL_DMATRIX *P, double *pi)
{
  double S = 0.;
  int    i,j;

  for (i = 0; i < P->m; i++)
    for (j = 0; j < P->n; j++)
      S += pi[j] * pi[i] * log(P->mx[i][j] / pi[j]);
  return S / eslCONST_LOG2;
}

/*****************************************************************
 * 4. Benchmark driver
 *****************************************************************/

#ifdef eslRATEMATRIX_BENCHMARK

/*
  without GSL:
  gcc -O2 -I. -L. -o benchmark -DeslRATEMATRIX_BENCHMARK esl_ratematrix.c -leasel -lm

  with GSL:
  gcc -g -Wall -I. -L. -o benchmark -DeslRATEMATRIX_BENCHMARK -DHAVE_LIBGSL esl_dmatrix.c esl_ratematrix.c -leasel -lgsl -lgslcblas -lm
 */

#include <esl_config.h>

#ifdef HAVE_LIBGSL
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#endif

#include "easel.h"
#include "esl_stopwatch.h"
#include "esl_dmatrix.h"
#include "esl_ratematrix.h"

int
main(void)
{
  ESL_STOPWATCH *w = NULL;
  ESL_DMATRIX *Q  = NULL;
  ESL_DMATRIX *P  = NULL;
  double       t = 5.0;
  int          esl_iterations = 100;
  int          i;
#ifdef HAVE_LIBGSL
  gsl_matrix  *Qg = NULL;
  gsl_matrix  *Pg = NULL;
  int          gsl_iterations = 100;
#endif

  w = esl_stopwatch_Create();
  Q = esl_dmatrix_Create(20, 20);
  P = esl_dmatrix_Create(20, 20);
  esl_rmx_SetWAG(Q, NULL);

  esl_stopwatch_Start(w);
  for (i = 0; i < esl_iterations; i++)
    esl_dmx_Exp(Q, t, P);
  esl_stopwatch_Stop(w);
  printf("Easel takes:   %g sec\n", w->user / (double) esl_iterations);

#ifdef HAVE_LIBGSL
  if (esl_dmx_MorphGSL(Q, &Qg)             != eslOK) esl_fatal("morph to gsl_matrix failed");
  if ((Pg = gsl_matrix_alloc(20, 20))      == NULL)  esl_fatal("gsl alloc failed");
  gsl_matrix_scale(Qg, t);

  esl_stopwatch_Start(w);
  for (i = 0; i < gsl_iterations; i++)
    gsl_linalg_exponential_ss(Qg, Pg, GSL_PREC_DOUBLE);
  esl_stopwatch_Stop(w);
  printf("  GSL takes:   %g sec\n", w->user / (double) gsl_iterations);

  gsl_matrix_free(Qg);
  gsl_matrix_free(Pg);
#endif /*HAVE_LIBGSL*/

  esl_dmatrix_Destroy(Q);
  esl_dmatrix_Destroy(P);
  esl_stopwatch_Destroy(w);
  return 0;
}

#endif /*eslRATEMATRIX_BENCHMARK*/

/*****************************************************************
 * 5. Regression test driver
 *****************************************************************/
#ifdef eslRATEMATRIX_REGRESSION
#ifdef HAVE_LIBGSL

/* This tests rate matrix exponentiation against the GSL's
 * undocumented implementation of a matrix exponential.
 */
/*
  gcc -g -Wall -I. -L. -o ratematrix_regression -DeslRATEMATRIX_REGRESSION -DHAVE_LIBGSL esl_dmatrix.c esl_ratematrix.c -leasel -lgsl -lgslcblas -lm
 */

#include "esl_config.h"

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>

#include "easel.h"
#include "esl_dmatrix.h"
#include "esl_ratematrix.h"

int
main(void)
{
  char errbuf[eslERRBUFSIZE];
  char *alphabet = "ACDEFGHIKLMNPQRSTVWY";
  ESL_DMATRIX *Q  = NULL;
  ESL_DMATRIX *P  = NULL;
  gsl_matrix  *Qg = NULL;
  gsl_matrix  *Pg = NULL;
  ESL_DMATRIX *Pge = NULL;
  double       t = 15.0;

  if ((Q = esl_dmatrix_Create(20, 20))     == NULL)  esl_fatal("malloc failed");
  if ((P = esl_dmatrix_Create(20, 20))     == NULL)  esl_fatal("malloc failed");

  if (esl_rmx_SetWAG(Q, NULL)              != eslOK) esl_fatal("_SetWAG() failed");
  if (esl_rmx_ValidateQ(Q, 0.0001, errbuf) != eslOK) esl_fatal("Q validation failed: %s", errbuf);

  if (esl_dmx_Exp(Q, t, P)                 != eslOK) esl_fatal("matrix exponentiation failed");
  if (esl_rmx_ValidateP(P, 0.0001, errbuf) != eslOK) esl_fatal("P validation failed: %s", errbuf);

  if (esl_dmx_MorphGSL(Q, &Qg)             != eslOK) esl_fatal("morph to gsl_matrix failed");
  if ((Pg = gsl_matrix_alloc(20, 20))      == NULL)  esl_fatal("gsl alloc failed");
  gsl_matrix_scale(Qg, t);
  if (gsl_linalg_exponential_ss(Qg, Pg, GSL_PREC_DOUBLE) != 0) esl_fatal("gsl's exponentiation failed");
  if (esl_dmx_UnmorphGSL(Pg, &Pge)         != eslOK) esl_fatal("morph from gsl_matrix failed");

  esl_dmatrix_Dump(stdout, P, alphabet, alphabet);

  if (esl_dmatrix_Compare(Pge, P, 0.00001) != eslOK) esl_fatal("whoops, different answers.");

  esl_dmatrix_Destroy(Q);
  esl_dmatrix_Destroy(P);
  esl_dmatrix_Destroy(Pge);
  gsl_matrix_free(Qg);
  gsl_matrix_free(Pg);
  return 0;
}
#else
  /* if we don't have GSL, then compile in a dummy main(), solely
   *  to quiet any tests that are verifying that all drivers compile
   *  and run. */
int main(void) { return 0; }
#endif /*HAVE_LIBGSL*/

#endif /*eslRATEMATRIX_REGRESSION*/

/*****************************************************************
 * 6. Unit tests.
 *****************************************************************/
#ifdef eslRATEMATRIX_TESTDRIVE

static void
utest_SetWAG(void)
{
  char errbuf[eslERRBUFSIZE];
  ESL_DMATRIX *Q = NULL;
  ESL_DMATRIX *P = NULL;
  double       t = 50.0;	/* sufficiently large to drive e^tQ to stationarity  */
  double       pi[20];
  int          i;

  if ((Q = esl_dmatrix_Create(20, 20))     == NULL)  esl_fatal("malloc failed");
  if ((P = esl_dmatrix_Create(20, 20))     == NULL)  esl_fatal("malloc failed");

  /* This tests that exponentiating WAG gives a stable conditional
   * probability matrix solution. (It doesn't particularly test that
   * WAG was set correctly, but how could we have screwed that up?)
   */
  if (esl_rmx_SetWAG(Q, NULL)              != eslOK) esl_fatal("_SetWAG() failed");
  if (esl_dmx_Exp(Q, t, P)                 != eslOK) esl_fatal("matrix exponentiation failed");
  if (esl_rmx_ValidateP(P, 1e-7, errbuf)   != eslOK) esl_fatal("P validation failed: %s", errbuf);
  if (esl_rmx_ValidateQ(Q, 1e-7, errbuf)   != eslOK) esl_fatal("Q validation failed: %s", errbuf);

  /* This tests setting WAG to different stationary pi's than default,
   * then tests that exponentiating to large t reaches those stationaries.
   */
  esl_vec_DSet(pi, 20, 0.05);
  if (esl_rmx_SetWAG(Q, pi)                != eslOK) esl_fatal("_SetWAG() failed");
  if (esl_dmx_Exp(Q, t, P)                 != eslOK) esl_fatal("matrix exponentiation failed");
  if (esl_rmx_ValidateP(P, 1e-7, errbuf)   != eslOK) esl_fatal("P validation failed: %s", errbuf);
  if (esl_rmx_ValidateQ(Q, 1e-7, errbuf)   != eslOK) esl_fatal("Q validation failed: %s", errbuf);
  for (i = 0; i < 20; i++)
    if (esl_vec_DCompare(P->mx[i], pi, 20, 1e-7) != eslOK) esl_fatal("P didn't converge to right pi's");

  esl_dmatrix_Destroy(Q);
  esl_dmatrix_Destroy(P);
  return;
}

#ifdef HAVE_LIBLAPACK
static void
utest_Diagonalization(void)
{
  ESL_DMATRIX *P      = NULL;
  ESL_DMATRIX *P2     = NULL;
  ESL_DMATRIX *C      = NULL;
  ESL_DMATRIX *D      = NULL;
  double      *lambda = NULL;		/* eigenvalues */
  ESL_DMATRIX *U      = NULL;		/* left eigenvectors */
  ESL_DMATRIX *Ui     = NULL;		/* inverse of U */
  int  i,j;

  /* Create a J/C probability matrix for t=1:
   *    1/4 + 3/4 e^{-4/3 at}
   *    1/4 - 1/4 e^{-4/3 at}
   */
  if ((P  = esl_dmatrix_Create(4, 4))    == NULL)  esl_fatal("malloc failed");
  if ((C  = esl_dmatrix_Create(4, 4))    == NULL)  esl_fatal("malloc failed");
  if ((Ui = esl_dmatrix_Create(4, 4))    == NULL)  esl_fatal("malloc failed");
  if ((D  = esl_dmatrix_Create(4, 4))    == NULL)  esl_fatal("malloc failed");
  if ((P2 = esl_dmatrix_Create(4, 4))    == NULL)  esl_fatal("malloc failed");
  for (i = 0; i < 4; i++)
    for (j = 0; j < 4; j++)
      if (i == j) P->mx[i][j] = 0.25 + 0.75 * exp(-4./3.);
      else        P->mx[i][j] = 0.25 - 0.25 * exp(-4./3.);

  /* Diagonalize it
   */
  if (esl_dmx_Diagonalize(P, &lambda, NULL, &U, NULL) != eslOK) esl_fatal("diagonalization failed");

  /* Calculate P^k by U [diag(lambda_i)]^k U^{-1}
   */
  esl_dmatrix_SetZero(D);
  for (i = 0; i < P->n; i++) D->mx[i][i] = lambda[i];
  esl_dmx_Invert(U, Ui);
  esl_dmx_Multiply(U, D,  C);
  esl_dmx_Multiply(C, Ui, P2);

  if (esl_dmatrix_Compare(P, P2, 1e-7) != eslOK) esl_fatal("diagonalization unit test failed");

  free(lambda);
  esl_dmatrix_Destroy(P2);
  esl_dmatrix_Destroy(Ui);
  esl_dmatrix_Destroy(U);
  esl_dmatrix_Destroy(D);
  esl_dmatrix_Destroy(C);
  esl_dmatrix_Destroy(P);
  return;
}
#endif /*HAVE_LIBLAPACK*/

#endif /*eslRATEMATRIX_TESTDRIVE*/

/*****************************************************************
 * 7. Test driver
 *****************************************************************/

#ifdef eslRATEMATRIX_TESTDRIVE
/* gcc -g -Wall -o test -I. -L. -DeslRATEMATRIX_TESTDRIVE esl_ratematrix.c -leasel -lm
 * ./test
 *
 * gcc -g -Wall -o test -I. -L. -DHAVE_LIBLAPACK -DeslRATEMATRIX_TESTDRIVE esl_ratematrix.c esl_dmatrix.c -leasel -llapack -lm
 */
#include "esl_config.h"

#include "easel.h"
#include "esl_ratematrix.h"

int
main(void)
{
  utest_SetWAG();
#ifdef HAVE_LIBLAPACK
  utest_Diagonalization();
#endif

  return 0;

}
#endif /*eslRATEMATRIX_TESTDRIVE*/

/*****************************************************************
 * 8. Example driver
 *****************************************************************/

/*****************************************************************
 * Easel - a library of C functions for biological sequence analysis
 * Version h3.1b2; February 2015
 * Copyright (C) 2015 Howard Hughes Medical Institute.
 * Other copyrights also apply. See the COPYRIGHT file for a full list.
 *
 * Easel is distributed under the Janelia Farm Software License, a BSD
 * license. See the LICENSE file for more details.
 *
 * SVN $Id: esl_ratematrix.c 849 2013-01-31 15:10:24Z eddys $
 * SVN $URL: https://svn.janelia.org/eddylab/eddys/easel/branches/hmmer/3.1/esl_ratematrix.c $
 *****************************************************************/


text:
/* spdgemm.c
 *
 * Copyright (C) 2014 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_spblas.h>
#include <gsl/gsl_errno.h>

/*
gsl_spblas_dgemm()
  Multiply two sparse matrices

Inputs: alpha - scalar factor
        A     - sparse matrix
        B     - sparse matrix
        C     - (output) C = alpha * A * B

Return: success or error

Notes:
1) based on CSparse routine cs_multiply
*/

int
gsl_spblas_dgemm(const double alpha, const gsl_spmatrix *A,
                 const gsl_spmatrix *B, gsl_spmatrix *C)
{
  if (A->size2 != B->size1 || A->size1 != C->size1 || B->size2 != C->size2)
    {
      GSL_ERROR("matrix dimensions do not match", GSL_EBADLEN);
    }
  else if (A->sptype != B->sptype || A->sptype != C->sptype)
    {
      GSL_ERROR("matrix storage formats do not match", GSL_EINVAL);
    }
  else if (!GSL_SPMATRIX_ISCCS(A))
    {
      GSL_ERROR("compressed column format required", GSL_EINVAL);
    }
  else
    {
      int status = GSL_SUCCESS;
      const size_t M = A->size1;
      const size_t N = B->size2;
      size_t *Bi = B->i;
      size_t *Bp = B->p;
      double *Bd = B->data;
      size_t *w = (size_t *) A->work; /* workspace of length M */
      double *x = (double *) C->work; /* workspace of length M */
      size_t *Cp, *Ci;
      double *Cd;
      size_t j, p;
      size_t nz = 0;

      if (C->nzmax < A->nz + B->nz)
        {
          status = gsl_spmatrix_realloc(A->nz + B->nz, C);
          if (status)
            {
              GSL_ERROR("unable to realloc matrix C", status);
            }
        }

      /* initialize workspace to 0 */
      for (j = 0; j < M; ++j)
        w[j] = 0;

      Cp = C->p;
      Ci = C->i;
      Cd = C->data;

      for (j = 0; j < N; ++j)
        {
          if (nz + M > C->nzmax)
            {
              status = gsl_spmatrix_realloc(2 * C->nzmax + M, C);
              if (status)
                {
                  GSL_ERROR("unable to realloc matrix C", status);
                }

              /* these pointers could have changed due to reallocation */
              Ci = C->i;
              Cd = C->data;
            }

          Cp[j] = nz; /* column j of C starts here */

          for (p = Bp[j]; p < Bp[j + 1]; ++p)
            {
              nz = gsl_spblas_scatter(A, Bi[p], Bd[p], w, x, j + 1, C, nz);
            }

          for (p = Cp[j]; p < nz; ++p)
            Cd[p] = x[Ci[p]];
        }

      Cp[N] = nz;
      C->nz = nz;

      /* scale by alpha */
      gsl_spmatrix_scale(C, alpha);

      return status;
    }
} /* gsl_spblas_dgemm() */

/*
gsl_spblas_scatter()

  Keep a running total x -> x + alpha*A(:,j) for adding matrices together in CCS,
which will eventually be stored in C(:,j)

  When a new non-zero element with row index i is found, update C->i with
the row index. C->data is updated only by the calling function after all
matrices have been added via this function.

Inputs: A     - sparse matrix m-by-n
        j     - column index
        alpha - scalar factor
        w     - keeps track which rows of column j have been added to C;
                initialize to 0 prior to first call
        x     - column vector of length m
        mark  -
        C     - output matrix whose jth column will be added to A(:,j)
        nz    - (input/output) number of non-zeros in matrix C

Notes:
1) This function is designed to be called successively when adding multiple
matrices together. Column j of C is stored contiguously as per CCS but not
necessarily in order - ie: the row indices C->i may not be in ascending order.

2) based on CSparse routine cs_scatter
*/

size_t
gsl_spblas_scatter(const gsl_spmatrix *A, const size_t j, const double alpha,
                   size_t *w, double *x, const size_t mark, gsl_spmatrix *C,
                   size_t nz)
{
  size_t p;
  size_t *Ai = A->i;
  size_t *Ap = A->p;
  double *Ad = A->data;
  size_t *Ci = C->i;

  for (p = Ap[j]; p < Ap[j + 1]; ++p)
    {
      size_t i = Ai[p];          /* A(i,j) is nonzero */

      if (w[i] < mark)           /* check if row i has been stored in column j yet */
        {
          w[i] = mark;           /* i is new entry in column j */
          Ci[nz++] = i;          /* add i to pattern of C(:,j) */
          x[i] = alpha * Ad[p];  /* x(i) = alpha * A(i,j) */
        }
      else                       /* this (i,j) exists in C from a previous call */
        {
          x[i] += alpha * Ad[p]; /* add alpha*A(i,j) to C(i,j) */
        }
    }

  return (nz) ;
} /* gsl_spblas_scatter() */

ext: c

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>
#include <sys/time.h>
#include "../ulmblas.h"
#define M 1024
#define K 128
#define N 1024

void
ULMBLAS(dgemm)(const enum Trans  transA,
               const enum Trans  transB,
               const int         m,
               const int         n,
               const int         k,
               const double      alpha,
               const double      *A,
               const int         ldA,
               const double      *B,
               const int         ldB,
               const double      beta,
               double            *C,
               const int         ldC);
int main(){
  // init
  double (*a)[K] = (double (*)[K])malloc(sizeof(double)*M*K);
  double (*b)[N] = (double (*)[N])malloc(sizeof(double)*N*K);
  double (*c)[N] = (double (*)[N])malloc(sizeof(double)*M*N);
  struct timeval stop, start;
  for(int i=0;i<M;i++){
    for(int j=0;j<K;j++){
      a[i][j] = 0.33*(float)i+0.29*(float)j;
    }
  }
  for(int i=0;i<K;i++){
    for(int j=0;j<N;j++){
      b[i][j] = 0.67*(float)i+0.98*(float)j;
    }
  }
  //C=A*B^T
  // openblas
  for(int i=0;i<M;i++){
    for(int j=0;j<N;j++){
      c[i][j] = 0;
    }
  }
  gettimeofday(&start, NULL);
  cblas_dgemm(CblasRowMajor,
              CblasNoTrans,
              CblasNoTrans,
              M,N,K,1,(double*)a,K,(double*)b,N,0,(double*)c,N);
  gettimeofday(&stop, NULL);
  printf("li took %lu get %lf\n", (stop.tv_sec - start.tv_sec) * 1000000 + stop.tv_usec - start.tv_usec, c[1][1]);
  // me
  /*
  for(int i=0;i<M;i++){
    for(int j=0;j<N;j++){
      c[i][j] = 0;
    }
  }
  gettimeofday(&start, NULL);
  for(int i=0;i<M;i++)
    for(int j=0;j<N;j++){
      double *cp = &c[i][j];
      double *ap = a[i];
      double *bp = b[j];
      for(int p=0;p<K;p++){
        *cp += *(ap++) * *(bp++);
      }
    }
  gettimeofday(&stop, NULL);
  printf("me took %lu get %lf\n", (stop.tv_sec - start.tv_sec) * 1000000 + stop.tv_usec - start.tv_usec, c[1][1]);
  */
  // ulm
  for(int i=0;i<M;i++){
    for(int j=0;j<N;j++){
      c[i][j] = 0;
    }
  }
  gettimeofday(&start, NULL);
  ULMBLAS(dgemm)(Trans,
         Trans,
         M,N,K,1,(double*)a,K,(double*)b,N,0,(double*)c,N);
  gettimeofday(&stop, NULL);
  printf("ul took %lu get %lf\n", (stop.tv_sec - start.tv_sec) * 1000000 + stop.tv_usec - start.tv_usec, c[1][1]);
  return 0;
}

max_stars_repo_name: hzhangxyz/ulmBLAS

max_issues_repo_name: hzhangxyz/ulmBLAS

max_forks_repo_name: hzhangxyz/ulmBLAS


text:
#include <gsl/gsl_linalg.h>

/* user defined functions */

ext: h

max_stars_repo_name: yipf/lua-graphics

max_issues_repo_name: yipf/lua-graphics

max_forks_repo_name: yipf/lua-graphics


text:
/**
 * @file feature.h
 * @brief Speech recognition front end.
 * @author John McDonough, Tobias Gehrig, Kenichi Kumatani, Friedrich Faubel
 */

#ifndef FEATURE_H
#define FEATURE_H

#include <math.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_vector_complex.h>
#include "matrix/gslmatrix.h"
#include "stream/stream.h"
#include "common/mlist.h"

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <errno.h>
#include <string.h>
#include <sys/types.h>
#include <netinet/in.h>
#include <sys/socket.h>
#include <sys/time.h>
#include <netdb.h>

#include <pthread.h>

#include "btk.h"

#ifdef HAVE_LIBFFTW3
#include <fftw3.h>
#endif /* #ifdef HAVE_LIBFFTW3 */

#include "feature/spectralestimator.h"

/* sflib and sndfile both define sf_perror so we put them into seperate namespaces */
namespace sndfile {
#include <sndfile.h>
}

void unpack_half_complex(gsl_vector_complex* tgt, const double* src);
void unpack_half_complex(gsl_vector_complex* tgt, const unsigned N2, const double* src, const unsigned N);
void pack_half_complex(double* tgt, const gsl_vector_complex* src, unsigned size = 0);

/**
* \defgroup AudioFeature Audio Feature Hierarchy
* This hierarchy of classes provides the capability for extracting audio
* features for use in automatic speech recognition.
*/
/*@{*/

/**
* \defgroup FileFeature File Feature.
*/
/*@{*/

// ----- definition for class `FileFeature' -----
//
class FileFeature : public VectorFloatFeatureStream {
 public:
  FileFeature(unsigned sz, const String& nm = "FileFeature") :
    VectorFloatFeatureStream(sz, nm), feature_(NULL) {}

  virtual ~FileFeature() {
    if (feature_ != NULL)
      gsl_matrix_float_free(feature_);
  }

  FileFeature& operator=(const FileFeature& f);

  virtual const gsl_vector_float* next(int frame_no = -5);

  unsigned size() const {
    if (GSL_MATRIX_NCOLS(feature_) == 0)
      throw j_error("Matrix not loaded yet.");
    return (unsigned) GSL_MATRIX_NCOLS(feature_);
  }

  void bload(const String& fileName, bool old = false);

  void copy(gsl_matrix_float* matrix);

 private:
  gsl_matrix_float*  feature_;
};

typedef Inherit<FileFeature, VectorFloatFeatureStreamPtr> FileFeaturePtr;

/*@}*/

/**
* \defgroup ConversionbitToShort Conversion bit Short
*/
/*@{*/

// ----- definition of 'Conversion24bit2Short' -----
//
class Conversion24bit2Short : public VectorShortFeatureStream {
 public:
 Conversion24bit2Short(VectorCharFeatureStreamPtr& src,
                       const String& nm = "Conversion24bit2Short") :
  VectorShortFeatureStream(src->size()/3, nm), src_(src) {};
  virtual void reset() { src_->reset(); VectorShortFeatureStream::reset(); }

  virtual const gsl_vector_short* next(int frame_no = -5);
 private:
  VectorCharFeatureStreamPtr src_;
};

typedef Inherit<Conversion24bit2Short, VectorShortFeatureStreamPtr> Conversion24bit2ShortPtr;

/*@}*/

/**
* \defgroup Conversion24bit2Float Conversion 24 bit 2 Float
*/
/*@{*/

// ----- definition of Conversion24bit2Float -----
//
class Conversion24bit2Float : public VectorFloatFeatureStream {
 public:
  Conversion24bit2Float(VectorCharFeatureStreamPtr& src,
                        const String& nm = "Conversion24bit2Float") :
    VectorFloatFeatureStream(src->size()/3, nm), src_(src) {};
  virtual void reset() { src_->reset(); VectorFloatFeatureStream::reset(); }

  virtual const gsl_vector_float* next(int frame_no = -5);
 private:
  VectorCharFeatureStreamPtr src_;
};

typedef Inherit<Conversion24bit2Float, VectorFloatFeatureStreamPtr> Conversion24bit2FloatPtr;

/**
* \defgroup SampleFeature Sample Feature
*/
/*@{*/

// ----- definition for class `SampleFeature' -----
//
class SampleFeature;
typedef Inherit<SampleFeature, VectorFloatFeatureStreamPtr> SampleFeaturePtr;
class SampleFeature : public VectorFloatFeatureStream {
 public:
  SampleFeature(const String& fn = "", unsigned blockLen = 320,
                unsigned shiftLen = 160, bool padZeros = false, const String& nm = "Sample");
  virtual ~SampleFeature();

  unsigned read(const String& fn, int format = 0, int samplerate = 16000,
                int chX = 1, int chN = 1, int cfrom = 0, int to = -1, int outsamplerate = -1, float norm = 0.0);

  void write(const String& fn, int format = sndfile::SF_FORMAT_WAV|sndfile::SF_FORMAT_PCM_16, int sampleRate = -1);

  void cut(unsigned cfrom, unsigned cto);

  void randomize(int startX, int endX, double sigma2);

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { cur_ = 0; VectorFloatFeatureStream::reset(); is_end_ = false; }

  void exit(){ reset(); throw jiterator_error("end of samples!");}

  const gsl_vector_float* data();

  const gsl_vector* dataDouble();

  void copySamples(SampleFeaturePtr& src, unsigned cfrom, unsigned to);

  unsigned samplesN() const { return ttlsamples_; }

  int getSampleRate() const { return samplerate_; }

  int getChanN() const { return nChan_; }

  void zeroMean();

  void addWhiteNoise( float snr );

  void setSamples(const gsl_vector* samples, unsigned sampleRate);

protected:
  float*              samples_;
  float               norm_;
  unsigned            ttlsamples_;
  const unsigned      shiftLen_;
  int                 samplerate_;
  int                 nChan_;
  int                 format_;
  unsigned            cur_;
  bool                pad_zeros_;
  gsl_vector_float*   copy_fsamples_;
  gsl_vector*         copy_dsamples_;

private:
  SampleFeature(const SampleFeature& s);
  SampleFeature& operator=(const SampleFeature& s);
};

// ----- definition for class `SampleFeatureRunon' -----
//
class SampleFeatureRunon : public SampleFeature {
 public:
  SampleFeatureRunon(const String& fn = "", unsigned blockLen = 320,
                     unsigned shiftLen = 160, bool padZeros = false, const String& nm = "Sample") :
    SampleFeature(fn, blockLen, shiftLen, padZeros, nm) { }

  virtual void reset() { VectorFloatFeatureStream::reset(); }

  virtual int frame_no() const { return (cur_ / shiftLen_) - 1; }

  virtual int frameN() const { return (ttlsamples_ / shiftLen_) - 1; }
};

typedef Inherit<SampleFeatureRunon, SampleFeaturePtr> SampleFeatureRunonPtr;

/*@}*/

/**
* \defgroup IterativeSampleFeature Iterative Sample Feature for the single channel data
*/
/*@{*/

// ----- definition for class `IterativeSingleChannelSampleFeature' -----
//
class IterativeSingleChannelSampleFeature : public VectorFloatFeatureStream {
 public:
 public:
  IterativeSingleChannelSampleFeature(unsigned blockLen = 320, const String& nm = "IterativeSingleChannelSampleFeature");
  virtual ~IterativeSingleChannelSampleFeature();

  void read(const String& fileName, int format = 0, int samplerate = 44100, int cfrom = 0, int cto = -1 );

  unsigned samplesN() const { return ttlsamples_; }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset();

 private:

  float*				samples_;
  sndfile::SNDFILE*			sndfile_;
  sndfile::SF_INFO			sfinfo_;
  unsigned				interval_;
  unsigned				blockN_;
  unsigned				sampleN_;
  unsigned				ttlsamples_;

  const unsigned			blockLen_;
  unsigned				cur_;
  bool					last_;
  int					cto_;
};

typedef Inherit<IterativeSingleChannelSampleFeature, VectorFloatFeatureStreamPtr> IterativeSingleChannelSampleFeaturePtr;

/*@}*/

/**
* \defgroup IterativeSampleFeature Iterative Sample Feature
*/
/*@{*/

// ----- definition for class `IterativeSampleFeature' -----
//
class IterativeSampleFeature : public VectorFloatFeatureStream {
 public:
  IterativeSampleFeature(unsigned chX, unsigned blockLen = 320, unsigned firstChanX = 0, const String& nm = "Iterative Sample");
  virtual ~IterativeSampleFeature();

  void read(const String& fileName, int format = 0, int samplerate = 44100, int chN = 1, int cfrom = 0, int cto = -1 );

  unsigned samplesN() const { return ttlsamples_; }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset();

  void changeFirstChannelID( unsigned firstChanX ){ firstChanX_ = firstChanX; }

 private:
  IterativeSampleFeature(const IterativeSampleFeature& s);
  IterativeSampleFeature& operator=(const IterativeSampleFeature& s);

  static float*					allSamples_;
  static sndfile::SNDFILE*			sndfile_;
  static sndfile::SF_INFO			sfinfo_;
  static unsigned				interval_;
  static unsigned				blockN_;
  static unsigned				sampleN_;
  static unsigned				allSampleN_;
  static unsigned				ttlsamples_;

  const unsigned				blockLen_;
  const unsigned				chanX_;
  unsigned				        firstChanX_;
  unsigned					cur_;
  bool						last_;
  int						cto_;
};

typedef Inherit<IterativeSampleFeature, VectorFloatFeatureStreamPtr> IterativeSampleFeaturePtr;

/*@}*/

/**
* \defgroup BlockSizeConversionFeature Block Size Conversion Feature
*/
/*@{*/

// ----- definition for class `BlockSizeConversionFeature' -----
//
class BlockSizeConversionFeature : public VectorFloatFeatureStream {
 public:
  BlockSizeConversionFeature(const VectorFloatFeatureStreamPtr& src,
                             unsigned blockLen = 320,
                             unsigned shiftLen = 160,
                             const String& nm = "BlockSizeConversionFeature");

  virtual void reset() { curin_ = curout_ = 0; src_frame_no_ = -1; src_->reset();  VectorFloatFeatureStream::reset(); }

  virtual const gsl_vector_float* next(int frame_no = -5);

 private:
  void inputLonger_();
  void outputLonger_();

  VectorFloatFeatureStreamPtr			src_;
  const unsigned				inputLen_;
  const unsigned				blockLen_;
  const unsigned				shiftLen_;
  const unsigned				overlapLen_;
  unsigned					curin_;
  unsigned					curout_;
  int                                           src_frame_no_;
  const gsl_vector_float*			src_feat_;
};

typedef Inherit<BlockSizeConversionFeature, VectorFloatFeatureStreamPtr> BlockSizeConversionFeaturePtr;

/*@}*/

/**
* \defgroup BlockSizeConversionFeatureShort Block Size Conversion Feature Short
*/
/*@{*/

// ----- definition for class `BlockSizeConversionFeatureShort' -----
//
class BlockSizeConversionFeatureShort : public VectorShortFeatureStream {
 public:
  BlockSizeConversionFeatureShort(VectorShortFeatureStreamPtr& src,
				  unsigned blockLen = 320,
				  unsigned shiftLen = 160, const String& nm = "Block Size Conversion");

  virtual void reset() { curin_ = curout_ = 0;  src_->reset();  VectorShortFeatureStream::reset(); }

  virtual const gsl_vector_short* next(int frame_no = -5);

 private:
  void inputLonger_();
  void outputLonger_();

  VectorShortFeatureStreamPtr			src_;
  const unsigned				inputLen_;
  const unsigned				blockLen_;
  const unsigned				shiftLen_;
  const unsigned				overlapLen_;
  unsigned					curin_;
  unsigned					curout_;
  const gsl_vector_short*			src_feat_;
};

typedef Inherit<BlockSizeConversionFeatureShort, VectorShortFeatureStreamPtr> BlockSizeConversionFeatureShortPtr;

/*@}*/

#ifdef SMARTFLOW

namespace sflib {
#include "sflib.h"
}

/**
* \defgroup SmartFlowFeature Smart Flow Feature
*/
/*@{*/

// ----- definition for class `SmartFlowFeature' -----
//
class SmartFlowFeature : public VectorShortFeatureStream {
 public:
  SmartFlowFeature(sflib::sf_flow_sync* sfflow, unsigned blockLen = 320,
                   unsigned shiftLen = 160, const String& nm = "SmartFloatFeature") :
    VectorShortFeatureStream(blockLen, nm),
    sfflow_(sfflow), blockLen_(blockLen), shiftLen_(shiftLen) { }

  virtual ~SmartFlowFeature() { }

  virtual const gsl_vector_short* next(int frame_no = -5);

  virtual void reset() { VectorShortFeatureStream::reset(); }

 private:
  sflib::sf_flow_sync* 				sfflow_;
  const unsigned				blockLen_;
  const unsigned				shiftLen_;
};

typedef Inherit<SmartFlowFeature, VectorShortFeatureStreamPtr> SmartFlowFeaturePtr;

#endif

/*@}*/

/**
* \defgroup PreemphasisFeature Preemphasis Feature
*/
/*@{*/

// ----- definition for class `PreemphasisFeature' -----
//
class PreemphasisFeature : public VectorFloatFeatureStream {
 public:
 PreemphasisFeature(const VectorFloatFeatureStreamPtr& samp, double mu = 0.95, const String& nm = "Preemphasis");
  virtual ~PreemphasisFeature() { }

  virtual const gsl_vector_float* next(int frame_no = -5);
  virtual void reset();
  void next_speaker();

#ifdef ENABLE_LEGACY_BTK_API
  void nextSpeaker(){ next_speaker(); }
#endif

 private:
  VectorFloatFeatureStreamPtr			samp_;
  float						prior_;
  const double					mu_;
};

typedef Inherit<PreemphasisFeature, VectorFloatFeatureStreamPtr> PreemphasisFeaturePtr;

/*@}*/

/**
* \defgroup HammingFeatureShort Hamming Feature Short
*/
/*@{*/

// ----- definition for class `HammingFeatureShort' -----
//
class HammingFeatureShort : public VectorFloatFeatureStream {
 public:
  HammingFeatureShort(const VectorShortFeatureStreamPtr& samp, const String& nm = "HammingShort");
  virtual ~HammingFeatureShort() { delete[] window_; }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { samp_->reset(); VectorFloatFeatureStream::reset(); }

 private:
  VectorShortFeatureStreamPtr			samp_;
  unsigned					windowLen_;
  double*					window_;
};

typedef Inherit<HammingFeatureShort, VectorFloatFeatureStreamPtr> HammingFeatureShortPtr;

/*@}*/

/**
* \defgroup HammingFeature Hamming Feature
*/
/*@{*/

// ----- definition for class `HammingFeature' -----
//
class HammingFeature : public VectorFloatFeatureStream {
 public:
  HammingFeature(const VectorFloatFeatureStreamPtr& samp, const String& nm = "Hamming");
  virtual ~HammingFeature() { delete[] window_; }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { samp_->reset(); VectorFloatFeatureStream::reset(); }

 private:
  VectorFloatFeatureStreamPtr			samp_;
  unsigned					windowLen_;
  double*					window_;
};

typedef Inherit<HammingFeature, VectorFloatFeatureStreamPtr> HammingFeaturePtr;

/*@}*/

/**
* \defgroup FFTFeature FFT Feature
*/
/*@{*/

// ----- definition for class `FFTFeature' -----
//
class FFTFeature : public VectorComplexFeatureStream {
 public:
  FFTFeature(const VectorFloatFeatureStreamPtr& samp, unsigned fftLen, const String& nm = "FFT");
  virtual ~FFTFeature();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset() { samp_->reset(); VectorComplexFeatureStream::reset(); }

  unsigned fftLen()    const { return fftLen_;    }
  unsigned windowLen() const { return windowLen_; }

  unsigned nBlocks()     const { return 4; }
  unsigned subsamplerate() const { return 2; }

#ifdef ENABLE_LEGACY_BTK_API
  unsigned subSampRate() { return subsamplerate(); }
#endif

 private:
  VectorFloatFeatureStreamPtr			samp_;
  unsigned					fftLen_;
  unsigned					windowLen_;
  double*					samples_;

#ifdef HAVE_LIBFFTW3
  fftw_plan					fftwPlan_;
  fftw_complex*					output_;
#endif
};

typedef Inherit<FFTFeature, VectorComplexFeatureStreamPtr> FFTFeaturePtr;

/**
* \defgroup SpectralPowerFeature Spectral Power Feature
*/
/*@{*/

// ----- definition for class `SpectralPowerFeature' -----
//
class SpectralPowerFloatFeature : public VectorFloatFeatureStream {
 public:
  SpectralPowerFloatFeature(const VectorComplexFeatureStreamPtr& fft, unsigned powN = 0, const String& nm = "PowerFloat");

  virtual ~SpectralPowerFloatFeature() { }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { fft_->reset(); VectorFloatFeatureStream::reset(); }

 private:
  VectorComplexFeatureStreamPtr    fft_;
};

typedef Inherit<SpectralPowerFloatFeature, VectorFloatFeatureStreamPtr> SpectralPowerFloatFeaturePtr;

/*@}*/

/**
* \defgroup SpectralPowerFeature Spectral Power Feature
*/
/*@{*/

// ----- definition for class `SpectralPowerFeature' -----
//
class SpectralPowerFeature : public VectorFeatureStream {
 public:
  SpectralPowerFeature(const VectorComplexFeatureStreamPtr& fft, unsigned powN = 0, const String& nm = "Power");

  virtual ~SpectralPowerFeature() { }

  virtual const gsl_vector* next(int frame_no = -5);

  virtual void reset() { fft_->reset(); VectorFeatureStream::reset(); }

 private:
  VectorComplexFeatureStreamPtr			fft_;
};

typedef Inherit<SpectralPowerFeature, VectorFeatureStreamPtr> SpectralPowerFeaturePtr;

/*@}*/

/**
* \defgroup SignalPowerFeature Signal Power Feature
*/
/*@{*/

// ----- definition for class `SignalPowerFeature' -----
//
static const int ADCRANGE = 65536;

class SignalPowerFeature : public VectorFloatFeatureStream {
 public:
  SignalPowerFeature(const VectorFloatFeatureStreamPtr& samp, const String& nm = "Signal Power") :
    VectorFloatFeatureStream(/* size= */ 1, nm), samp_(samp), range_(float(ADCRANGE) * float(ADCRANGE) / 4.0) { }

  virtual ~SignalPowerFeature() { }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { samp_->reset(); VectorFloatFeatureStream::reset(); }

 private:
  VectorFloatFeatureStreamPtr			samp_;
  double					range_;
};

typedef Inherit<SignalPowerFeature, VectorFloatFeatureStreamPtr> SignalPowerFeaturePtr;

/*@}*/

/**
* \defgroup ALogFeature A-Log Feature
*/
/*@{*/

// ----- definition for class `ALogFeature' -----
//
class ALogFeature : public VectorFloatFeatureStream {
 public:
  ALogFeature(const VectorFloatFeatureStreamPtr& samp, double m = 1.0, double a = 4.0,
              bool runon = false, const String& nm = "ALogPower") :
    VectorFloatFeatureStream(/* size= */ 1, nm), samp_(samp), m_(m), a_(a),
      min_(HUGE_VAL), max_(-HUGE_VAL), minMaxFound_(false), runon_(runon) { }

  virtual ~ALogFeature() { }

  virtual const gsl_vector_float* next(int frame_no = -5);
  virtual void reset();
  void next_speaker() { min_ = HUGE_VAL; max_ = -HUGE_VAL; minMaxFound_ = false; }

#ifdef ENABLE_LEGACY_BTK_API
  void nextSpeaker(){ next_speaker(); }
#endif

 private:
  void find_min_max_(const gsl_vector_float* block);

  VectorFloatFeatureStreamPtr			samp_;
  double					m_;
  double					a_;
  double					min_;
  double					max_;
  bool						minMaxFound_;
  bool						runon_;
};

typedef Inherit<ALogFeature, VectorFloatFeatureStreamPtr> ALogFeaturePtr;

/*@}*/

/**
* \defgroup NormalizeFeature Normalize Feature
*/
/*@{*/

// ----- definition for class `NormalizeFeature' -----
//
class NormalizeFeature : public VectorFloatFeatureStream {
 public:
  NormalizeFeature(const VectorFloatFeatureStreamPtr& samp, double min = 0.0, double max = 1.0,
                   bool runon = false, const String& nm = "Normalize");

  virtual const gsl_vector_float* next(int frame_no = -5);
  virtual void reset();
  void next_speaker() { xmin_ = HUGE_VAL; xmax_ = -HUGE_VAL; minMaxFound_ = false; }

#ifdef ENABLE_LEGACY_BTK_API
  void nextSpeaker(){ next_speaker(); }
#endif
 private:
  void find_min_max_(const gsl_vector_float* block);

  VectorFloatFeatureStreamPtr			samp_;
  double					min_;
  double					max_;
  double					range_;

  double					xmin_;
  double					xmax_;
  bool						minMaxFound_;
  bool						runon_;
};

typedef Inherit<NormalizeFeature, VectorFloatFeatureStreamPtr> NormalizeFeaturePtr;

/*@}*/

/**
* \defgroup ThresholdFeature Threshold Feature
*/
/*@{*/

// ----- definition for class `ThresholdFeature' -----
//
class ThresholdFeature : public VectorFloatFeatureStream {
 public:
  ThresholdFeature(const VectorFloatFeatureStreamPtr& samp, double value = 0.0, double thresh = 1.0,
		   const String& mode = "upper", const String& nm = "Threshold");

  virtual ~ThresholdFeature() { }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { samp_->reset(); VectorFloatFeatureStream::reset(); }

 private:
  VectorFloatFeatureStreamPtr			samp_;
  double					value_;
  double					thresh_;
  int						compare_;
};

typedef Inherit<ThresholdFeature, VectorFloatFeatureStreamPtr> ThresholdFeaturePtr;

/*@}*/

/**
* \defgroup SpectralResamplingFeature Spectral Resampling Feature
*/
/*@{*/

// ----- definition for class `SpectralResamplingFeature' -----
//
class SpectralResamplingFeature : public VectorFeatureStream {
  static const double SampleRatio;
 public:
  SpectralResamplingFeature(const VectorFeatureStreamPtr& src, double ratio = SampleRatio, unsigned len = 0,
			    const String& nm = "Resampling");

  virtual ~SpectralResamplingFeature();

  virtual const gsl_vector* next(int frame_no = -5);

  virtual void reset() { src_->reset(); VectorFeatureStream::reset(); }

 private:
  VectorFeatureStreamPtr			src_;
  const double					ratio_;
};

typedef Inherit<SpectralResamplingFeature, VectorFeatureStreamPtr> SpectralResamplingFeaturePtr;

/*@}*/

/**
* \defgroup SamplerateConversionFeature Samplerate Conversion Feature
*/
/*@{*/

// ----- definition for class `SamplerateConversionFeature' -----
//
#ifdef SRCONV

#include <samplerate.h>

class SamplerateConversionFeature : public VectorFloatFeatureStream {
 public:
  // ratio : Equal to input_sample_rate / output_sample_rate.
  SamplerateConversionFeature(const VectorFloatFeatureStreamPtr& src, unsigned sourcerate = 22050, unsigned destrate = 16000,
			      unsigned len = 0, const String& method = "fastest", const String& nm = "SamplerateConversion");
  virtual ~SamplerateConversionFeature();

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset();

 private:
  VectorFloatFeatureStreamPtr			src_;
  SRCSTATE_*                                    state_;
  SRC_DATA                                      data_;
  int                                           error_;

  unsigned					dataInSamplesN_;
  unsigned					dataOutStartX_;
  unsigned					dataOutSamplesN_;
};

typedef Inherit<SamplerateConversionFeature, VectorFloatFeatureStreamPtr> SamplerateConversionFeaturePtr;

#endif

/*@}*/

/**
* \defgroup VTLNFeature VTLN Feature
*/
/*@{*/

// ----- definition for class `VTLNFeature' -----
//
// -------------------------------------------
// Piecewise linear: Y = X/Ratio  for X < edge
// -------------------------------------------
class VTLNFeature : public VectorFeatureStream {
 public:
  VTLNFeature(const VectorFeatureStreamPtr& pow,
	      unsigned coeffN = 0, double ratio = 1.0, double edge = 1.0, int version = 1,
	      const String& nm = "VTLN");
  virtual ~VTLNFeature() {  }

  virtual const gsl_vector* next(int frame_no = -5);

  virtual void reset() { pow_->reset(); VectorFeatureStream::reset(); }

  // specify the warp factor
  void warp(double w) { ratio_ = w; }

  void matrix(gsl_matrix* mat) const;

private:
  virtual const gsl_vector* nextFF(int frame_no);
  virtual const gsl_vector* nextOrg(int frame_no);

 private:
  VectorFeatureStreamPtr pow_;
  double                 ratio_;
  const double           edge_;
  const int              version_;
  gsl_vector*            auxV_;
};

typedef Inherit<VTLNFeature, VectorFeatureStreamPtr> VTLNFeaturePtr;

/*@}*/

/**
* \defgroup MelFeature Mel Feature
*/
/*@{*/

// ----- definition for class 'MelFeature' -----
//
class MelFeature : public VectorFeatureStream {
  class SparseMatrix_ {
  public:
    SparseMatrix_(unsigned m, unsigned n, unsigned version);
    ~SparseMatrix_();

    void melScale(int powN,  float rate, float low, float up, int filterN);
    void melScaleOrg(int powN,  float rate, float low, float up, int filterN);
    void melScaleFF(int powN,  float rate, float low, float up, int filterN);
    gsl_vector* fmatrixBMulotOrg( gsl_vector* C, const gsl_vector* A) const;
    gsl_vector* fmatrixBMulotFF( gsl_vector* C, const gsl_vector* A) const;
    void fmatrixBMulot( gsl_vector* C, const gsl_vector* A) const;
    void readBuffer(const String& fb);

    void matrix(gsl_matrix* mat) const;

  private:
    void alloc_(unsigned m, unsigned n);
    void dealloc_();

    float mel_(float hz);
    float hertz_(float m);

    float**					data_;
    unsigned					m_;
    unsigned					n_;
    unsigned*					offset_;// offset
    unsigned*					coefN_; // number of coefficients
    float					rate_;  // sampling rate in Hz
    int version_; // SparseMatrix_ version number, 1:Org, 2:Friedich's changes
  };

 public:
  MelFeature(const VectorFeatureStreamPtr& mag, int powN = 0,
             float rate = 16000.0, float low = 0.0, float up = 0.0,
             unsigned filterN = 30, unsigned version = 1, const String& nm = "MelFFT");

  virtual ~MelFeature();

  virtual const gsl_vector* next(int frame_no = -5);

  virtual void reset() { mag_->reset(); VectorFeatureStream::reset(); }

  void read(const String& fileName);

  void matrix(gsl_matrix* mat) const { mel_.matrix(mat); }

 private:
  // gsl_vector* _fmatrixBMulot(gsl_vector* C, const gsl_vector* A, FBMatrix* B) const;

  const unsigned				nmel_;
  const unsigned				powN_;
  VectorFeatureStreamPtr			mag_;
  SparseMatrix_					mel_;
};

typedef Inherit<MelFeature, VectorFeatureStreamPtr> MelFeaturePtr;

// ----- definition for class 'SphinxMelFeature' -----
//
class SphinxMelFeature : public VectorFeatureStream {
  class Boundary_ {
  public:
    Boundary_(unsigned min_k, unsigned max_k)
      : min_k_(min_k), max_k_(max_k) { }

    unsigned					min_k_;
    unsigned					max_k_;
  };
  typedef vector<Boundary_>			Boundaries_;

 public:
  SphinxMelFeature(const VectorFeatureStreamPtr& mag, unsigned fftN = 512, unsigned powerN = 0,
                   float sampleRate = 16000.0, float lowerF = 0.0, float upperF = 0.0,
                   unsigned filterN = 30, const String& nm = "SphinxMelFilterBank");

  virtual ~SphinxMelFeature();

  virtual const gsl_vector* next(int frame_no = -5);

  virtual void reset() { mag_->reset(); VectorFeatureStream::reset(); }

  void read(const String& fileName);

 private:
  static double melFrequency_(double frequency);
  static double melInverseFrequency_(double frequency);

  const unsigned				fftN_;
  const unsigned				filterN_;
  const unsigned				powN_;
  const double					samplerate_;
  VectorFeatureStreamPtr			mag_;
  gsl_matrix*					filters_;
  Boundaries_					boundaries_;
};

typedef Inherit<SphinxMelFeature, VectorFeatureStreamPtr> SphinxMelFeaturePtr;

/*@}*/

/**
* \defgroup LogFeature Log Feature
*/
/*@{*/

// ----- definition for class `LogFeature' -----
//
class LogFeature : public VectorFloatFeatureStream {
 public:
  LogFeature(const VectorFeatureStreamPtr& mel, double m = 1.0, double a = 1.0,
             bool sphinxFlooring = false, const String& nm = "LogMel");
  virtual ~LogFeature() { }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { mel_->reset(); VectorFloatFeatureStream::reset(); }

 private:
  const unsigned                                nmel_;
  VectorFeatureStreamPtr                        mel_;
  const double                                  m_;
  const double                                  a_;
  const bool                                    SphinxFlooring_;
};

typedef Inherit<LogFeature, VectorFloatFeatureStreamPtr> LogFeaturePtr;

// ----- definition for class `FloatToDoubleConversionFeature' -----
//
class FloatToDoubleConversionFeature : public VectorFeatureStream {
 public:
  FloatToDoubleConversionFeature(const VectorFloatFeatureStreamPtr& src, const String& nm = "FloatToDoubleConversion") : VectorFeatureStream(src->size(), nm), src_(src) {};

  virtual ~FloatToDoubleConversionFeature() { }

  virtual const gsl_vector* next(int frame_no = -5);

  virtual void reset() { src_->reset(); VectorFeatureStream::reset(); }

 private:
  VectorFloatFeatureStreamPtr     src_;

};

typedef Inherit<FloatToDoubleConversionFeature, VectorFeatureStreamPtr> FloatToDoubleConversionFeaturePtr;

/*@}*/

/**
* \defgroup CepstralFeature Cepstral Feature
*/
/*@{*/

// ----- definition for class `CepstralFeature' -----
//
// type:
//   0 =
//   1 = Type 2 DCT
//   2 = Sphinx Legacy
class CepstralFeature : public VectorFloatFeatureStream {
 public:
  CepstralFeature(const VectorFloatFeatureStreamPtr& mel, unsigned ncep = 13,
                  int type = 1, const String& nm = "Cepstral");

  virtual ~CepstralFeature() { gsl_matrix_float_free(cos_); }

  virtual const gsl_vector_float* next(int frame_no = -5);

  gsl_matrix* matrix() const;

  virtual void reset() { mel_->reset(); VectorFloatFeatureStream::reset(); }

 private:
  void sphinxLegacy_();

  gsl_matrix_float*				cos_;
  VectorFloatFeatureStreamPtr			mel_;
};

typedef Inherit<CepstralFeature, VectorFloatFeatureStreamPtr> CepstralFeaturePtr;

/*@}*/

/**
* \defgroup MeanSubtractionFeature Mean Subtraction Feature
*/
/*@{*/

// ----- definition for class `MeanSubtractionFeature' -----
//
class MeanSubtractionFeature : public VectorFloatFeatureStream {
 public:
  MeanSubtractionFeature(const VectorFloatFeatureStreamPtr& src, const VectorFloatFeatureStreamPtr& weight = NULL, double devNormFactor = 0.0, bool runon = false, const String& nm = "Mean Subtraction");

  virtual ~MeanSubtractionFeature();

  virtual const gsl_vector_float* next(int frame_no = -5);
  virtual void reset();
  const gsl_vector_float* mean() const { return mean_; }
  void write(const String& fileName, bool variance = false) const;
  void next_speaker();

#ifdef ENABLE_LEGACY_BTK_API
  void nextSpeaker();
#endif

 private:
  const gsl_vector_float* nextRunon_(int frame_no);
  const gsl_vector_float* nextBatch_(int frame_no);
  void calcMeanVariance_();
  void normalize_(const gsl_vector_float* srcVec);

  static const float				variance_floor_;
  static const float				before_wgt_;
  static const float				after_wgt;
  static const unsigned				framesN2change_;

  VectorFloatFeatureStreamPtr			src_;
  VectorFloatFeatureStreamPtr			wgt_;

  gsl_vector_float*				mean_;
  gsl_vector_float*				var_;
  const double					devNormFactor_;
  unsigned					framesN_;

  bool						runon_;
  bool						mean_var_found_;
};

typedef Inherit<MeanSubtractionFeature, VectorFloatFeatureStreamPtr> MeanSubtractionFeaturePtr;

/*@}*/

/**
* \defgroup FileMeanSubtractionFeature File Mean Subtraction Feature
*/
/*@{*/

// ----- definition for class `FileMeanSubtractionFeature' -----
//
class FileMeanSubtractionFeature : public VectorFloatFeatureStream {
 public:
  FileMeanSubtractionFeature(const VectorFloatFeatureStreamPtr& src,
			     double devNormFactor = 0.0, const String& nm = "File Mean Subtraction");

  virtual ~FileMeanSubtractionFeature();

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset();

  void read(const String& fileName, bool variance = false);

 private:
  static const float				variance_floor_;

  VectorFloatFeatureStreamPtr			src_;
  gsl_vector_float*				mean_;
  gsl_vector_float*				variance_;
  const double					devNormFactor_;
};

typedef Inherit<FileMeanSubtractionFeature, VectorFloatFeatureStreamPtr> FileMeanSubtractionFeaturePtr;

/*@}*/

/**
* \defgroup AdjacentFeature Adjacent Feature
*/
/*@{*/

// ----- definition for class `AdjacentFeature' -----
//
class AdjacentFeature : public VectorFloatFeatureStream {
 public:
  AdjacentFeature(const VectorFloatFeatureStreamPtr& single, unsigned delta = 5,
		  const String& nm = "Adjacent");

  virtual ~AdjacentFeature() { }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset();

 private:
  void buffer_next_frame_(int frame_no);

  const unsigned				delta_;
  VectorFloatFeatureStreamPtr			single_;
  const unsigned				singleSize_;
  const unsigned				plen_;
  unsigned					framesPadded_;
};

typedef Inherit<AdjacentFeature, VectorFloatFeatureStreamPtr> AdjacentFeaturePtr;

/*@}*/

/**
* \defgroup LinearTransformFeature Linear Transform Feature
*/
/*@{*/

// ----- definition for class `LinearTransformFeature' -----
//
class LinearTransformFeature : public VectorFloatFeatureStream {
 public:
#if 0
  LinearTransformFeature(const VectorFloatFeatureStreamPtr& src,
			 gsl_matrix_float* mat = NULL, unsigned sz = 0, const String& nm = "Transform");
#else
  LinearTransformFeature(const VectorFloatFeatureStreamPtr& src, unsigned sz = 0, const String& nm = "Transform");
#endif

  virtual ~LinearTransformFeature() { gsl_matrix_float_free(trans_); }

  virtual const gsl_vector_float* next(int frame_no = -5);

  gsl_matrix_float* matrix() const;

  void load(const String& fileName, bool old = false);

  void identity();

  virtual void reset() { src_->reset(); VectorFloatFeatureStream::reset(); }

 protected:
  VectorFloatFeatureStreamPtr			src_;
  gsl_matrix_float*				trans_;
};

typedef Inherit<LinearTransformFeature, VectorFloatFeatureStreamPtr> LinearTransformFeaturePtr;

/*@}*/

/**
* \defgroup StorageFeature Storage Feature
*/
/*@{*/

// ----- definition for class `StorageFeature' -----
//
class StorageFeature : public VectorFloatFeatureStream {
  typedef vector<gsl_vector_float*>	_StorageVector;
  static const int MaxFrames;
 public:
  StorageFeature(const VectorFloatFeatureStreamPtr& src, const String& nm = "Storage");

  virtual ~StorageFeature();

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { src_->reset(); VectorFloatFeatureStream::reset(); }

  void write(const String& fileName, bool plainText = false) const;
  void read(const String& fileName);
  int evaluate();

 private:
  VectorFloatFeatureStreamPtr			src_;
  _StorageVector				frames_;
};

typedef Inherit<StorageFeature, VectorFloatFeatureStreamPtr> StorageFeaturePtr;

/**
* \defgroup StaticStorageFeature Static Storage Feature
*/
/*@{*/

// ----- definition for class `StaticStorageFeature' -----
//
class StaticStorageFeature : public VectorFloatFeatureStream {
  typedef vector<gsl_vector_float*> _StorageVector;
  static const int MaxFrames;
 public:
  StaticStorageFeature(unsigned dim, const String& nm = "Static Storage");

  virtual ~StaticStorageFeature();

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { VectorFloatFeatureStream::reset(); }

  //void write(const String& fileName) const;
  void read(const String& fileName);
  int evaluate();
  unsigned currentNFrames() const { return frame_no_; };

 private:
  //VectorFloatFeatureStreamPtr     src_;
  _StorageVector        frames_;
  int                  framesN_;
};

typedef Inherit<StaticStorageFeature, VectorFloatFeatureStreamPtr> StaticStorageFeaturePtr;

/*@}*/

/**
* \defgroup CircularStorageFeature Circular Storage Feature
*/
/*@{*/

// ----- definition for class `CircularStorageFeature' -----
//
class CircularStorageFeature : public VectorFloatFeatureStream {
  typedef vector<gsl_vector_float*>	_StorageVector;
  static const int MaxFrames;
 public:
  CircularStorageFeature(const VectorFloatFeatureStreamPtr& src, unsigned framesN = 3, const String& nm = "CircularStorage");

  virtual ~CircularStorageFeature();

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset();

 private:
  unsigned get_index_(int diff) const;

  VectorFloatFeatureStreamPtr			src_;
  const unsigned				framesN_;
  _StorageVector				frames_;
  unsigned					pointerX_;
};

typedef Inherit<CircularStorageFeature, VectorFloatFeatureStreamPtr> CircularStorageFeaturePtr;

/*@}*/

/**
* \defgroup FilterFeature Filter Feature
*/
/*@{*/

// ----- definition for class `FilterFeature' -----
//
class FilterFeature : public VectorFloatFeatureStream {
  class Buffer_ {
  public:
    Buffer_(unsigned len, unsigned nsamp)
      : len_(len), nsamp_(nsamp), offset_(int((nsamp_ - 1) / 2)),
      zero_(0), samples_(new gsl_vector_float*[nsamp_])
    {
      assert (nsamp_ % 2 == 1);
      for (unsigned i = 0; i < nsamp_; i++)
        samples_[i] = gsl_vector_float_calloc(len_);
    }
    ~Buffer_()
    {
      for (unsigned i = 0; i < nsamp_; i++)
        gsl_vector_float_free(samples_[i]);
      delete[] samples_;
    }

    const gsl_vector_float* sample(unsigned timeX) const {
      return samples_[index_(timeX)];
    }

    const double sample(int timeX, unsigned binX) const {
      unsigned idx = index_(timeX);
      const gsl_vector_float* vec = samples_[idx];
      return gsl_vector_float_get(vec, binX);
    }

    void nextSample(const gsl_vector_float* s = NULL) {
      zero_ = (zero_ + 1) % nsamp_;
      gsl_vector_float* nextBlock = samples_[(zero_ + offset_) % nsamp_];

      if (s == NULL) {
        gsl_vector_float_set_zero(nextBlock);
      } else {
        assert( s->size == len_ );
        gsl_vector_float_memcpy(nextBlock, s);
      }
    }

    void zero() {
      for (unsigned i = 0; i < nsamp_; i++)
        gsl_vector_float_set_zero(samples_[i]);
      zero_ = nsamp_ - offset_;
    }

    void print() const {
      for (int i = -offset_; i <= offset_; i++)
        printf("        %4d", i);
      printf("\n     --------------------------------------------------------------------------------\n");
      for (unsigned l = 0; l < len_; l++) {
        for (int i = -offset_; i <= offset_; i++)
          printf("  %10.4f", sample(i, l));
        printf("\n");
      }
    }

  private:
    unsigned index_(int idx) const {
      assert ( abs(idx) <= offset_);
      unsigned ret = (zero_ + nsamp_ + idx) % nsamp_;
      return ret;
    }

    const unsigned				len_;
    const unsigned				nsamp_;
    int						offset_;
    unsigned					zero_; // index of most recent sample
    gsl_vector_float**				samples_;
  };

 public:
  FilterFeature(const VectorFloatFeatureStreamPtr& src, gsl_vector* coeffA, const String& nm = "Filter");
  virtual ~FilterFeature();

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset();

  void printBuffer() const { buffer_.print(); }

 private:
  void buffer_next_frame_(int frame_no);

  VectorFloatFeatureStreamPtr			src_;
  unsigned					lenA_;
  gsl_vector*					coeffA_;
  int						offset_;
  Buffer_					buffer_;
  unsigned					framesPadded_;
};

typedef Inherit<FilterFeature, VectorFloatFeatureStreamPtr> FilterFeaturePtr;

/*@}*/

/**
* \defgroup MergeFeature Merge Feature
*/
/*@{*/

// ----- definition for class `MergeFeature' -----
//
class MergeFeature : public VectorFloatFeatureStream {
  typedef list<VectorFloatFeatureStreamPtr>	FeatureList_;
  typedef FeatureList_::iterator		FeatureListIterator_;
  typedef FeatureList_::const_iterator		FeatureListConstIterator_;
 public:
  MergeFeature(VectorFloatFeatureStreamPtr& stat,
	       VectorFloatFeatureStreamPtr& delta,
	       VectorFloatFeatureStreamPtr& deltaDelta,
	       const String& nm = "Merge");

  virtual ~MergeFeature() { }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset();

 private:
  FeatureList_					flist_;
};

typedef Inherit<MergeFeature, VectorFloatFeatureStreamPtr> MergeFeaturePtr;

/**
* \defgroup MultiModalFeature
*/
/*@{*/

// ----- definition for class `MultiModalFeature MergeFeature' -----
//
class MultiModalFeature : public VectorFloatFeatureStream {
  typedef list<VectorFloatFeatureStreamPtr>	FeatureList_;
  typedef FeatureList_::iterator		FeatureListIterator_;
  typedef FeatureList_::const_iterator		FeatureListConstIterator_;
 public:
  MultiModalFeature(unsigned nModality, unsigned totalVecSize, const String& nm = "Multi");

  ~MultiModalFeature();

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset();

  void addModalFeature( VectorFloatFeatureStreamPtr &feature, unsigned samplePeriodinNanoSec=1 );

 private:
  unsigned     *samplePeriods_;  /* sample period (in nano sec.) */
  unsigned     minSamplePeriod_;
  unsigned     nModality_;
  unsigned     curr_vecsize_;
  FeatureList_ flist_;
};

typedef Inherit<MultiModalFeature, VectorFloatFeatureStreamPtr> MultiModalFeaturePtr;

/*@}*/

/**
* \defgroup FeatureSet Feature Set
*/
/*@{*/

// ----- definition for class `FeatureSet' -----
//
class FeatureSet {
  typedef List <VectorFloatFeatureStreamPtr>	List_;
 public:
  FeatureSet(const String& nm = "FeatureSet") :
    name_(nm), list_(nm) { }

  const String& name() const { return name_; }

  void add(VectorFloatFeatureStreamPtr& feat) { list_.add(feat->name(), feat); }
  VectorFloatFeatureStreamPtr& feature(const String& nm) { return list_[nm]; }

 private:
  const String		name_;
  List_			list_;
};

typedef refcount_ptr<FeatureSet>	FeatureSetPtr;

/*@}*/

void writeGSLMatrix(const String& fileName, const gsl_matrix* mat);

#ifdef JACK
#include <vector>
#include <jack/jack.h>
#include <jack/ringbuffer.h>

typedef struct {
  jack_port_t *port;
  jack_ringbuffer_t *buffer;
  unsigned buffersize;
  unsigned overrun;
  bool can_process;
} jack_channel_t;

/**
* \defgroup Jack Jack Object
*/
/*@{*/

class Jack {
 public:
  Jack(const String& nm);
  ~Jack();
  jack_channel_t* addPort(unsigned buffersize, const String& connection, const String& nm);
  void start(void) { can_capture = true; };
  unsigned getSampleRate() { return (unsigned)jack_get_sample_rate(client); }

 private:
  int process_callback (jack_nframes_t nframes);

  static int _process_callback(jack_nframes_t nframes, void *arg) {
    return static_cast<Jack *> (arg)->process_callback (nframes);
  }

  void shutdown_callback (void);
  static void _shutdown_callback(void *arg) {
    static_cast<Jack *> (arg)->shutdown_callback();
  }

  jack_client_t*				client;
  volatile bool					can_capture;
  volatile bool					can_process;
  vector<jack_channel_t*>			channel;
};

typedef refcount_ptr<Jack>			JackPtr;

/*@}*/

/**
* \defgroup JackFeature Jack Feature
*/
/*@{*/

// ----- definition for class `JackFeature' -----
//
class JackFeature;
typedef Inherit<JackFeature, VectorFloatFeatureStreamPtr> JackFeaturePtr;
class JackFeature : public VectorFloatFeatureStream {
 public:
  JackFeature(JackPtr& jack, unsigned blockLen, unsigned buffersize,
              const String& connection, const String& nm);

  virtual ~JackFeature() { }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { VectorFloatFeatureStream::reset(); }

 private:
  JackPtr					jack_;
  jack_channel_t*				channel;
};

#endif

/*@}*/

/**
* \defgroup ZeroCrossingRateHammingFeature Zero Crossing Rate Hamming Feature
*/
/*@{*/

// ----- definition for class `ZeroCrossingRateHammingFeature' -----
//
class ZeroCrossingRateHammingFeature : public VectorFloatFeatureStream {
 public:
  ZeroCrossingRateHammingFeature(const VectorFloatFeatureStreamPtr& samp, const String& nm = "Zero Crossing Rate Hamming");
  virtual ~ZeroCrossingRateHammingFeature() { delete[] window_; }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { samp_->reset(); VectorFloatFeatureStream::reset(); }

 private:
  VectorFloatFeatureStreamPtr			samp_;
  unsigned					windowLen_;
  double*					window_;
};

typedef Inherit<ZeroCrossingRateHammingFeature, VectorFloatFeatureStreamPtr> ZeroCrossingRateHammingFeaturePtr;

/*@}*/

/**
* \defgroup YINPitchFeature YIN Pitch Feature
*/
/*@{*/

// ----- definition for class `YINPitchFeature' -----
//
class YINPitchFeature : public VectorFloatFeatureStream {
 public:
  YINPitchFeature(const VectorFloatFeatureStreamPtr& samp, unsigned samplerate = 16000, float threshold = 0.5, const String& nm = "YIN Pitch");
  virtual ~YINPitchFeature() { }

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { samp_->reset(); VectorFloatFeatureStream::reset(); }

 private:
  float	_getPitch(const gsl_vector_float *input, gsl_vector_float *yin, float tol);
  VectorFloatFeatureStreamPtr			samp_;
  unsigned sr_;
  float	tr_;
};

typedef Inherit<YINPitchFeature, VectorFloatFeatureStreamPtr> YINPitchFeaturePtr;

/*@}*/

/**
* \defgroup SpikeFilter Spike Filter
*/
/*@{*/

// ----- definition for class `SpikeFilter' -----
//
class SpikeFilter : public VectorFloatFeatureStream {
 public:
  SpikeFilter(VectorFloatFeatureStreamPtr& src, unsigned tapN = 3, const String& nm = "Spike Filter");
  virtual ~SpikeFilter();

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset() { src_->reset(); VectorFloatFeatureStream::reset(); }

private:
  VectorFloatFeatureStreamPtr		src_;
  const unsigned			adcN_;
  const unsigned			queueN_;
  float*				queue_;
  const unsigned			windowN_;
  float*				window_;
};

typedef Inherit<SpikeFilter, VectorFloatFeatureStreamPtr> SpikeFilterPtr;

/*@}*/

/**
* \defgroup SpikeFilter2 Spike Filter 2
*/
/*@{*/

// ----- definition for class `SpikeFilter2' -----
//
class SpikeFilter2 : public VectorFloatFeatureStream {
 public:
  SpikeFilter2(VectorFloatFeatureStreamPtr& src,
	       unsigned width = 3, float maxslope = 7000.0, float startslope = 100.0, float thresh = 15.0, float alpha = 0.2, unsigned verbose = 1,
	       const String& nm = "Spike Filter 2");

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset();

  unsigned spikesN() const { return count_; }

private:
  VectorFloatFeatureStreamPtr		src_;
  const unsigned 			adcN_;
  const unsigned			width_;
  const float				maxslope_;
  const float				startslope_;
  const float				thresh_;
  float 				alpha_;
  float 				beta_;
  float					meanslope_;
  unsigned				count_;
  const unsigned			verbose_;
};

typedef Inherit<SpikeFilter2, VectorFloatFeatureStreamPtr> SpikeFilter2Ptr;

/*@}*/

namespace sndfile {
#include <sndfile.h>

/**
* \defgroup SoundFile Sound File
*/
/*@{*/

// ----- definition for class `SoundFile' -----
//
class SoundFile {
 public:
  SoundFile(const String& fn,
            int mode = sndfile::SFM_RDWR,
            int format = sndfile::SF_FORMAT_WAV|sndfile::SF_FORMAT_PCM_16,
            int samplerate = 16000,
            int channels = 1,
            bool normalize = false);
  ~SoundFile() { sf_close(sndfile_); }
  sf_count_t frames() const { return sfinfo_.frames; }
  int samplerate() const { return sfinfo_.samplerate; }
  int channels() const { return sfinfo_.channels; }
  int format() const { return sfinfo_.format; }
  int sections() const { return sfinfo_.sections; }
  int seekable() const { return sfinfo_.seekable; }
  sf_count_t readf(float *ptr, sf_count_t frames)
    { return sf_readf_float(sndfile_, ptr, frames); }
  sf_count_t writef(float *ptr, sf_count_t frames)
    { return sf_writef_float(sndfile_, ptr, frames); }
  sf_count_t read(float *ptr, sf_count_t items)
    { return sf_read_float(sndfile_, ptr, items); }
  sf_count_t write(float *ptr, sf_count_t items)
    { return sf_write_float(sndfile_, ptr, items); }
  sf_count_t seek(sf_count_t frames, int whence = SEEK_SET)
    { return sf_seek(sndfile_, frames, whence); }
 private:
  SNDFILE* sndfile_;
  SF_INFO sfinfo_;
};
}
typedef refcount_ptr<sndfile::SoundFile>	SoundFilePtr;

/*@}*/

/**
* \defgroup DirectSampleFeature Direct Sample Feature
*/
/*@{*/

// ----- definition for class `DirectSampleFeature' -----
//
class DirectSampleFeature;
typedef Inherit<DirectSampleFeature, VectorFloatFeatureStreamPtr> DirectSampleFeaturePtr;
class DirectSampleFeature : public VectorFloatFeatureStream {
 public:
  DirectSampleFeature(const SoundFilePtr &sndfile,
                      unsigned blockLen = 320,
                      unsigned start = 0,
                      unsigned end = (unsigned)-1,
                      const String& nm = "DirectSample");
  virtual ~DirectSampleFeature() {}
  virtual const gsl_vector_float* next(int frame_no = -5);
  int sampleRate() const { return sndfile_->samplerate(); }
  int channels() const { return sndfile_->channels(); }
  void setRegion(unsigned start = 0, unsigned end = (unsigned)-1) {
    start_ = start;
    end_ = end;
  }
  virtual void reset() {
    sndfile_->seek(start_, SEEK_SET);
    cur_ = 0;
    VectorFloatFeatureStream::reset();
  }
 private:
  SoundFilePtr sndfile_;
  unsigned blockLen_;
  unsigned start_;
  unsigned end_;
  unsigned cur_;
};

/*@}*/

/**
* \defgroup DirectSampleOutputFeature Direct Sample Output Feature
*/
/*@{*/

// ----- definition for class `DirectSampleOutputFeature' -----
//
class DirectSampleOutputFeature;
typedef Inherit<DirectSampleOutputFeature, VectorFloatFeatureStreamPtr> DirectSampleOutputFeaturePtr;
class DirectSampleOutputFeature : public VectorFloatFeatureStream {
 public:
  DirectSampleOutputFeature(const VectorFloatFeatureStreamPtr& src,
			    const SoundFilePtr &sndfile,
			    const String& nm = "DirectSampleOutput");
  virtual ~DirectSampleOutputFeature() {}
  virtual const gsl_vector_float* next(int frame_no = -5);
  virtual void reset() { src_->reset(); sndfile_->seek(0, SEEK_SET); VectorFloatFeatureStream::reset(); }

 private:
  VectorFloatFeatureStreamPtr src_;
  SoundFilePtr sndfile_;
  unsigned blockLen_;
};

/*@}*/

/**
* \defgroup ChannelExtractionFeature Channel Extraction Feature
*/
/*@{*/

// ----- definition for class `ChannelExtractionFeature' -----
//
class ChannelExtractionFeature;
typedef Inherit<ChannelExtractionFeature, VectorFloatFeatureStreamPtr> ChannelExtractionFeaturePtr;
class ChannelExtractionFeature : public VectorFloatFeatureStream {
 public:
  ChannelExtractionFeature(const VectorFloatFeatureStreamPtr& src,
                           unsigned chX = 0,
                           unsigned chN = 1,
                           const String& nm = "ChannelExtraction")
    : VectorFloatFeatureStream(src->size()/chN, nm), src_(src), chX_(chX), chN_(chN)
    {
      assert(chX < chN);
      assert((src->size() % chN) == 0);
    }
  virtual ~ChannelExtractionFeature() {}
  virtual const gsl_vector_float* next(int frame_no = -5);
  virtual void reset() { src_->reset(); VectorFloatFeatureStream::reset(); }
 private:
  VectorFloatFeatureStreamPtr src_;
  unsigned chX_;
  unsigned chN_;
};

/*@}*/

/**
* \defgroup SignalInterferenceFeature Signal Interference Feature
*/
/*@{*/

// ----- definition for class SignalInterferenceFeature -----
//
class SignalInterferenceFeature;
typedef Inherit<SignalInterferenceFeature, VectorFloatFeatureStreamPtr> SignalInterferenceFeaturePtr;
class SignalInterferenceFeature : public VectorFloatFeatureStream{
public:
  SignalInterferenceFeature(VectorFloatFeatureStreamPtr& signal, VectorFloatFeatureStreamPtr& interference, double dBInterference = 0.0, unsigned blockLen = 512, const String& nm = "Signal Interference");

  virtual ~SignalInterferenceFeature() {};
  virtual const gsl_vector_float* next(int frame_no = -5);
  virtual void reset() { signal_->reset(); interference_->reset(); VectorFloatFeatureStream::reset(); }

private:
  VectorFloatFeatureStreamPtr   	signal_;
  VectorFloatFeatureStreamPtr		interference_;
  const double				level_;
};

/*@}*/

/**
* \defgroup AmplificationFeature Amplification Feature
*/
/*@{*/

// ----- definition for class `AmplificationFeature' -----
//
class AmplificationFeature;
typedef Inherit<AmplificationFeature, VectorFloatFeatureStreamPtr> AmplificationFeaturePtr;
class AmplificationFeature : public VectorFloatFeatureStream {
 public:
 AmplificationFeature(const VectorFloatFeatureStreamPtr& src,
                      double amplify = 1.0,
                      const String& nm = "Amplification")
   : VectorFloatFeatureStream(src->size(), nm), src_(src), amplify_(amplify)
  {}
  virtual ~AmplificationFeature() {}
  virtual const gsl_vector_float* next(int frame_no = -5);
  virtual void reset() { src_->reset(); VectorFloatFeatureStream::reset(); }

 private:
  VectorFloatFeatureStreamPtr src_;
  double amplify_;
};

// ----- definition for class `WriteSoundFile' -----
//
class WriteSoundFile {
public:
  WriteSoundFile(const String& fn, int sampleRate, int nChan=1, int format = sndfile::SF_FORMAT_WAV|sndfile::SF_FORMAT_PCM_32);
  ~WriteSoundFile();
  int write( gsl_vector *vector );
  int writeInt( gsl_vector *vector );
  int writeShort( gsl_vector *vector );
  int writeFloat( gsl_vector *vector );

private:
  sndfile::SNDFILE* sndfile_;
  sndfile::SF_INFO sfinfo_;
};

typedef refcount_ptr<WriteSoundFile> WriteSoundFilePtr;

#endif


text:
#ifndef ControlByRow_h
#define ControlByRow_h
/**
 * @file
 * $Revision: 1.2 $
 * $Date: 2007/01/30 22:12:21 $
 *
 *   Unless noted otherwise, the portions of Isis written by the USGS are
 *   public domain. See individual third-party library and package descriptions
 *   for intellectual property information, user agreements, and related
 *   information.
 *
 *   Although Isis has been used by the USGS, no warranty, expressed or
 *   implied, is made by the USGS as to the accuracy and functioning of such
 *   software and related material nor shall the fact of distribution
 *   constitute any such warranty, and no responsibility is assumed by the
 *   USGS in connection therewith.
 *
 *   For additional information, launch
 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html
 *   in a browser or see the Privacy & Disclaimers page on the Isis website,
 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on
 *   http://www.usgs.gov/privacy.html.
 */

#include <vector>
#include "ControlMeasure.h"
#include "ControlMeasureLogData.h"
#include "Statistics.h"
#include "CollectorMap.h"
#include "IString.h"
#include <gsl/gsl_math.h>

namespace Isis {

  /**
   * @brief Container for point collection
   */
  struct PointData {
    ControlMeasure refPoint;
    ControlMeasure chpPoint;
  };

  /**
   * @brief Less than test for Control point group
   *
   * This function tests the reference line numbers and returns true if the first
   * point reference line is less than the second.
   *
   * @param p1  First PointData set to compare
   * @param p2  Second PointData set to compare
   *
   * @return bool If first point reference line is less than the second
   */
  inline bool PointLess(const PointData &p1, const PointData &p2) {
    return (p1.refPoint.GetLine() < p2.refPoint.GetLine());
  }

  /**
   * @brief Equality test for Control point group
   *
   * This function tests the reference line numbers for equality and returns true
   * if the line references are equivalent, according to an approximation using an
   * epsilon of 1.0e-6.
   *
   * @param p1  First PointData set to compare
   * @param p2  Second PointData set to compare
   *
   * @return bool If the reference point lines are (approximately) equivalent
   */
  inline bool PointEqual(const PointData &p1, const PointData &p2) {
    return (gsl_fcmp(p1.refPoint.GetLine(), p2.refPoint.GetLine(), 1.0E-6) == 0);
  }

  /**
   * @brief Collector for Control points within the same row for analysis
   *
   * This class is designed to be used as a Functor object collecting control net
   * file and collapsing all column measures into on row.  This is primarily used
   * for analysis of coregistration results with one or more columns specified
   * in the search/pattern chip strategy.
   *
   * @author ????-??-?? Unknown
   *
   * @internal
   */
  class ControlByRow {
    public:
      /**
       * @brief Structure to return control point statistics for row
       *
       * This structure contains the row statistics of merged control points.
       * This will eventually be used to compute the spline interpolations for
       * line and sample offsets.
       *
       */
      struct RowPoint {
        RowPoint() : refLine(0.0), refSamp(0.0), chpLine(0.0), chpSamp(0.0),
          total(0), count(0) { }
        double refLine;    //!<  Reference line (row)
        double refSamp;    //!<  Reference sample
        double chpLine;    //!<  Registered line
        double chpSamp;    //!<  Registered sample
        int total;         //!<  Total points in row
        int count;         //!<  Valid points found

        Statistics rSStats;
        Statistics cLStats;
        Statistics cSStats;
        Statistics cLOffset;
        Statistics cSOffset;
        Statistics GOFStats;
      };

    public:
      /**
       * @brief Default constructor
       */
      ControlByRow()  {
        _minGOF = DBL_MIN;
        _maxGOF = DBL_MAX;
      }

      /**
       * @brief Constructor that sets the maximum goodness of fit tolerance
       * @param maxGOF Value that specifies the maximum goodness of fit which is
       *               typically never expected to exceed 1.0 for a good fit.
       */
      ControlByRow(double maxGOF)  {
        _minGOF = DBL_MIN;
        _maxGOF = maxGOF;
      }

      /**
       * @brief Constructor that sets the maximum goodness of fit tolerance
       * @param minGOF Value of minimum goodness of fit.  Allows user selectable
       *               adjustment to coregistration minimum tolerance
       * @param maxGOF Value that specifies the maximum goodness of fit which is
       *               typically never expected to exceed 1.0 for a good fit.
       */
      ControlByRow(double minGOF, double maxGOF)  {
        _minGOF = minGOF;
        _maxGOF = maxGOF;
      }

      /**
       * @brief Destructor
       */
      ~ControlByRow() { }

      /**
       * @brief Determines the number of points (rows) found valid
       *
       * The number returned is really the number of unique rows of coregistration
       * chips gleened from the control net.
       *
       * @return unsigned int Row/line count
       */
      inline unsigned int size() const {
        return (_rowList.size());
      }

      /**
       * @brief Set the minimum acceptable goodness of fit value
       *
       * This sets the minimum (absolute) value used to gleen valid points from
       * the control net data.
       *
       * @param minGOF Minimum goodness of fit tolerance
       */
      void setMinGOF(double minGOF) {
        _minGOF = minGOF;
      }

      /**
       * @brief Set the maximum acceptable goodness of fit value
       *
       * This sets the maximum (absolute) value used to gleen valid points from
       * the control net data.  This is intended to use to exclude wild points
       * that exceed the level of reasonable tolerance.  This is typically 1.0 for
       * most coregistration algorithms.
       *
       * @param maxGOF Maximum goodness of fit tolerance
       */
      void setMaxGOF(double maxGOF) {
        _maxGOF = maxGOF;
      }

      /**
       * @brief Operator used to add a point to the data set
       *
       * This method provides support for the STL for_each algorithm.  This allows
       * this class to be used as a functor object.
       *
       * @param p Point to tested for acceptance in the list
       */
      void operator()(const PointData &p) {
        addPoint(p);
        return;
      }

      /**
       * @brief Formal method of adding a control point to the data set
       *
       * This method will add the provided point to the collection of rows (or
       * lines of points).
       *
       * @param p  Point to add to the list
       */
      void addPoint(const PointData &p) {
        if(_rowList.exists(p.refPoint.GetLine())) {
          PointList &r = _rowList.get(p.refPoint.GetLine());
          r.push_back(p);
        }
        else {
          PointList pl;
          pl.push_back(p);
          _rowList.add(p.refPoint.GetLine(), pl);
        }
        return;
      }

      /**
       * @brief Returns a point at the ith location
       *
       * Traverses the list of points after computing the merge statistics for the
       * ith row.
       *
       * @param i Index of point to return
       *
       * @return RowPoint  Structure containing the (statistically) merged row
       */
      RowPoint operator[](int i) const {
        try {
          return (computeStats(_rowList.getNth(i)));
        }
        catch(IException &oor) {
          QString msg = "Requested value (" + toString(i) +  ") not found";
          throw IException(oor, IException::User, msg, _FILEINFO_);
        }
      }

    private:
      typedef std::vector<PointData> PointList;    //!< Composite list
      typedef CollectorMap<double, PointList, RobustFloatCompare> CNetRow;  /**
                                  * Nifty templated collector class works just
                                  * nicely for merging rows of correlations
                                                                         */
      double            _minGOF;  //!<  Minimum acceptable goodness of fit
      double            _maxGOF;  //!<  Maximum acceptable goodness of fit
      CNetRow           _rowList; //!<  Collection of merged rows/lines

      /**
       * @brief Convenience method for adding point to Statistics class
       *
       * The interface to the Isis Statistics class requires the value to be added
       * by address.  This method expedites the addition of values from say a
       * function or method.
       *
       * @param value  Value to add to Statistics class
       * @param stats  Statistics class to add the data point to
       */
      inline void addToStats(double value, Statistics &stats) const {
        stats.AddData(&value, 1);
        return;
      }

      /**
       * @brief  All important method that computes statistics for a row
       *
       * This method computes the statistics for a potentially merged row of
       * coregistration chips.  It applies the minimum and maximum goodness of fit
       * tolerance checks and add valid points to each statistical component of
       * the merge.
       *
       * @param cols  List of column chip registrations
       *
       * @return RowPoint Structure containing the statistics for row/line of
       *         merged registration columns
       */
      RowPoint computeStats(const std::vector<PointData> &cols) const {
        RowPoint rp;
        rp.refLine = cols[0].refPoint.GetLine();
        for(unsigned int i = 0; i < cols.size() ; i++) {
          double regGOF =
            cols[i].chpPoint.GetLogData(ControlMeasureLogData::GoodnessOfFit)
              .GetNumericalValue();
          if(fabs(regGOF) > _maxGOF) continue;
          if(fabs(regGOF) < _minGOF) continue;
          (rp.count)++;
          addToStats(cols[i].refPoint.GetSample(), rp.rSStats);
          addToStats(cols[i].chpPoint.GetLine(), rp.cLStats);
          addToStats(cols[i].chpPoint.GetSample(), rp.cSStats);
          addToStats(cols[i].chpPoint.GetLineResidual(), rp.cLOffset);
          addToStats(cols[i].chpPoint.GetSampleResidual(), rp.cSOffset);
          addToStats(regGOF, rp.GOFStats);
        }

        rp.total = cols.size();
        rp.refSamp = rp.rSStats.Average();
        rp.chpLine = rp.cLStats.Average();
        rp.chpSamp = rp.cSStats.Average();
        return (rp);
      }

  };

}  // namespace Isis

#endif