sdmodel.c

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

/*
 * Before calling generate sequence, the global GSL random number generator must
 * be initialized by calling gsl_rng_init().
 *
 */
//returns the extended sequence struct with state matrix
sequence_t *sdmodel_generate_sequences(sdmodel* mo, int seed, int global_len, long seq_number, int Tmax) {
# define CUR_PROC "sdmodel_generate_sequences"

  /* An end state is characterized by not having out-going transition. */
  unsigned long tm; /* Time seed */
  sequence_t *sq = NULL;
  int state, n, i, j, m,  reject_os, reject_tmax, badseq, trans_class;
  double p, sum, osum = 0.0;
  int len = global_len, up = 0, stillbadseq = 0, reject_os_tmp = 0;
  int obsLength = 0;
  double dummy = 0.0;
  int silent_len = 0, badSilentStates = 0;
  int	lastStateSilent = 0;
  int matchcount = 0;

  sq = sequence_calloc(seq_number);
  
  if (!sq) { mes_proc(); goto STOP; }
  if (len <= 0)
    /* A specific length of the sequences isn't given. As a model should have
       an end state, the konstant MAX_SEQ_LEN is used. */
    len = (int)MAX_SEQ_LEN;
  
  if (seed > 0) {
    gsl_rng_set(RNG,seed);
  } else
    {
      tm = rand();
      gsl_rng_set(RNG, tm);
      fprintf(stderr, "# using GSL rng '%s' seed=%ld\n", gsl_rng_name(RNG), tm);      
    }

  n = 0;
  reject_os = reject_tmax = 0;

  while (n < seq_number) 
    {
      /* Test: A new seed for each sequence */
      /*   gsl_rng_timeseed(RNG); */
      stillbadseq = badseq = 0;
      matchcount = 0;
      if(!m_calloc(sq->seq[n], len)) {mes_proc(); goto STOP;}
      if(!m_calloc(sq->states[n], len)) {mes_proc(); goto STOP;}

      /* Get a random initial state i */
      p = gsl_rng_uniform(RNG);
      sum = 0.0;
      for (i = 0; i < mo->N; i++) {
	sum += mo->s[i].pi;
	if (sum >= p)
	  break;
      }
      /* assert( !mo->silent[i] ); */

      if ( mo->model_type == kSilentStates ) { 	
	
	if ( !mo->silent[i] ) { /* fDte emits */
	  lastStateSilent = 0;
	  silent_len      = 0;				 
	  /* Get a random initial output m */
	  p = gsl_rng_uniform(RNG);
	  sum = 0.0;   
	  for (m = 0; m < mo->M; m++) {
	    sum += mo->s[i].b[m];
	    if (sum >= p)
	      break;
	  }
	  sq->seq[n][0] = m;
	  sq->states[n][0] = i;
	  if (mo->s[i].countme){
	    matchcount++;
	  }
	  state = 1;
	} else { /* silent state: we do nothing, no output */
	  state = 0;
	  lastStateSilent = 1;
	  silent_len      = 1;	
	}
      } else {
	/* Get a random initial output m */
	p = gsl_rng_uniform(RNG);
	sum = 0.0;   
	for (m = 0; m < mo->M; m++) {
	  sum += mo->s[i].b[m];
	  if (sum >= p)
	    break;
	}
	sq->seq[n][0] = m;
	sq->states[n][0] = i;
	if (mo->s[i].countme){
	  matchcount++;
	}
	state = 1;
      }
      
      obsLength = 0;

      /* class NOT dependent on observable but on path!!!! */
      while (state < len && 
	     obsLength < Tmax ) {
	/* Get a new state */

	if (mo->cos>1)
	  trans_class = mo->get_class(sq->seq[n],state);
	else
	  trans_class = 0;
	p = gsl_rng_uniform(RNG);
	sum = 0.0;   
	for (j = 0; j < mo->s[i].out_states; j++) {
	  sum += mo->s[i].out_a[trans_class][j];   
	  if (sum >= p) break;
	}
	
	if (sum == 0.0) {
	  if (mo->s[i].out_states > 0) {
	    /* Repudiate the sequence, if all smo->s[i].out_a[class][.] == 0,
	       that is, class "class" isn't used in the original data:
	       go out of the while-loop, n should not be counted. */
	    /* printf("Zustand %d, class %d, len %d out_states %d \n", i, class,
	       state, smo->s[i].out_states); */
	    badseq = 1;
	    /* break; */
	    
	    /* Try: If the class is "empty", try the neighbour class;
	       first, sweep down to zero; if still no success, sweep up to
	       COS - 1. If still no success --> Repudiate the sequence. */
	    if (trans_class > 0 && up == 0) {
	      //trans_class--;
	      continue;
	    } else {
	      if (trans_class < mo->cos - 1) {
		//trans_class++;          // as it is not dependent on time!
		up = 1;
		continue;
	      } else {
		stillbadseq = 1;
		break;
	      }
	    }
	  } else {
	    /* An end state is reached, get out of the while-loop */
	    break;
	  }
	}
	
	i = mo->s[i].out_id[j];

	if (mo->model_type == kSilentStates && mo->silent[i]) { /* Get a silent state i */
	  silent_len++;
	  if (silent_len >= Tmax) {
	    badSilentStates = 1;
	    break;         /* reject this sequence */ 
	  } else {
	    badSilentStates = 0;
	    lastStateSilent = 1;
	  }
	} else {
	  lastStateSilent = 0;
	  silent_len  = 0;
	  /* Get a random output m from state i */
	  p = gsl_rng_uniform(RNG);
	  sum = 0.0;   
	  for (m = 0; m < mo->M; m++) {
	    sum += mo->s[i].b[m];
	    if (sum >= p)
	      break;
	  }
	  sq->seq[n][state] = m;
	  sq->states[n][state] = i;
	  if (mo->s[i].countme){
	    matchcount++;
	  }
	  state++;
	}			
	 
	/* Decide the class for the next step */
	//trans_class = mo->get_class(sq->seq[n],state);
	up = 0;
	obsLength++;
	
      } /* while (state < len) , global_len depends on the data */  

next_sequence:;   

      if (badseq) {
	reject_os_tmp++;
      }
      
      if (stillbadseq) {
	reject_os++;
	m_free(sq->seq[n]);
	m_free(sq->states[n]);
	/*  printf("cl %d, s %d, %d\n", class, i, n); */
      } else {
	if (badSilentStates) {
	  reject_tmax++;
	  m_free(sq->seq[n]);
	  m_free(sq->states[n]);
	} else {
	  if ( obsLength > Tmax ) { 
	    reject_tmax++;
	    m_free(sq->seq[n]);
	    m_free(sq->seq[n]);
	  } else {
	    if (state < len)
	      {
		if(m_realloc(sq->seq[n], state)) {mes_proc(); goto STOP;}
		if(m_realloc(sq->states[n], state)) {mes_proc(); goto STOP;}
	      }
	    sq->seq_len[n] = state;
	    /* sq->seq_label[n] = label; */
	    /* vector_d_print(stdout, sq->seq[n], sq->seq_len[n]," "," ",""); */
	    n++;
	  }
	}
      }
	

      /*    printf("reject_os %d, reject_tmax %d\n", reject_os, reject_tmax); */
      if (reject_os > 10000) {
	mes_prot("Reached max. no. of rejections\n");
	break;
      }
      if (!(n%1000)) printf("%d Seqs. generated\n", n);
    } /* n-loop, while (n < seq_number) */
	  

  if (reject_os > 0) printf("%d sequences rejected (os)!\n", reject_os);
  if (reject_os_tmp > 0) printf("%d sequences changed class\n", 
				reject_os_tmp - reject_os);	
  if (reject_tmax > 0) printf("%d sequences rejected (Tmax, %d)!\n", 
			      reject_tmax, Tmax);


  return(sq);
 STOP:

  sequence_free(&sq);
  return(NULL);
# undef CUR_PROC
} /* data */

