ghmm-0.9-rc3/ghmm/matrix.c

/*******************************************************************************
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version __VERSION__, see http://ghmm.org
*
*       Filename: ghmm/ghmm/matrix.c
*       Authors:  Bernhard Knab, Benjamin Georgi
*
*       Copyright (C) 1998-2004 Alexander Schliep
*       Copyright (C) 1998-2001 ZAIK/ZPR, Universitaet zu Koeln
*       Copyright (C) 2002-2004 Max-Planck-Institut fuer Molekulare Genetik,
*                               Berlin
*
*       Contact: schliep@ghmm.org
*
*       This library is free software; you can redistribute it and/or
*       modify it under the terms of the GNU Library General Public
*       License as published by the Free Software Foundation; either
*       version 2 of the License, or (at your option) any later version.
*
*       This library 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
*       Library General Public License for more details.
*
*       You should have received a copy of the GNU Library General Public
*       License along with this library; if not, write to the Free
*       Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*
*       This file is version $Revision: 2267 $
*                       from $Date: 2009-04-24 11:01:58 -0400 (Fri, 24 Apr 2009) $
*             last change by $Author: grunau $.
*
*******************************************************************************/

#ifdef HAVE_CONFIG_H
#  include "../config.h"
#endif

#include <math.h>
#include <float.h>

#include "ghmmconfig.h"

#include "matrix.h"
#include "vector.h"
#include "model.h"
#include "rng.h"
#include "randvar.h"
#include "mes.h"
#include "ghmm_internals.h"

#include "obsolete.h"

#ifdef GHMM_OBSOLETE
static void lrdecomp (int dim, double **a, double *p);
static void lyequalsb (double **a, double *b, double *p, int dim, double *y);
static void ltranspxequalsy (double **a, double *y, double *p, int dim,
                             double *x);

/*============================================================================*/

int ighmm_cmatrix_read (scanner_t * s, double **matrix, int max_zeile,
                   int max_spalte)
{
#define CUR_PROC "ighmm_cmatrix_read"
  int len = 0, zeile = 0;
  ighmm_scanner_consume (s, '{');
  if (s->err)
    return (-1);
  while (!s->eof && !s->err && s->c - '}') {
    if (zeile >= max_zeile) {
      ighmm_scanner_error (s, "too many rows in matrix");
      return (-1);
    }
    /* Memory allocation takes place in scanner_get_double_earray */
    matrix[zeile] = scanner_get_double_earray (s, &len);
    if (len != max_spalte) {
      ighmm_scanner_error (s, "wrong number of elements in matrix");
      return (-1);
    }
    ighmm_scanner_consume (s, ';');
    if (s->err) {
      ighmm_scanner_error (s, "missing ';' or wrong number of columns");
      return (-1);
    }
    zeile++;
  }
  ighmm_scanner_consume (s, '}');
  if (s->err)
    return (-1);
  if (zeile < max_zeile) {
    ighmm_scanner_error (s, "rows missing in matrix");
    return (-1);
  }
  return (0);
#undef CUR_PROC
}                               /* ighmm_cmatrix_read */

/*============================================================================*/
#endif /* GHHM_OBSOLETE */

#ifdef GHMM_UNSUPPORTED
int ighmm_dmatrix_read (scanner_t * s, int **matrix, int max_zeile, int max_spalte)
{
#define CUR_PROC "ighmm_dmatrix_read"
  int len = 0, zeile = 0;
  ighmm_scanner_consume (s, '{');
  if (s->err)
    return (-1);
  while (!s->eof && !s->err && s->c - '}') {
    if (zeile >= max_zeile) {
      ighmm_scanner_error (s, "too many rows in matrix");
      return (-1);
    }
    /* Memory allocation takes place in scanner_get_int_array */
    matrix[zeile] = scanner_get_int_array (s, &len);
    if (len != max_spalte) {
      ighmm_scanner_error (s, "wrong number of elements in matrix");
      return (-1);
    }
    ighmm_scanner_consume (s, ';');
    if (s->err) {
      ighmm_scanner_error (s, "missing ';' or wrong number of columns");
      return (-1);
    }
    zeile++;
  }
  ighmm_scanner_consume (s, '}');
  if (s->err)
    return (-1);
  if (zeile < max_zeile) {
    ighmm_scanner_error (s, "rows missing in matrix");
    return (-1);
  }
  return (0);
#undef CUR_PROC
}                               /* ighmm_dmatrix_read */
#endif /* GHMM_UNSUPPORTED */

