vector.c

来自「General Hidden Markov Model Library 一个通用」· C语言 代码 · 共 206 行

/*******************************************************************************
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version 0.8_beta1, see http://ghmm.org
*
*       Filename: ghmm/ghmm/vector.c
*       Authors:  Bernhard Knab
*
*       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: 1588 $
*                       from $Date: 2006-01-11 15:15:50 +0100 (Wed, 11 Jan 2006) $
*             last change by $Author: grunau $.
*
*******************************************************************************/

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

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

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

#include "mes.h"
#include "vector.h"
#include "rng.h"
#include "mprintf.h"
#include "ghmm_internals.h"
#include "obsolete.h"

/*============================================================================*/
/* Scales the elements of a vector to have the sum 1 */
/* PROBLEM: Entries can get very small and be rounded to 0 */
int ighmm_cvector_normalize (double *v, int len)
{
#define CUR_PROC "ighmm_cvector_normalize"
  int i;
  double sum = 0.0;
  char * estr;

  for (i = 0; i < len; i++)
    sum += v[i];
  if (i>0 && sum<GHMM_EPS_PREC) {
    estr = ighmm_mprintf(NULL, 0, "Can't normalize vector. Sum smaller than %g\n"
			, GHMM_EPS_PREC);
    GHMM_LOG(LWARN, estr);
    m_free(estr);
    return (-1);
  }
  for (i = 0; i < len; i++)
    v[i] /= sum;
  return 0;
#undef CUR_PROC
}                               /* ighmm_cvector_normalize */

#ifdef GHMM_OBSOLETE
/*============================================================================*/
void ighmm_cvector_const_values (double *v, int len, double c)
{
  int i;
  for (i = 0; i < len; i++)
    v[i] = c;
}                               /* ighmm_cvector_const_values */


/*============================================================================*/
void ighmm_cvector_const_preserve_struct (double *v, int len, double c)
{
  int i;
  for (i = 0; i < len; i++) {
    if (v[i] != 0.0)
      v[i] = c;
  }
}                               /* ighmm_cvector_const_preserve_struct */


/*============================================================================*/
void ighmm_cvector_random_values (double *v, int len)
{
  int i;
  for (i = 0; i < len; i++)
    v[i] = GHMM_RNG_UNIFORM (RNG);
}                               /* ighmm_cvector_random_values */


/*============================================================================*/
void ighmm_cvector_random_preserve_struct (double *v, int len)
{
  int i;
  for (i = 0; i < len; i++) {
    if (v[i] != 0.0)
      v[i] = GHMM_RNG_UNIFORM (RNG);
  }
}                               /* ighmm_cvector_random_preserve_struct */
#endif /* GHMM_OBSOLETE */

/*============================================================================*/
void ighmm_cvector_print (FILE * file, double *vector, int len,
                     char *tab, char *separator, char *ending)
{
  int j;

  fprintf (file, "%s", tab);
  if (len > 0)
    fprintf (file, "%6.2f", vector[0]);

  for (j = 1; j < len; j++)
    fprintf (file, "%s %6.2f", separator, vector[j]);
  fprintf (file, "%s\n", ending);
}                               /* ighmm_cvector_print */

/*============================================================================*/
void ighmm_cvector_print_prec (FILE * file, double *vector, int len, int width,
                          int prec, char *tab, char *separator, char *ending)
{
  int j;
  if (len > 0)
    fprintf (file, "%s%*.*f", tab, width, prec, vector[0]);
  for (j = 1; j < len; j++)
    fprintf (file, "%s %*.*f", separator, width, prec, vector[j]);
  fprintf (file, "%s\n", ending);
}                               /* ighmm_cvector_print */

/*============================================================================*/
void ighmm_dvector_print (FILE * file, int *vector, int len,
                     char *tab, char *separator, char *ending)
{
  int j;
  fprintf (file, "%s", tab);
  if (len > 0)
    fprintf (file, "%3d", vector[0]);
  for (j = 1; j < len; j++)
    fprintf (file, "%s %3d", separator, vector[j]);
  fprintf (file, "%s\n", ending);
}                               /* ighmm_dvector_print */

#ifdef GHMM_OBSOLETE
/*============================================================================*/
int ighmm_cvector_mat_times_vec (double **A, double *x, int n, int m, double *v)
{
#define CUR_PROC "ighmm_cvector_mat_times_vec"
  int i, j;
  for (i = 0; i < n; i++) {
    v[i] = 0.0;
    for (j = 0; j < m; j++)
      v[i] += A[i][j] * x[j];
  }
  return 0;
#undef CUR_PROC
}                               /* ighmm_cvector_mat_times_vec */
#endif /* GHMM_OBSOLETE */


/*============================================================================*/
double ighmm_cvector_log_sum (double *a, int length) {
#define CUR_PROC "ighmm_cvector_log_sum"
  int i;
  double max = 1.0;
  int argmax = 0;
  double result;

  /* find maximum value in a: */
  for (i = 0; i < length; i++)
    if (max == 1.0 || (a[i] > max && a[i] != 1.0)) {
      max = a[i];
      argmax = i;
    }

  /* calculate max+log(1+sum[i!=argmax; exp(a[i]-max)])  */
  result = 1.0;
  for (i = 0; i < length; i++) {
    if (a[i] != 1.0 && i != argmax)
      result += exp (a[i] - max);
  }
  result = log (result);
  result += max;
  return result;
#undef CUR_PROC
}