smodel.c

来自「General Hidden Markov Model Library 一个通用」· C语言 代码 · 共 1,896 行 · 第 1/4 页

           smo->M, smo->N, (int) smo->s[0].density[0], smo->cos);
  fprintf (file, "\tprior = %.3f;\n", smo->prior);
  fprintf (file, "\tPi = vector {\n");
  ghmm_cmodel_Pi_print (file, smo, "\t", ",", ";");
  fprintf (file, "\t};\n");
  fprintf (file, "\tfix_state = vector {\n");
  ghmm_cmodel_fix_print (file, smo, "\t", ",", ";");
  fprintf (file, "\t};\n");
  /* Attention: assumption is: A_k are all the same */
  fprintf (file, "\tA = matrix {\n");
  ghmm_cmodel_Ak_print (file, smo, 0, "\t", ",", ";");
  fprintf (file, "\t};\n");
  /***/
  fprintf (file, "\tC = matrix {\n");
  ghmm_cmodel_C_print (file, smo, "\t", ",", ";");
  fprintf (file, "\t};\n\tMue = matrix {\n");
  ghmm_cmodel_Mue_print (file, smo, "\t", ",", ";");
  fprintf (file, "\t};\n\tU = matrix {\n");
  ghmm_cmodel_U_print (file, smo, "\t", ",", ";");
  fprintf (file, "\t};\n");
  fprintf (file, "};\n\n");
}                               /* ghmm_cmodel_print */


/*============================================================================*/
double ghmm_cmodel_calc_cmbm (ghmm_cmodel * smo, int state, int m, double omega)
{
  double bm = 0.0;
  switch (smo->s[state].density[m]) {
  case normal:
    bm = ighmm_rand_normal_density (omega, smo->s[state].mue[m],
                                 smo->s[state].u[m]);
    break;
  case normal_right:
    bm = ighmm_rand_normal_density_trunc (omega, smo->s[state].mue[m],
                                     smo->s[state].u[m],smo->s[state].a[m]);
    break;
  case normal_left:
    bm = ighmm_rand_normal_density_trunc (-omega, -smo->s[state].mue[m],
                                     smo->s[state].u[m],-smo->s[state].a[m]);
    break;
  case normal_approx:
    bm = ighmm_rand_normal_density_approx (omega,
                                        smo->s[state].mue[m],
                                        smo->s[state].u[m]);
    break;
  case uniform:
    bm = ighmm_rand_uniform_density (omega, 
                                  smo->s[state].mue[m],
                                  smo->s[state].u[m]);
    break;
  default:
    ighmm_mes (MES_WIN, "Warning: density function not specified!\n");
  }
  if (bm == -1) {
    ighmm_mes (MES_WIN, "Warning: density function returns -1!\n");
    bm = 0.0;
  }
  return (smo->s[state].c[m] * bm);
}                               /* ghmm_cmodel_calc_cmbm */


/*============================================================================*/
/* PDF(omega) in a given state */
double ghmm_cmodel_calc_b (ghmm_cmodel * smo, int state, double omega)
{
  int m;
  double b = 0.0;
  for (m = 0; m < smo->s[state].M; m++)
    b += ghmm_cmodel_calc_cmbm (smo, state, m, omega);
  return (b);
}                               /* ghmm_cmodel_calc_b */