/*=======================================================================
generate sequences by calling generate_sequences_ext
========================================*/

/*
sequence_t *sdmodel_generate_sequences(sdmodel* mo, int seed, int global_len,
				     long seq_number, int Tmax) {
  sequence_ext_t *sq = sdmodel_generate_sequences_ext(mo,seed,global_len,seq_number,Tmax);
  printf("jaja\n\n");
  return sequence_unext(&sq);
}
*/
/*============================================================================*/
/* Some outputs */
/*============================================================================*/

void sdmodel_states_print(FILE *file, sdmodel *mo) {
  int i, j;
  fprintf(file, "Modelparameters: \n M = %d \t N = %d\n", 
	 mo->M, mo->N);
  for (i = 0; i < mo->N; i++) {
    fprintf(file, "\nState %d \n PI = %.3f \n out_states = %d \n in_states = %d \n",
	   i, mo->s[i].pi, mo->s[i].out_states, mo->s[i].in_states);
    fprintf(file, " Output probability:\t");
    for (j = 0; j < mo->M; j++)
      fprintf(file, "%.3f \t", mo->s[i].b[j]);
    fprintf(file, "\n Transition probability \n");
    fprintf(file, "  Out states (Id, a):\t");
    for (j = 0; j < mo->s[i].out_states; j++)
      fprintf(file, "(%d, %.3f) \t", mo->s[i].out_id[j], mo->s[i].out_a[j]);
    fprintf(file, "\n");
    fprintf(file, "  In states (Id, a):\t");
    for (j = 0; j < mo->s[i].in_states; j++)
      fprintf(file, "(%d, %.3f) \t", mo->s[i].in_id[j], mo->s[i].in_a[j]);
    fprintf(file, "\n");
  }
} /* model_states_print */

