ghmm_sequences.cpp

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

/** 
    Copies all entries from one sequence in a source array to a target array.
    No memory allocation here.
    @param target double sequence target
    @param source double sequence source
    @param t_num position in target array
    @param s_num position in source array
*/
//void GHMM_Sequences::sequence_d_copy_all(sequence_d_t *target, long t_num, 
//					sequence_d_t *source, long s_num) {
//}

/** Log-Likelihood function in a mixture model:
    (mathe mode?)
    $\sum_k w^k \log( \sum_c (\alpha_c p(O^k | \lambda_c)))$
    @param smo pointer to array of smodels
    @param smo\_number number of models
    @param sqd sequence struct
    @param like log likelihood
*/
//int GHMM_Sequences::sequence_d_mix_like(smodel **smo, int  smo_number, sequence_d_t *sqd, double *like) {
//}


void GHMM_Sequences::clean() {

  if (c_i_sequences)
    sequence_free(&c_i_sequences);
  if (c_d_sequences)
    sequence_d_free(&c_d_sequences);

  c_i_sequences = NULL;
  c_d_sequences = NULL;

  if (own_alphabet)
    SAFE_DELETE(alphabet);
  
  alphabet     = NULL;
  own_alphabet = false;

  clean_cpp();
}


int* GHMM_Sequences::getIntSequence(int index) const {
  if (!c_i_sequences) {
    fprintf(stderr,"GHMM_Sequences::getIntSequence(int) object does not contain int sequences.\n");
    exit(1);
  }

  return c_i_sequences->seq[index];
}


double* GHMM_Sequences::getDoubleSequence(int index) const {
  if (!c_d_sequences) {
    fprintf(stderr,"GHMM_Sequences::getDoubleSequence(int) object does not contain double sequences.\n");
    exit(1);
  }

  return c_d_sequences->seq[index];
}


unsigned int GHMM_Sequences::getLength(int index) const {
  if (c_d_sequences)
    return c_d_sequences->seq_len[index];

  if (c_i_sequences)
    return c_i_sequences->seq_len[index];

  fprintf(stderr,"GHMM_Sequences::getLength(int) object does not contain int sequences.\n");
  exit(1);

  return 0;
}


XMLIO_Element* GHMM_Sequences::XMLIO_startTag(const string& my_tag, XMLIO_Attributes& my_attributes) {
  /* discrete models */
  if (c_i_sequences) {
    if (my_tag == "alphabet") {
      alphabet     = new GHMM_Alphabet();
      own_alphabet = true;
      
      return alphabet;
    }

    if (my_tag == "sequence") {
      GHMM_Sequence* sequence = new GHMM_Sequence(alphabet,0,last_weight);
      sequences.push_back(sequence);
      
      last_weight = 1;
      
      return sequence;
    }
  }

  /* continuous models */
  if (c_d_sequences) {
    if (my_tag == "sequence") {
      GHMM_Sequence* sequence = new GHMM_Sequence(GHMM_DOUBLE,0,last_weight);
      sequences.push_back(sequence);
      
      last_weight = 1;
      
      return sequence;
    }
  }

  fprintf(stderr,"unexpected tag: %s in <sequences>\n",my_tag.c_str());
  exit(1);

  return NULL;
}


void GHMM_Sequences::XMLIO_finishedReading() {
  if (sequence_type == GHMM_INT)
    c_i_sequences = sequence_calloc(0);

  if (sequence_type == GHMM_DOUBLE)
    c_d_sequences = sequence_d_calloc(0);

  unsigned int i;
  for (i = 0; i < sequences.size(); ++i)
    add(sequences[i]);

  clean_cpp();
}


void GHMM_Sequences::clean_cpp() {
  unsigned int i;
  for (i = 0; i < sequences.size(); ++i)
    SAFE_DELETE(sequences[i]);
  sequences.clear();
}


void GHMM_Sequences::read(const string& filename) {
  int seq_number;
  int i;

  /* first clean up old data. */
  clean();

  if (sequence_type == GHMM_DOUBLE) {
    sequence_d_t** seq = sequence_d_read(filename.c_str(),&seq_number);

    if (seq_number > 0)
      copyFromSequences(seq[0]);

    for (i = 0; i < seq_number; ++i)
      sequence_d_free(&seq[i]);

    free(seq);
  }
  else {
    sequence_t** seq = sequence_read(filename.c_str(),&seq_number);

    if (seq_number > 0)
      copyFromSequences(seq[0]);

    for (i = 0; i < seq_number; ++i)
      sequence_free(&seq[i]);

    free(seq);
  }
}