/*============================================================================*/
double ghmm_cmodel_prob_distance (ghmm_cmodel * cm0, ghmm_cmodel * cm, int maxT,
                             int symmetric, int verbose)
{
#define CUR_PROC "ghmm_cmodel_prob_distance"

#define STEPS 100

  double p0, p;
  double d = 0.0;
  double *d1;
  ghmm_cseq *seq0 = NULL;
  ghmm_cseq *tmp = NULL;
  ghmm_cmodel *smo1, *smo2;
  int i, t, a, k, n;
  int true_len;
  int true_number;
  int left_to_right = 0;
  long total, index;
  int step_width = 0;
  int steps = 1;

  if (verbose) {                /* If we are doing it verbosely we want to have STEPS steps */
    step_width = maxT / STEPS;
    steps = STEPS;
  }
  else                          /* else just one */
    step_width = maxT;

  ARRAY_CALLOC (d1, steps);

  smo1 = cm0;
  smo2 = cm;

  for (k = 0; k < 2; k++) {

    seq0 = ghmm_cmodel_generate_sequences (smo1, 0, maxT + 1, 1, 0, maxT + 1);

    /*ghmm_cseq_print(stdout,seq0,0);*/

    if (seq0->seq_len[0] < maxT) {      /* There is an absorbing state */

      /* For now check that Pi puts all weight on state */
      /*
         t = 0;
         for (i = 0; i < smo1->N; i++) {
         if (smo1->s[i].pi > 0.001)
         t++;
         }    
         if (t > 1) {
         GHMM_LOG(LCONVERTED, "No proper left-to-right model. Multiple start states");
         goto STOP;
         } */

      left_to_right = 1;
      total = seq0->seq_len[0];

      while (total <= maxT) {

        /* create a additional sequences at once */
        a = (maxT - total) / (total / seq0->seq_number) + 1;
        /* printf("total=%d generating %d", total, a); */
        tmp = ghmm_cmodel_generate_sequences (smo1, 0, 0, a, 0, maxT + 1);
        ghmm_cseq_add (seq0, tmp);
        ghmm_cseq_free (&tmp);

        total = 0;
        for (i = 0; i < seq0->seq_number; i++)
          total += seq0->seq_len[i];
      }
    }

    if (left_to_right) {

      for (t = step_width, i = 0; t <= maxT; t += step_width, i++) {

        index = 0;
        total = seq0->seq_len[0];

        /* Determine how many sequences we need to get a total of t
           and adjust length of last sequence to obtain total of 
           exactly t */

        while (total < t) {
          index++;
          total += seq0->seq_len[index];
        }

        true_len = seq0->seq_len[index];
        true_number = seq0->seq_number;

        if ((total - t) > 0)
          seq0->seq_len[index] = total - t;
        seq0->seq_number = index + 1;

        if (ghmm_cmodel_likelihood (smo1, seq0, &p0) == -1) {
          /* error! */
          GHMM_LOG(LCONVERTED, "seq0 can't be build from smo1!");
          goto STOP;


        }
        n = ghmm_cmodel_likelihood (smo2, seq0, &p); /* ==-1: KEINE Seq. erzeugbar */
        if (n < seq0->seq_number) {
          ighmm_mes (MES_WIN,
               "problem: some seqences in seq0 can't be build from smo2\n");
          /* what shall we do now? */
          goto STOP;
        }

        /*printf("1/%d *(%f - %f)\n",t,p0,p);*/

        d = 1.0 / t * (p0 - p);

        if (symmetric) {
          if (k == 0)
            /* save d */
            d1[i] = d;
          else {
            /* calculate d */
            d = 0.5 * (d1[i] + d);
          }
        }

        if (verbose && (!symmetric || k == 1))
          printf ("%d\t%f\t%f\t%f\n", t, p0, p, d);

        seq0->seq_len[index] = true_len;
        seq0->seq_number = true_number;
      }
    }

    else {                      /* no left to right model */

      true_len = seq0->seq_len[0];

      for (t = step_width, i = 0; t <= maxT; t += step_width, i++) {
        seq0->seq_len[0] = t;

        if (ghmm_cmodel_likelihood (smo1, seq0, &p0) == -1) {
          /* error case */
          GHMM_LOG(LCONVERTED, "seq0 can't be build from smo1!");
          goto STOP;
        }
        n = ghmm_cmodel_likelihood (smo2, seq0, &p);/*== -1: KEINE Seq. erzeugbar*/
        if (n < seq0->seq_number) {
          ighmm_mes (MES_WIN,
               "problem: some sequences in seq0 can't be build from smo2\n");
          /* what shall we do now? */
          goto STOP;
        }

        /*printf("1/%d *(%f - %f)\n",t,p0,p);*/

        d = 1.0 / t * (p0 - p);

        if (symmetric) {
          if (k == 0)
            /* save d */
            d1[i] = d;
          else {
            /* calculate d */
            d = 0.5 * (d1[i] + d);
          }
        }

        if (verbose && (!symmetric || k == 1))
          printf ("%d\t%f\t%f\t%f\n", t, p0, p, d);

      }
      seq0->seq_len[0] = true_len;
    }

    if (symmetric) {
      ghmm_cseq_free (&seq0);
      smo1 = cm;
      smo2 = cm0;
    }
    else
      break;

  }                             /* k = 1,2 */

  ghmm_cseq_free (&seq0);

  return d;

STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ghmm_cseq_free (&seq0);
  return (-1.0);
#undef CUR_PROC
}


