sequence.c

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

  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];
#ifdef GHMM_OBSOLETE
    target->seq_label[i] = old_seq_label[i];
#endif /* GHMM_OBSOLETE */
    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++) {
    ARRAY_CALLOC (target->seq[i + old_seq_number], source->seq_len[i]);

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

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

    /* ghmm_dseq_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];
#ifdef GHMM_OBSOLETE
    target->seq_label[i + old_seq_number] = source->seq_label[i];
#endif /* GHMM_OBSOLETE */
    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);
#ifdef GHMM_OBSOLETE
  m_free (old_seq_label);
#endif /* GHMM_OBSOLETE */
  m_free (old_seq_id);
  m_free (old_seq_w);
  res = 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return res;
#undef CUR_PROC
}


/*============================================================================*/
int ghmm_cseq_add (ghmm_cseq * target, ghmm_cseq * source)
{
#define CUR_PROC "ghmm_cseq_add"

  int res = -1;
  double **old_seq = target->seq;
  int *old_seq_len = target->seq_len;
#ifdef GHMM_OBSOLETE
  long *old_seq_label = target->seq_label;
#endif /* GHMM_OBSOLETE */
  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;

  ARRAY_CALLOC (target->seq, target->seq_number);
  ARRAY_CALLOC (target->seq_len, target->seq_number);
#ifdef GHMM_OBSOLETE
  ARRAY_CALLOC (target->seq_label, target->seq_number);
#endif /* GHMM_OBSOLETE */
  ARRAY_CALLOC (target->seq_id, target->seq_number);
  ARRAY_CALLOC (target->seq_w, target->seq_number);

  for (i = 0; i < old_seq_number; i++) {
    target->seq[i] = old_seq[i];
    target->seq_len[i] = old_seq_len[i];
#ifdef GHMM_OBSOLETE
    target->seq_label[i] = old_seq_label[i];
#endif /* GHMM_OBSOLETE */
    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++) {
    ARRAY_CALLOC (target->seq[i + old_seq_number], source->seq_len[i]);

    ghmm_cseq_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];
#ifdef GHMM_OBSOLETE
    target->seq_label[i + old_seq_number] = source->seq_label[i];
#endif /* GHMM_OBSOLETE */
    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);
#ifdef GHMM_OBSOLETE
  m_free (old_seq_label);
#endif /* GHMM_OBSOLETE */
  m_free (old_seq_id);
  m_free (old_seq_w);
  res = 0;

STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return res;
#undef CUR_PROC
}

/*============================================================================*/
int ghmm_dseq_check (ghmm_dseq * sq, int max_symb)
{
#define CUR_PROC "ghmm_dseq_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 =
          ighmm_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);
        GHMM_LOG(LCONVERTED, str);
        m_free (str);
        return (-1);
      }
    }
  }
  return 0;
#undef CUR_PROC
}                               /* ghmm_dseq_check */

/*============================================================================*/
int ghmm_dseq_best_model (ghmm_dmodel ** 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++) {
    ghmm_dmodel_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
}                               /* ghmm_dseq_best_model */


/*============================================================================*/
void ghmm_dseq_print (ghmm_dseq * sq, FILE * file)
{
  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]);
#ifdef GHMM_OBSOLETE
    if (sq->seq_label[i] != -1)
      fprintf (file, "\t<%ld>", sq->seq_label[i]);
#endif /* GHMM_OBSOLETE */
    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");
}                               /* ghmm_dseq_print */

/*============================================================================*/
 /**/ void ghmm_dseq_print_xml (ghmm_dseq * sq, FILE * file)
{
  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, ">");
#ifdef GHMM_OBSOLETE
    if (sq->seq_label[i] != -1)
      fprintf (file, "<Label>%ld</Label>", sq->seq_label[i]);
#endif /* GHMM_OBSOLETE */
    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");
}                               /* ghmm_dseq_print_xml */

void ghmm_cseq_print_xml (FILE * file, ghmm_cseq * 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, ">");
#ifdef GHMM_OBSOLETE
    if (sq->seq_label[i] != -1)
      fprintf (file, "<Label>%ld</Label>", sq->seq_label[i]);
#endif /* GHMM_OBSOLETE */
    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");
}                               /* ghmm_dseq_print_xml */

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

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

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

void ghmm_cseq_gnu_print (ghmm_cseq * sqd, FILE * file)
{
  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");
  }
}

