/*******************************************************************************
*
* This file is part of the General Hidden Markov Model Library,
* GHMM version __VERSION__, see http://ghmm.org
*
* Filename: ghmm/ghmm/smixturehmm.c
* Authors: Bernd Wichern
*
* 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: schliep@ghmm.org
*
* 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: 2267 $
* from $Date: 2009-04-24 11:01:58 -0400 (Fri, 24 Apr 2009) $
* last change by $Author: grunau $.
*
*******************************************************************************/
#ifdef HAVE_CONFIG_H
# include "../config.h"
#endif
#include "ghmmconfig.h"
#ifdef GHMM_OBSOLETE
#include <stdio.h>
#include <float.h>
#include <math.h>
#include <time.h>
#include "ghmm.h"
#include "mprintf.h"
#include "mes.h"
#include "matrix.h"
#include "smixturehmm.h"
#include "smodel.h"
#include "sreestimate.h"
#include "sfoba.h"
#include "rng.h"
#include "sequence.h"
#include "smap_classify.h"
#include "ghmm_internals.h"
/* used in main and deactivated other functions */
#if 0
double total_train_w = 0.0;
double total_test_w = 0.0;
#endif
#if 0 /* no main */
/*============================================================================*/
int main (int argc, char *argv[])
{
#define CUR_PROC "smixturehmm_main"
int exitcode = -1;
ghmm_cseq *sqd = NULL, *sqd_train = NULL, *sqd_test = NULL;
ghmm_cseq **sqd_dummy = NULL;
ghmm_cmodel **smo_initial = NULL, **smo = NULL;
int i, k, iter, iterations, field_number, min_smo, max_smo, smo_number,
total_smo_number, print_models, mode, free_parameter;
long errors_train = 0, errors_test = 0;
char *str;
double train_likelihood, test_likelihood, train_ratio, bic;
double *avg_comp_like = NULL;
double **cp = NULL; /* component probabilities for all train-seqs. */
FILE *outfile = NULL, *likefile = NULL, *smofile = NULL;
char filename[256], likefilename[256];
ghmm_rng_init ();
if (argc == 10 || argc == 11) {
sprintf (filename, "%s.smo", argv[3]);
if (!(smofile = ighmm_mes_fopen (filename, "wt"))) {
GHMM_LOG_QUEUED(LCONVERTED);
goto STOP;
}
sprintf (likefilename, "%s.like", argv[3]);
if (!(likefile = ighmm_mes_fopen (likefilename, "at"))) {
GHMM_LOG_QUEUED(LCONVERTED);
goto STOP;
}
sprintf (filename, "%s", argv[3]);
if (!(outfile = ighmm_mes_fopen (filename, "at"))) {
GHMM_LOG_QUEUED(LCONVERTED);
goto STOP;
}
iterations = atoi (argv[4]);
min_smo = atoi (argv[5]);
max_smo = atoi (argv[6]);
printf ("Iter %d, min model no %d, max model no %d\n", iterations,
min_smo, max_smo);
if (max_smo < min_smo) {
printf (" Max. Model No. < Min. Model No.! \n");
goto STOP;
}
train_ratio = atof (argv[7]);
if (train_ratio > 1 || train_ratio < 0) {
printf ("Fehler Train_Ratio (argv[7]): %f\n", train_ratio);
goto STOP;
}
print_models = atoi (argv[8]);
mode = atoi (argv[9]);
if (mode < 1 || mode > 5) {
printf ("Error: initial mode out of range\n");
goto STOP;
}
printf ("TRAIN RATIO %.2f\n", train_ratio);
if (argc == 11)
GHMM_RNG_SET (RNG, atoi (argv[10]));
else {
ghmm_rng_timeseed (RNG);
}
}
else {
printf ("Insufficient arguments. Usage: \n");
printf
("smixturehmm [Seq.File] [Model File] [Out File] [No. of runs] \n");
printf
("[Min Model No.] [Max Model No.] [train ratio] [print models] [mode] <seed>\n");
goto STOP;
}
ghmm_smixturehmm_print_header (likefile, argv, 1);
ghmm_smixturehmm_print_header (outfile, argv, 0);
/* Read Seqs. and Initial Models only once */
sqd_dummy = ghmm_cseq_read (argv[1], &field_number);
printf ("Length first Seq: %d\n", sqd_dummy[0]->seq_len[0]);
if (!sqd_dummy) {
GHMM_LOG_QUEUED(LCONVERTED);
goto STOP;
}
sqd = sqd_dummy[0];
if (field_number > 1)
GHMM_LOG(LCONVERTED,
"Warning: Seq. File contains multiple Seq. Fields; use only the first one\n");
smo_initial = ghmm_cmodel_read (argv[2], &total_smo_number);
if (!smo_initial) {
GHMM_LOG_QUEUED(LCONVERTED);
goto STOP;
}
ARRAY_CALLOC (smo, total_smo_number);
if (total_smo_number < max_smo) {
str = ighmm_mprintf (NULL, 0,
"Need %d initial models, read only %d from model file\n",
max_smo, total_smo_number);
GHMM_LOG(LCONVERTED, str);
m_free (str);
goto STOP;
}
if (smo_initial[0]->density == 0) {
printf ("normal\n");
fprintf (outfile, "normal\n");
}
else {
printf ("normal_pos\n");
fprintf (outfile, "normal_pos\n");
}
/*-------------- Main loop over different no. of models-------------------------*/
for (smo_number = min_smo; smo_number <= max_smo; smo_number++) {
/* different partitions of test und train seqs. */
for (iter = 0; iter < iterations; iter++) {
if (!outfile)
if (!(outfile = ighmm_mes_fopen (filename, "at"))) {
GHMM_LOG_QUEUED(LCONVERTED);
goto STOP;
}
if (!likefile)
if (!(likefile = ighmm_mes_fopen (likefilename, "at"))) {
GHMM_LOG_QUEUED(LCONVERTED);
goto STOP;
}
fprintf (outfile,
"*********** Model No. %d, CV Iter %d ****************\n",
smo_number, iter + 1);
/*==================INIT========================================*/
/* divide into train and test seqs. */
ARRAY_CALLOC (sqd_train, 1);
ARRAY_CALLOC (sqd_test, 1);
if (ghmm_cseq_partition (sqd, sqd_train, sqd_test, train_ratio) == -1) {
str =
ighmm_mprintf (NULL, 0, "Error partitioning seqs, (model %d, iter %d)\n",
smo_number, iter);
GHMM_LOG(LCONVERTED, str);
m_free (str);
if (sqd_train->seq == NULL) {
str =
ighmm_mprintf (NULL, 0,
"Error: empty train seqs, (model %d, iter %d)\n",
smo_number, iter);
GHMM_LOG(LCONVERTED, str);
m_free (str);
}
goto STOP;
}
printf ("%ld seqs divided into %ld train seqs and %ld test seqs\n",
sqd->seq_number, sqd_train->seq_number, sqd_test->seq_number);
/* copy appropriate no. of initial models */
for (k = 0; k < smo_number; k++) {
smo[k] = ghmm_cmodel_copy (smo_initial[k]);
if (!smo[k]) {
str =
ighmm_mprintf (NULL, 0, "Error copying models, (model %d, iter %d)\n",
smo_number, iter);
GHMM_LOG(LCONVERTED, str);
m_free (str);
goto STOP;
}
}
/* find out no. of free parameters; smo_number - 1 = component weights */
free_parameter = ghmm_cmodel_count_free_parameter (smo, smo_number)
+ smo_number - 1;
cp = ighmm_cmatrix_alloc (sqd_train->seq_number, smo_number);
if (!cp) {
GHMM_LOG_QUEUED(LCONVERTED);
goto STOP;
}
/* Initial values for component probs for all seqs. and model priors
before the actual clustering starts */
if (ghmm_smixturehmm_init (cp, sqd_train, smo, smo_number, mode) == -1) {
str =
ighmm_mprintf (NULL, 0,
"Error in initialization comp. prob (model %d, iter %d)\n",
smo_number, iter);
GHMM_LOG(LCONVERTED, str);
m_free (str);
goto STOP;
}
/*==================END INIT=====================================*/
/* clustering */
if (ghmm_smixturehmm_cluster (outfile, cp, sqd_train, smo, smo_number) == -1) {
str = ighmm_mprintf (NULL, 0, "Error in clustering, (model %d, iter %d)\n",
smo_number, iter);
GHMM_LOG(LCONVERTED, str);
m_free (str);
goto STOP;
}
/* TEST */
/*
if (ghmm_smixturehmm_calc_cp(cp, sqd_train, smo, smo_number) == -1) {
str = ighmm_mprintf(NULL, 0, "Error after clustering, (model %d, CV Iter %d)\n",
smo_number, iter + 1);
GHMM_LOG(LCONVERTED, str); m_free(str); goto STOP;
}
fprintf(outfile, "Component Probs. after Clustering:\n");
ighmm_cmatrix_print(outfile, cp, sqd_train->seq_number, smo_number, " ", ",", ";");
*/
/* End TEST */
/*================ OUTPUT=====================================*/
/* print trained models */
if (print_models) {
fprintf (smofile,
"# ****************** Models: %d, CV-Iter %d *****************\n",
smo_number, iter);
for (k = 0; k < smo_number; k++)
ghmm_cmodel_print (smofile, smo[k]);
}
avg_comp_like = ghmm_smixturehmm_avg_like (cp, sqd_train, smo, smo_number);
if (avg_comp_like == NULL) {
GHMM_LOG(LCONVERTED, "Error calculating avg_like \n");
goto STOP;
}
fprintf (outfile, "\nTrain Set:\nModel Priors \t Likel. per Seq\n");
for (k = 0; k < smo_number; k++)
fprintf (outfile, "%.4f\t%.4f\t%.4f\n", smo[k]->prior,
avg_comp_like[k]);
errors_train = ghmm_cseq_mix_like
(smo, smo_number, sqd_train, &train_likelihood);
fprintf (outfile,
"Likelihood Train Set %.4f (%ld seqs); (per seq %.4f)\n",
train_likelihood, sqd_train->seq_number,
train_likelihood / total_train_w);
/* likelihood on test set */
if (sqd_test != NULL) {
errors_test = ghmm_cseq_mix_like
(smo, smo_number, sqd_test, &test_likelihood);
fprintf (outfile,
"Likelihood Test Set %.4f (%ld seqs); (per normseq ",
test_likelihood, sqd_test->seq_number);
if (total_test_w > 0)
fprintf (outfile, " %.4f)\n", test_likelihood / total_test_w);
else
fprintf (outfile, " ---)\n");
}
/* Bayes Information Criterion, see: Kass, Raftery: Bayes Faktors. */
bic =
2 * train_likelihood - free_parameter * log (sqd_train->seq_number);
fprintf (likefile, "%d\t%d\t%ld\t%ld\t%.4f\t%.4f\t%.4f\t", smo_number,
iter, sqd_train->seq_number, sqd_test->seq_number,
train_likelihood, test_likelihood,
train_likelihood / total_train_w);
if (total_test_w > 0)
fprintf (likefile, " %.4f\t", test_likelihood / total_test_w);
else
fprintf (likefile, " ---\t");
fprintf (likefile, "%ld\t%ld\t%.4f\n", errors_train, errors_test, bic);
/*================ FREE=====================================*/
for (k = 0; k < smo_number; k++)
ghmm_cmodel_free (&(smo[k]));
ighmm_cmatrix_free (&cp, sqd_train->seq_number);
ghmm_cseq_free (&sqd_train);
ghmm_cseq_free (&sqd_test);
m_free (avg_comp_like);
if (outfile)
fclose (outfile); /* save results, open to append at the beginning of the loop */
if (likefile)
fclose (likefile);
likefile = NULL;
outfile = NULL;
} /* for (iterations) */
} /* for (smo_number) */
if (smofile)
fclose (smofile);
exitcode = 0;
/*------------------------------------------------------------------------*/
STOP: /* Label STOP from ARRAY_[CM]ALLOC */
ighmm_mes (MES_WIN, "\n(%2.2T): Program finished with exitcode %d.\n", exitcode);
ighmm_mes_exit ();
return (exitcode);
# undef CUR_PROC
} /* main */
#endif /* nomain */
/*============================================================================*/
/* Original version, without attempt to avoid lokal traps */
int ghmm_smixturehmm_cluster (FILE * outfile, double **cp, ghmm_cseq * sqd,
ghmm_cmodel ** smo, int smo_number)
{
#define CUR_PROC "ghmm_smixturehmm_cluster"
int i, k, iter = 0;
double likelihood, old_likelihood, delta_likelihood = 1000000.0;
double total_train_w = 0.0;
char *str;
double *save_w;
double model_weight, log_p, sum = 0.0;
ghmm_cmodel_baum_welch_context *smo_sqd; /* this structure is used by sreestimate() */
ARRAY_CALLOC (smo_sqd, 1);
/* smo_sqd->max_iter = MAX_ITER_BW; */
smo_sqd->max_iter = 10;
smo_sqd->eps = GHMM_EPS_ITER_BW;
smo_sqd->logp = &log_p;
smo_sqd->sqd = sqd;
ARRAY_CALLOC (save_w, sqd->seq_number);
for (i = 0; i < sqd->seq_number; i++) {
save_w[i] = sqd->seq_w[i];
total_train_w += save_w[i];
}
for (k = 0; k < smo_number; k++) {
sum = 0.0;
for (i = 0; i < sqd->seq_number; i++)
sum += cp[i][k] * sqd->seq_w[i];
smo[k]->prior = sum / total_train_w;
}
ghmm_cseq_mix_like (smo, smo_number, sqd, &old_likelihood);
printf ("Initial Likelihood %.4f\n", old_likelihood);
fprintf (outfile, "Initial Likelihood %.4f\n", old_likelihood);
while (((-1) * delta_likelihood / old_likelihood) > 0.001 && iter < 75) {
iter++;
for (k = 0; k < smo_number; k++) {
printf ("Model %d\n", k);
smo_sqd->smo = smo[k];
model_weight = 0.0;
for (i = 0; i < sqd->seq_number; i++) {
/* combine seq_w with appropriate comp. prob. */
sqd->seq_w[i] = save_w[i] * cp[i][k];
model_weight += sqd->seq_w[i];
}
if (ghmm_cmodel_baum_welch (smo_sqd) == -1) {
str = ighmm_mprintf (NULL, 0, "Error iteration %d, model %d\n", iter, k);
GHMM_LOG(LCONVERTED, str);
m_free (str);
goto STOP;
}
smo[k]->prior = model_weight / total_train_w;
} /* for (k ...) */
/* reset weigths for smixturehmm_like */
for (i = 0; i < sqd->seq_number; i++)
sqd->seq_w[i] = save_w[i];
ghmm_cseq_mix_like (smo, smo_number, sqd, &likelihood);
if (ghmm_smixturehmm_calc_cp (cp, sqd, smo, smo_number, &total_train_w) == -1) {
str = ighmm_mprintf (NULL, 0, "Error iteration %d\n", iter);
GHMM_LOG(LCONVERTED, str);
m_free (str);
goto STOP;
}
printf ("Iter %d, likelihood: %.4f\n", iter, likelihood);
fprintf (outfile, "BIter %d, likelihood: %.4f\n", iter, likelihood);
delta_likelihood = likelihood - old_likelihood;
old_likelihood = likelihood;
} /* while delta_likelihood */
m_free (smo_sqd);
m_free (save_w);
return 0;
STOP: /* Label STOP from ARRAY_[CM]ALLOC */
m_free (smo_sqd);
m_free (save_w);
return -1;
#undef CUR_PROC
} /* ghmm_smixturehmm_cluster */
/*============================================================================*/
/* Initial component probs. (cp) for each sequence;
dimension cp: cp[sqd->seq_number] [smo_number]
and (depending on cp): model priors for all models
*/
int ghmm_smixturehmm_init (double **cp, ghmm_cseq * sqd, ghmm_cmodel ** smo,
int smo_number, int mode)
{
#define CUR_PROC "ghmm_smixturehmm_init"
int i, j;
double p;
double *result;
char *str;
int bm;
for (i = 0; i < sqd->seq_number; i++)
for (j = 0; j < smo_number; j++)
cp[i][j] = 0.0;
if (mode < 1 || mode > 5) {
GHMM_LOG(LCONVERTED, "Error: initial mode out of range\n");
goto STOP;
}
/* 1. strict random partition; cp = 0/1 */
if (mode == 1) {
for (i = 0; i < sqd->seq_number; i++) {
p = GHMM_RNG_UNIFORM (RNG);
j = (int) floor (smo_number * p); /* ??? */
if (j < 0 || j >= smo_number) {
GHMM_LOG(LCONVERTED, "Error: initial model out of range\n");
goto STOP;
}
cp[i][j] = 1.0;
}
}
/* 2. ghmm_smap_bayes from initial models */
else if (mode == 2) {
for (i = 0; i < sqd->seq_number; i++)
if (ghmm_smap_bayes (smo, cp[i], smo_number, sqd->seq[i], sqd->seq_len[i]) ==
-1) {
str =
ighmm_mprintf (NULL, 0, "Can't determine comp. prob for seq ID %.0f \n",
sqd->seq_id[i]);
GHMM_LOG(LCONVERTED, str);
m_free (str); /* goto STOP; */
}
}
/* another possibility: make partition with best model from smap_bayes... */
/* 3. cp = 1 for best model, cp = 0 for other models */
else if (mode == 3) {
ARRAY_CALLOC (result, smo_number);
for (i = 0; i < sqd->seq_number; i++) {
bm = ghmm_smap_bayes (smo, result, smo_number, sqd->seq[i], sqd->seq_len[i]);
if (bm == -1) {
str = ighmm_mprintf (NULL, 0,
"Can't determine comp. prob for seq ID %.0f \n",
sqd->seq_id[i]);
GHMM_LOG(LCONVERTED, str);
m_free (str); /* goto STOP; */
}
cp[i][bm] = 1.0;
}
m_free (result);
}
/* mode == 4 used to be kmeans labels */
/* 5. no start partition == equal cp for each model */
else if (mode == 5) {
for (i = 0; i < sqd->seq_number; i++)
for (j = 0; j < smo_number; j++)
cp[i][j] = 1 / (double) smo_number;
}
else {
printf ("Unknown Init Mode %d \n", mode);
return -1;
}
return 0;
STOP: /* Label STOP from ARRAY_[CM]ALLOC */
return -1;
#undef CUR_PROC
} /* smixturehmm_compprob_init */
/*============================================================================*/
/* currently not activated */
# if 0
int smixturehmm_calc_priors (double **cp, ghmm_cseq * sqd, ghmm_cmodel ** smo,
int smo_number)
{
#define CUR_PROC "smixturehmm_calc_priors"
int i, k;
double sum;
if (total_train_w == 0) {
GHMM_LOG(LCONVERTED, "total_train_w == 0!\n");
goto STOP;
}
for (k = 0; k < smo_number; k++) {
sum = 0.0;
for (i = 0; i < sqd->seq_number; i++)
sum += cp[i][k] * sqd->seq_w[i];
smo[k]->prior = sum / total_train_w;
printf ("\nPriors[%d] is %.6f\n", k, smo[k]->prior);
}
return 0;
STOP: /* Label STOP from ARRAY_[CM]ALLOC */
return -1;
#undef CUR_PROC
} /* smixturehmm_calc_priors */
#endif
/*============================================================================*/
/* calculate component probs for all sequences and all components */
/* also recalculate total_train_w; if a seq has cp = 0 don't add its weigth to
total_train_w, otherwise errors in calculating model priors occur
*/
int ghmm_smixturehmm_calc_cp (double **cp, ghmm_cseq * sqd, ghmm_cmodel ** smo,
int smo_number, double *total_train_w)
{
#define CUR_PROC "ghmm_smixturehmm_calc_cp"
int i;
char *str;
double errorseqs = 0.0;
*total_train_w = 0.0;
for (i = 0; i < sqd->seq_number; i++)
if (ghmm_smap_bayes (smo, cp[i], smo_number, sqd->seq[i], sqd->seq_len[i]) ==
-1) {
/* all cp[i] [ . ] are set to zero; seq. will be ignored for reestimation!!! */
str = ighmm_mprintf (NULL, 0,
"Warning[%d]: Can't determine comp. prob for seq ID %.0f\n",
i, sqd->seq_id[i]);
GHMM_LOG(LCONVERTED, str);
m_free (str);
errorseqs++;
if (errorseqs > 0.1 * (double) sqd->seq_number) {
printf ("errorseqs %.1f, max false %.1f\n", errorseqs,
0.1 * (double) sqd->seq_number);
GHMM_LOG(LCONVERTED, "max. no of errors from ghmm_smap_bayes exceeded\n");
goto STOP;
}
}
else
*total_train_w += sqd->seq_w[i];
return 0;
STOP: /* Label STOP from ARRAY_[CM]ALLOC */
return -1;
#undef CUR_PROC
} /* ghmm_smixturehmm_calc_cp */
/*============================================================================*/
/* Is currently not used. Use this function in calc_cp and smixturehmm_like
later on (saves half of the ghmm_cmodel_logp() calls).
Danger: Is it neccessary to take seq_w into account?
*/
void ghmm_smixture_calc_logp (double **logp, int **error, ghmm_cseq * sqd,
ghmm_cmodel ** smo, int smo_number)
{
int i, k;
for (i = 0; i < sqd->seq_number; i++)
for (k = 0; k < smo_number; k++) {
if (ghmm_cmodel_logp (smo[k], sqd->seq[i], sqd->seq_len[i], &(logp[i][k])) ==
-1)
error[i][k] = 1;
else
error[i][k] = 0;
}
}
/*============================================================================*/
/* flag == 1 --> Header for .like-file, else --> other file */
void ghmm_smixturehmm_print_header (FILE * file, char *argv[], int flag)
{
time_t zeit;
int mode = atoi (argv[9]);
time (&zeit);
fprintf (file,
"\n************************************************************************\n");
fprintf (file, "Date: %ssmixturehmm:\n", ctime ((const time_t *) &zeit));
fprintf (file, "Seq. File\t%s\nInit-model File\t%s\nInit-Mode\t", argv[1],
argv[2]);
switch (mode) {
case 1:
fprintf (file, "%d SP_ZUF (random start partition)\n", mode);
break;
case 2:
fprintf (file, "%d SP_VERT (distr. accord smap_bayes)\n", mode);
break;
case 3:
fprintf (file, "%d SP_BEST (best model)\n", mode);
break;
case 4:
fprintf (file, "%d SP_KM (partition from k-means)\n", mode);
break;
case 5:
fprintf (file, "%d NO_SP (no start partition)\n", mode);
break;
default:
fprintf (file, "???\n");
break;
}
fprintf (file, "Train Ratio\t %.4f\n\n", atof (argv[7]));
if (flag == 1) /* only in .like-File */
fprintf (file,
"smo no.\tCV Iter\t SeqTrain\tSeqTest\tLikeTrain\tLikeTest\tavrgTrain\tavrgTest\tErrorTrain\tErrorTest\tBIC\n");
}
/*============================================================================*/
/* calculates for each model the average likelihood per seq. weight:
avg_like[j] = sum_i (cp[i][j] * seq_w[i] * log_p(i | j) ) / sum_i (cp[i][j] * seq_w[i]).
Usefull for identifying models with poor likelihood
*/
double *ghmm_smixturehmm_avg_like (double **cp, ghmm_cseq * sqd,
ghmm_cmodel ** smo, int smo_number)
{
#define CUR_PROC "ghmm_smixturehmm_avg_like"
double *avg_like = NULL;
int i, k;
double num = 0.0, denom = 0.0, log_p = 0.0;
ARRAY_CALLOC (avg_like, smo_number);
for (k = 0; k < smo_number; k++) {
num = denom = 0.0;
for (i = 0; i < sqd->seq_number; i++) {
if (ghmm_cmodel_logp (smo[k], sqd->seq[i], sqd->seq_len[i], &log_p) != -1) {
num += cp[i][k] * sqd->seq_w[i] * log_p;
denom += cp[i][k] * sqd->seq_w[i];
}
}
if (denom > 0)
avg_like[k] = num / denom;
else
avg_like[k] = -1;
} /* for models k ... */
return avg_like;
STOP: /* Label STOP from ARRAY_[CM]ALLOC */
return NULL;
#undef CUR_PROC
} /* ghmm_smixturehmm_avg_like */
#endif /* GHMM_OBSOLETE */