/*============================================================================*/

void sdmodel_Ak_print(FILE *file, sdmodel *mo, int k, char *tab, char *separator, 
		   char *ending) {
  int i, j, out_state;
  for (i = 0; i < mo->N; i++) {
    out_state = 0;
    fprintf(file, "%s", tab);
    if (mo->s[i].out_states > 0 && 
	mo->s[i].out_id[out_state] == 0) {
	fprintf(file, "%.2f", mo->s[i].out_a[k][out_state]);
	out_state++;
    }
    else fprintf(file, "0.00");
    for (j = 1; j < mo->N; j++) {
      if (mo->s[i].out_states > out_state && 
	  mo->s[i].out_id[out_state] == j) {
	fprintf(file, "%s %.2f", separator, mo->s[i].out_a[k][out_state]);
	out_state++;
      }
      else fprintf(file, "%s 0.00", separator);
    }
    fprintf(file, "%s\n", ending);
  }
} /* model_A_print */

/*============================================================================*/

void sdmodel_B_print(FILE *file, sdmodel *mo, char *tab, char *separator, 
		   char *ending) {
  int i, j;
  for (i = 0; i < mo->N; i++) {
    fprintf(file, "%s", tab);
    fprintf(file, "%.2f", mo->s[i].b[0]);
    for (j = 1; j < mo->M; j++)
      fprintf(file, "%s %.2f", separator, mo->s[i].b[j]);
    fprintf(file, "%s\n", ending);
  }
} /* model_B_print */

/*============================================================================*/

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


void model_to_sdmodel(const model *mo, const sdmodel *smo, int klass) {
#define CUR_PROC "model_to_sdmodel"
  int i,j,m, nachf, vorg;

  for (i = 0; i < mo->N; i++) {
    nachf = mo->s[i].out_states;
    vorg  = mo->s[i].in_states;
    /* Copy the values */
    for (j = 0; j < nachf; j++) {
      smo->s[i].out_a[klass][j]  = mo->s[i].out_a[j];
      smo->s[i].out_id[j]        = mo->s[i].out_id[j];
    }
    for (j = 0; j < vorg; j++) {
      smo->s[i].in_a[klass][j]  = mo->s[i].in_a[j];
      smo->s[i].in_id[j]        = mo->s[i].in_id[j];
    }
    for (m = 0; m < mo->M; m++)
      smo->s[i].b[m] = mo->s[i].b[m];
    smo->s[i].pi = mo->s[i].pi;
    smo->s[i].out_states = nachf;
    smo->s[i].in_states = vorg;
  }
  smo->prior = mo->prior;
}