/*============================================================================*/
double ghmm_cmodel_calc_cmBm (ghmm_cmodel * smo, int state, int m, double omega)
{
  double Bm = 0.0;
  switch (smo->s[state].density[m]) {
  case normal_approx:
  case normal:
    Bm = ighmm_rand_normal_cdf (omega, smo->s[state].mue[m], smo->s[state].u[m]);
    break;
  case normal_right:
    Bm = ighmm_rand_normal_right_cdf (omega, smo->s[state].mue[m],
                              smo->s[state].u[m],smo->s[state].a[m]);
    break;
  case normal_left:
    Bm = (ighmm_rand_normal_right_cdf (omega, -smo->s[state].mue[m], smo->s[state].u[m],-smo->s[state].a[m]));
    break;
  case uniform:
    Bm = (ighmm_rand_uniform_cdf (omega, smo->s[state].mue[m], smo->s[state].u[m]));
    break;
  default:
    ighmm_mes (MES_WIN, "Warning: density function not specified!\n");
  }
  if (Bm == -1) {
    ighmm_mes (MES_WIN, "Warning: density function returns -1!\n");
    Bm = 0.0;
  }
  return (smo->s[state].c[m] * Bm);
}                               /* smodel_calc_Bm */


/*============================================================================*/
/* CDF(omega) in a given state */
double ghmm_cmodel_calc_B (ghmm_cmodel * smo, int state, double omega)
{
  int m;
  double B = 0.0;
  for (m = 0; m < smo->s[state].M; m++)
    B += ghmm_cmodel_calc_cmBm (smo, state, m, omega);
  return (B);
}                               /* ghmm_cmodel_calc_B */

/*============================================================================*/
/* What is the dimension of the modell ( = dimension of the parameter vektor) ?
   count the number of free parameters in a field of models; used for calc. BIC
   Only those parameters, that can be changed during  training.
   mixture coeff from smix and priors are not counted! 
*/
int ghmm_cmodel_count_free_parameter (ghmm_cmodel ** smo, int smo_number)
{
  int i, k;
  int pi_counted = 0, cnt = 0;

  for (k = 0; k < smo_number; k++) {
    pi_counted = 0;
    /* for states */
    for (i = 0; i < smo[k]->N; i++) {
      if (smo[k]->s[i].out_states > 1)
        /* multipl. with COS correct ??? */
        cnt += smo[k]->cos * (smo[k]->s[i].out_states - 1);
      if (smo[k]->s[i].pi != 0 && smo[k]->s[i].pi != 1) {
        pi_counted = 1;
        cnt++;
      }
      if (!smo[k]->s[i].fix) {
        if (smo[k]->s[i].M == 1)
          cnt += 2;             /* mu, sigma */
        else
          cnt += (3 * smo[k]->s[i].M);       /* c, mu, sigma */
      }
    }                           /* for (i ..) */
    if (pi_counted)
      cnt--;                    /* due to condition: sum(pi) = 1 */
    if (smo[k]->s[0].M > 1)
      cnt--;                    /* due to condition: sum(c) = 1 */
  }

  return cnt;
}

