matrix.c

一个通用的隐性马尔可夫C代码库 开发环境:C语言 简要说明:这是一个通用的隐性马尔可夫C代码库

/*============================================================================*/
int matrix_d_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;
} /* matrix_d_notzero_columns */

/*============================================================================*/
int matrix_d_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 matrix_d_normalize(double **matrix, int rows, int cols) {
#define CUR_PROC "matrix_d_normalize"
  int i;
  for (i = 0; i < rows; i++)
    if (vector_normalize(matrix[i], cols) == -1)
      mes(MES_WIN, "WARNING: sum row[%d] == 0!\n", i);
      /* return (-1); */
  return 0;
#undef CUR_PROC
} /* matrix_d_normalize */

/*============================================================================*/
void matrix_d_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 + gsl_rng_uniform(RNG)*interval;
} /* matrix_d_random_values */

/*============================================================================*/
/* Fixed value for final state */
void matrix_d_random_const_values(double **matrix, int rows, int cols,
				  double min, double max, double c) {
  int i, j;
  double interval;
  if (rows < 1) {
    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 + gsl_rng_uniform(RNG)*interval;
  for (j = 0; j < cols; j++)
    matrix[rows - 1][j] = c;
} /* matrix_d_random_const_values */


/*============================================================================*/
void matrix_d_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;
} /* matrix_d_const_values */

/*============================================================================*/
void matrix_d_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] = gsl_rng_uniform(RNG);
      else
	matrix[i][j] = 0.0;
} /* matrix_d_random_values */

/*============================================================================*/
void matrix_d_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;
} /* matrix_d_left_right_strict */

/*============================================================================*/
int matrix_d_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) {mes_prot("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) */
    if (!m_calloc(mean, rows)) {mes_proc(); goto STOP;}
    for (i = 0; i < rows; i++)
      mean[i] = gsl_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] = randvar_normal_density((double)j, mean[i], u);
      if (matrix[i][j] == -1) {mes_proc(); 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:
  return(res);
# undef CUR_PROC
} /* matrix_gaussrows_values */

/*============================================================================*/
void matrix_d_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;
    }
} /* matrix_d_const_preserve_struct */

/*============================================================================*/
void matrix_d_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] = gsl_rng_uniform(RNG);
    }
} /* matrix_d_random_preserve_struct */

/*============================================================================*/
void matrix_d_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];
} /* matrix_d_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) 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 matrix_cholesky(double **a, double *b, int dim, double *x) {
#define CUR_PROC "matrix_cholesky"
  int res = -1;
  double *p, *y;  
  if (!m_calloc(p, dim)) {mes_proc(); goto STOP;}
  if (!m_calloc(y, dim)) {mes_proc(); goto STOP;}
  lrdecomp(dim,a,p);
  lyequalsb(a,b,p,dim,y);
  ltranspxequalsy(a,y,p,dim,x);
  res = 0;
STOP:
  return(res);
#undef CUR_PROC
}

/* Finds the determinant of a symetric, positiv definit matrix. */
int matrix_det_symposdef(double **a, int dim, double *det) {
#define CUR_PROC "matrix_det_symposdef"
  int res = -1;
  int i;
  double *p, r;  
  if (!m_calloc(p, dim)) {mes_proc(); goto STOP;}
  lrdecomp(dim,a,p);
  *det = 1.0;
  for (i = 0; i < dim; i++) {
    r = 1/p[i];
    *det *= (r*r);
  }
  res = 0;
STOP:
  return(res);
#undef CUR_PROC
}

/*============================================================================*/
/* Copies a matrix, the memory allocations must to be done outside! */
void matrix_d_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 ];
}