model *sdmodel_to_model(const sdmodel *mo, int kclass) {
#define CUR_PROC "sdmodel_to_model"
  /*
   * Set the pointer appropriately
   */
  int i, j, nachf, vorg, m;
  model *m2 = NULL;
  if(!m_calloc(m2, 1)) {mes_proc(); goto STOP;}
  if (!m_calloc(m2->s, mo->N)) {mes_proc(); goto STOP;}
  for (i = 0; i < mo->N; i++) {
    nachf = mo->s[i].out_states;
    vorg = mo->s[i].in_states;
    if(!m_calloc(m2->s[i].out_id, nachf)) {mes_proc(); goto STOP;}
    if(!m_calloc(m2->s[i].out_a, nachf)) {mes_proc(); goto STOP;}
    if(!m_calloc(m2->s[i].in_id, vorg)) {mes_proc(); goto STOP;}
    if(!m_calloc(m2->s[i].in_a, vorg)) {mes_proc(); goto STOP;}
    if(!m_calloc(m2->s[i].b, mo->M)) {mes_proc(); goto STOP;}
    /* Copy the values */
    for (j = 0; j < nachf; j++) {
      m2->s[i].out_a[j]  = mo->s[i].out_a[kclass][j];
      m2->s[i].out_id[j] = mo->s[i].out_id[j];
    }
    for (j = 0; j < vorg; j++) {
      m2->s[i].in_a[j]  = mo->s[i].in_a[kclass][j];
      m2->s[i].in_id[j] = mo->s[i].in_id[j];
    }
    for (m = 0; m < mo->M; m++)
      m2->s[i].b[m] = mo->s[i].b[m];
    m2->s[i].pi = mo->s[i].pi;
    m2->s[i].out_states = nachf;
    m2->s[i].in_states = vorg;
  }
  m2->N = mo->N;
  m2->M = mo->M;
  m2->prior = mo->prior;

  m2->model_type=mo->model_type;
  if ( mo->model_type == kSilentStates ) {
    assert( mo->silent != NULL );
    if(!m_calloc(m2->silent, mo->N)) {mes_proc(); goto STOP;}
    for(i=0; i < m2->N; i++) {
      m2->silent[i]=mo->silent[i];
    }
    m2->topo_order_length=mo->topo_order_length;
    if(!m_calloc(m2->topo_order, mo->topo_order_length)) {mes_proc(); goto STOP;}
    for(i=0; i < m2->topo_order_length; i++) {
      m2->topo_order[i]=mo->topo_order[i];
    }
  }
  return(m2);
STOP:
  m_free(m2->silent);
  m_free(m2->topo_order);
  model_free(&m2);
  return(NULL);
}

/*============================================================================*/

/*============================================================================*/

/*state* state_copy(state *my_state) {
  state* new_state = (state*) malloc(sizeof(state));

  state_copy_to(my_state,new_state);

  return new_state;

  }*/ /* state_copy */

/*============================================================================*/

/*void state_copy_to(state *source, state* dest) {
  dest->pi         = source->pi;
  dest->out_states = source->out_states;
  dest->in_states  = source->in_states;
  dest->fix        = source->fix;

  dest->b          = malloc(xxx);
  memcpy(dest->b,source->b,xxx);

  dest->out_id     = malloc(xxx);
  memcpy(dest->out_id,source->out_id,xxx);

  dest->in_id      = malloc(xxx);
  memcpy(dest->in_id,source->in_id,xxx);

  dest->out_a      = malloc(xxx);
  memcpy(dest->out_a,source->out_a,xxx);

  dest->in_a       = malloc(xxx);
  memcpy(dest->in_a,source->in_a,xxx);
  }*/ /* state_copy_to */

		  	  
/*===================== E n d   o f  f i l e  "model.c"       ===============*/