matrix.c

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

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  double NxN matrix to be checked
  @param  matrix dimension N (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 */