/*============================================================================*/

/* allocation of matrices with fixed dimensions  */
double **ighmm_cmatrix_stat_alloc (int n, int m)
{
#define CUR_PROC "ighmm_cmatrix_stat_alloc"
  int i;
  double **A;
  double *tmp;

  /* if (!(A = ighmm_calloc (n * sizeof (*A) + n * m * sizeof (**A)))) { */


  if (!(A = ighmm_calloc (n * sizeof (double*) + n * m * sizeof (double)))) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }

  tmp = (double *) (A + n);
  for (i = 0; i < n; i++) {
    A[i] = tmp;
    tmp += m;
  }
  return A;
STOP:
  ighmm_cmatrix_stat_free (&A);
  return NULL;
#undef CUR_PROC
}


int ighmm_cmatrix_stat_free (double ***matrix)
{
#define CUR_PROC "ighmm_cmatrix_stat_free"
  mes_check_ptr (matrix, return (-1));
  if (!*matrix)
    return (0);
  free (*matrix);
  return 0;
#undef CUR_PROC
}

/*============================================================================*/

/* allocation of matrices with fixed dimensions  */
int **ighmm_dmatrix_stat_alloc(int n, int m)
{
#define CUR_PROC "ighmm_dmatrix_stat_alloc"
  int i;
  int **A;
  int *tmp;
  if (!(A = ighmm_calloc (n * sizeof(int*) + n * m * sizeof (int)))) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }

  tmp = (int *) (A + n);
  for (i = 0; i < n; i++) {
    A[i] = tmp;
    tmp += m;
  }
  return A;
STOP:
  ighmm_dmatrix_stat_free (&A);
  return NULL;
#undef CUR_PROC
}


int ighmm_dmatrix_stat_free (int ***matrix)
{
#define CUR_PROC "ighmm_dmatrix_stat_free"
  mes_check_ptr (matrix, return (-1));
  if (!*matrix)
    return (0);
  free (*matrix);
  return 0;
#undef CUR_PROC
}


/*============================================================================*/


double **ighmm_cmatrix_alloc (int zeilen, int spalten)
{
#define CUR_PROC "ighmm_cmatrix_alloc"
  double **matrix;
  int i;

  /*printf("*** ighmm_cmatrix_alloc %d zeilen, %d spalten:\n",zeilen, spalten);*/

  ARRAY_CALLOC (matrix, zeilen);
  for (i = 0; i < zeilen; i++)
    ARRAY_CALLOC (matrix[i], spalten);
  return matrix;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ighmm_cmatrix_free (&matrix, zeilen);
  return NULL;
#undef CUR_PROC
}                               /* ighmm_cmatrix_alloc */


double *** ighmm_cmatrix_3d_alloc(int i, int j, int k) {
#define CUR_PROC "ighmm_cmatrix_3d_alloc"
  double *** matrix;
  int a, b;

  /* printf("*** ighmm_cmatrix_alloc %d zeilen, %d spalten:\n",zeilen, spalten); */

  ARRAY_CALLOC (matrix, i);
  for (a = 0; a < i; a++) {
    ARRAY_CALLOC (matrix[a], j);
    for (b=0; b<j; b++) {
      ARRAY_CALLOC (matrix[a][b], k);
    }
  }
  return matrix;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ighmm_cmatrix_3d_free(&matrix, i, j);
  return NULL;
#undef CUR_PROC
} /* ighmm_cmatrix_3d_alloc */

