smodel.c

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

void smodel_C_print(FILE *file, smodel *smo, char *tab, char *separator, 
		    char *ending) {
  int i, j;
  for (i = 0; i < smo->N; i++) {
    fprintf(file, "%s", tab);
    fprintf(file, "%.4f", smo->s[i].c[0]);
    for (j = 1; j < smo->M; j++)
      fprintf(file, "%s %.4f", separator, smo->s[i].c[j]);
    fprintf(file, "%s\n", ending);
  }
} /* smodel_C_print */


/*============================================================================*/
void smodel_Mue_print(FILE *file, smodel *smo, char *tab, char *separator, 
		      char *ending) {
  int i, j;
  for (i = 0; i < smo->N; i++) {
    fprintf(file, "%s", tab);
    fprintf(file, "%.4f", smo->s[i].mue[0]);
    for (j = 1; j < smo->M; j++)
      fprintf(file, "%s %.4f", separator, smo->s[i].mue[j]);
    fprintf(file, "%s\n", ending);
  }
} /* smodel_Mue_print */


/*============================================================================*/
void smodel_U_print(FILE *file, smodel *smo, char *tab, char *separator, 
			char *ending) {
  /* attention: choose precision big enough to allow printing of  
     EPS_U in const.h */
  int i, j;
  for (i = 0; i < smo->N; i++) {
    fprintf(file, "%s", tab);
    fprintf(file, "%.4f", smo->s[i].u[0]);
    for (j = 1; j < smo->M; j++)
      fprintf(file, "%s %.4f", separator, smo->s[i].u[j]);
    fprintf(file, "%s\n", ending);
  }
} /* smodel_U_print */


/*============================================================================*/
void smodel_Pi_print(FILE *file, smodel *smo, char *tab, char *separator, 
		     char *ending) {
  int i;
  fprintf(file, "%s%.4f", tab, smo->s[0].pi);
  for (i = 1; i < smo->N; i++)
    fprintf(file, "%s %.4f", separator, smo->s[i].pi);
  fprintf(file, "%s\n", ending);
} /* smodel_Pi_print */

/*============================================================================*/
void smodel_fix_print(FILE *file, smodel *smo, char *tab, char *separator, 
		     char *ending) {
  int i;
  fprintf(file, "%s%d", tab, smo->s[0].fix);
  for (i = 1; i < smo->N; i++)
    fprintf(file, "%s %d", separator, smo->s[i].fix);
  fprintf(file, "%s\n", ending);
} 


/*============================================================================*/
void smodel_print(FILE *file, smodel *smo) {
  int k;
  fprintf(file, "SHMM = {\n\tM = %d;\n\tN = %d;\n\tdensity = %d;\n\tcos = %d;\n", 
	  smo->M, smo->N, (int)smo->density, smo->cos);
  fprintf(file, "\tprior = %.5f;\n", smo->prior);
  fprintf(file, "\tPi = vector {\n");
  smodel_Pi_print(file, smo, "\t", ",", ";");
  fprintf(file, "\t};\n");
  fprintf(file, "\tfix_state = vector {\n");
  smodel_fix_print(file, smo, "\t", ",", ";");
  fprintf(file, "\t};\n");
  for (k = 0; k < smo->cos; k++) {
    fprintf(file, "\tAk_%d = matrix {\n", k);
    smodel_Ak_print(file, smo, k, "\t", ",", ";");
    fprintf(file, "\t};\n");
  }
  fprintf(file, "\tC = matrix {\n");
  smodel_C_print(file, smo,"\t", ",", ";");
  fprintf(file, "\t};\n\tMue = matrix {\n");
  smodel_Mue_print(file, smo,"\t", ",", ";");
  fprintf(file, "\t};\n\tU = matrix {\n");
  smodel_U_print(file, smo,"\t", ",", ";");
  fprintf(file, "\t};\n");
  fprintf(file, "};\n\n");
} /* smodel_print */

/*============================================================================*/
/* needed for hmm_input: only one A (=Ak_1 = Ak_2...) is written */
void smodel_print_oneA(FILE *file, smodel *smo) {
  fprintf(file, "SHMM = {\n\tM = %d;\n\tN = %d;\n\tdensity = %d;\n\tcos = %d;\n", 
	  smo->M, smo->N, (int)smo->density,smo->cos);
  fprintf(file, "\tprior = %.3f;\n", smo->prior);
  fprintf(file, "\tPi = vector {\n");
  smodel_Pi_print(file, smo, "\t", ",", ";");
  fprintf(file, "\t};\n");
  fprintf(file, "\tfix_state = vector {\n");
  smodel_fix_print(file, smo, "\t", ",", ";");
  fprintf(file, "\t};\n");
  /* Attention: assumption is: A_k are all the same */
  fprintf(file, "\tA = matrix {\n");
  smodel_Ak_print(file, smo, 0, "\t", ",", ";");
  fprintf(file, "\t};\n");
  /***/
  fprintf(file, "\tC = matrix {\n");
  smodel_C_print(file, smo,"\t", ",", ";");
  fprintf(file, "\t};\n\tMue = matrix {\n");
  smodel_Mue_print(file, smo,"\t", ",", ";");
  fprintf(file, "\t};\n\tU = matrix {\n");
  smodel_U_print(file, smo,"\t", ",", ";");
  fprintf(file, "\t};\n");
  fprintf(file, "};\n\n");
} /* smodel_print */