/*============================================================================*/
/* interval (a,b) with ~ B(a) < 0.01, B(b) > 0.99 */
void ghmm_cmodel_get_interval_B (ghmm_cmodel * smo, int state, double *a, double *b)
{
  int m;
  double mue, delta, max, min;
  *a = DBL_MAX;
  *b = -DBL_MAX;
  for (m = 0; m < smo->s[state].M; m++) {
    switch (smo->s[state].density[m]) {
    case normal:
    case normal_approx:
      mue = smo->s[state].mue[m];
      delta = 3 * sqrt (smo->s[state].u[m]);
      if (*a > mue - delta)
        *a = floor (mue - delta);
      if (*b < mue + delta)
        *b = ceil (mue + delta);
    case normal_right:
    case normal_left:
      if (smo->s[state].density[state] == normal_right && *a < smo->s[state].a[m])
        *a = smo->s[state].a[m];
      if (smo->s[state].density[state] == normal_left && *b > smo->s[state].a[m])
        *b = smo->s[state].a[m] ;
    break;
    case uniform:
      max = smo->s[state].mue[m];
      min = smo->s[state].u[m];
      *a = floor ((0.01 * (max - min))+min);
      *b = ceil ((0.99 * (max - min))+min);
    break;      
    default:
      ighmm_mes (MES_WIN, "Warning: density function not specified!\n");
    }
  }

  return;
}                               /* ghmm_cmodel_get_interval_B */


double ghmm_cmodel_get_transition(ghmm_cmodel* smo, int i, int j, int c)
{
# define CUR_PROC "ghmm_dmodel_get_transition"
  int out;

  if (smo->s && smo->s[i].out_a && smo->s[j].in_a) {
    for (out=0; out < smo->s[i].out_states; out++) {
      if (smo->s[i].out_id[out] == j)
        return smo->s[i].out_a[c][out];
    }
  }
  return 0.0;
# undef CUR_PROC
}   /* ghmm_cmodel_get_transition */

void ghmm_cmodel_set_transition (ghmm_cmodel* smo, int i, int j, int c, double prob)
{
# define CUR_PROC "ghmm_cmodel_set_transition"
  int in, out;

  if (smo->s && smo->s[i].out_a && smo->s[j].in_a) {
    for (out = 0; out < smo->s[i].out_states; out++) {
      if (smo->s[i].out_id[out] == j) {
        smo->s[i].out_a[c][out] = prob;
        /* fprintf(stderr, CUR_PROC": State %d, %d, = %f\n", i, j,
                prob); */
        break;
      }
    }

    for (in = 0; in < smo->s[j].in_states; in++) {
      if (smo->s[j].in_id[in] == i) {
        smo->s[j].in_a[c][in] = prob;
        break;
      }
    }
  }
# undef CUR_PROC
}   /* ghmm_cmodel_set_transition */


/*============================================================================*/
double ghmm_cmodel_ifunc (ghmm_cmodel * smo, int state, double c, double x)
{
  return (fabs (ghmm_cmodel_calc_B (smo, state, x) - c));
}

/*============================ TEST =========================================*/
#ifdef XXX
int smodel_test_callback(int pos){
   char* ModuleName = "class_change";
   char* FunctionName = "getClass";
   int class;
   PyObject *pName, *pModule, *pDict, *pFunc, *pArgs, *pValue;
    
   /*Py_Initialize();*/      /* Init Python Interpreter*/
   
   /*PyRun_SimpleString("import sys\n");*/
   /*PyRun_SimpleString("sys.stdout.write('Hello from an embedded Python Script\\n')\n"); */

   
   printf("C: Importing Python module ... ");
   pName = PyString_FromString(ModuleName);
   pModule = PyImport_Import(pName);       /* Import module*/
   pDict = PyModule_GetDict(pModule);
   printf("done.\n");    
    
   printf("C: Calling Python with value %d\n",pos);
   pFunc = PyDict_GetItemString(pDict, FunctionName);

   pArgs = PyTuple_New(1);
   pValue = PyInt_FromLong(pos);

   PyTuple_SetItem(pArgs, 0, pValue); 
   pValue = PyObject_CallObject(pFunc, pArgs); /* Calling Python */

   /* parsing the result from Python to C data type
   class = PyInt_AsLong(pValue);
   printf("C: The returned class is %d\n",class); */
     
   /*Py_Finalize();         */
   
   return class; 
 
}
#endif