/*============================================================================*/
void ghmm_cseq_print (ghmm_cseq * sqd, FILE * file, int discrete)
{
  int i, j;
  fprintf (file, "SEQD = {\n\tO = {\n");
  for (i = 0; i < sqd->seq_number; i++) {
    if (sqd->seq_id[i] != -1.0)
      fprintf (file, "\t(%10.0f)", sqd->seq_id[i]);
#ifdef GHMM_OBSOLETE
    if (sqd->seq_label[i] != -1)
      fprintf (file, "\t<%ld>", sqd->seq_label[i]);
#endif /* GHMM_OBSOLETE */
    if (sqd->seq_w[i] != 1)
      fprintf (file, "\t|%.0f|", sqd->seq_w[i]);
    fprintf (file, "\t");
    if (sqd->seq_len[i] > 0) {
      if (discrete)
        fprintf (file, "%3.0f", sqd->seq[i][0]);
      else {
        if (sqd->seq[i][0] > 500)
          fprintf (file, "%8.0f", sqd->seq[i][0]);
        else
          fprintf (file, "%8.2f", sqd->seq[i][0]);
      }
      for (j = 1; j < sqd->seq_len[i]; j++) {
        if (discrete)
          fprintf (file, ", %3.0f", sqd->seq[i][j]);
        else {
          if (sqd->seq[i][j] > 500)
            fprintf (file, ", %8.0f", sqd->seq[i][j]);
          else
            fprintf (file, ", %8.2f", sqd->seq[i][j]);
        }
      }
    }
    fprintf (file, ";\n");
  }
  fprintf (file, "\t};\n};\n\n");
}                               /* ghmm_dseq_print */

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

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

/*============================================================================*/
void ghmm_dseq_clean (ghmm_dseq * sq)
{
#define CUR_PROC "ghmm_dseq_clean"
  /* keep data, only delete references */
  m_free (sq->seq);
  m_free (sq->seq_len);
#ifdef GHMM_OBSOLETE
  m_free (sq->seq_label);
#endif /* GHMM_OBSOLETE */
  m_free (sq->seq_id);
  m_free (sq->seq_w);
  /* m_free(sq->states);*/
  sq->seq_number = 0;
  sq->total_w = 0.0;
#undef CUR_PROC
}                               /* ghmm_dseq_clean */

/*============================================================================*/
void ghmm_cseq_clean (ghmm_cseq * sqd)
{
#define CUR_PROC "ghmm_cseq_clean"
  /* keep data, only delete references */
  m_free (sqd->seq);
  m_free (sqd->seq_len);
#ifdef GHMM_OBSOLETE
  m_free (sqd->seq_label);
#endif /* GHMM_OBSOLETE */
  m_free (sqd->seq_id);
  m_free (sqd->seq_w);
  sqd->seq_number = 0;
  sqd->total_w = 0.0;
#undef CUR_PROC
}                               /* ghmm_cseq_clean */

/*============================================================================*/
int ghmm_dseq_free (ghmm_dseq ** sq) {
# define CUR_PROC "ghmm_dseq_free"
  mes_check_ptr(sq, return -1);
  if (!*sq)
    return -1;

  /* ighmm_dmatrix_free also takes care of (*sq)->seq */
  if ((*sq)->flags & kBlockAllocation)
    free((*sq)->seq);
  else if (ighmm_dmatrix_free(&(*sq)->seq, (*sq)->seq_number) == -1)
    GHMM_LOG(LWARN, "Error in ghmm_dseq_free!");

  m_free((*sq)->seq_len);
#ifdef GHMM_OBSOLETE
  m_free((*sq)->seq_label);
#endif /* GHMM_OBSOLETE */
  m_free((*sq)->seq_id);
  m_free((*sq)->seq_w);

  if ((*sq)->states)
    ighmm_dmatrix_free(&(*sq)->states, (*sq)->seq_number);

  if ((*sq)->state_labels) {
    ighmm_dmatrix_free (&(*sq)->state_labels, (*sq)->seq_number);
    m_free((*sq)->state_labels_len);
  }

  m_free(*sq);
  return 0;
# undef CUR_PROC
}                               /* ghmm_dseq_free */


/*============================================================================*/
int ghmm_cseq_free (ghmm_cseq ** sqd)
{
# define CUR_PROC "ghmm_cseq_free"
  mes_check_ptr(sqd, return -1);
  if (!*sqd)
    return -1;

  ighmm_cmatrix_free(&(*sqd)->seq, (*sqd)->seq_number);
  m_free((*sqd)->seq_len);
#ifdef GHMM_OBSOLETE
  m_free((*sqd)->seq_label);
#endif /* GHMM_OBSOLETE */
  m_free((*sqd)->seq_id);
  m_free((*sqd)->seq_w);
  m_free(*sqd);
  return 0;
# undef CUR_PROC
}                               /* ghmm_cseq_free */

/*============================================================================*/
ghmm_cseq *ghmm_cseq_create_from_dseq (const ghmm_dseq * sq)
{
# define CUR_PROC "ghmm_cseq_create_from_dseq"
  int j, i;
  ghmm_cseq *sqd = NULL;     /* target seq. array */
  if (!(sqd = ghmm_cseq_calloc (sq->seq_number))) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  for (j = 0; j < sq->seq_number; j++) {
    ARRAY_CALLOC (sqd->seq[j], sq->seq_len[j]);
    for (i = 0; i < sq->seq_len[j]; i++)
      sqd->seq[j][i] = (double) (sq->seq[j][i]);
    sqd->seq_len[j] = sq->seq_len[j];
#ifdef GHMM_OBSOLETE