/** gets a pointer on a 3d matrix and rows and cols **/
int ighmm_cmatrix_3d_free(double **** matrix, int i, int j) {
# define CUR_PROC "ighmm_cmatrix_3d_free"
  int a,b;
  mes_check_ptr(matrix, return(-1));
  if ( !*matrix) return(0);
  for (a = i - 1; a >=  0; a--) {
    for (b=j-1; b>=0; b--)
      m_free((*matrix)[a][b]);
    m_free((*matrix)[a]);
  }
  m_free(*matrix);
  return (0);
# undef CUR_PROC
} /* ighmm_cmatrix_3d_free */

int ighmm_cmatrix_free(double ***matrix, long rows) {
# define CUR_PROC "ighmm_cmatrix_free"
  long i;
  mes_check_ptr(matrix, return (-1));

  if (*matrix) {
    for (i=0; i<rows; i++)
      m_free((*matrix)[i]);
    m_free(*matrix);
  }
  return (0);
# undef CUR_PROC
}                               /* ighmm_cmatrix_free */


/*============================================================================*/


double **ighmm_cmatrix_alloc_copy (int zeilen, int spalten, double **copymatrix)
{
#define CUR_PROC "ighmm_cmatrix_alloc_copy"
  double **matrix;
  int i, j;
  ARRAY_CALLOC (matrix, zeilen);
  for (i = 0; i < zeilen; i++) {
    ARRAY_CALLOC (matrix[i], spalten);
    for (j = 0; j < spalten; j++)
      matrix[i][j] = copymatrix[i][j];
  }
  return matrix;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ighmm_cmatrix_free (&matrix, zeilen);
  return NULL;
#undef CUR_PROC
}                               /* ighmm_cmatrix_alloc_copy */

/*============================================================================*/

int **ighmm_dmatrix_alloc (int zeilen, int spalten)
{
#define CUR_PROC "ighmm_dmatrix_alloc"
  int **matrix;
  int i;
  ARRAY_CALLOC (matrix, zeilen);
  for (i = 0; i < zeilen; i++)
    ARRAY_CALLOC (matrix[i], spalten);
  return matrix;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ighmm_dmatrix_free (&matrix, zeilen);
  return NULL;
#undef CUR_PROC
}                               /* ighmm_dmatrix_alloc */

/*============================================================================*/

int ighmm_dmatrix_free (int ***matrix, long rows)
{
# define CUR_PROC "ighmm_dmatrix_free"
  long i;
  mes_check_ptr (matrix, return (-1));
  if (*matrix) {
    for (i=0; i<rows; i++)
      m_free((*matrix)[i]);
    m_free(*matrix);
  }
  return (0);
# undef CUR_PROC
}                               /* ighmm_dmatrix_free */

/*============================================================================*/
/*============================================================================*/

int*** ighmm_dmatrix_3d_alloc(int zeilen, int spalten, int hoehe) {
#define CUR_PROC "ighmm_dmatrix_alloc"
  int ***matrix;
  int i, j;
  ARRAY_CALLOC (matrix, zeilen);
  for (i = 0; i < zeilen; i++) {
    ARRAY_CALLOC (matrix[i], spalten);
    for (j=0; j < spalten; j++)
      ARRAY_CALLOC (matrix[i][j], hoehe);
  }
  return matrix;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ighmm_dmatrix_3d_free(&matrix, zeilen, spalten);
  return NULL;
#undef CUR_PROC
} /* ighmm_dmatrix_3d_alloc */

/*============================================================================*/

int ighmm_dmatrix_3d_free(int **** matrix, int zeilen, int spalten) {
# define CUR_PROC "ighmm_dmatrix_free"
  long i, j;
  mes_check_ptr(matrix, return(-1));
  if ( !*matrix ) return(0);
  for (i = zeilen - 1; i >= 0; i--) {
    for (j=spalten - 1; j >= 0; j--)
      m_free((*matrix)[i][j]);
    m_free((*matrix)[i]);
  }
  m_free(*matrix);
  return (0);
# undef CUR_PROC
} /* ighmm_dmatrix_3d_free */