void GHMM_Sequences::copyFromSequences(sequence_d_t* seq) {
  int i;
  int j;

  clean();

  c_d_sequences = sequence_d_calloc(seq->seq_number);

  for (i = 0; i < seq->seq_number; ++i) {
    c_d_sequences->seq[i] = (double*) malloc(sizeof(double) * seq->seq_len[i]);
    for (j = 0; j < seq->seq_len[i]; ++j)
      c_d_sequences->seq[i][j] = seq->seq[i][j];

    c_d_sequences->seq_len[i]   = seq->seq_len[i];
    c_d_sequences->seq_label[i] = seq->seq_label[i];
    c_d_sequences->seq_id[i]    = seq->seq_id[i];
    c_d_sequences->seq_w[i]     = seq->seq_w[i];
  }
  c_d_sequences->seq_number = seq->seq_number;
  c_d_sequences->total_w    = seq->total_w;
}


void GHMM_Sequences::copyFromSequences(sequence_t* seq) {
  int i;
  int j;

  clean();

  c_i_sequences = sequence_calloc(seq->seq_number);

  for (i = 0; i < seq->seq_number; ++i) {
    c_i_sequences->seq[i] = (int*) malloc(sizeof(double) * seq->seq_len[i]);
    for (j = 0; j < seq->seq_len[i]; ++j)
      c_d_sequences->seq[i][j] = seq->seq[i][j];

    c_i_sequences->seq_len[i]   = seq->seq_len[i];
    c_i_sequences->seq_label[i] = seq->seq_label[i];
    c_i_sequences->seq_id[i]    = seq->seq_id[i];
    c_i_sequences->seq_w[i]     = seq->seq_w[i];
  }
  c_i_sequences->seq_number = seq->seq_number;
  c_i_sequences->total_w    = seq->total_w;
}


/** Generates all possible integer sequence of lenght n from an alphabet with
    M letters. Use lexicographical ordering. Memory allocation here.
    @param n     length of sequences
    @param M     size of alphabet
*/
void GHMM_Sequences::lexWords(int n, int M) {
  clean();

  c_i_sequences = sequence_lexWords(n,M);
}


void GHMM_Sequences::init() {
  tag               = "sequences";
  c_i_sequences     = NULL;
  c_d_sequences     = NULL;
  alphabet          = NULL;
  own_alphabet      = false;
  last_weight       = 1;
  xmlio_indent_type = XMLIO_INDENT_BOTH;
}


void GHMM_Sequences::init_INT(GHMM_Alphabet* my_alphabet, sequence_t* my_c_i_sequences) {
  init();

  sequence_type      = GHMM_INT;
  alphabet           = my_alphabet;
  c_i_sequences      = my_c_i_sequences;
  attributes["type"] = "discrete";

  /* always initialize sequence data structure. */
  if (!c_i_sequences)
    c_i_sequences = sequence_calloc(0);
}


void GHMM_Sequences::init_DOUBLE(sequence_d_t* my_c_d_sequences) {
  init();

  sequence_type      = GHMM_DOUBLE;
  c_d_sequences      = my_c_d_sequences;
  attributes["type"] = "continuous";

  /* always initialize sequence data structure. */
  if (!c_d_sequences)
    c_d_sequences = sequence_d_calloc(0);
}


const int GHMM_Sequences::XMLIO_writeContent(XMLIO_Document& writer) {
  unsigned int i;
  int result = 0;

  writer.changeIndent(2);
  result += writer.writeEndl();
  result += writer.writeElement(alphabet);
  result += writer.writeEndl();

  for (i = 0; i < getNumberOfSequences(); ++i) {
    result += writer.writeIndent();
    result += writer.writeStartTag("sequence");
    result += writer.writef("%s",getSequence(i).c_str());
    result += writer.writeEndTag("sequence");
    result += writer.writeEndl();
  }

  return result;
}


unsigned int GHMM_Sequences::getNumberOfSequences() const {
  if (c_i_sequences)
    return c_i_sequences->seq_number;

  if (c_d_sequences)
    return c_d_sequences->seq_number;

  return 0;
}


string GHMM_Sequences::getSequence(int index) const {
  string seq;
  unsigned int i;

  if (c_i_sequences) {
    for (i = 0; i < getLength(index); ++i)
      if (alphabet)
	seq += alphabet->getSymbol(c_i_sequences->seq[index][i]);
      else {
	if (i > 0) seq += " ";
	seq += GHMM_Toolkit::toString(c_i_sequences->seq[index][i]);
      }
  }

  if (c_d_sequences) {
    for (i = 0; i < getLength(index); ++i) {
      if (i > 0) seq += " ";
      seq += GHMM_Toolkit::toString(c_i_sequences->seq[index][i]);
    }
  }
  
  return seq;
}