randvar.c

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

/*******************************************************************************
  author       : Bernhard Knab
  filename     : ghmm/ghmm/randvar.c
  created      : TIME: 16:40:03     DATE: Wed 17. February 1999
  $Id: randvar.c,v 1.17 2002/03/05 15:37:37 pipenb Exp $

Copyright (C) 1998-2001, ZAIK/ZPR, Universit則 zu K鱨n

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


*******************************************************************************/

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

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

#include <math.h>
#ifdef HAVE_LIBPTHREAD
# include <pthread.h>
#endif /* HAVE_LIBPTHREAD */
#include "mes.h"
#include "mprintf.h"
#include "randvar.h"
#include "rng.h"
#include "float.h"
#include "const.h"

#ifdef DO_WITH_GSL

# include <gsl/gsl_math.h>
# include <gsl/gsl_sf_erf.h>
# include <gsl/gsl_randist.h>

/* missing functions in gsl-0.7 */

#define EVAL_RESULT(fn) \
   gsl_sf_result result; \
   int status = fn; \
   if (status == GSL_EDOM) { \
     return GSL_NAN; \
   } else if (status != GSL_SUCCESS) { \
     GSL_ERROR(#fn, status); \
   } ; \
   return result.val;

# ifndef HAVE_GSL_SF_ERF
double gsl_sf_erfc(double x)
{
  EVAL_RESULT(gsl_sf_erfc_e(x, &result));
}
# endif

# ifndef HAVE_GSL_SF_ERFC
double gsl_sf_erf(double x)
{
  EVAL_RESULT(gsl_sf_erf_e(x, &result));
}
# endif

#else

#include "scanner.h"

#endif /* DO_WITH_GSL */

/* A list of already calculated values of the density function of a 
   N(0,1)-distribution, with x in [0.00, 19.99] */
#define PDFLEN 2000
#define X_STEP_PDF 0.01  /* step size */
#define X_FAKT_PDF 100   /* equivalent to step size */
static double pdf_stdnormal[PDFLEN];
static int pdf_stdnormal_exists = 0;

/* A list of already calulated values PHI of the Gauss distribution is
   read in, x in [-9.999, 0] */
#define X_STEP_PHI 0.001  /* step size */
#define X_FAKT_PHI 1000   /* equivalent to step size*/
static int PHI_len = 0;
#if 0 /* needed for randvar_get_xPHIxgleichPHIy() */
static double x_PHI_xy = -1.0;
#endif
static double x_PHI_1 = -1.0;


#ifndef DO_WITH_GSL
static double* PHI = NULL;

/*----------------------------------------------------------------------------*/
static int randvar_read_PHI() {
# define CUR_PROC "randvar_read_PHI"
  int res = -1;
  char filename[] = PHI_DATA_FILE;
  #warning "PHI_DATA_FILE depreciated!"
  scanner_t *s = NULL;
  s = scanner_alloc(filename); if(!s) {mes_proc(); goto STOP;}
  scanner_get_name(s);
  scanner_consume(s, '='); if(s->err) goto STOP; 
  if (!strcmp(s->id, "PHI")) {
    scanner_consume(s, '{'); if (s->err) goto STOP;
    PHI = scanner_get_double_earray(s, &PHI_len); if (s->err) goto STOP;
    scanner_consume(s, ';'); if (s->err) goto STOP;
    scanner_consume(s, '}'); if (s->err) goto STOP;
    scanner_consume(s, ';'); if (s->err) goto STOP;
  }
  else {
    scanner_error(s, "unknown identifier"); goto STOP;
  }
  /* printf("%.4f\n", PHI[PHI_len-1]); */

  res = 0;;
STOP:
  scanner_free(&s);
  return res;
# undef CUR_PROC
} /* randvar_read_PHI */


/*----------------------------------------------------------------------------*/
static int randvar_init_PHI() {
# define CUR_PROC "randvar_init_PHI"
#ifdef HAVE_LIBPTHREAD
  static pthread_mutex_t lock;
#endif /* HAVE_LIBPTHREAD */
  /* Read PHI in */
  if (!PHI_len) {
#ifdef HAVE_LIBPTHREAD
    pthread_mutex_lock(&lock); /* Put on a lock, because the clustering is parallel */
#endif /* HAVE_LIBPTHREAD */
    if (randvar_read_PHI() == -1) {mes_proc(); goto STOP;};
#ifdef HAVE_LIBPTHREAD
    pthread_mutex_unlock(&lock); /* Take the lock off */
#endif /* HAVE_LIBPTHREAD */
  }
  return 0;
STOP:
  return(-1);
# undef CUR_PROC
} /* randvar_init_PHI */


#endif /* DO_WITH_GSL */


/*============================================================================*/
/* needed by pmue_interpol */

double randvar_get_xfaktphi() {
  return X_FAKT_PHI;
}

double randvar_get_xstepphi() {
  return X_STEP_PHI;
}

double randvar_get_philen() {
#ifdef DO_WITH_GSL
  return PHI_len;
#else
  return randvar_get_xPHIless1()/X_STEP_PHI;
#endif
}


/*============================================================================*/
double randvar_get_PHI(double x) {
# define CUR_PROC "randvar_get_PHI"

#ifdef DO_WITH_GSL
  return (gsl_sf_erf(x*M_SQRT1_2)+1.0)/2.0;
#else
  int i;
  double phi_x;

  if (randvar_init_PHI() == -1) {mes_proc(); goto STOP;}

  /* Linear interpolation (Alternative: Round off with i=m_int(fabs(x)*X_FAKT))*/
  i = (int)(fabs(x) * X_FAKT_PHI); 
  if (i >= PHI_len-1) {
    i = PHI_len-1;
    phi_x = PHI[i];
  }
  else
    phi_x = PHI[i] + (fabs(x) - i*X_STEP_PHI) * (PHI[i+1]-PHI[i])/X_STEP_PHI;
  /* NOTA BENE: PHI is tabulated for negative values! */
  if (x > 0.0)
    return(1.0 - phi_x);
  else
    return(phi_x);
#endif /* DO_WITH_GSL */
  return(-1.0);
# undef CUR_PROC
} /* randvar_get_PHI */


/*============================================================================*/
/* When is PHI[x,0,1] == 1? */
double randvar_get_xPHIless1() {
# define CUR_PROC "randvar_get_xPHIless1"
#ifdef DO_WITH_GSL
  if (x_PHI_1 == -1)
    {
      double low,up,half;
      low=0;
      up=100;
      while (up-low>0.001)
	{
	  half=(low+up)/2.0;
	  if (randvar_get_PHI(half)<1.0)
	    low=half;
	  else
	    up=half;
	}
      x_PHI_1=low;
    }
  return(x_PHI_1);
#else
  double x;
  int i;
  if (x_PHI_1 == -1) {
    if (randvar_init_PHI() == -1) {mes_proc(); goto STOP;}    
    /* The last value of the table is 1 */
    for (x = (PHI_len-1)*X_STEP_PHI, i=PHI_len-1; i > 0; x -= X_STEP_PHI, i--)
      if (randvar_get_PHI(-x) > 0.0)
	break;
    /* Modification: x exactly between 2 sampling points! */
    x_PHI_1 = x - (double)X_STEP_PHI/2.0;
  }
  return(x_PHI_1);
#endif
  return(-1.0);
# undef CUR_PROC
}

#if 0 /* is not in use */
/*============================================================================*/
/* When is PHI[x,0,1] ==  PHI[y,0,1] for consecutive x and y ?*/
double randvar_get_xPHIxgleichPHIy() {
# define CUR_PROC "randvar_get_xPHIxgleichPHIy"
  double x,y;
  int i;
  if (x_PHI_xy == -1) {
    if (randvar_init_PHI() == -1) {mes_proc(); goto STOP;}     
    y = -1.0;
    for (x = 0.0, i=0; i < PHI_len; x += X_STEP_PHI, i++) {
      if (randvar_get_PHI(-x) == randvar_get_PHI(-y))
	break;
      y = x;
    }
    x_PHI_xy = y;
  }
  return(x_PHI_xy);
STOP:
  return(-1.0);
# undef CUR_PROC
}
#endif /* 0 */

/*============================================================================*/
double randvar_get_1overa(double x, double mean, double u) {
  /* Calulates 1/a(x, mean, u), with a = the integral from x til \infty over
     the Gauss density function */
# define CUR_PROC "randvar_get_1overa"

#ifdef DO_WITH_GSL
  double erfc_value;
#else
  int i;
  double y, z, phi_z, a;
#endif

  if (u <= 0.0) {
    mes_prot("u <= 0.0 not allowed\n"); goto STOP;
  }

#ifdef DO_WITH_GSL
  /* int gsl_sf_erfc (double x) 
     erfc(x) = 1 - erf(x) = 2/\sqrt(\pi) \int_x^\infty \exp(-t^2)
  */
  erfc_value=gsl_sf_erfc((x-mean)/sqrt(u*2));
  if (erfc_value <= DBL_MIN)
    {
      mes(MES_WIN,
	  "a ~= 0.0 critical! (mue = %.2f, u =%.2f)\n",
	  mean, u);
      return(erfc_value);
    }
  else
    return(2.0/erfc_value);
#else

  if (randvar_init_PHI() == -1) {mes_proc(); goto STOP;};

  y = 1/sqrt(u);
  z = (x - mean) * y; 
 /* Linear interpolation (Alternative: Round off with i=m_int(fabs(z)*X_FAKT))*/
  i = (int)(fabs(z)*X_FAKT_PHI); 
 
  if (i >= PHI_len-1) { 
    i = PHI_len-2;
    /* Originally:
       i = PHI_len-1; but then, the last value in the table is zero! */
    phi_z = PHI[i];
  }
  else
    phi_z = PHI[i] + (fabs(z)-i*X_STEP_PHI) * (PHI[i+1]-PHI[i])/X_STEP_PHI;
  /* NOTA BENE: PHI is tabulated for negative values! */
  if (z > 0.0) {
    if (phi_z ==  0) {
      mes_proc(); goto STOP;
    }
    else
      a = 1 / phi_z; /* PHI between 0.5 and 1 */ /*??? between 0.5 and 0 ! */   
  }
  else {