#ifdef GHMM_OBSOLETE
/*============================================================================*/
void ighmm_cmatrix_print(FILE *file, double **matrix, int zeilen, int spalten,
		    char *tab, char *separator, char *ending) {
  int i;
  for (i = 0; i < zeilen; i++)
    ighmm_cvector_print (file, matrix[i], spalten, tab, separator, ending);
}                               /* ighmm_cmatrix_print */
#endif /* GHMM_OBSOLETE */

#ifdef GHMM_UNSUPPORTED
/*============================================================================*/
void ighmm_cmatrix_print_prec (FILE * file, double **matrix, int zeilen,
                          int spalten, int width, int prec, char *tab,
                          char *separator, char *ending)
{
  int i;
  for (i = 0; i < zeilen; i++)
    ighmm_cvector_print_prec (file, matrix[i], spalten, width, prec,
                         tab, separator, ending);
}                               /* ighmm_cmatrix_print_prec */

/*============================================================================*/
void ighmm_dmatrix_print (FILE * file, int **matrix, int zeilen, int spalten,
                     char *tab, char *separator, char *ending)
{
  int i;
  for (i = 0; i < zeilen; i++)
    ighmm_dvector_print (file, matrix[i], spalten, tab, separator, ending);
}                               /* ighmm_dmatrix_print */
#endif /* GHMM_UNSUPPORTED */

#ifdef GHMM_OBSOLETE
/*============================================================================*/
int ighmm_cmatrix_notzero_columns (double **matrix, int row, int max_col)
{
  int i, count = 0;
  for (i = 0; i < max_col; i++)
    if (matrix[row][i])
      count++;
  return count;
}                               /* ighmm_cmatrix_notzero_columns */

/*============================================================================*/
int ighmm_cmatrix_notzero_rows (double **matrix, int col, int max_row)
{
  int i, count = 0;
  for (i = 0; i < max_row; i++)
    if (matrix[i][col])
      count++;
  return count;
}                               /* matrix_d__notzero_rows */

/*============================================================================*/
/* Scales the row vectors of a matrix to have the sum 1 */
int ighmm_cmatrix_normalize (double **matrix, int rows, int cols)
{
#define CUR_PROC "ighmm_cmatrix_normalize"
  int i;
  for (i = 0; i < rows; i++)
    if (ighmm_cvector_normalize (matrix[i], cols) == -1)
      ighmm_mes (MES_WIN, "WARNING: sum row[%d] == 0!\n", i);
  /* return (-1); */
  return 0;
#undef CUR_PROC
}                               /* ighmm_cmatrix_normalize */

/*============================================================================*/
void ighmm_cmatrix_random_values (double **matrix, int rows, int cols,
                             double min, double max)
{
  int i, j;
  double interval;
  if (max < min) {
    min = 0.0;
    max = 1.0;
  }
  interval = max - min;
  for (i = 0; i < rows; i++)
    for (j = 0; j < cols; j++)
      matrix[i][j] = min + GHMM_RNG_UNIFORM (RNG) * interval;
}                               /* ighmm_cmatrix_random_values */

/*============================================================================*/
/* Fixed value for final state */
void ighmm_cmatrix_random_const_values (double **matrix, int rows, int cols,
                                   double min, double max, double c)
{
  int i, j;
  double interval;
  if (rows < 1) {
    ighmm_mes (MES_WIN, "WARNING: rows = %d not allowed\n", rows);
    return;
  }
  if (max < min) {
    min = 0.0;
    max = 1.0;
  }
  interval = max - min;
  for (i = 0; i < rows - 1; i++)
    for (j = 0; j < cols; j++)
      matrix[i][j] = min + GHMM_RNG_UNIFORM (RNG) * interval;
  for (j = 0; j < cols; j++)
    matrix[rows - 1][j] = c;
}                               /* ighmm_cmatrix_random_const_values */


