matrix.c

General Hidden Markov Model Library 一个通用的隐马尔科夫模型的C代码库

/*******************************************************************************
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version 0.8_beta1, 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: [e-mail]
*
*       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: 1650 $
*                       from $Date: 2006-07-20 17:57:27 +0200 (Thu, 20 Jul 2006) $
*             last change by $Author: grunau $.
*
*******************************************************************************/

#ifdef WIN32
#  include "win_config.h"
#endif

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

#include <math.h>
#include <float.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;