/*============================================================================*/
double smodel_calc_cmbm(smodel *smo, int state, int m, double omega) {
  double bm = 0.0;
  switch (smo->density) {
  case normal: bm = randvar_normal_density(omega, smo->s[state].mue[m], 
					   smo->s[state].u[m]);
    break;
  case normal_pos: bm = randvar_normal_density_pos(omega,smo->s[state].mue[m],
						   smo->s[state].u[m]);
    break;
  case normal_approx: bm = randvar_normal_density_approx(omega, 
							 smo->s[state].mue[m],
							 smo->s[state].u[m]);
    break;
  default: mes(MES_WIN, "Warning: density function not specified!\n"); 
  }
  if (bm == -1) {
    mes(MES_WIN, "Warning: density function returns -1!\n"); 
    bm = 0.0;
  }
  return(smo->s[state].c[m] * bm);  
} /* smodel_calc_cmbm */


/*============================================================================*/
/* PDF(omega) in a given state */
double smodel_calc_b(smodel *smo, int state, double omega) {
  int m;
  double b = 0.0;
  for (m = 0; m < smo->M; m++)
    b += smodel_calc_cmbm(smo, state, m, omega);
  return(b);
} /* smodel_calc_b */


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

#define STEPS 100

  double p0, p;
  double d = 0.0;
  double *d1;
  sequence_d_t *seq0 = NULL;
  sequence_d_t *tmp = NULL;
  smodel *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;

  if( !m_calloc(d1, steps) ) {mes_proc();goto STOP;} 

  smo1 = cm0;
  smo2 = cm;
 
  for (k = 0; k < 2; k++) {
    
    seq0 = smodel_generate_sequences(smo1, 0, maxT+1, 1, 0, maxT+1);
    
	//sequence_d_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) {
	mes_prot("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 = smodel_generate_sequences(smo1, 0, 0, a, 0, maxT+1);
	sequence_d_add(seq0, tmp);
	sequence_d_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 (smodel_likelihood(smo1, seq0, &p0) == -1) {
	  /* error! */
	  mes_prot("seq0 can't be build from smo1!");
	  goto STOP;
	
	
	}
	n = smodel_likelihood(smo2, seq0, &p); /* ==-1: KEINE Seq. erzeugbar*/
	if (n < seq0->seq_number) {
	  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 (smodel_likelihood(smo1, seq0, &p0) == -1) {
	  /* error case */
	  mes_prot("seq0 can't be build from smo1!");
	  goto STOP;
	}
	n = smodel_likelihood(smo2, seq0, &p); /*== -1: KEINE Seq. erzeugbar*/
	if (n < seq0->seq_number) {
	  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) {
      sequence_d_free(&seq0);
      smo1 = cm ;
      smo2 = cm0;
    }
    else
      break;
   
  } /* k = 1,2 */ 

  sequence_d_free(&seq0);
  
  return d;

STOP:
  sequence_d_free(&seq0);
  return(-1.0);
#undef CUR_PROC
}


/*============================================================================*/
double smodel_calc_cmBm(smodel *smo, int state, int m, double omega) {
  double Bm = 0.0;
  switch (smo->density) {
  case normal_approx: 
  case normal:
    Bm = randvar_normal_cdf(omega, smo->s[state].mue[m], smo->s[state].u[m]);
    break;
  case normal_pos: 
    Bm = randvar_normal_pos_cdf(omega,smo->s[state].mue[m],smo->s[state].u[m]);
    break;
  default: mes(MES_WIN, "Warning: density function not specified!\n"); 
  }
  if (Bm == -1) {
    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 smodel_calc_B(smodel *smo, int state, double omega) {
  int m;
  double B = 0.0;
  for (m = 0; m < smo->M; m++)
    B += smodel_calc_cmBm(smo, state, m, omega);
  return(B);
} /* smodel_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 smodel_count_free_parameter(smodel **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]->M == 1) cnt += 2; /* mu, sigma */
	else cnt += (3 * smo[k]->M); /* c, mu, sigma */
      }
    } /* for (i ..) */
    if (pi_counted) cnt--; /* due to condition: sum(pi) = 1 */
    if (smo[k]->M > 1) cnt--; /* due to condition: sum(c) = 1 */
  }

  return cnt;
}

/*============================================================================*/
/* interval (a,b) with ~ B(a) < 0.01, B(b) > 0.99 */  
void smodel_get_interval_B(smodel *smo, int state, double *a, double *b) {
  int m;
  double mue, delta;
  switch (smo->density) {
  case normal:
  case normal_approx: 
  case normal_pos: 
    *a = DBL_MAX;
    *b = -DBL_MAX;
    for (m = 0; m < smo->M; m++) {
      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);
    }
    break;
  default: mes(MES_WIN, "Warning: density function not specified!\n"); 
  }
  if (smo->density == normal_pos && *a < 0.0) 
    *a = 0.0;
  return;
} /* smodel_get_interval_B */

/*============================================================================*/
double smodel_ifunc(smodel *smo, int state, double c, double x) {
  return( fabs( smodel_calc_B(smo, state, x) - c));
}