/*============================================================================*/
void ighmm_cmatrix_const_values (double **matrix, int rows, int cols, double c)
{
  int i, j;
  for (i = 0; i < rows; i++)
    for (j = 0; j < cols; j++)
      matrix[i][j] = c;
}                               /* ighmm_cmatrix_const_values */

/*============================================================================*/
void ighmm_cmatrix_random_left_right (double **matrix, int rows, int cols)
{
  int i, j;
  for (i = 0; i < rows; i++)
    for (j = 0; j < cols; j++)
      if (j == i || j == i + 1)
        matrix[i][j] = GHMM_RNG_UNIFORM (RNG);
      else
        matrix[i][j] = 0.0;
}                               /* ighmm_cmatrix_random_values */

/*============================================================================*/
void ighmm_cmatrix_left_right_strict (double **matrix, int rows, int cols)
{
  int i, j;
  for (i = 0; i < rows; i++)
    for (j = 0; j < cols; j++)
      if (j == i + 1)
        matrix[i][j] = 1.0;
      else
        matrix[i][j] = 0.0;
}                               /* ighmm_cmatrix_left_right_strict */

/*============================================================================*/
int ighmm_cmatrix_gaussrows_values (double **matrix, int rows, int cols,
                               double **mue, double u)
{
# define CUR_PROC "matrix_gaussrows_values"
  int res = -1;
  double *mean;
  int i, j;
  if (u <= 0.0) {
    GHMM_LOG(LCONVERTED, "sigma^2 <= 0.0 not allowed\n");
    goto STOP;
  }
  if (*mue == NULL) {
    /* for each row, a random mean value mean[i] in (0, cols-1) */
    ARRAY_CALLOC (mean, rows);
    for (i = 0; i < rows; i++)
      mean[i] = GHMM_RNG_UNIFORM (RNG) * (cols - 1);
    /* for (i = 0; i < rows; i++) printf("%6.4f ", mean[i]); printf("\n"); */
    *mue = mean;
  }
  else
    /* Check, if the mean value is on the right interval */
    mean = *mue;
  for (i = 0; i < rows; i++) {
    /* Gauss-distribution around the mean value for each state. */
    for (j = 0; j < cols; j++) {
      matrix[i][j] = ighmm_rand_normal_density ((double) j, mean[i], u);
      if (matrix[i][j] == -1) {
        GHMM_LOG_QUEUED(LCONVERTED);
        goto STOP;
      }
      /* To avoid zero: (Cheap version!!) */
      if (matrix[i][j] < 0.0001)
        matrix[i][j] = 0.0001;  /* The Output has only 4 significant digits! */
    }
  }
  res = 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return (res);
# undef CUR_PROC
}                               /* matrix_gaussrows_values */

/*============================================================================*/
void ighmm_cmatrix_const_preserve_struct (double **matrix, int rows, int cols,
                                     double c)
{
  int i, j;
  for (i = 0; i < rows; i++)
    for (j = 0; j < cols; j++) {
      if (matrix[i][j] != 0)
        matrix[i][j] = c;
    }
}                               /* ighmm_cmatrix_const_preserve_struct */

/*============================================================================*/
void ighmm_cmatrix_random_preserve_struct (double **matrix, int rows, int cols)
{
  int i, j;
  for (i = 0; i < rows; i++)
    for (j = 0; j < cols; j++) {
      if (matrix[i][j] != 0)
        matrix[i][j] = GHMM_RNG_UNIFORM (RNG);
    }
}                               /* ighmm_cmatrix_random_preserve_struct */

/*============================================================================*/
void ighmm_cmatrix_transpose (double **A, int rows, int cols, double **A_T)
{
  int i, j;
  for (i = 0; i < rows; i++)
    for (j = 0; j < cols; j++)
      A_T[j][i] = A[i][j];
}                               /* ighmm_cmatrix_transpose */


