sequence.c

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

    char *str = mprintf(NULL, 0,
			"Number of sequences %ld exceeds possible range", 
			seq_number); 
    mes_prot(str);
    m_free(str);
    goto STOP;
  }
  if(!m_calloc(sqd, 1)) {mes_proc(); goto STOP;}
  if(!m_calloc(sqd->seq, seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(sqd->seq_len, seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(sqd->seq_label, seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(sqd->seq_id, seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(sqd->seq_w, seq_number)) {mes_proc(); goto STOP;}
  sqd->seq_number = seq_number;
  
  sqd-> total_w = 0.0;
  for (i = 0; i < seq_number; i++) {
    sqd->seq_label[i] = -1;
    sqd->seq_id[i] = -1.0;
    sqd->seq_w[i] = 1;    
  }
    
  return sqd;
 STOP:
  sequence_d_free(&sqd);
  return NULL;
#undef CUR_PROC
} /* sequence_d_calloc */

/*============================================================================*/
sequence_t *sequence_calloc(long seq_number) {
#define CUR_PROC "sequence_calloc"
  int i;
  sequence_t *sq = NULL;
  if (seq_number > MAX_SEQ_NUMBER) {
    char *str = mprintf(NULL, 0,
			"Number of sequences %ld exceeds possible range", 
			seq_number); 
    mes_prot(str);
    m_free(str);
    goto STOP;
  }
  if(!m_calloc(sq, 1)) {mes_proc(); goto STOP;}
  if(!m_calloc(sq->seq, seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(sq->states, seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(sq->seq_len, seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(sq->seq_label, seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(sq->seq_id, seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(sq->seq_w, seq_number)) {mes_proc(); goto STOP;}
  sq->seq_number = seq_number;
  for (i = 0; i < seq_number; i++) {
    sq->seq_label[i] = -1;
    sq->seq_id[i] = -1.0;
    sq->seq_w[i] = 1;
  }
  return sq;
 STOP:
  sequence_free(&sq);
  return NULL;
#undef CUR_PROC
} /* sequence_calloc */


/*============================================================================*/
sequence_t *sequence_lexWords(int n, int M) {
# define CUR_PROC "sequence_lexWords"
  
  sequence_t *sq = NULL;
  long seq_number, cnt = 0;
  int j = n - 1;
  int i;
  int *seq;
  if ((n < 0) || (M <= 0)) { mes_proc(); goto STOP; }
  seq_number = (long) pow((double) M, (double) n);
  sq = sequence_calloc(seq_number);
  if (!sq) { mes_proc(); goto STOP; }
  for (i = 0; i < seq_number; i++) {
    if (!m_calloc(sq->seq[i], n)) { mes_proc(); goto STOP; }
    sq->seq_len[i] = n;
    sq->seq_id[i] = i;
  }

  if (!m_calloc(seq, n)) { mes_proc(); goto STOP; }
  while (!(j < 0)) {
    sequence_copy(sq->seq[cnt], seq, n);
    j = n - 1;
    while (seq[j] == M - 1) {
      seq[j] = 0;
      j--;
    }
    seq[j]++;
    cnt++;
  }
  m_free(seq);
  return sq;
STOP:
  sequence_free(&sq);
  return NULL;
# undef CUR_PROC
} /* sequence_lewWords */

/*============================================================================*/
int sequence_max_symbol(sequence_t *sq) {
  long i, j;
  int max_symb = -1;
  for (i = 0; i < sq->seq_number; i++)
    for (j = 0; j < sq->seq_len[i]; j++) {
      if (sq->seq[i][j] > max_symb)
	max_symb = sq->seq[i][j];
    }
  return max_symb;
} /* sequence_max_symbol */

/*============================================================================*/
void sequence_copy(int *target, int *source, int len) {
  int i;
  for (i = 0; i < len; i ++)
    target[i] = source[i];
} /* sequence_copy */


/*============================================================================*/
void sequence_d_copy(double *target, double *source, int len) {
  int i;
  for (i = 0; i < len; i ++)
    target[i] = source[i];
} /* sequence_copy */


/*============================================================================*/
int sequence_add(sequence_t *target, sequence_t *source) {
#define CUR_PROC "sequence_add"

  int res = -1;
  int **old_seq       = target->seq;
  int **old_seq_st    = target->states;
  int *old_seq_len    = target->seq_len;
  long *old_seq_label = target->seq_label;
  double *old_seq_id = target->seq_id;
  double *old_seq_w = target->seq_w;
  long old_seq_number = target->seq_number;
  long i;

  target->seq_number = old_seq_number + source->seq_number;
  target->total_w += source->total_w;

  if(!m_calloc(target->seq, target->seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(target->states, target->seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(target->seq_len, target->seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(target->seq_label, target->seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(target->seq_id, target->seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(target->seq_w, target->seq_number)) {mes_proc(); goto STOP;}

  for (i = 0; i < old_seq_number; i++) {
    target->seq[i] = old_seq[i];
    target->states[i] = old_seq_st[i];
    target->seq_len[i] = old_seq_len[i];
    target->seq_label[i] = old_seq_label[i];
    target->seq_id[i] = old_seq_id[i];
    target->seq_w[i] = old_seq_w[i];
  }

  for (i = 0; i < (target->seq_number - old_seq_number); i++) {
    if(!m_calloc(target->seq[i+old_seq_number], source->seq_len[i])) 
      {mes_proc(); goto STOP;}

    sequence_copy(target->seq[i+old_seq_number], source->seq[i], 
		  source->seq_len[i]);

    if(!m_calloc(target->states[i+old_seq_number], source->seq_len[i])) 
      {mes_proc(); goto STOP;}

    /*sequence_copy(target->states[i+old_seq_number], source->states[i], 
      source->seq_len[i]);*/

    target->seq_len[i+old_seq_number] = source->seq_len[i];
    target->seq_label[i+old_seq_number] = source->seq_label[i];
    target->seq_id[i+old_seq_number] = source->seq_id[i];
    target->seq_w[i+old_seq_number] = source->seq_w[i];
  } 
 
  /*m_free(old_seq);
    m_free(old_seq_st);*/
  m_free(old_seq_len);
  m_free(old_seq_label);
  m_free(old_seq_id);
  m_free(old_seq_w);
  res = 0;
 STOP:
  return res;
#undef CUR_PROC
}


/*============================================================================*/
int sequence_d_add(sequence_d_t *target, sequence_d_t *source) {
#define CUR_PROC "sequence_d_add"

  int res = -1;
  double **old_seq    = target->seq;
  int *old_seq_len    = target->seq_len;
  long *old_seq_label = target->seq_label;
  double *old_seq_id  = target->seq_id;
  double *old_seq_w  = target->seq_w;
  long old_seq_number = target->seq_number;  
  long i;

  target->seq_number = old_seq_number + source->seq_number;
  target->total_w += source->total_w;

  if(!m_calloc(target->seq, target->seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(target->seq_len, target->seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(target->seq_label, target->seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(target->seq_id, target->seq_number)) {mes_proc(); goto STOP;}
  if(!m_calloc(target->seq_w, target->seq_number)) {mes_proc(); goto STOP;}

  for (i = 0; i < old_seq_number; i++) {
    target->seq[i] = old_seq[i];
    target->seq_len[i] = old_seq_len[i];
    target->seq_label[i] = old_seq_label[i];
    target->seq_id[i] = old_seq_id[i];
    target->seq_w[i] = old_seq_w[i];
  }

  for (i = 0; i < (target->seq_number - old_seq_number); i++) {
    if(!m_calloc(target->seq[i+old_seq_number], source->seq_len[i])) 
      {mes_proc(); goto STOP;}

    sequence_d_copy(target->seq[i+old_seq_number], source->seq[i], 
		    source->seq_len[i]);
    target->seq_len[i+old_seq_number] = source->seq_len[i];
    target->seq_label[i+old_seq_number] = source->seq_label[i];
    target->seq_id[i+old_seq_number] = source->seq_id[i];
    target->seq_w[i+old_seq_number] = source->seq_w[i];
  } 
 
  m_free(old_seq);
  m_free(old_seq_len);
  m_free(old_seq_label);
  m_free(old_seq_id);
  m_free(old_seq_w);
  res = 0;

 STOP:
  return res;
#undef CUR_PROC
}

/*============================================================================*/
int sequence_check(sequence_t *sq, int max_symb) {
#define CUR_PROC "sequence_check"
  int i, j;
  for (j = 0; j < sq->seq_number; j++) {
    for (i = 0; i < sq->seq_len[j]; i++) {
      if ((sq->seq[j][i] >= max_symb) || (sq->seq[j][i] < 0)) {
	char *str =  
	  mprintf(NULL, 0, "Wrong symbol \'%d\' in sequence %d at Pos. %d;\
                            Should be within [0..%d]\n",
		  sq->seq[j][i], j + 1, i + 1, max_symb - 1);
	mes_prot(str);	
	m_free(str); return (-1);
      }
    }
  }
  return 0;
#undef CUR_PROC
} /* sequence_check */

/*============================================================================*/
int sequence_best_model(model **mo, int model_number, int *sequence, 
			int seq_len, double *log_p) {
# define CUR_PROC "seqence_best_model"
  double log_ptmp;
  int model_index, i;
  *log_p = -DBL_MAX;
  model_index = -1;
  for (i = 0; i < model_number; i++) {    
    foba_logp(mo[i], sequence, seq_len, &log_ptmp);
    if (log_ptmp != +1 && log_ptmp > *log_p) {
      *log_p = log_ptmp;
      model_index = i;
    }
  }
  if (*log_p == -DBL_MAX)
    *log_p = +1;
  return (model_index);
# undef CUR_PROC
} /* sequence_best_model */


/*============================================================================*/
void sequence_print(FILE *file, sequence_t *sq) {
  int i, j;
  fprintf(file, "SEQ = {\n\tO = {\n");
  for (i = 0; i < sq->seq_number; i++) {
    if (sq->seq_id[i] != -1.0)
      fprintf(file, "\t(%10.0f)", sq->seq_id[i]);
    if (sq->seq_label[i] != -1)
      fprintf(file, "\t<%ld>", sq->seq_label[i]);
    if (sq->seq_w[i] != 1)
      fprintf(file, "\t|%.0f|", sq->seq_w[i]);
    fprintf(file,"\t");
    if (sq->seq_len[i] > 0) {
      fprintf(file, "%d", sq->seq[i][0]);
      for (j = 1; j < sq->seq_len[i]; j++)
	fprintf(file, ", %d", sq->seq[i][j]);
    }
    fprintf(file, ";\n");
  }
  fprintf(file, "\t};\n};\n\n");
} /* sequence_print */

/*============================================================================*/
/**/
void sequence_print_xml(FILE *file, sequence_t *sq) {
  int i, j;
  /* coding missing */
  fprintf(file, "<Sequences type=\"int\" >\n");
  fprintf(file, " <DiscretePD>\n");
  for (i = 0; i < sq->seq_number; i++) {
    fprintf(file, "  %.0f <Sequence", sq->seq_w[i]);
    if (sq->seq_id[i] != -1.0)
      fprintf(file, " id=\"seq%f\" ", sq->seq_id[i]);
    fprintf(file,">");
    if (sq->seq_label[i] != -1)
      fprintf(file, "<Label>%ld</Label>", sq->seq_label[i]);
    if (sq->seq_len[i] > 0) {
      fprintf(file, "<!-- Length: %d -->", sq->seq_len[i]);
      for (j = 0; j < sq->seq_len[i]; j++)
	fprintf(file, " %d", sq->seq[i][j]);
    }
    fprintf(file, "  </Sequence>\n");
  }
  fprintf(file, " </DiscretePD>\n");
  fprintf(file, "</Sequences>\n");
} /* sequence_print_xml */

void sequence_d_print_xml(FILE *file, sequence_d_t *sq) {
  int i, j;
  /* coding missing */
  fprintf(file, "<Sequences type=\"int\" >\n");
  fprintf(file, " <DiscretePD>\n");
  for (i = 0; i < sq->seq_number; i++) {
    fprintf(file, "  %.0f <Sequence", sq->seq_w[i]);
    if (sq->seq_id[i] != -1.0)
      fprintf(file, " id=\"seq%f\" ", sq->seq_id[i]);
    fprintf(file,">");
    if (sq->seq_label[i] != -1)
      fprintf(file, "<Label>%ld</Label>", sq->seq_label[i]);
    if (sq->seq_len[i] > 0) {
      fprintf(file, "<!-- Length: %d -->", sq->seq_len[i]);
      for (j = 0; j < sq->seq_len[i]; j++)
	fprintf(file, " %f", sq->seq[i][j]);
    }
    fprintf(file, "  </Sequence>\n");
  }
  fprintf(file, " </DiscretePD>\n");
  fprintf(file, "</Sequences>\n");
} /* sequence_print_xml */

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

void sequence_mathematica_print(FILE *file, sequence_t *sq, char *name) {
  int i;
  fprintf(file, "%s = {\n", name);
  for (i = 0; i < sq->seq_number - 1; i++)
    vector_i_print(file, sq->seq[i], sq->seq_len[i], 
      "{", ",", "},");
  /* no comma after last seq. */
  vector_i_print(file, sq->seq[sq->seq_number - 1], 
		 sq->seq_len[sq->seq_number - 1], 
		 "{", ",", "}"); 
  fprintf(file, "};\n");
} /* sequence_d_mathematica_print */

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

void sequence_d_gnu_print(FILE *file, sequence_d_t *sqd) {
  int i, j;
  for (i = 0; i < sqd->seq_number; i++) {
    for (j = 0; j < sqd->seq_len[i]; j++)
      fprintf(file, "%.8f\n", sqd->seq[i][j]);
    fprintf(file, "\n\n");
  }
}