gql.c

General Hidden Markov Model Library 一个通用的隐马尔科夫模型的C代码库

/*******************************************************************************
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version 0.8_beta1, see http://ghmm.org
*
*       Filename: gql.c
*       Authors:   Wasinee Rungsarityotin, Benjamin Georgi
*
*       Copyright (C) 1998-2004 Alexander Schliep
*       Copyright (C) 1998-2001 ZAIK/ZPR, Universitaet zu Koeln
*       Copyright (C) 2002-2004 Max-Planck-Institut fuer Molekulare Genetik,
*                               Berlin
*
*       Contact: [e-mail]
*
*       This library is free software; you can redistribute it and/or
*       modify it under the terms of the GNU Library General Public
*       License as published by the Free Software Foundation; either
*       version 2 of the License, or (at your option) any later version.
*
*       This library is distributed in the hope that it will be useful,
*       but WITHOUT ANY WARRANTY; without even the implied warranty of
*       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
*       Library General Public License for more details.
*
*       You should have received a copy of the GNU Library General Public
*       License along with this library; if not, write to the Free
*       Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*
*       This file is version $Revision: 1773 $
*                       from $Date: 2006-11-07 13:44:18 +0100 (Tue, 07 Nov 2006) $
*             last change by $Author: grunau $.
*
*******************************************************************************/

/* XXX FIXME: breaks out of tree build of ghmmwrapper */
#include "../config.h"

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <assert.h>

#include <ghmm/rng.h>
#include <ghmm/sequence.h>
#include <ghmm/smodel.h>
#include <ghmm/sfoba.h>
#include <ghmm/obsolete.h>
/*#include <ghmm++/GHMM_convertXMLtoC.h>*/

static double *log_p_pt;

static int seq_rank_fp(const void *a1, const void *a2) {
  
  int i, j;

  i = *((int*)a1);
  j = *((int*)a2);

  if ( log_p_pt[i] == log_p_pt[j] ) return 0;
  else
    if ( log_p_pt[i] < log_p_pt[j] ) return -1;
    else
      return 1;
}


ghmm_cmodel *smodel_alloc_fill(int N, int M, int cos, double prior, int density) {
#define CUR_PROC "smodel_alloc_fill"
  int i;
  ghmm_cmodel *smo=NULL;
  if (!(smo = malloc(sizeof(ghmm_cmodel)))) {goto STOP;}  
  smo->M   = M;
  smo->N   = N;
  smo->cos = cos;
  smo->prior = prior;
  if (!(smo->s = calloc(sizeof(*(smo->s)), (smo->N))))
    {goto STOP;}

  for(i=0; i < smo->N; i++) {
    ghmm_cstate_alloc(&smo->s[i], smo->M, smo->N, smo->N, cos);
  }
  return smo;
STOP:
  fprintf (stderr, "smodel_alloc_fill(1): out of memory\n");
  return NULL;
#undef CUR_PROC
}

void smodel_set_pivector(ghmm_cmodel *smo, int i, double piv) {
  if (smo->s != NULL) {
    smo->s[i].pi = piv;
  }
}

void smodel_set_fixvector(ghmm_cmodel *smo, int i, double fixv) {
  if (smo->s != NULL) {
    smo->s[i].fix = fixv;
  }
}

void smodel_set_transition(ghmm_cmodel *smo, int i, int j, int cos, double prob) {
  int in, out;
  if (cos >= smo->cos) {
    fprintf(stderr, "smodel_set_transition(cos): cos > state->cos\n");
    exit(-1);
  }
  if (smo->s && smo->s[i].out_a && smo->s[j].in_a) {
    for(out=0; out < smo->s[i].out_states; out++) {
      if ( smo->s[i].out_id[out] == j ) {
	smo->s[i].out_a[cos][out] = prob;
	fprintf(stderr, "smodel_set_transition(0):State %d, %d, = %f\n", i, j, prob);
	break;
      }
    }

    for(in=0; in < smo->s[j].in_states; in++) {
      if ( smo->s[j].in_id[in] == i ) {
	smo->s[j].in_a[cos][in] = prob;
	break;
      }
    }
  }
}


double smodel_get_transition(ghmm_cmodel *smo, int i, int j, int cos) {
  int out;
  if (cos >= smo->cos) {
    fprintf(stderr, "smodel_get_transition(0): cos > state->cos\n");
    exit(-1);
  }
  if (smo->s && smo->s[i].out_a && smo->s[j].in_a) {
    for(out=0; out < smo->s[i].out_states; out++) {
      if ( smo->s[i].out_id[out] == j ) {
	return 	smo->s[i].out_a[cos][out];
      }
    }
  }
  fprintf(stderr, "smodel_get_transition(1): data structure not initialized\n");
  return -1.0;
}

void smodel_set_mean(ghmm_cmodel *smo, int i, double *mu) {
  int m;
  if (smo->s != NULL) {
    for(m = 0; m < smo->M; m++)
      smo->s[i].mue[m] = mu[m];
  }
}

void smodel_set_variance(ghmm_cmodel *smo, int i, double *variance) {
  int m;
  if (smo->s != NULL) {
    for(m = 0; m < smo->M; m++) {
      smo->s[i].u[m] = variance[m];
      assert( smo->s[i].u[m] > 0.0 );
    }
  }
}


int smodel_sorted_individual_likelihoods(ghmm_cmodel *smo, ghmm_cseq *sqd, double *log_ps, int *seq_rank) {
  int matched, res;
  double log_p_i;
  int i;
  res=-1;
  matched=0;
  for ( i = 0; i < sqd->seq_number; i++) {
    seq_rank[i] = i;
    if (ghmm_cmodel_logp(smo, sqd->seq[i], sqd->seq_len[i], &log_p_i) != -1) { 
      log_ps[i] = log_p_i;
      matched++;
    }
    else  {
      /* Test: very small log score for sequence cannot be produced */
      log_ps[i] = -DBL_MAX;
    }
  }

  res=matched;
  if (matched == 0) { 
    fprintf(stderr, "smodel_likelihood: NO sequence can be build.\n"); 
  }
  else
    qsort(seq_rank, sqd->seq_number, sizeof(int), (int (*)(const void *, const void *))seq_rank_fp);

  /* return number of "matched" sequences */
  return res;
}