/*----------------------------------------------------------------------------*/
/* Decomposes a positiv definit, symmetric matrix A in L*L-T, that is except for
   the diagonal, L (= lower triangular marix) is stored in A. The Diagonal is
   stored in a vector p.
*/
static void lrdecomp (int dim, double **a, double *p)
{
  int k, i, j;
  double x;
  for (i = 0; i < dim; i++) {
    for (j = i; j < dim; j++) {
      x = a[i][j];
      for (k = i - 1; k >= 0; k--)
        x = x - a[j][k] * a[i][k];
      if (i == j) {
        if (x < DBL_MIN)
          ighmm_mes (MES_WIN, "FEHLER: Pivotel.<=0!");
        p[i] = 1 / sqrt (x);
      }
      else
        a[j][i] = x * p[i];
    }
  }
}

/*----------------------------------------------------------------------------*/
/* Solves L*y=b, L is a lower triangular matrix stored in A, p=1/diagonal elements. */
static void lyequalsb (double **a, double *b, double *p, int dim, double *y)
{
  int k, j;
  for (k = 0; k < dim; k++) {
    y[k] = b[k];
    for (j = 0; j < k; j++)
      y[k] = y[k] - a[k][j] * y[j];
    y[k] = y[k] * p[k];
  }
}

/*----------------------------------------------------------------------------*/
/* Solves L-T*x=y, L-T an upper triangular matrix, BUT saved in A as L! */
static void ltranspxequalsy (double **a, double *y, double *p, int dim,
                             double *x)
{
  int k, j;
  for (k = dim - 1; k >= 0; k--) {
    x[k] = y[k];
    for (j = k + 1; j < dim; j++)
      x[k] = x[k] - a[j][k] * x[j];
    x[k] = x[k] * p[k];
  }
}


/*============================================================================*/
/* Solves a linear equation system for a symmetric, positiv definite matrix. */
int ighmm_cmatrix_cholesky (double **a, double *b, int dim, double *x)
{
#define CUR_PROC "ighmm_cmatrix_cholesky"
  int res = -1;
  double *p, *y;
  ARRAY_CALLOC (p, dim);
  ARRAY_CALLOC (y, dim);
  lrdecomp (dim, a, p);
  lyequalsb (a, b, p, dim, y);
  ltranspxequalsy (a, y, p, dim, x);
  res = 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return (res);
#undef CUR_PROC
}

/* Finds the determinant of a symetric, positiv definit matrix. */
int ighmm_cmatrix_det_symposdef (double **a, int dim, double *det)
{
#define CUR_PROC "ighmm_cmatrix_det_symposdef"
  int res = -1;
  int i;
  double *p, r;
  ARRAY_CALLOC (p, dim);
  lrdecomp (dim, a, p);
  *det = 1.0;
  for (i = 0; i < dim; i++) {
    r = 1 / p[i];
    *det *= (r * r);
  }
  res = 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return (res);
#undef CUR_PROC
}

/*============================================================================*/
/* Copies a matrix, the memory allocations must to be done outside! */
void ighmm_cmatrix_copy (double **src, double **target, int rows, int cols)
{
  int i, j;

  for (i = 0; i < rows; i++)
    for (j = 0; j < cols; j++)
      target[i][j] = src[i][j];
}

/*============================================================================*/
/**
  Checks whether a quadratic double matrix is stochastic
  @return 0/1 flag for true/false
  @param  matrix double NxN matrix to be checked
  @param  N      number of rows/columns (matrix must be quadaratic)
  */
int ighmm_cmatrix_check_stochasticity (double **matrix, int N)
{
  int i, j;
  double row_sum;
  int stochastic = 1;

  for (i = 0; i < N; i++) {
    row_sum = 0.0;
    for (j = 0; j < N; j++) {
      row_sum = row_sum + matrix[i][j];
    }
    if (row_sum != 1.0) {
      stochastic = 0;
      break;
    }
  }
  return (stochastic);
}
#endif /* GHMM_OBSOLETE */