482-sphinx3 Report

**482-sphinx3**
Exec Time (%)	Codelet Name	Error (%)
15.9	subvq_subvq_init_270	100.0
1.0	subvq_subvq_maha_precomp_107	0.4
24.0	mdef_mdef_init_689	100.0
1.5	dict2pid_dict2pid_build_297	1.6
1.9	mdef_sseq_compress_537	5.3
13.0	cont_mgau_mgau_precomp_476	2.2

Exec Time : 1.834394e+09 cycles

Matching : 14.9%

                /* ====================================================================
 * Copyright (c) 1999-2001 Carnegie Mellon University.  All rights
 * reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer. 
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 * This work was supported in part by funding from the Defense Advanced 
 * Research Projects Agency and the National Science Foundation of the 
 * United States of America, and the CMU Sphinx Speech Consortium.
 *
 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * ====================================================================
 *
 */
/*
 * dict2pid.c -- Triphones for dictionary
 * 
 * **********************************************
 * CMU ARPA Speech Project
 *
 * Copyright (c) 1999 Carnegie Mellon University.
 * ALL RIGHTS RESERVED.
 * **********************************************
 * 
 * HISTORY
 * 
 * 14-Sep-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Added dict2pid_comsseq2sen_active().
 * 
 * 04-May-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Started.
 */


#include "dict2pid.h"
#include "logs3.h"


/*
 * Build a glist of triphone senone-sequence IDs (ssids) derivable from [b][r] at the word
 * begin position.  If no triphone found in mdef, include the ssid for basephone b.
 * Return the generated glist.
 */
static glist_t ldiph_comsseq (mdef_t *mdef, int32 b, int32 r)
{
    int32 l, p, ssid;
    glist_t g;
    
    g = NULL;
    for (l = 0; l < mdef_n_ciphone(mdef); l++) {
	p = mdef_phone_id (mdef, (s3cipid_t)b, (s3cipid_t)l, (s3cipid_t)r, WORD_POSN_BEGIN);
	
	if (IS_S3PID(p)) {
	    ssid = mdef_pid2ssid(mdef, p);
	    if (! glist_chkdup_int32 (g, ssid))
		g = glist_add_int32 (g, ssid);
	}
    }
    if (! g)
	g = glist_add_int32 (g, mdef_pid2ssid(mdef, b));
    
    return g;
}


/*
 * Build a glist of triphone senone-sequence IDs (ssids) derivable from [r][b] at the word
 * end position.  If no triphone found in mdef, include the ssid for basephone b.
 * Return the generated glist.
 */
static glist_t rdiph_comsseq (mdef_t *mdef, int32 b, int32 l)
{
    int32 r, p, ssid;
    glist_t g;
    
    g = NULL;
    for (r = 0; r < mdef_n_ciphone(mdef); r++) {
	p = mdef_phone_id (mdef, (s3cipid_t)b, (s3cipid_t)l, (s3cipid_t)r, WORD_POSN_END);
	
	if (IS_S3PID(p)) {
	    ssid = mdef_pid2ssid(mdef, p);
	    if (! glist_chkdup_int32 (g, ssid))
		g = glist_add_int32 (g, ssid);
	}
    }
    if (! g)
	g = glist_add_int32 (g, mdef_pid2ssid(mdef, b));
    
    return g;
}


/*
 * Build a glist of triphone senone-sequence IDs (ssids) derivable from [b] as a single
 * phone word.  If no triphone found in mdef, include the ssid for basephone b.
 * Return the generated glist.
 */
static glist_t single_comsseq (mdef_t *mdef, int32 b)
{
    int32 l, r, p, ssid;
    glist_t g;
    
    g = NULL;
    for (l = 0; l < mdef_n_ciphone(mdef); l++) {
	for (r = 0; r < mdef_n_ciphone(mdef); r++) {
	    p = mdef_phone_id (mdef, (s3cipid_t)b, (s3cipid_t)l, (s3cipid_t)r, WORD_POSN_SINGLE);
	    
	    if (IS_S3PID(p)) {
		ssid = mdef_pid2ssid(mdef, p);
		if (! glist_chkdup_int32 (g, ssid))
		    g = glist_add_int32 (g, ssid);
	    }
	}
    }
    if (! g)
	g = glist_add_int32 (g, mdef_pid2ssid(mdef, b));
    
    return g;
}


/*
 * Build a glist of triphone senone-sequence IDs (ssids) derivable from [b] as a single
 * phone word, with a given left context l.  If no triphone found in mdef, include the ssid
 * for basephone b.  Return the generated glist.
 */
static glist_t single_lc_comsseq (mdef_t *mdef, int32 b, int32 l)
{
    int32 r, p, ssid;
    glist_t g;
    
    g = NULL;
    for (r = 0; r < mdef_n_ciphone(mdef); r++) {
	p = mdef_phone_id (mdef, (s3cipid_t)b, (s3cipid_t)l, (s3cipid_t)r, WORD_POSN_SINGLE);
	
	if (IS_S3PID(p)) {
	    ssid = mdef_pid2ssid(mdef, p);
	    if (! glist_chkdup_int32 (g, ssid))
		g = glist_add_int32 (g, ssid);
	}
    }
    if (! g)
	g = glist_add_int32 (g, mdef_pid2ssid(mdef, b));
    
    return g;
}


/*
 * Convert the glist of ssids to a composite sseq id.  Return the composite ID.
 */
static s3ssid_t ssidlist2comsseq (glist_t g, mdef_t *mdef, dict2pid_t *dict2pid,
				  hash_table_t *hs,	/* For composite states */
				  hash_table_t *hp)	/* For composite senone seq */
{
    int32 i, j, n, s, ssid;
    s3senid_t **sen;
    s3senid_t *comsenid;
    gnode_t *gn;
    
    n = glist_count (g);
    if (n <= 0)
	E_FATAL("Panic: length(ssidlist)= %d\n", n);
    
    /* Space for list of senones for each state, derived from the given glist */
    sen = (s3senid_t **) ckd_calloc (mdef_n_emit_state (mdef), sizeof(s3senid_t *));
    for (i = 0; i < mdef_n_emit_state (mdef); i++) {
	sen[i] = (s3senid_t *) ckd_calloc (n+1, sizeof(s3senid_t));
	sen[i][0] = BAD_S3SENID;	/* Sentinel */
    }
    /* Space for composite senone ID for each state position */
    comsenid = (s3senid_t *) ckd_calloc (mdef_n_emit_state (mdef), sizeof(s3senid_t));
    
    for (gn = g; gn; gn = gnode_next(gn)) {
	ssid = gnode_int32 (gn);
	
	/* Expand ssid into individual states (senones); insert in sen[][] if not present */
	for (i = 0; i < mdef_n_emit_state (mdef); i++) {
	    s = mdef->sseq[ssid][i];
	    
	    for (j = 0; (IS_S3SENID(sen[i][j])) && (sen[i][j] != s); j++);
	    if (NOT_S3SENID(sen[i][j])) {
		sen[i][j] = s;
		sen[i][j+1] = BAD_S3SENID;
	    }
	}
    }
    
    /* Convert senones list for each state position into composite state */
    for (i = 0; i < mdef_n_emit_state (mdef); i++) {
	for (j = 0; IS_S3SENID(sen[i][j]); j++);
	assert (j > 0);
	
	j = hash_enter_bkey (hs, (char *)(sen[i]), j*sizeof(s3senid_t), dict2pid->n_comstate);
	if (j == dict2pid->n_comstate)
	    dict2pid->n_comstate++;	/* New composite state */
	else
	    ckd_free ((void *) sen[i]);
	
	comsenid[i] = j;
    }
    ckd_free (sen);
    
    /* Convert sequence of composite senids to composite sseq ID */
    j = hash_enter_bkey (hp, (char *)comsenid, mdef->n_emit_state * sizeof(s3senid_t),
			 dict2pid->n_comsseq);
    if (j == dict2pid->n_comsseq) {
	dict2pid->n_comsseq++;
	if (dict2pid->n_comsseq >= MAX_S3SENID)
	    E_FATAL("#Composite sseq limit(%d) reached; increase MAX_S3SENID\n",
		    dict2pid->n_comsseq);
    } else
	ckd_free ((void *) comsenid);
    
    return ((s3ssid_t)j);
}


/* RAH 4.16.01 This code has several leaks that must be fixed */
dict2pid_t *dict2pid_build (mdef_t *mdef, dict_t *dict)
{
    dict2pid_t *dict2pid;
    s3ssid_t *internal, **ldiph, **rdiph, *single;
    int32 pronlen;
    hash_table_t *hs, *hp;
    glist_t g;
    gnode_t *gn;
    s3senid_t *sen;
    hash_entry_t *he;
    int32 *cslen;
    int32 i, j, b, l, r, w, n, p;
    
    E_INFO("Building PID tables for dictionary\n");

    dict2pid = (dict2pid_t *) ckd_calloc (1, sizeof(dict2pid_t));
    dict2pid->internal = (s3ssid_t **) ckd_calloc (dict_size(dict), sizeof(s3ssid_t *));
    dict2pid->ldiph_lc = (s3ssid_t ***) ckd_calloc_3d (mdef->n_ciphone,
						       mdef->n_ciphone,
						       mdef->n_ciphone,
						       sizeof(s3ssid_t));
    dict2pid->single_lc = (s3ssid_t **) ckd_calloc_2d (mdef->n_ciphone,
						       mdef->n_ciphone,
						       sizeof(s3ssid_t));
    dict2pid->n_comstate = 0;
    dict2pid->n_comsseq = 0;
    
    hs = hash_new (mdef->n_ciphone * mdef->n_ciphone * mdef->n_emit_state, HASH_CASE_YES);
    hp = hash_new (mdef->n_ciphone * mdef->n_ciphone, HASH_CASE_YES);
    
    for (w = 0, n = 0; w < dict_size(dict); w++) {
	pronlen = dict_pronlen(dict, w);
	if (pronlen < 0)
	    E_FATAL("Pronunciation-length(%s)= %d\n", dict_wordstr(dict, w), pronlen);
	n += pronlen;
    }

    internal = (s3ssid_t *) ckd_calloc (n, sizeof(s3ssid_t));
    
    /* Temporary */
    ldiph = (s3ssid_t **) ckd_calloc_2d (mdef->n_ciphone, mdef->n_ciphone, sizeof(s3ssid_t));
    rdiph = (s3ssid_t **) ckd_calloc_2d (mdef->n_ciphone, mdef->n_ciphone, sizeof(s3ssid_t));
    single = (s3ssid_t *) ckd_calloc (mdef->n_ciphone, sizeof(s3ssid_t));
    for (b = 0; b < mdef->n_ciphone; b++) {
	for (l = 0; l < mdef->n_ciphone; l++) {
	    for (r = 0; r < mdef->n_ciphone; r++)
		dict2pid->ldiph_lc[b][r][l] = BAD_S3SSID;
	    
	    dict2pid->single_lc[b][l] = BAD_S3SSID;
	    
	    ldiph[b][l] = BAD_S3SSID;
	    rdiph[b][l] = BAD_S3SSID;
	}
	single[b] = BAD_S3SSID;
    }
    
    for (w = 0; w < dict_size(dict); w++) {
	dict2pid->internal[w] = internal;
	pronlen = dict_pronlen(dict,w);
	
	if (pronlen >= 2) {
	    b = dict_pron(dict, w, 0);
	    r = dict_pron(dict, w, 1);
	    if (NOT_S3SSID(ldiph[b][r])) {
		g = ldiph_comsseq(mdef, b, r);
		ldiph[b][r] = ssidlist2comsseq (g, mdef, dict2pid, hs, hp);
		glist_free (g);
		
		for (l = 0; l < mdef_n_ciphone(mdef); l++) {
		    p = mdef_phone_id_nearest (mdef, (s3cipid_t)b, (s3cipid_t)l, (s3cipid_t)r, WORD_POSN_BEGIN);
		    dict2pid->ldiph_lc[b][r][l] = mdef_pid2ssid(mdef, p);
		}
	    }
	    internal[0] = ldiph[b][r];
	    
	    for (i = 1; i < pronlen-1; i++) {
		l = b;
		b = r;
		r = dict_pron(dict, w, i+1);
		
		p = mdef_phone_id_nearest(mdef, (s3cipid_t)b, (s3cipid_t)l, (s3cipid_t)r, WORD_POSN_INTERNAL);
		internal[i] = mdef_pid2ssid(mdef, p);
	    }
	    
	    l = b;
	    b = r;
	    if (NOT_S3SSID(rdiph[b][l])) {
		g = rdiph_comsseq(mdef, b, l);
		rdiph[b][l] = ssidlist2comsseq (g, mdef, dict2pid, hs, hp);
		glist_free (g);
	    }
	    internal[pronlen-1] = rdiph[b][l];
	} else if (pronlen == 1) {
	    b = dict_pron(dict, w, 0);
	    if (NOT_S3SSID(single[b])) {
		g = single_comsseq(mdef, b);
		single[b] = ssidlist2comsseq (g, mdef, dict2pid, hs, hp);
		glist_free (g);
		
		for (l = 0; l < mdef_n_ciphone(mdef); l++) {
		    g = single_lc_comsseq(mdef, b, l);
		    dict2pid->single_lc[b][l] = ssidlist2comsseq (g, mdef, dict2pid, hs, hp);
		    glist_free (g);
		}
	    }
	    internal[0] = single[b];
	}
	
	internal += pronlen;
    }
    
    ckd_free_2d ((void **) ldiph);
    ckd_free_2d ((void **) rdiph);
    ckd_free ((void *) single);
    
    /* Allocate space for composite state table */
    cslen = (int32 *) ckd_calloc (dict2pid->n_comstate, sizeof(int32));
    g = hash_tolist(hs, &n);
    assert (n == dict2pid->n_comstate);
    n = 0;
    for (gn = g; gn; gn = gnode_next(gn)) {
	he = (hash_entry_t *) gnode_ptr (gn);
	sen = (s3senid_t *) hash_entry_key(he);
	for (i = 0; IS_S3SENID(sen[i]); i++);
	
	cslen[hash_entry_val(he)] = i+1;	/* +1 for terminating sentinel */
	
	n += (i+1);
    }
    dict2pid->comstate = (s3senid_t **) ckd_calloc (dict2pid->n_comstate, sizeof(s3senid_t *));
    sen = (s3senid_t *) ckd_calloc (n, sizeof(s3senid_t));
    for (i = 0; i < dict2pid->n_comstate; i++) {
	dict2pid->comstate[i] = sen;
	sen += cslen[i];
    }
    
    /* Build composite state table from hash table hs */
    for (gn = g; gn; gn = gnode_next(gn)) {
	he = (hash_entry_t *) gnode_ptr (gn);
	sen = (s3senid_t *) hash_entry_key(he);
	i = hash_entry_val(he);
	
	for (j = 0; j < cslen[i]; j++)
	    dict2pid->comstate[i][j] = sen[j];
	assert (sen[j-1] == BAD_S3SENID);

	ckd_free ((void *)sen);
    }
    ckd_free (cslen);
    glist_free (g);
    hash_free (hs);
    
    /* Allocate space for composite sseq table */
    dict2pid->comsseq = (s3senid_t **) ckd_calloc (dict2pid->n_comsseq, sizeof(s3senid_t *));
    g = hash_tolist (hp, &n);
    assert (n == dict2pid->n_comsseq);
    
    /* Build composite sseq table */
    for (gn = g; gn; gn = gnode_next(gn)) {
	he = (hash_entry_t *) gnode_ptr (gn);
	i = hash_entry_val(he);
	dict2pid->comsseq[i] = (s3senid_t *) hash_entry_key(he);
    }
    glist_free (g);
    hash_free (hp);
    
    /* Weight for each composite state */
    dict2pid->comwt = (int32 *) ckd_calloc (dict2pid->n_comstate, sizeof(int32));
    for (i = 0; i < dict2pid->n_comstate; i++) {
	sen = dict2pid->comstate[i];
	
	for (j = 0; IS_S3SENID(sen[j]); j++);
#if 0
	/* if comstate i has N states, its weight= (1/N^2) (Major Hack!!) */
	dict2pid->comwt[i] = - (logs3 ((float64)j) << 1);
#else
	/* if comstate i has N states, its weight= 1/N */
	dict2pid->comwt[i] = - logs3 ((float64)j);
#endif
    }
    
    E_INFO("%d composite states; %d composite sseq\n",
	   dict2pid->n_comstate, dict2pid->n_comsseq);
    
    return dict2pid;
}


void dict2pid_comsenscr (dict2pid_t *d2p, int32 *senscr, int32 *comsenscr)
{
    int32 i, j;
    int32 best;
    s3senid_t *comstate, k;
    
    for (i = 0; i < d2p->n_comstate; i++) {
	comstate = d2p->comstate[i];
	
	best = senscr[comstate[0]];
	for (j = 1;; j++) {
	    k = comstate[j];
	    if (NOT_S3SENID(k))
		break;
	    if (best < senscr[k])
		best = senscr[k];
	}
	
	comsenscr[i] = best + d2p->comwt[i];
    }
}


void dict2pid_comsseq2sen_active (dict2pid_t *d2p, mdef_t *mdef, int32 *comssid, int32 *sen)
{
    int32 ss, cs, i, j;
    s3senid_t *csp, *sp;	/* Composite state pointer */
    
    for (ss = 0; ss < d2p->n_comsseq; ss++) {
	if (comssid[ss]) {
	    csp = d2p->comsseq[ss];

	    for (i = 0; i < mdef_n_emit_state(mdef); i++) {
		cs = csp[i];
		sp = d2p->comstate[cs];
		
		for (j = 0; IS_S3SENID(sp[j]); j++)
		    sen[sp[j]] = 1;
	    }
	}
    }
}


void dict2pid_dump (FILE *fp, dict2pid_t *d2p, mdef_t *mdef, dict_t *dict)
{
    int32 w, p, pronlen;
    int32 i, j, b, l, r;
    
    fprintf (fp, "# INTERNAL (wd comssid ssid ssid ... ssid comssid)\n");
    for (w = 0; w < dict_size(dict); w++) {
	fprintf (fp, "%30s ", dict_wordstr(dict, w));
	
	pronlen = dict_pronlen(dict, w);
	for (p = 0; p < pronlen; p++)
	    fprintf (fp, " %5d", d2p->internal[w][p]);
	fprintf (fp, "\n");
    }
    fprintf (fp, "#\n");
    
    fprintf (fp, "# LDIPH_LC (b r l ssid)\n");
    for (b = 0; b < mdef_n_ciphone(mdef); b++) {
	for (r = 0; r < mdef_n_ciphone(mdef); r++) {
	    for (l = 0; l < mdef_n_ciphone(mdef); l++) {
		if (IS_S3SSID(d2p->ldiph_lc[b][r][l]))
		    fprintf (fp, "%6s %6s %6s %5d\n",
			     mdef_ciphone_str (mdef, (s3cipid_t)b),
			     mdef_ciphone_str (mdef, (s3cipid_t)r),
			     mdef_ciphone_str (mdef, (s3cipid_t)l),
			     d2p->ldiph_lc[b][r][l]); /* RAH, ldiph_lc is returning an int32, %d expects an int16 */
	    }
	}
    }
    fprintf (fp, "#\n");
    
    fprintf (fp, "# SINGLE_LC (b l comssid)\n");
    for (b = 0; b < mdef_n_ciphone(mdef); b++) {
	for (l = 0; l < mdef_n_ciphone(mdef); l++) {
	    if (IS_S3SSID(d2p->single_lc[b][l]))
		fprintf (fp, "%6s %6s %5d\n",
			 mdef_ciphone_str (mdef, (s3cipid_t)b),
			 mdef_ciphone_str (mdef, (s3cipid_t)l),
			 d2p->single_lc[b][l]);	/* RAH, single_lc is returning an int32, %d expects an int16 */
	}
    }
    fprintf (fp, "#\n");
    
    fprintf (fp, "# SSEQ %d (senid senid ...)\n", mdef->n_sseq);
    for (i = 0; i < mdef->n_sseq; i++) {
	fprintf (fp, "%5d ", i);
	for (j = 0; j < mdef_n_emit_state(mdef); j++)
	    fprintf (fp, " %5d", mdef->sseq[i][j]);
	fprintf (fp, "\n");
    }
    fprintf (fp, "#\n");
    
    fprintf (fp, "# COMSSEQ %d (comstate comstate ...)\n", d2p->n_comsseq);
    for (i = 0; i < d2p->n_comsseq; i++) {
	fprintf (fp, "%5d ", i);
	for (j = 0; j < mdef_n_emit_state(mdef); j++)
	    fprintf (fp, " %5d", d2p->comsseq[i][j]);
	fprintf (fp, "\n");
    }
    fprintf (fp, "#\n");
    
    fprintf (fp, "# COMSTATE %d (senid senid ...)\n", d2p->n_comstate);
    for (i = 0; i < d2p->n_comstate; i++) {
	fprintf (fp, "%5d ", i);
	for (j = 0; IS_S3SENID(d2p->comstate[i][j]); j++)
	    fprintf (fp, " %5d", d2p->comstate[i][j]);
	fprintf (fp, "\n");
    }
    fprintf (fp, "#\n");
    fprintf (fp, "# END\n");
    
    fflush (fp);
}

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	2.916783e+07	2.864011e+07	1.8

                /* ====================================================================
 * Copyright (c) 1999-2001 Carnegie Mellon University.  All rights
 * reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer. 
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 * This work was supported in part by funding from the Defense Advanced 
 * Research Projects Agency and the National Science Foundation of the 
 * United States of America, and the CMU Sphinx Speech Consortium.
 *
 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * ====================================================================
 *
 */
/*
 * cont_mgau.c -- Mixture Gaussians for continuous HMM models.
 *
 * **********************************************
 * CMU ARPA Speech Project
 *
 * Copyright (c) 1997 Carnegie Mellon University.
 * ALL RIGHTS RESERVED.
 * **********************************************
 *
 * HISTORY
 * 
 * 26-May-2004  ARCHAN (archan@cs.cmu.edu)
 *              Incorporate code for numerous fixes. 
 * 20.Apr.2001  RAH (rhoughton@mediasite.com, ricky.houghton@cs.cmu.edu)
 *              Added mgau_free to free memory allocated by mgau_init()
 * 15-Dec-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Added mgau_var_nzvec_floor().
 * 
 * 28-Mar-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Started.
 */


#include "s3types.h"
#include "bio.h"
#include "vector.h"
#include "logs3.h"
#include "cont_mgau.h"


#define MGAU_PARAM_VERSION	"1.0"	/* Sphinx-3 file format version for mean/var */
#define MGAU_MIXW_VERSION	"1.0"	/* Sphinx-3 file format version for mixw */
#define MGAU_MEAN		1
#define MGAU_VAR		2

#ifdef SPEC_CPU
long considered = 0;
#endif

/*
 * Sphinx-3 model mean and var files have the same format.  Use this routine for reading
 * either one.
 */
static int32 mgau_file_read(mgau_model_t *g, char *file_name, int32 type)
{
    char tmp;
    FILE *fp;
    int32 i, k, n;
    int32 n_mgau;
    int32 n_feat;
    int32 n_density;
    int32 *veclen;
    int32 blk;
    int32 byteswap, chksum_present;
    float32 *buf, **pbuf;
    char **argname, **argval;
    uint32 chksum;
    float64 f;
    
    E_INFO("Reading mixture gaussian file '%s'\n", file_name);
    
    fp = myfopen (file_name, "rb");
    
    /* Read header, including argument-value info and 32-bit byteorder magic */
    if (bio_readhdr (fp, &argname, &argval, &byteswap) < 0)
	E_FATAL("bio_readhdr(%s) failed\n", file_name);
    
    /* Parse argument-value list */
    chksum_present = 0;
    for (i = 0; argname[i]; i++) {
	if (strcmp (argname[i], "version") == 0) {
	    if (strcmp(argval[i], MGAU_PARAM_VERSION) != 0)
		E_WARN("Version mismatch(%s): %s, expecting %s\n",
		       file_name, argval[i], MGAU_PARAM_VERSION);
	} else if (strcmp (argname[i], "chksum0") == 0) {
	    chksum_present = 1;	/* Ignore the associated value */
	}
    }
    bio_hdrarg_free (argname, argval);
    argname = argval = NULL;
    
    chksum = 0;
    
    /* #Codebooks */
    if (bio_fread (&n_mgau, sizeof(int32), 1, fp, byteswap, &chksum) != 1)
	E_FATAL("fread(%s) (#codebooks) failed\n", file_name);
    if (n_mgau >= MAX_S3MGAUID) {
	E_FATAL("%s: #Mixture Gaussians (%d) exceeds limit(%d) enforced by MGAUID type\n",
		file_name, n_mgau, MAX_S3MGAUID);
    }
    
    /* #Features/codebook */
    if (bio_fread (&n_feat, sizeof(int32), 1, fp, byteswap, &chksum) != 1)
	E_FATAL("fread(%s) (#features) failed\n", file_name);
    
    if(g->gau_type==CONTHMM){
      if (n_feat != 1)
	E_FATAL("#Features streams(%d) != 1 in continuous HMM\n", n_feat);
    }else if (g->gau_type==SEMIHMM){
      if (n_feat != 4)
	E_FATAL("#Features streams(%d) != 1 in semi-continuous HMM\n", n_feat);
    }
    
    /* #Gaussian densities/feature in each codebook */
    if (bio_fread (&n_density, sizeof(int32), 1, fp, byteswap, &chksum) != 1)
	E_FATAL("fread(%s) (#density/codebook) failed\n", file_name);
    
    /* Vector length of feature stream */

    veclen = ckd_calloc(n_feat, sizeof(uint32));

    if (bio_fread (veclen, sizeof(int32), n_feat, fp, byteswap, &chksum) != n_feat)
	E_FATAL("fread(%s) (feature-lengths) failed\n", file_name);

    for (i = 0, blk = 0; i < n_feat; i++)
	blk += veclen[i];

    /*    if (bio_fread (&veclen, sizeof(int32), 1, fp, byteswap, &chksum) != 1)
	  E_FATAL("fread(%s) (feature vector-length) failed\n", file_name);*/
    
    /* #Floats to follow; for the ENTIRE SET of CODEBOOKS */
    if (bio_fread (&n, sizeof(int32), 1, fp, byteswap, &chksum) != 1)
	E_FATAL("fread(%s) (total #floats) failed\n", file_name);

    if (n != n_mgau * n_density * blk) {
	E_FATAL("%s: #float32s(%d) doesn't match dimensions: %d x %d x %d\n",
		file_name, n, n_mgau, n_density, blk);
    }
    
    if(g->gau_type==SEMIHMM){
      E_FATAL("Currently S2 semi-continous HMM is not supported\n");
    }

    if (type == MGAU_MEAN) {
	/* Allocate memory for mixture gaussian densities */
	g->n_mgau = n_mgau;
	g->max_comp = n_density;
	g->veclen = blk;
	g->mgau = (mgau_t *) ckd_calloc (n_mgau, sizeof(mgau_t));
	
	buf = (float32 *) ckd_calloc (n, sizeof(float));
	pbuf = (float32 **) ckd_calloc (n_mgau * n_density, sizeof(float32 *));


	for (i = 0; i < n_mgau; i++) {
	    g->mgau[i].n_comp = n_density;
	    g->mgau[i].mean = pbuf;
	    
	    for (k = 0; k < n_density; k++) {
		g->mgau[i].mean[k] = buf;
		buf += blk;

	    }
	    pbuf += n_density;
	}
	
	buf = g->mgau[0].mean[0];	/* Restore buf to original value */
    } else {
	assert (type == MGAU_VAR);
	
	if (g->n_mgau != n_mgau)
	    E_FATAL("#Mixtures(%d) doesn't match that of means(%d)\n", n_mgau, g->n_mgau);
	if (g->max_comp != n_density)
	    E_FATAL("#Components(%d) doesn't match that of means(%d)\n", n_density, g->max_comp);
	if (g->veclen != blk)
	    E_FATAL("#Vector length(%d) doesn't match that of means(%d)\n", blk, g->veclen);
	
	buf = (float32 *) ckd_calloc (n, sizeof(float32));
	pbuf = (float32 **) ckd_calloc (n_mgau * n_density, sizeof(float32 *));


	for (i = 0; i < n_mgau; i++) {
	    if (g->mgau[i].n_comp != n_density)
		E_FATAL("Mixture %d: #Components(%d) doesn't match that of means(%d)\n",
			i, n_density, g->mgau[i].n_comp);
	    
	    g->mgau[i].var = pbuf;
	    
	    for (k = 0; k < n_density; k++) {
		g->mgau[i].var[k] = buf;
		buf += blk;

		
	    }
	    pbuf += n_density;
	}
	
	buf = (float32 *) ckd_calloc (n_mgau * n_density, sizeof(float32));

	for (i = 0; i < n_mgau; i++) {
	    g->mgau[i].lrd = buf;
	    buf += n_density;
	}
	
	buf = g->mgau[0].var[0];	/* Restore buf to original value */
    }
    
    /* Read mixture gaussian densities data */
    if (bio_fread (buf, sizeof(float32), n, fp, byteswap, &chksum) != n)
	E_FATAL("fread(%s) (densitydata) failed\n", file_name);

    f = log_to_logs3_factor();    

    if (chksum_present)
	bio_verify_chksum (fp, byteswap, chksum);

    if (fread (&tmp, 1, 1, fp) == 1)
	E_FATAL("%s: More data than expected\n", file_name);
    
    fclose(fp);
    
    E_INFO("%d mixture Gaussians, %d components, %d streams, veclen %d\n", n_mgau, n_density, n_feat, blk);
    
    return 0;
}


static int32 mgau_mixw_read(mgau_model_t *g, char *file_name, float64 mixwfloor)
{
    char **argname, **argval;
    char eofchk;
    FILE *fp;
    int32 byteswap, chksum_present;
    uint32 chksum;
    int32 *buf;
    float32 *pdf;
    int32 i, j, n;
    int32 n_mgau;
    int32 n_feat;
    int32 n_comp;
    int32 n_err;
    
    E_INFO("Reading mixture weights file '%s'\n", file_name);
    
    fp = myfopen (file_name, "rb");
    
    /* Read header, including argument-value info and 32-bit byteorder magic */
    if (bio_readhdr (fp, &argname, &argval, &byteswap) < 0)
	E_FATAL("bio_readhdr(%s) failed\n", file_name);
    
    /* Parse argument-value list */
    chksum_present = 0;
    for (i = 0; argname[i]; i++) {
	if (strcmp (argname[i], "version") == 0) {
	    if (strcmp(argval[i], MGAU_MIXW_VERSION) != 0)
		E_WARN("Version mismatch(%s): %s, expecting %s\n",
			file_name, argval[i], MGAU_MIXW_VERSION);
	} else if (strcmp (argname[i], "chksum0") == 0) {
	    chksum_present = 1;	/* Ignore the associated value */
	}
    }
    bio_hdrarg_free (argname, argval);
    argname = argval = NULL;

    chksum = 0;

    /* Read #senones, #features, #codewords, arraysize */
    if ((bio_fread (&n_mgau, sizeof(int32), 1, fp, byteswap, &chksum) != 1) ||
	(bio_fread (&n_feat, sizeof(int32), 1, fp, byteswap, &chksum) != 1) ||
	(bio_fread (&n_comp, sizeof(int32), 1, fp, byteswap, &chksum) != 1) ||
	(bio_fread (&n, sizeof(int32), 1, fp, byteswap, &chksum) != 1)) {
	E_FATAL("bio_fread(%s) (arraysize) failed\n", file_name);
    }
    if(g->gau_type==CONTHMM){
      if (n_feat != 1)
	E_FATAL("#Features streams(%d) != 1 in continuous HMM\n", n_feat);
    }else if (g->gau_type==SEMIHMM){
      if (n_feat != 4)
	E_FATAL("#Features streams(%d) != 4 in semi-continuous HMM\n", n_feat);
    }else{
	E_FATAL("How can this happen? Someone must have moved this part of the code somewhere! Not my fault! ARCHAN at 20040504 :-)\n");
    }

    if (n != n_mgau * n_comp) {
	E_FATAL("%s: #float32s(%d) doesn't match header dimensions: %d x %d\n",
		file_name, i, n_mgau, n_comp);
    }
    if (n_mgau != g->n_mgau)
	E_FATAL("%s: #Mixture Gaussians(%d) doesn't match mean/var parameters(%d)\n",
		n_mgau, g->n_mgau);
    
    buf = (int32 *) ckd_calloc (n_mgau * n_comp, sizeof(int32));
    for (i = 0; i < n_mgau; i++) {
	if (n_comp != mgau_n_comp(g,i))
	    E_FATAL("Mixture %d: #Weights(%d) doesn't match #Gaussian components(%d)\n",
		    i, n_comp, mgau_n_comp(g,i));
	
	g->mgau[i].mixw = buf;
	buf += g->mgau[i].n_comp;
    }
    
    /* Temporary structure to read in floats before conversion to (int32) logs3 */
    pdf = (float32 *) ckd_calloc (n_comp, sizeof(float32));
    
    /* Read mixw data, normalize, floor, convert to logs3 */
    n_err = 0;
    for (i = 0; i < n_mgau; i++) {
	if (bio_fread((void *)pdf, sizeof(float32), n_comp, fp, byteswap, &chksum) != n_comp)
	    E_FATAL("bio_fread(%s) (arraydata) failed\n", file_name);
	
	/* Normalize and floor */
	if (vector_is_zero (pdf, n_comp)) {
	    n_err++;
	    for (j = 0; j < n_comp; j++)
		mgau_mixw(g,i,j) = S3_LOGPROB_ZERO;
	} else {
	    vector_nz_floor (pdf, n_comp, mixwfloor);
	    vector_sum_norm (pdf, n_comp);
	    for (j = 0; j < n_comp; j++)
		mgau_mixw(g,i,j) = (pdf[j] != 0.0) ? logs3(pdf[j]) : S3_LOGPROB_ZERO;
	}
    }
//SPEC if (n_err > 0)
//SPEC	E_ERROR("Weight normalization failed for %d senones\n", n_err);

    ckd_free (pdf);

    if (chksum_present)
	bio_verify_chksum (fp, byteswap, chksum);
    
    if (fread (&eofchk, 1, 1, fp) == 1)
	E_FATAL("More data than expected in %s\n", file_name);

    fclose(fp);

    E_INFO("Read %d x %d mixture weights\n", n_mgau, n_comp);
    
    return 0;
}


/*
 * Compact each mixture Gaussian in the given model by removing any uninitialized components.
 * A component is considered to be uninitialized if its variance is the 0 vector.  Compact by
 * copying the data rather than moving pointers.  Otherwise, malloc pointers could get
 * corrupted.
 */
static void mgau_uninit_compact (mgau_model_t *g)
{
    int32 m, c, c2, n, nm;
    
    E_INFO("Removing uninitialized Gaussian densities\n");
    
    n = 0;
    nm = 0;
    for (m = 0; m < mgau_n_mgau(g); m++) {
	for (c = 0, c2 = 0; c < mgau_n_comp(g,m); c++) {
	    if (! vector_is_zero (mgau_var(g,m,c), mgau_veclen(g))) {
		if (c2 != c) {
		    memcpy (mgau_mean(g,m,c2), mgau_mean(g,m,c),
			    mgau_veclen(g) * sizeof(float32));
		    memcpy (mgau_var(g,m,c2), mgau_var(g,m,c),
			    mgau_veclen(g) * sizeof(float32));
		    mgau_mixw(g,m,c2) = mgau_mixw(g,m,c);
		}
		c2++;
	    } else {
		n++;
	    }
	}
	mgau_n_comp(g,m) = c2;
	if (c2 == 0) {
	    fprintf (stdout, " %d", m);
	    //SPEC fflush (stderr);
	    nm++;
	}
    }
    if (nm > 0)
	fprintf (stdout, "\n");
    
    if ((nm > 0) || (n > 0))
	E_INFO ("%d densities removed (%d mixtures removed entirely)\n", n, nm);
}


static void mgau_var_floor (mgau_model_t *g, float64 floor)
{
  int32 m, c, i, n;
  
  E_INFO("Applying variance floor\n");
  n = 0;
  for (m = 0; m < mgau_n_mgau(g); m++) {
    for (c = 0; c < mgau_n_comp(g,m); c++) {
      for (i = 0; i < mgau_veclen(g); i++) {
	if (g->mgau[m].var[c][i] < floor) {
	  g->mgau[m].var[c][i] = (float32) floor;
	  n++;
	}
      }
    }
  }
  E_INFO("%d variance values floored\n", n);
}


int32 mgau_var_nzvec_floor (mgau_model_t *g, float64 floor)
{
  int32 m, c, i, n, l;
  float32 *var;
  
  E_INFO("Applying variance floor to non-zero variance vectors\n");
  
  l = mgau_veclen(g);
  
  n = 0;
  for (m = 0; m < mgau_n_mgau(g); m++) {
    for (c = 0; c < mgau_n_comp(g,m); c++) {
      var = g->mgau[m].var[c];
      
      if (! vector_is_zero (var, l)) {
	for (i = 0; i < l; i++) {
	  if (var[i] < floor) {
	    var[i] = (float32) floor;
	    n++;
	  }
	}
      }
    }
  }
  
  E_INFO("%d variance values floored\n", n);
  
  return n;
}


/*
 * Some of the Mahalanobis distance computation (between Gaussian density means and given
 * vectors) can be carried out in advance.  (See comment in .h file.)
 */
static int32 mgau_precomp (mgau_model_t *g)
{
    int32 m, c, i;
    float64 lrd;
    float64 f;

    f = log_to_logs3_factor();            
    E_INFO("Precomputing Mahalanobis distance invariants\n");

    for (m = 0; m < mgau_n_mgau(g); m++) {
	for (c = 0; c < mgau_n_comp(g,m); c++) {
	    lrd = 0.0;
	    for (i = 0; i < mgau_veclen(g); i++) {
		lrd += log(g->mgau[m].var[c][i]);
		
		/* Precompute this part of the exponential */
		g->mgau[m].var[c][i] = (float32) (1.0 / (g->mgau[m].var[c][i] * 2.0));
	    }
	    
	    lrd += mgau_veclen(g) * log(2.0 * PI);	/* (2pi)^velen */
	    mgau_lrd(g,m,c) = (float32)(-0.5 * lrd);	/* Reciprocal, sqrt */
	}
    }
    
    return 0;
}


/* At the moment, S3 models have the same #means in each codebook and 1 var/mean */
mgau_model_t *mgau_init (char *meanfile, char *varfile, float64 varfloor,
			 char *mixwfile, float64 mixwfloor,
			 int32 precomp,char* senmgau)
{
    mgau_model_t *g;
    /*    int32 c_id,f_id; */

    assert (meanfile != NULL);
    assert (varfile != NULL);
    assert (varfloor >= 0.0);
    assert (mixwfile != NULL);
    assert (mixwfloor >= 0.0);
    
    g = (mgau_model_t *) ckd_calloc (1, sizeof(mgau_model_t));

    if(strcmp(senmgau,".cont.") == 0) {
      g->gau_type=CONTHMM;
    }else if(strcmp(senmgau,".semi.") == 0){
      g->gau_type=SEMIHMM;
    }else{
      E_FATAL("Feature should be either .semi. or .cont.");
    }
    
    /* Read means and (diagonal) variances for all mixture gaussians */
    mgau_file_read (g, meanfile, MGAU_MEAN);
    mgau_file_read (g, varfile, MGAU_VAR);
    mgau_mixw_read (g, mixwfile, mixwfloor);
    
    mgau_uninit_compact (g);		/* Delete uninitialized components */
    
    if (varfloor > 0.0)
	mgau_var_floor (g, varfloor);	/* Variance floor after above compaction */
    
    if (precomp)
	mgau_precomp (g);		/* Precompute Mahalanobis distance invariants */
    
    g->distfloor = logs3_to_log (S3_LOGPROB_ZERO);	/* Floor for Mahalanobis distance values */

    
    return g;
}


int32 mgau_comp_eval (mgau_model_t *g, int32 s, float32 *x, int32 *score)
{
    mgau_t *mgau;
    int32 veclen;
    float32 *m, *v;
    float64 dval, diff, f;
    int32 bs;
    int32 i, c;
    
    veclen = mgau_veclen(g);
    mgau = &(g->mgau[s]);
    f = log_to_logs3_factor();

    bs = MAX_NEG_INT32;
    for (c = 0; c < mgau->n_comp; c++) {
	m = mgau->mean[c];
	v = mgau->var[c];
	dval = mgau->lrd[c];
	
	for (i = 0; i < veclen; i++) {
	    diff = x[i] - m[i];
	    dval -= diff * diff * v[i];
	}
	
	if (dval < g->distfloor)
	    dval = g->distfloor;
	
	score[c] = (int32) (f * dval);
	if (score[c] > bs)
	    bs = score[c];
    }
    
    return bs;
}


int32 mgau_eval (mgau_model_t *g, int32 m, int32 *active, float32 *x)
{
    mgau_t *mgau;
    int32 veclen, score;
    /*tmpscore;*/
    float32 *m1, *m2, *v1, *v2;
    float64 dval1, dval2, diff1, diff2, f;


    int32 i, j, c;
    
    veclen = mgau_veclen(g);
    mgau = &(g->mgau[m]);
    f = log_to_logs3_factor();
    score = S3_LOGPROB_ZERO;


    if (! active) {	/* No short list; use all */

     
#if 1
	for (c = 0; c < mgau->n_comp-1; c += 2) {	/* Interleave 2 components for speed */
	    m1 = mgau->mean[c];
	    m2 = mgau->mean[c+1];
	    v1 = mgau->var[c];
	    v2 = mgau->var[c+1];
	    dval1 = mgau->lrd[c];
	    dval2 = mgau->lrd[c+1];
	    
	    for (i = 0; i < veclen; i++) {
		diff1 = x[i] - m1[i];
		dval1 -= diff1 * diff1 * v1[i];
		diff2 = x[i] - m2[i];
		dval2 -= diff2 * diff2 * v2[i];
		/*		E_INFO("x %10f m1 %10f m2 %10f v1 %10f, v2 %10f\n",x[i],m1[i],m2[i],v1[i],v2[i]);
				E_INFO("diff1 %10f,dval1 %10f, diff2 %10f, dval2 %10f\n",diff1,dval1,diff2,dval2);*/
	    }
	    
	    if (dval1 < g->distfloor)	/* Floor */
		dval1 = g->distfloor;
	    if (dval2 < g->distfloor)
		dval2 = g->distfloor;
	    
	    score = logs3_add (score, (int32)(f * dval1) + mgau->mixw[c]);
	    score = logs3_add (score, (int32)(f * dval2) + mgau->mixw[c+1]);
	}
	
	/* Remaining iteration if n_mean odd */
	if (c < mgau->n_comp) {
	    m1 = mgau->mean[c];
	    v1 = mgau->var[c];
	    dval1 = mgau->lrd[c];
	    
	    for (i = 0; i < veclen; i++) {
		diff1 = x[i] - m1[i];
		dval1 -= diff1 * diff1 * v1[i];
	    }
	    
	    if (dval1 < g->distfloor)
		dval1 = g->distfloor;
	    
	    score = logs3_add (score, (int32)(f * dval1) + mgau->mixw[c]);
	}

#endif

    } else {
	for (j = 0; active[j] >= 0; j++) {
#ifdef SPEC_CPU
            considered++;
#endif
	    c = active[j];
	    
	    m1 = mgau->mean[c];
	    v1 = mgau->var[c];
	    dval1 = mgau->lrd[c];
	    
	    for (i = 0; i < veclen; i++) {
		diff1 = x[i] - m1[i];
		dval1 -= diff1 * diff1 * v1[i];
	    }
	    
	    if (dval1 < g->distfloor)
		dval1 = g->distfloor;
	    
	    score = logs3_add (score, (int32)(f * dval1) + mgau->mixw[c]);
	}
    }

    
    if(score == S3_LOGPROB_ZERO){
      /*      E_INFO("Warning!! Score is S3_LOGPROB_ZERO\n");*/
    }
    return score;
}

/* RAH, free memory allocated in mgau_init
   I've not verified that this function catches all of the leaks, just most of them.
 */
void mgau_free (mgau_model_t *g)
{
  //  int i,j;
  if (g) {
    /* Free memory allocated for the mean structure*/
    //    for (i=0;i<g->n_mgau;i++) {
    //      for (j=0;j<g->max_comp;j++) {
    //	if (g->mgau[i].mean[j]) 
    //	  ckd_free ((void *) g->mgau[i].mean[j]);
    //      }
    //    }
    if (g->mgau[0].mean) 
      ckd_free ((void *) g->mgau[0].mean);

    /* Free memory allocated for the var structure*/
//    for (i=0;i<g->n_mgau;i++) {
//      for (j=0;j<g->max_comp;j++) {
//	if (g->mgau[i].var[j]) 
//	  ckd_free ((void *) g->mgau[i].var[j]);
//      }
//    }
    if (g->mgau[0].var) 
      ckd_free ((void *) g->mgau[0].var);
    if (g->mgau[0].lrd) 
      ckd_free ((void *) g->mgau[0].lrd);

    /* Free memory allocated for the mixture weights*/
//    for (i=0;i<g->n_mgau;i++) 
//      if (g->mgau[i].mixw[0]) 
//	ckd_free ((void *) g->mgau[i].mixw[0]);
    
    if (g->mgau[0].mixw) 
      ckd_free ((void *) g->mgau[0].mixw);
    
    if (g->mgau)
      ckd_free ((void *) g->mgau);
    ckd_free ((void *) g);
  }
}

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	2.452782e+08	2.395778e+08	2.3

                /* ====================================================================
 * Copyright (c) 1999-2001 Carnegie Mellon University.  All rights
 * reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer. 
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 * This work was supported in part by funding from the Defense Advanced 
 * Research Projects Agency and the National Science Foundation of the 
 * United States of America, and the CMU Sphinx Speech Consortium.
 *
 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * ====================================================================
 *
 */
/*
 * mdef.c -- HMM model definition: base (CI) phones and triphones
 *
 * **********************************************
 * CMU ARPA Speech Project
 *
 * Copyright (c) 1999 Carnegie Mellon University.
 * ALL RIGHTS RESERVED.
 * **********************************************
 * 
 * HISTORY
 * 
 * 19.Apr-2001  Ricky Houghton, added code for free allocated memory
 *
 * 14-Oct-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon
 * 		Added mdef_sseq2sen_active().
 * 
 * 06-May-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon
 * 		In mdef_phone_id(), added backing off to silence phone context from filler
 * 		context if original triphone not found.
 * 
 * 30-Apr-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon
 * 		Added senone-sequence id (ssid) to phone_t and appropriate functions to
 * 		maintain it.  Instead, moved state sequence info to mdef_t.
 * 
 * 13-Jul-96	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Added mdef_phone_str().
 * 
 * 01-Jan-96	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Allowed mdef_phone_id_nearest to return base phone id if either
 * 		left or right context (or both) is undefined.
 * 
 * 01-Jan-96	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Created.
 */


/*
 * Major assumptions:
 *   All phones have same #states, same topology.
 *   Every phone has exactly one non-emitting, final state--the last one.
 *   CI phones must appear first in model definition file.
 */


#include "mdef.h"


#define MODEL_DEF_VERSION	"0.3"


void mdef_dump (FILE *fp, mdef_t *m)
{
    int32 i, j;
    int32 ssid;
    char buf[1024];
    
    fprintf (fp, "%d ciphone\n", m->n_ciphone);
    fprintf (fp, "%d phone\n", m->n_phone);
    fprintf (fp, "%d emitstate\n", m->n_emit_state);
    fprintf (fp, "%d cisen\n", m->n_ci_sen);
    fprintf (fp, "%d sen\n", m->n_sen);
    fprintf (fp, "%d tmat\n", m->n_tmat);
    
    for (i = 0; i < m->n_phone; i++) {
	mdef_phone_str (m, i, buf);
	ssid = m->phone[i].ssid;
	
	fprintf (fp, "%3d %5d", m->phone[i].tmat, ssid);
	for (j = 0; j < m->n_emit_state; j++)
	    fprintf (fp, " %5d", m->sseq[ssid][j]);
	fprintf (fp, "\t");
	for (j = 0; j < m->n_emit_state; j++)
	    fprintf (fp, " %3d", m->cd2cisen[m->sseq[ssid][j]]);
	fprintf (fp, "\t%s\n", buf);
    }
    
    fflush (fp);
}


#if 0
int32 mdef_hmm_cmp (mdef_t *m, s3pid_t p1, s3pid_t p2)
{
    int32 i;
    
    if (m->phone[p1].tmat != m->phone[p2].tmat)
	return -1;
    
    for (i = 0; i < m->n_emit_state; i++)
	if (m->phone[p1].state[i] != m->phone[p2].state[i])
	    return -1;
    
    return 0;
}
#endif


static void ciphone_add (mdef_t *m, char *ci, s3pid_t p)
{
    assert (p < m->n_ciphone);

    m->ciphone[p].name = (char *) ckd_salloc (ci); /* freed in mdef_free */
    if (hash_enter (m->ciphone_ht, m->ciphone[p].name, p) != p)
	E_FATAL("hash_enter(%s) failed; duplicate CIphone?\n", m->ciphone[p].name);
}


static ph_lc_t *find_ph_lc (ph_lc_t *lclist, s3cipid_t lc)
{
    ph_lc_t *lcptr;

    for (lcptr = lclist; lcptr && (lcptr->lc != lc); lcptr = lcptr->next);
    return lcptr;
}


static ph_rc_t *find_ph_rc (ph_rc_t *rclist, s3cipid_t rc)
{
    ph_rc_t *rcptr;

    for (rcptr = rclist; rcptr && (rcptr->rc != rc); rcptr = rcptr->next);
    return rcptr;
}


static void triphone_add (mdef_t *m,
			  s3cipid_t ci, s3cipid_t lc, s3cipid_t rc, word_posn_t wpos,
			  s3pid_t p)
{
    ph_lc_t *lcptr;
    ph_rc_t *rcptr;
    
    assert (p < m->n_phone);

    /* Fill in phone[p] information (state and tmat mappings added later) */
    m->phone[p].ci = ci;
    m->phone[p].lc = lc;
    m->phone[p].rc = rc;
    m->phone[p].wpos = wpos;

    /* Create <ci,lc,rc,wpos> -> p mapping if not a CI phone */
    if (p >= m->n_ciphone) {
	if ((lcptr = find_ph_lc (m->wpos_ci_lclist[wpos][(int)ci], lc)) 
	    == NULL)
	{
	lcptr = (ph_lc_t *) ckd_calloc (1, sizeof(ph_lc_t)); /* freed at mdef_free, I believe */
	    lcptr->lc = lc;
	    lcptr->next = m->wpos_ci_lclist[wpos][(int)ci];
	    m->wpos_ci_lclist[wpos][(int)ci] = lcptr; /* This is what needs to be freed */
	}
	if ((rcptr = find_ph_rc (lcptr->rclist, rc)) != NULL) {
	    char buf[4096];
	    
	    mdef_phone_str (m, rcptr->pid, buf);
	    E_FATAL("Duplicate triphone: %s\n", buf);
	}
	
      rcptr = (ph_rc_t *) ckd_calloc (1, sizeof(ph_rc_t)); /* freed in mdef_free, I believe */
	rcptr->rc = rc;
	rcptr->pid = p;
	rcptr->next = lcptr->rclist;
	lcptr->rclist = rcptr;
    }
}


s3cipid_t mdef_ciphone_id (mdef_t *m, char *ci)
{
    int32 id;
    
    assert (m);
    assert (ci);
    
    if (hash_lookup (m->ciphone_ht, ci, &id) < 0)
	return (BAD_S3CIPID);
    return ((s3cipid_t) id);
}


const char *mdef_ciphone_str (mdef_t *m, s3cipid_t id)
{
    assert (m);
    assert ((id >= 0) && (id < m->n_ciphone));
    
    return (m->ciphone[(int)id].name);
}


int32 mdef_phone_str (mdef_t *m, s3pid_t pid, char *buf)
{
    char *wpos_name;
    
    assert (m);
    assert ((pid >= 0) && (pid < m->n_phone));
    wpos_name = WPOS_NAME;
    
    buf[0] = '\0';
    if (pid < m->n_ciphone)
	sprintf (buf, "%s", mdef_ciphone_str (m, (s3cipid_t) pid));
    else {
	sprintf (buf, "(%s,%s,%s,%c)",
		 mdef_ciphone_str(m, m->phone[pid].ci),
		 mdef_ciphone_str(m, m->phone[pid].lc),
		 mdef_ciphone_str(m, m->phone[pid].rc),
		 wpos_name[m->phone[pid].wpos]);
    }
    return 0;
}


s3pid_t mdef_phone_id (mdef_t *m, 
		       s3cipid_t ci, s3cipid_t lc, s3cipid_t rc, word_posn_t wpos)
{
    ph_lc_t *lcptr;
    ph_rc_t *rcptr;
    s3cipid_t newl, newr;
    
    assert (m);
    assert ((ci >= 0) && (ci < m->n_ciphone));
    assert ((lc >= 0) && (lc < m->n_ciphone));
    assert ((rc >= 0) && (rc < m->n_ciphone));
    assert ((wpos >= 0) && (wpos < N_WORD_POSN));

    if (((lcptr = find_ph_lc (m->wpos_ci_lclist[wpos][(int)ci], lc)) == NULL) ||
	((rcptr = find_ph_rc (lcptr->rclist, rc)) == NULL)) {
	/* Not found; backoff to silence context if non-silence filler context */
	if (NOT_S3CIPID(m->sil))
	    return BAD_S3PID;
	
	newl = m->ciphone[(int)lc].filler ? m->sil : lc;
	newr = m->ciphone[(int)rc].filler ? m->sil : rc;
	if ((newl == lc) && (newr == rc))
	    return BAD_S3PID;
	
	return (mdef_phone_id (m, ci, newl, newr, wpos));
    }
    
    return (rcptr->pid);
}


s3pid_t mdef_phone_id_nearest (mdef_t *m, 
			       s3cipid_t b, s3cipid_t l, s3cipid_t r, word_posn_t pos)
{
    word_posn_t tmppos;
    s3pid_t p;
    s3cipid_t newl, newr;
    char *wpos_name;
    
    assert (m);
    assert ((b >= 0) && (b < m->n_ciphone));
    assert ((pos >= 0) && (pos < N_WORD_POSN));

    if ((NOT_S3CIPID(l)) || (NOT_S3CIPID(r)))
	return ((s3pid_t) b);
    
    assert ((l >= 0) && (l < m->n_ciphone));
    assert ((r >= 0) && (r < m->n_ciphone));
    
    p = mdef_phone_id (m, b, l, r, pos);
    if (IS_S3PID(p))
	return p;
    
    /* Exact triphone not found; backoff to other word positions */
    for (tmppos = 0; tmppos < N_WORD_POSN; tmppos++) {
	if (tmppos != pos) {
	    p = mdef_phone_id (m, b, l, r, tmppos);
	    if (IS_S3PID(p))
		return p;
	}
    }
    
    /* Nothing yet; backoff to silence phone if non-silence filler context */
    if (IS_S3CIPID(m->sil)) {
	newl = m->ciphone[(int)l].filler ? m->sil : l;
	newr = m->ciphone[(int)r].filler ? m->sil : r;
	if ((newl != l) || (newr != r)) {
	    p = mdef_phone_id (m, b, newl, newr, pos);
	    if (IS_S3PID(p))
		return p;
	    
	    for (tmppos = 0; tmppos < N_WORD_POSN; tmppos++) {
		if (tmppos != pos) {
		    p = mdef_phone_id (m, b, newl, newr, tmppos);
		    if (IS_S3PID(p))
			return p;
		}
	    }
	}
    }
    
    /* Nothing yet; backoff to base phone */
    if ((m->n_phone > m->n_ciphone) && (! m->ciphone[(int)b].filler)) {
	wpos_name = WPOS_NAME;
#if 0
	E_WARN("Triphone(%s,%s,%s,%c) not found; backing off to CIphone\n",
	       mdef_ciphone_str(m, b),
	       mdef_ciphone_str(m, l),
	       mdef_ciphone_str(m, r),
	       wpos_name[pos]);
#endif
    }
    return ((s3pid_t) b);
}


int32 mdef_phone_components (mdef_t *m,
			     s3pid_t p,
			     s3cipid_t *b,
			     s3cipid_t *l,
			     s3cipid_t *r,
			     word_posn_t *pos)
{
    assert (m);
    assert ((p >= 0) && (p < m->n_phone));

    *b = m->phone[p].ci;
    *l = m->phone[p].lc;
    *r = m->phone[p].rc;
    *pos = m->phone[p].wpos;

    return 0;
}


int32 mdef_is_ciphone (mdef_t *m, s3pid_t p)
{
    assert (m);
    assert ((p >= 0) && (p < m->n_phone));
    
    return ((p < m->n_ciphone) ? 1 : 0);
}


/* Parse tmat and state->senone mappings for phone p and fill in structure */
static void parse_tmat_senmap (mdef_t *m, char *line, int32 off, s3pid_t p)
{
    int32 wlen, n, s;
    char word[1024], *lp;

    lp = line + off;
    
    /* Read transition matrix id */
    if ((sscanf (lp, "%d%n", &n, &wlen) != 1) || (n < 0))
        E_FATAL("Missing or bad transition matrix id: %s\n", line);
    m->phone[p].tmat = n;
    if (m->n_tmat <= n)
	E_FATAL("tmat-id(%d) > #tmat in header(%d): %s\n", n, m->n_tmat, line);
    lp += wlen;
    
    /* Read senone mappings for each emitting state */
    for (n = 0; n < m->n_emit_state; n++) {
        if ((sscanf (lp, "%d%n", &s, &wlen) != 1) || (s < 0))
	    E_FATAL("Missing or bad state[%d]->senone mapping: %s\n", n, line);
	if ((p < m->n_ciphone) && (m->n_ci_sen <= s))
	    E_FATAL("CI-senone-id(%d) > #CI-senones(%d): %s\n", s, m->n_ci_sen, line);
	if (m->n_sen <= s)
	    E_FATAL("Senone-id(%d) > #senones(%d): %s\n", s, m->n_sen, line);

	m->sseq[p][n] = s;
	lp += wlen;
    }

    /* Check for the last non-emitting state N */
    if ((sscanf (lp, "%s%n", word, &wlen) != 1) || (strcmp (word, "N") != 0))
        E_FATAL("Missing non-emitting state spec: %s\n", line);
    lp += wlen;

    /* Check for end of line */
    if (sscanf (lp, "%s%n", word, &wlen) == 1)
        E_FATAL("Non-empty beyond non-emitting final state: %s\n", line);
}


static void parse_base_line (mdef_t *m, char *line, s3pid_t p)
{
    int32 wlen, n;
    char word[1024], *lp;
    s3cipid_t ci;

    lp = line;
    
    /* Read base phone name */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing base phone name: %s\n", line);
    lp += wlen;
    
    /* Make sure it's not a duplicate */
    ci = mdef_ciphone_id (m, word);
    if (IS_S3CIPID(ci))
        E_FATAL("Duplicate base phone: %s\n", line);

    /* Add ciphone to ciphone table with id p */
    ciphone_add (m, word, p);
    ci = (s3cipid_t) p;

    /* Read and skip "-" for lc, rc, wpos */
    for (n = 0; n < 3; n++) {
	if ((sscanf (lp, "%s%n", word, &wlen) != 1) || (strcmp (word, "-") != 0))
	    E_FATAL("Bad context info for base phone: %s\n", line);
	lp += wlen;
    }
    
    /* Read filler attribute, if present */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing filler atribute field: %s\n", line);
    lp += wlen;
    if (strcmp (word, "filler") == 0)
        m->ciphone[(int)ci].filler = 1;
    else if (strcmp (word, "n/a") == 0)
        m->ciphone[(int)ci].filler = 0;
    else
        E_FATAL("Bad filler attribute field: %s\n", line);

    triphone_add (m, ci, BAD_S3CIPID, BAD_S3CIPID, WORD_POSN_UNDEFINED, p);

    /* Parse remainder of line: transition matrix and state->senone mappings */
    parse_tmat_senmap (m, line, lp-line, p);
}


static void parse_tri_line (mdef_t *m, char *line, s3pid_t p)
{
    int32 wlen;
    char word[1024], *lp;
    s3cipid_t ci, lc, rc;
    word_posn_t wpos = WORD_POSN_BEGIN;

    lp = line;
    
    /* Read base phone name */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing base phone name: %s\n", line);
    lp += wlen;

    ci = mdef_ciphone_id (m, word);
    if (NOT_S3CIPID(ci))
        E_FATAL("Unknown base phone: %s\n", line);

    /* Read lc */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing left context: %s\n", line);
    lp += wlen;
    lc = mdef_ciphone_id (m, word);
    if (NOT_S3CIPID(lc))
        E_FATAL("Unknown left context: %s\n", line);

    /* Read rc */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing right context: %s\n", line);
    lp += wlen;
    rc = mdef_ciphone_id (m, word);
    if (NOT_S3CIPID(rc))
        E_FATAL("Unknown right  context: %s\n", line);
    
    /* Read tripone word-position within word */
    if ((sscanf (lp, "%s%n", word, &wlen) != 1) || (word[1] != '\0'))
        E_FATAL("Missing or bad word-position spec: %s\n", line);
    lp += wlen;
    switch (word[0]) {
    case 'b': wpos = WORD_POSN_BEGIN; break;
    case 'e': wpos = WORD_POSN_END; break;
    case 's': wpos = WORD_POSN_SINGLE; break;
    case 'i': wpos = WORD_POSN_INTERNAL; break;
    default: E_FATAL("Bad word-position spec: %s\n", line);
    }

    /* Read filler attribute, if present.  Must match base phone attribute */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing filler attribute field: %s\n", line);
    lp += wlen;
    if (((strcmp (word, "filler") == 0) && (m->ciphone[(int)ci].filler)) ||
	((strcmp (word, "n/a") == 0) && (! m->ciphone[(int)ci].filler))) {
	/* Everything is fine */
    } else
        E_FATAL("Bad filler attribute field: %s\n", line);
    
    triphone_add (m, ci, lc, rc, wpos, p);

    /* Parse remainder of line: transition matrix and state->senone mappings */
    parse_tmat_senmap (m, line, lp-line, p);
}


static void sseq_compress (mdef_t *m)
{
    hash_table_t *h;
    s3senid_t **sseq;
    int32 n_sseq;
    int32 p, j, k;
    glist_t g;
    gnode_t *gn;
    hash_entry_t *he;
    
    k = m->n_emit_state * sizeof(s3senid_t);
    
    h = hash_new (m->n_phone, HASH_CASE_YES);
    n_sseq = 0;
    
    /* Identify unique senone-sequence IDs.  BUG: tmat-id not being considered!! */
    for (p = 0; p < m->n_phone; p++) {
	/* Add senone sequence to hash table */
	if ((j = hash_enter_bkey (h, (char *)(m->sseq[p]), k, n_sseq)) == n_sseq)
	    n_sseq++;
	
	m->phone[p].ssid = j;
    }
    
    /* Generate compacted sseq table */
    sseq = (s3senid_t **) ckd_calloc_2d (n_sseq, m->n_emit_state, sizeof(s3senid_t));/* freed in mdef_free() */
    
    g = hash_tolist (h, &j);
    assert (j == n_sseq);
    
    for (gn = g; gn; gn = gnode_next(gn)) {
	he = (hash_entry_t *) gnode_ptr (gn);
	j = hash_entry_val(he);
	memcpy (sseq[j], hash_entry_key(he), k);
    }
    glist_free (g);
    
    /* Free the old, temporary senone sequence table, replace with compacted one */
    ckd_free_2d ((void **) m->sseq);
    m->sseq = sseq;
    m->n_sseq = n_sseq;
    
    hash_free (h);
}


static int32 noncomment_line(char *line, int32 size, FILE *fp)
{
    while (fgets (line, size, fp) != NULL) {
        if (line[0] != '#')
	    return 0;
    }
    return -1;
}


/*
 * Initialize phones (ci and triphones) and state->senone mappings from .mdef file.
 */
mdef_t *mdef_init (char *mdeffile)
{
    FILE *fp;
    int32 n_ci, n_tri, n_map, n;
    char tag[1024], buf[1024];
    s3senid_t **senmap;
    s3pid_t p;
    int32 s, ci, cd;
    mdef_t *m;
    int32 *cdsen_start, *cdsen_end;

    if (! mdeffile)
	E_FATAL("No mdef-file\n");

    E_INFO("Reading model definition: %s\n", mdeffile);

    m = (mdef_t *) ckd_calloc (1, sizeof(mdef_t)); /* freed in mdef_free */
    
    if ((fp = fopen(mdeffile, "r")) == NULL)
        E_FATAL_SYSTEM("fopen(%s,r) failed\n", mdeffile);

    if (noncomment_line(buf, sizeof(buf), fp) < 0)
        E_FATAL("Empty file: %s\n", mdeffile);

    if (strncmp(buf, MODEL_DEF_VERSION, strlen(MODEL_DEF_VERSION)) != 0)
        E_FATAL("Version error: Expecing %s, but read %s\n", MODEL_DEF_VERSION, buf);

    /* Read #base phones, #triphones, #senone mappings defined in header */
    n_ci = -1;
    n_tri = -1;
    n_map = -1;
    m->n_ci_sen = -1;
    m->n_sen = -1;
    m->n_tmat = -1;
    do {
	if (noncomment_line(buf, sizeof(buf), fp) < 0)
	    E_FATAL("Incomplete header\n");

	if ((sscanf(buf, "%d %s", &n, tag) != 2) || (n < 0))
	    E_FATAL("Error in header: %s\n", buf);

	if (strcmp(tag, "n_base") == 0)
	    n_ci = n;
	else if (strcmp(tag, "n_tri") == 0)
	    n_tri = n;
	else if (strcmp(tag, "n_state_map") == 0)
	    n_map = n;
	else if (strcmp(tag, "n_tied_ci_state") == 0)
	    m->n_ci_sen = n;
	else if (strcmp(tag, "n_tied_state") == 0)
	    m->n_sen = n;
	else if (strcmp(tag, "n_tied_tmat") == 0)
	    m->n_tmat = n;
	else
	    E_FATAL("Unknown header line: %s\n", buf);
    } while ((n_ci < 0) || (n_tri < 0) || (n_map < 0) ||
	     (m->n_ci_sen < 0) || (m->n_sen < 0) || (m->n_tmat < 0));

    if ((n_ci == 0) || (m->n_ci_sen == 0) || (m->n_tmat == 0) || (m->n_ci_sen > m->n_sen))
        E_FATAL("%s: Error in header\n", mdeffile);
    
    /* Check typesize limits */
    if (n_ci >= MAX_S3CIPID)
	E_FATAL("%s: #CI phones (%d) exceeds limit (%d)\n", mdeffile, n_ci, MAX_S3CIPID);
    if (n_ci + n_tri >= MAX_S3PID)
	E_FATAL("%s: #Phones (%d) exceeds limit (%d)\n", mdeffile, n_ci+n_tri, MAX_S3PID);
    if (m->n_sen >= MAX_S3SENID)
	E_FATAL("%s: #senones (%d) exceeds limit (%d)\n", mdeffile, m->n_sen, MAX_S3SENID);
    if (m->n_tmat >= MAX_S3TMATID)
	E_FATAL("%s: #tmats (%d) exceeds limit (%d)\n", mdeffile, m->n_tmat, MAX_S3TMATID);
    
    m->n_emit_state = (n_map / (n_ci+n_tri)) - 1;
    if ((m->n_emit_state+1) * (n_ci+n_tri) != n_map)
        E_FATAL("Header error: n_state_map not a multiple of n_ci*n_tri\n");

    /* Initialize ciphone info */
    m->n_ciphone = n_ci;
    m->ciphone_ht = hash_new (n_ci, 1);	/* With case-insensitive string names */ /* freed in mdef_free */
    m->ciphone = (ciphone_t *) ckd_calloc (n_ci, sizeof(ciphone_t)); /* freed in mdef_free */
    /* RAH, let's null the pointers so that we can reliably deallocate them */
    /*    for (i=0;i<m->n_ciphone;i++) { */ /*  */
    /*m->ciphone[i].name = NULL; */ /*  */
    /*} */ /*  */

    /* Initialize phones info (ciphones + triphones) */
    m->n_phone = n_ci + n_tri;
    m->phone = (phone_t *) ckd_calloc (m->n_phone, sizeof(phone_t)); /* freed in mdef_free */

    /* Allocate space for state->senone map for each phone */
    senmap = (s3senid_t **) ckd_calloc_2d (m->n_phone, m->n_emit_state, sizeof(s3senid_t));/* freed in mdef_free */
    m->sseq = senmap;	/* TEMPORARY; until it is compressed into just the unique ones */
    
    /* Allocate initial space for <ci,lc,rc,wpos> -> pid mapping */
    m->wpos_ci_lclist = (ph_lc_t ***) ckd_calloc_2d (N_WORD_POSN, m->n_ciphone, sizeof(ph_lc_t *)); /* freed in mdef_free */

    /*
     * Read base phones and triphones.  They'll simply be assigned a running sequence
     * number as their "phone-id".  If the phone-id < n_ci, it's a ciphone.
     */

    /* Read base phones */
    for (p = 0; p < n_ci; p++) {
        if (noncomment_line(buf, sizeof(buf), fp) < 0)
	    E_FATAL("Premature EOF reading CIphone %d\n", p);
        parse_base_line (m, buf, p);
    }
    m->sil = mdef_ciphone_id (m, S3_SILENCE_CIPHONE);
    
    /* Read triphones, if any */
    for (; p < m->n_phone; p++) {
        if (noncomment_line(buf, sizeof(buf), fp) < 0)
	    E_FATAL("Premature EOF reading phone %d\n", p);
        parse_tri_line (m, buf, p);
    }

    if (noncomment_line(buf, sizeof(buf), fp) >= 0)
	E_ERROR("Non-empty file beyond expected #phones (%d)\n", m->n_phone);

    /* Build CD senones to CI senones map */
    if (m->n_ciphone * m->n_emit_state != m->n_ci_sen)
	E_FATAL("#CI-senones(%d) != #CI-phone(%d) x #emitting-states(%d)\n",
		m->n_ci_sen, m->n_ciphone, m->n_emit_state);
    m->cd2cisen = (s3senid_t *) ckd_calloc (m->n_sen, sizeof(s3senid_t)); /* freed in mdef_free */

    m->sen2cimap = (s3cipid_t *) ckd_calloc (m->n_sen, sizeof(s3cipid_t)); /* freed in mdef_free */

    for (s = 0; s < m->n_sen; s++)
	m->sen2cimap[s] = BAD_S3CIPID;
    for (s = 0; s < m->n_ci_sen; s++) {		/* CI senones */
	m->cd2cisen[s] = (s3senid_t) s;
	m->sen2cimap[s] = s / m->n_emit_state;
    }
    for (p = n_ci; p < m->n_phone; p++) {	/* CD senones */
	for (s = 0; s < m->n_emit_state; s++) {
	    cd = m->sseq[p][s];
	    ci = m->sseq[(int)m->phone[p].ci][s];
	    m->cd2cisen[cd] = (s3senid_t) ci;
	    m->sen2cimap[cd] = m->phone[p].ci;
	}
    }
    
    /*
     * Count #senones (CI+CD) for each CI phone.
     * HACK!!  For handling holes in senone-CIphone mappings.  Does not work if holes
     * are present at the beginning or end of senones for a given CIphone.
     */
    cdsen_start = (int32 *) ckd_calloc (m->n_ciphone, sizeof(int32)); /* freed locally */

    cdsen_end = (int32 *) ckd_calloc (m->n_ciphone, sizeof(int32)); /* freed locally */

    for (s = m->n_ci_sen; s < m->n_sen; s++) {
	if (NOT_S3CIPID(m->sen2cimap[s]))
	    continue;
	
	if (! cdsen_start[(int)m->sen2cimap[s]])
	    cdsen_start[(int)m->sen2cimap[s]] = s;
	cdsen_end[(int)m->sen2cimap[s]] = s;
    }

    /* Fill up holes */
    for (s = m->n_ci_sen; s < m->n_sen; s++) {
	if (IS_S3CIPID(m->sen2cimap[s]))
	    continue;

	/* Check if properly inside the observed ranges above */
	for (p = 0; p < m->n_ciphone; p++) {
	    if ((s > cdsen_start[p]) && (s < cdsen_end[p]))
		break;
	}
	if (p >= m->n_ciphone)
	    E_FATAL("Unreferenced senone %d; cannot determine parent CIphone\n", s);
	m->sen2cimap[s] = p;
    }

    /* Build #CD-senones for each CIphone */
    m->ciphone2n_cd_sen = (int32 *) ckd_calloc (m->n_ciphone, sizeof(int32));/* freed mdef_free */
    n = 0;
    for (p = 0; p < m->n_ciphone; p++) {
	if (cdsen_start[p] > 0) {
	    m->ciphone2n_cd_sen[p] = cdsen_end[p] - cdsen_start[p] + 1;
	    n += m->ciphone2n_cd_sen[p];
	}
    }
    n += m->n_ci_sen;
    assert (n == m->n_sen);

    ckd_free (cdsen_start);
    ckd_free (cdsen_end);

    sseq_compress (m);
    
    E_INFO("%d CI-phone, %d CD-phone, %d emitstate/phone, %d CI-sen, %d Sen, %d Sen-Seq\n",
	   m->n_ciphone, m->n_phone - m->n_ciphone, m->n_emit_state,
	   m->n_ci_sen, m->n_sen, m->n_sseq);
    
    fclose (fp);
    
    return m;
}


void mdef_sseq2sen_active (mdef_t *mdef, int32 *sseq, int32 *sen)
{
    int32 ss, i;
    s3senid_t *sp;
    
    for (ss = 0; ss < mdef_n_sseq(mdef); ss++) {
	if (sseq[ss]) {
	    sp = mdef->sseq[ss];
	    for (i = 0; i < mdef_n_emit_state(mdef); i++)
		sen[sp[i]] = 1;
	}
    }
}

/* RAH 4.23.01, Need to step down the ->next list to see if there are
   any more things to free
 */



/* RAH 4.19.01, Attempt to free memory that was allocated within this module
   I have not verified that all the memory has been freed. I've taken only a 
   reasonable effort for now.
   RAH 4.24.01 - verified that all memory is released.
 */
void mdef_free_recursive_lc (ph_lc_t *lc)
{
  if (lc == NULL) return;

  if (lc->rclist) 
    mdef_free_recursive_rc (lc->rclist);

  if (lc->next) 
    mdef_free_recursive_lc (lc->next);

  ckd_free ((void *) lc);
}

void mdef_free_recursive_rc (ph_rc_t *rc)
{
  if (rc == NULL) return;

    if (rc->next) 
      mdef_free_recursive_rc (rc->next);

    ckd_free ((void *) rc);
}


/* RAH, Free memory that was allocated in mdef_init 
   Rational purify shows that no leaks exist
 */
   
void mdef_free (mdef_t *m)
{
  int i,j;

  if (m) { 
    if (m->ciphone2n_cd_sen)
      ckd_free    ((void *)m->ciphone2n_cd_sen);
    if (m->sen2cimap)
      ckd_free    ((void *)m->sen2cimap);
    if (m->cd2cisen)
      ckd_free    ((void *)m->cd2cisen);

    /* RAH, go down the ->next list and delete all the pieces */
    for (i=0;i<N_WORD_POSN;i++)
      for (j=0;j<m->n_ciphone;j++) 
	if (m->wpos_ci_lclist[i][j]) {
	  mdef_free_recursive_lc (m->wpos_ci_lclist[i][j]->next);
	  mdef_free_recursive_rc (m->wpos_ci_lclist[i][j]->rclist);
	}
    
    for (i=0;i<N_WORD_POSN;i++)
      for (j=0;j<m->n_ciphone;j++) 
	if (m->wpos_ci_lclist[i][j])  
	  ckd_free ((void *) m->wpos_ci_lclist[i][j]);

    
    if (m->wpos_ci_lclist)
      ckd_free_2d ((void *)m->wpos_ci_lclist);
    if (m->sseq) 
      ckd_free_2d ((void *)m->sseq);
    /* Free phone context */
    if (m->phone) 
      ckd_free    ((void *)m->phone);    
    if (m->ciphone_ht)
      hash_free (m->ciphone_ht);

    for (i=0;i<m->n_ciphone;i++) {
      if (m->ciphone[i].name) 
	ckd_free    ((void *)m->ciphone[i].name);
    }
    if (m->ciphone) 
      ckd_free    ((void *)m->ciphone);
    
    ckd_free    ((void *)m);
  }
}

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
26255	### 1	80375.7622496213	7.266000e+03	9.680000e+02	86.6
48356	### 2	54447.8108987715	7.807000e+03	1.004000e+03	87.1

                /* ====================================================================
 * Copyright (c) 1999-2001 Carnegie Mellon University.  All rights
 * reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer. 
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 * This work was supported in part by funding from the Defense Advanced 
 * Research Projects Agency and the National Science Foundation of the 
 * United States of America, and the CMU Sphinx Speech Consortium.
 *
 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * ====================================================================
 *
 */
/*
 * subvq.c
 * 
 * **********************************************
 * CMU ARPA Speech Project
 *
 * Copyright (c) 1999 Carnegie Mellon University.
 * ALL RIGHTS RESERVED.
 * **********************************************
 * 
 * HISTORY
 * 
 * 20.Apr.2001  RAH (rhoughton@mediasite.com, ricky.houghton@cs.cmu.edu)
 *              Updated subvq_free () to free allocated memory
 * 
 * 17-Dec-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Added handling of a single sub-vector in subvq_mgau_shortlist().
 * 
 * 15-Dec-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Changes due to moving subvq_t.{frm_sen_eval,frm_gau_eval} to cont_mgau.h.
 * 
 * 14-Dec-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Added subvq_t.{frm_sen_eval,frm_gau_eval}.  Changed subvq_frame_eval to
 * 		return the normalization factor.
 * 
 * 06-Dec-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Added subvq_subvec_eval_logs3().
 * 
 * 14-Oct-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Changed ci_active flags input to sen_active in subvq_frame_eval().
 * 
 * 21-Jul-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Bugfix in subvq_init() that used g for allocating working areas, even
 *		though g wasn't defined.
 * 
 * 20-Jul-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Added subvq_gautbl_eval_logs3() and used it in subvq_frame_eval().
 * 
 * 12-Mar-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Started.
 */


#include "subvq.h"
/* #include "cmd_ln_args.h"	*/ /* RAH, added so we can allow for -vqeval parameter */
#include "s3types.h"

/* RAH, 5.8.01, VQ_EVAL determines how many vectors are used to
 * compute the shortlist, for now this value is only relevant when n_sv =3.
 * Setting it to 1 means that only the CEP values are estimated, 2 means that 
 * CEP and delta values are estimated, 3 means all three are estimated.
 * Note, we must adjust the beam widths as we muck around with these.
 */
static int VQ_EVAL = 3;		


/*
 * Precompute variances/(covariance-matrix-determinants) to simplify Mahalanobis distance
 * calculation.  Also, calculate 1/(det) for the original codebooks, based on the VQ vars.
 * (See comment in libmisc/vector.h.)
 */
static void subvq_maha_precomp (subvq_t *svq, float64 floor)
{
    int32 s;
    float32 *lrd;
    vector_gautbl_t *gautbl;
    
    E_INFO("Precomputing Mahalanobis distance invariants\n");
    
    lrd = (float32 *) ckd_calloc (svq->n_sv * svq->vqsize, sizeof(float32));
    
    for (s = 0; s < svq->n_sv; s++) {
	gautbl = &(svq->gautbl[s]);

	vector_gautbl_var_floor (gautbl, floor);
	
	gautbl->lrd = lrd;
	lrd += svq->vqsize;
	vector_gautbl_maha_precomp (gautbl);
    }
}


static void subvq_map_compact (subvq_t *vq, mgau_model_t *g)
{
    int32 r, c, c2, s;
    
    if (g) {
	if ((g->n_mgau != vq->origsize.r) || (g->max_comp != vq->origsize.c))
	    E_FATAL("Model size conflict: %d x %d (SubVQ) vs %d x %d (Original)\n",
		    vq->origsize.r, vq->origsize.c, g->n_mgau, g->max_comp);
    }
    
    /*
     * Compress map entries by removing invalid ones.  NOTE: There's not much of a consistency
     * check to ensure that the entries remaining do map correctly on to the valid entries in
     * the original (parent) mixture Gaussian model g.  The only check we have is to verify
     * the NUMBER of valid entries (components) in each mixture.
     */
    for (r = 0; r < vq->origsize.r; r++) {
	for (c = 0, c2 = 0; c < vq->origsize.c; c++) {
	    if (vq->map[r][c][0] < 0) {
		/* All ought to be < 0 */
		for (s = 1; s < vq->n_sv; s++) {
		    if (vq->map[r][c][s] >= 0)
			E_FATAL("Partially undefined map[%d][%d]\n", r, c);
		}
	    } else {
		if (c2 != c) {
		    for (s = 0; s < vq->n_sv; s++) {
			if (vq->map[r][c][s] < 0)
			    E_FATAL("Partially undefined map[%d][%d]\n", r, c);
			vq->map[r][c2][s] = vq->map[r][c][s];
		    }
		}
		c2++;
	    }
	}
	
	if (g && (c2 != mgau_n_comp (g, r)))
	    E_FATAL("Mixture %d: #Valid components conflict: %d (SubVQ) vs %d (Original)\n",
		    r, c2, mgau_n_comp(g,r));
	
	/* Invalidate the remaining map entries, for good measure */
	for (; c2 < vq->origsize.c; c2++) {
	    for (s = 0; s < vq->n_sv; s++)
		vq->map[r][c2][s] = -1;
	}
    }
}


/*
 * At this point, a map entries is an index within the subvector; i.e., map[r][c][s] points to
 * a subvq codeword within subvector s.  In evaluating a complete density using it subvq
 * component scores, these maps are used, with 2 lookups for each sub-vector.  To reduce the
 * lookups, linearize the map entries by viewing the subvq scores as a linear array rather than
 * a 2-D array.
 */
static void subvq_map_linearize (subvq_t *vq)
{
    int32 r, c, s;
    
    for (r = 0; r < vq->origsize.r; r++) {
	for (c = 0; (c < vq->origsize.c) && (vq->map[r][c][0] >= 0); c++) {
	    for (s = 0; s < vq->n_sv; s++)
		vq->map[r][c][s] = (s * vq->vqsize) + vq->map[r][c][s];
	}
    }
}


subvq_t *subvq_init (char *file, float64 varfloor, int32 max_sv, mgau_model_t *g)
{
    FILE *fp;
    char line[16384];
    int32 n_sv;		/* #Subvectors in file, as opposed to that used */
    int32 s, k, l, n, r, c;
    char *strp;
    subvq_t *vq;
    
    VQ_EVAL = cmd_ln_int32 ("-vqeval");	/* RAH, Currently only works when n_sv = 3, values computed but ignored in other cases */

    E_INFO("Loading Mixture Gaussian sub-VQ file '%s' (vq_eval: %d)\n", file,VQ_EVAL);
    
    vq = (subvq_t *) ckd_calloc (1, sizeof(subvq_t));
    
    fp = myfopen(file, "r");
    
    /* Read until "Sub-vectors" */
    for (;;) {
	if (fgets (line, sizeof(line), fp) == NULL)
	    E_FATAL("Failed to read VQParam header\n");
	if (sscanf (line, "VQParam %d %d -> %d %d",
		    &(vq->origsize.r), &(vq->origsize.c), &(vq->n_sv), &(vq->vqsize)) == 4)
	    break;
    }
    
    if (g) {
	if ((g->n_mgau != vq->origsize.r) || (g->max_comp != vq->origsize.c))
	    E_FATAL("Model size conflict: %d x %d (SubVQ) vs %d x %d (Original)\n",
		    vq->origsize.r, vq->origsize.c, g->n_mgau, g->max_comp);
    }
    
    if (max_sv < 0)
      max_sv = vq->n_sv;
    if (max_sv < vq->n_sv)
      E_INFO("Using %d subvectors out of %d\n", max_sv, vq->n_sv);
    else if (max_sv > vq->n_sv) {
      E_WARN("#Subvectors specified(%d) > available(%d); using latter\n", max_sv, vq->n_sv);
      max_sv = vq->n_sv;
    }
    
    n_sv = vq->n_sv;
    vq->n_sv = max_sv;
    if (vq->n_sv < VQ_EVAL)	/* RAH, 5.9.01, sanity check to make sure VQ_EVAL isn't higher than the n_sv */
      VQ_EVAL = vq->n_sv;
    vq->featdim = (int32 **) ckd_calloc (vq->n_sv, sizeof(int32 *));
    vq->gautbl = (vector_gautbl_t *) ckd_calloc (vq->n_sv, sizeof(vector_gautbl_t));
    vq->map = (int32 ***) ckd_calloc_3d (vq->origsize.r, vq->origsize.c, vq->n_sv,
					 sizeof(int32));
    
    /* Read subvector sizes and feature dimension maps */
    for (s = 0; s < n_sv; s++) {
	if ((fgets (line, sizeof(line), fp) == NULL) ||
	    (sscanf (line, "Subvector %d length %d%n", &k, &l, &n) != 2) || (k != s))
	    E_FATAL("Error reading length(subvector %d)\n", s);
	
	if (s < vq->n_sv) {
	    vq->gautbl[s].veclen = l;
	    vq->featdim[s] = (int32 *) ckd_calloc (vq->gautbl[s].veclen, sizeof(int32));
	
	    for (strp = line+n, c = 0; c < vq->gautbl[s].veclen; c++) {
		if (sscanf (strp, "%d%n", &(vq->featdim[s][c]), &n) != 1)
		    E_FATAL("Error reading subvector(%d).featdim(%d)\n", s, c);
		strp += n;
	    }
	    
	    vector_gautbl_alloc (&(vq->gautbl[s]), vq->vqsize, vq->gautbl[s].veclen);
	}
    }
    
    /* Echo info for sanity check */
    E_INFO("Original #codebooks(states)/codewords: %d x %d\n", vq->origsize.r, vq->origsize.c);
    E_INFO("Subvectors: %d, VQsize: %d\n", vq->n_sv, vq->vqsize);
    for (s = 0; s < vq->n_sv; s++) {
	E_INFO("SV %d feature dims(%d): ", s, vq->gautbl[s].veclen);
	for (c = 0; c < vq->gautbl[s].veclen; c++)
	    fprintf (stdout, " %2d", vq->featdim[s][c]);
	fprintf (stdout, "\n");
    }
    //SPEC fflush (stderr);
    
    /* Read VQ codebooks and maps for each subvector */
    for (s = 0; s < n_sv; s++) {
	E_INFO("Reading subvq %d%s\n", s, (s < vq->n_sv) ? "" : " (not used)");
	
	E_INFO("Reading codebook\n");
	if ((fgets (line, sizeof(line), fp) == NULL) ||
	    (sscanf (line, "Codebook %d", &k) != 1) || (k != s))
	    E_FATAL("Error reading codebook header\n", s);
	
	for (r = 0; r < vq->vqsize; r++) {
	    if (fgets (line, sizeof(line), fp) == NULL)
		E_FATAL("Error reading row(%d)\n", r);
	    
	    if (s >= vq->n_sv)
		continue;
	    
	    for (strp = line, c = 0; c < vq->gautbl[s].veclen; c++) {
		if (sscanf (strp, "%f %f%n",
			    &(vq->gautbl[s].mean[r][c]), &(vq->gautbl[s].var[r][c]), &k) != 2)
		    E_FATAL("Error reading row(%d) col(%d)\n", r, c);
		strp += k;
	    }
	}
	
	E_INFO("Reading map\n");
	if ((fgets (line, sizeof(line), fp) == NULL) ||
	    (sscanf (line, "Map %d", &k) != 1) || (k != s))
	    E_FATAL("Error reading map header\n", s);
	
	for (r = 0; r < vq->origsize.r; r++) {
	    if (fgets (line, sizeof(line), fp) == NULL)
		E_FATAL("Error reading row(%d)\n", r);
	    
	    if (s >= vq->n_sv)
		continue;
	    
	    for (strp = line, c = 0; c < vq->origsize.c; c++) {
		if (sscanf (strp, "%d%n", &(vq->map[r][c][s]), &k) != 1)
		    E_FATAL("Error reading row(%d) col(%d)\n", r, c);
		strp += k;
	    }
	}
	
	fflush (stdout);
    }
    
    if ((fscanf (fp, "%s", line) != 1) || (strcmp (line, "End") != 0))
	E_FATAL("Error reading 'End' token\n");
    
    fclose (fp);
    
    subvq_maha_precomp (vq, varfloor);
    subvq_map_compact (vq, g);
    subvq_map_linearize (vq);
    
    n = 0;
    for (s = 0; s < n_sv; s++) {
	if (vq->gautbl[s].veclen > n)
	    n = vq->gautbl[s].veclen;
    }
    assert (n > 0);
    vq->subvec = (float32 *) ckd_calloc (n, sizeof(float32));
    vq->vqdist = (int32 **) ckd_calloc_2d (vq->n_sv, vq->vqsize, sizeof(int32));
    vq->gauscore = (int32 *) ckd_calloc (vq->origsize.c, sizeof(int32));
    vq->mgau_sl = (int32 *) ckd_calloc (vq->origsize.c + 1, sizeof(int32));
    
    return vq;
}




/*
 * Based on previously computed subvq scores (Mahalanobis distances), determine the active
 * components in the given mixture (using the vq->map).
 * Return value: #Candidates in the returned shortlist.
 */
int32 subvq_mgau_shortlist (subvq_t *vq,
				   int32 m,	/* In: Mixture index */
				   int32 n,	/* In: #Components in specified mixture */
				   int32 beam)	/* In: Threshold to select active components */
{
    int32 *gauscore;
    int32 *map;
    int32 i, v, bv, th, nc;
    int32 *sl;
    int32 *vqdist;
    int32 sv_id;
    
    vqdist = vq->vqdist[0];	/* Since map is linearized for efficiency, must also
				   look at vqdist[][] as vqdist[] */
    gauscore = vq->gauscore;
    sl = vq->mgau_sl;
    /* Special case when vq->n_sv == 3; for speed */
    map = vq->map[m][0];
    bv = MAX_NEG_INT32;
    
    switch (vq->n_sv) {
    case 3:
	for (i = 0; i < n; i++) {
	  if (VQ_EVAL == 1) {
	    v = (int32) vqdist[*map];/* If we are not weighting the cep values, we need to adjust the subvqbeam */
	    map += 3;
	  } else {
	    /* RAH, we are ignoring the delta-delta, scoring the delta twice, strangely this works better than weighting the scores  */
	    /* I believe it has to do with the beam widths */
	    if (VQ_EVAL == 2) {
	    v = vqdist[*(map++)];
	      v += 2 * vqdist[*map]; /* RAH Count delta twice, we can keep the same subvqbeam as vq_eval = 3 if we double the delta*/
	      map += 2;
	    } else {
	      v = vqdist[*(map++)];/* Standard way */
	      v += vqdist[*(map++)]; /*  */
	      v += vqdist[*(map++)]; /*  */
	    }
	  }

	  gauscore[i] = v;
	    
	  if (bv < v)
	    bv = v;
	}
	break;
    case 2:
	for (i = 0; i < n; i++) {
	    v = vqdist[*(map++)];
	    v += vqdist[*(map++)];
	    gauscore[i] = v;
	   
	    if (bv < v)
		bv = v;
	}
	break;
    case 1:
      for (i = 0; i < n; i++) {
	v = vqdist[*(map++)];
	gauscore[i] = v;
	
	if (bv < v)
	  bv = v;
      }
      break;
    default:
      for (i = 0; i < n; i++) {
	v=0;
	for(sv_id =0 ; sv_id<vq->n_sv; sv_id++){
	  v+=vqdist[*(map++)];
	}
	gauscore[i] = v;
	
	if (bv < v)
	  bv = v;
      }
    }
    
    th = bv + beam;
    nc = 0;
    for (i = 0; i < n; i++) {
      if (gauscore[i] >= th)
	sl[nc++] = i;
    }
    sl[nc] = -1;
    
    return nc;
}


void subvq_subvec_eval_logs3 (subvq_t *vq, float32 *feat, int32 s)
{
    int32 i;
    int32 *featdim;
    
    /* Extract subvector from feat */
    featdim = vq->featdim[s];
    for (i = 0; i < vq->gautbl[s].veclen; i++)
	vq->subvec[i] = feat[featdim[i]];
    
    /* Evaluate distances between extracted subvector and corresponding codebook */
    vector_gautbl_eval_logs3(&(vq->gautbl[s]), 0, vq->vqsize, vq->subvec, vq->vqdist[s]);
}


void subvq_gautbl_eval_logs3 (subvq_t *vq, float32 *feat)
{
    int32 s, i;
    int32 *featdim;
    
    for (s = 0; s < vq->n_sv; s++) {
	/* Extract subvector from feat */
	featdim = vq->featdim[s];
	for (i = 0; i < vq->gautbl[s].veclen; i++)
	    vq->subvec[i] = feat[featdim[i]];
	
	/* Evaluate distances between extracted subvector and corresponding codebook */
    /* RAH, only evaluate the first VQ_EVAL set of features */
    if (s < VQ_EVAL) 
	vector_gautbl_eval_logs3(&(vq->gautbl[s]), 0, vq->vqsize, vq->subvec, vq->vqdist[s]);
    }
}


int32 subvq_frame_eval (subvq_t *vq, mgau_model_t *g, int32 beam, float32 *feat,
			int32 *sen_active, int32 *senscr)
{
  int32 s;
  int32 best, ns, ng;
  
  best = MAX_NEG_INT32;
  ns = 0;
  ng = 0;
  if (! vq) {
    /* No subvq model, use the original (SLOW!!) */
    for (s = 0; s < g->n_mgau; s++) {
      if ((! sen_active) || sen_active[s]) {
	senscr[s] = mgau_eval (g, s, NULL, feat);
	if (best < senscr[s])
	  best = senscr[s];
	ns++;
	ng += mgau_n_comp (g, s);
      } else
	senscr[s] = S3_LOGPROB_ZERO;
    }
  } else {
    /* Evaluate subvq model for given feature vector */
    subvq_gautbl_eval_logs3 (vq, feat);
    
    /* Find mixture component shortlists using subvq scores, and evaluate senones */
    for (s = 0; s < g->n_mgau; s++) {
      if ((! sen_active) || sen_active[s]) {
	ng += subvq_mgau_shortlist (vq, s, mgau_n_comp(g,s), beam);
	
	senscr[s] = mgau_eval (g, s, vq->mgau_sl, feat);
	if (best < senscr[s])
	  best = senscr[s];
	
	ns++;
      } else
	senscr[s] = S3_LOGPROB_ZERO;
    }
  }
    
    /* Normalize senone scores */
    for (s = 0; s < g->n_mgau; s++)
	senscr[s] -= best;
    
    g->frm_sen_eval = ns;
    g->frm_gau_eval = ng;
    
    return best;
}

int32 subvq_mgau_eval (mgau_model_t *g, subvq_t *vq, int32 m, int32 n, int32 *active)
{
    mgau_t *mgau;
    int32 *map;
    int32 i,j, v, sv_id;
    int32 c;
    int32 *vqdist;
    int32 score;
	int32 last_active;

    float64 f;
    f = log_to_logs3_factor();

    vqdist = vq->vqdist[0];	
    score = S3_LOGPROB_ZERO;
    mgau = &(g->mgau[m]);
    map = vq->map[m][0];

    if(!active){
      for (i = 0; i < n; i++) {
	v=0;
	for(sv_id =0 ; sv_id<vq->n_sv; sv_id++){
	  v+=vqdist[*(map++)];
	}
	score = logs3_add (score, v + mgau->mixw[i]);
	/*v=vqdist[*(map++)];
	  v+=vqdist[*(map++)];
	  v+=vqdist[*(map++)];*/
      }
    }else{

      last_active=0;
      for (i = 0; active[i] >=0; i++) {
	c=active[i];
      }
      for (i = 0; active[i] >=0; i++) {

	/*	E_INFO("Value of c %d\n",c);
		E_INFO("Value of last_active %d\n",last_active);*/
	c=active[i];
	map+=(c-last_active)*vq->n_sv;
	v=0;
	for(sv_id =0 ; sv_id<vq->n_sv; sv_id++){
	  v+=vqdist[*(map++)];
	}
	/*v=vqdist[*(map++)];
	v+=vqdist[*(map++)];
	v+=vqdist[*(map++)];*/

	last_active=c+1;

	score = logs3_add (score, v + mgau->mixw[i]);
      }
    }

    if(score == S3_LOGPROB_ZERO){
      E_INFO("Warning!! Score is S3_LOGPROB_ZERO\n");
    }

    return score;

}


/* RAH, free memory allocated by subvq_init() */
void subvq_free (subvq_t *s)
{
  int i;

  if (s) {
    
    for (i=0;i<s->n_sv;i++) {
      //      vector_gautbl_free (&(s->gautbl[i]));
      if (s->featdim[i]) ckd_free ((void *) s->featdim[i]);
    }


    if (s->featdim) 
      ckd_free ((void *) s->featdim);

    /* Free gaussian table */
    if (s->gautbl) 
      ckd_free ((void *)s->gautbl);      


    /* Free map */
    if (s->map)
      ckd_free_3d ((void ***) s->map);

    if (s->subvec) 
      ckd_free ((void *) s->subvec);

    if (s->vqdist) 
      ckd_free_2d ((void **) s->vqdist);

    if (s->gauscore) 
      ckd_free ((void *) s->gauscore);

    if (s->mgau_sl) 
      ckd_free ((void *) s->mgau_sl);

	
    ckd_free ((void *)s);


  }
}

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	1.856825e+07	1.861962e+07	0.2

                /* ====================================================================
 * Copyright (c) 1999-2001 Carnegie Mellon University.  All rights
 * reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer. 
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 * This work was supported in part by funding from the Defense Advanced 
 * Research Projects Agency and the National Science Foundation of the 
 * United States of America, and the CMU Sphinx Speech Consortium.
 *
 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * ====================================================================
 *
 */
/*
 * vector.c
 *
 * **********************************************
 * CMU ARPA Speech Project
 *
 * Copyright (c) 1997 Carnegie Mellon University.
 * ALL RIGHTS RESERVED.
 * **********************************************
 * 
 * HISTORY
 * 
 * 10-Mar-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Added vector_accum(), vector_vqlabel(), and vector_vqgen().
 * 
 * 09-Mar-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Added vector_is_zero(), vector_cmp(), and vector_dist_eucl().
 * 		Changed the name vector_dist_eval to vector_dist_maha.
 * 
 * 07-Oct-98	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Added distance computation related functions.
 * 
 * 12-Nov-95	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Copied from Eric Thayer.
 */


#include "vector.h"
#include "logs3.h"
#include "s3types.h"

#ifdef SPEC_CPU
   #include "specrand.h"
#else
   #if (WIN32)
   /* RAH #include <random.h> */
   #include <time.h>		/* RAH */
   #endif
#endif 



float64 vector_sum_norm (float32 *vec, int32 len)
{
    float64 sum, f;
    int32 i;
    
    sum = 0.0;
    for (i = 0; i < len; i++)
	sum += vec[i];

    if (sum != 0.0) {
	f = 1.0 / sum;
	for (i = 0; i < len; i++)
	  /*	  vec[i] *=  f;		*/  /*  */
	  /* The compiler was complaining about the above line, to make sure we don't lose accuracy, use this  */
	  vec[i] = (float32) ((float64) vec[i] * (float64)  f);		
    }
    
    return sum;
}


void vector_floor (float32 *vec, int32 len, float64 flr)
{
    int32 i;

    for (i = 0; i < len; i++)
	if (vec[i] < flr)
	    vec[i] = (float32)flr;
}


void vector_nz_floor (float32 *vec, int32 len, float64 flr)
{
    int32 i;

    for (i = 0; i < len; i++)
	if ((vec[i] != 0.0) && (vec[i] < flr))
	    vec[i] = (float32)flr;
}


void vector_print(FILE *fp, vector_t v, int32 dim)
{
    int32 i;
    
    for (i = 0; i < dim; i++)
	fprintf (fp, " %11.4e", v[i]);
    fprintf (fp, "\n");
    fflush (fp);
}


int32 vector_is_zero (float32 *vec, int32 len)
{
    int32 i;
    
    for (i = 0; (i < len) && (vec[i] == 0.0); i++);
    return (i == len);	/* TRUE iff all mean values are 0.0 */
}


int32 vector_maxcomp_int32 (int32 *val, int32 len)
{
    int32 i, bi;
    
    bi = 0;
    for (i = 1; i < len; i++) {
	if (val[i] > val[bi])
	    bi = i;
    }
    return bi;
}


int32 vector_mincomp_int32 (int32 *val, int32 len)
{
    int32 i, bi;
    
    bi = 0;
    for (i = 1; i < len; i++) {
	if (val[i] < val[bi])
	    bi = i;
    }
    return bi;
}


int32 vector_maxcomp_float32 (float32 *val, int32 len)
{
    int32 i, bi;
    
    bi = 0;
    for (i = 1; i < len; i++) {
	if (val[i] > val[bi])
	    bi = i;
    }
    return bi;
}


int32 vector_mincomp_float32 (float32 *val, int32 len)
{
    int32 i, bi;
    
    bi = 0;
    for (i = 1; i < len; i++) {
	if (val[i] < val[bi])
	    bi = i;
    }
    return bi;
}


void vector_accum (float32 *dst, float32 *src, int32 len)
{
    int32 i;
    
    for (i = 0; i < len; i++)
	dst[i] += src[i];
}


int32 vector_cmp (float32 *v1, float32 *v2, int32 len)
{
    int32 i;
    
    for (i = 0; i < len; i++) {
	if (v1[i] < v2[i])
	    return -1;
	if (v1[i] > v2[i])
	    return 1;
    }
    
    return 0;
}


int32 vector_mean (float32 *mean, float32 **data, int32 n_vec, int32 n_dim)
{
    int32 i, j;
    float64 f;
    
    assert ((n_vec > 0) && (n_dim > 0));
    
    for (i = 0; i < n_dim; i++)
	mean[i] = 0.0;
    
    for (i = 0; i < n_vec; i++) {
	for (j = 0; j < n_dim; j++)
	    mean[j] += data[i][j];
    }
    
    f = 1.0/(float64)n_vec;
    for (i = 0; i < n_dim; i++)
	mean[i] *= (float32)f;
    
    return 0;
}


float64 vector_dist_eucl (float32 *v1, float32 *v2, int32 len)
{
    float64 d;
    int32 i;
    
    d = 0.0;
    for (i = 0; i < len; i++)
	d += (v1[i] - v2[i]) * (v1[i] - v2[i]);
    
    return d;
}


float64 vector_maha_precomp (float32 *var, int32 len)
{
    float64 det;
    int32 i;
    
    for (det = (float64)0.0, i = 0; i < len; i++) {	/* log(1.0/prod(var[i])) */
	det -= (float64) (log(var[i]));
	var[i] = (float32)(1.0 / (var[i] * 2.0));
    }
    det -= log(2.0 * PI) * len;
    
    return (det * 0.5);	/* sqrt */
}


float64 vector_dist_maha (float32 *vec, float32 *mean, float32 *varinv, float64 loginvdet,
			  int32 len)
{
    float64 dist, diff;
    int32 i;
    
    dist = loginvdet;
    for (i = 0; i < len; i++) {
	diff = (vec[i] - mean[i]);
	dist -= diff * diff * varinv[i];
    }
    
    return dist;
}


int32 vector_vqlabel (float32 *vec, float32 **mean, int32 rows, int32 cols, float64 *sqerr)
{
    int32 i, besti;
    float64 d, bestd;
    
    assert ((rows > 0) && (cols > 0));
    
    bestd = vector_dist_eucl (mean[0], vec, cols);
    besti = 0;
    
    for (i = 1; i < rows; i++) {
	d = vector_dist_eucl (mean[i], vec, cols);
	if (bestd > d) {
	    bestd = d;
	    besti = i;
	}
    }
    
    if (sqerr)
	*sqerr = bestd;
    
    return besti;
}


float64 vector_vqgen (float32 **data, int32 rows, int32 cols, int32 vqrows,
		      float64 epsilon, int32 maxiter,
		      float32 **mean, int32 *map)
{
    int32 i, j, r, it;
#   ifdef SPEC_CPU
       int seed = 1;
#   else
       static uint32 seed = 1;
#   endif

    float64 sqerr, prev_sqerr=0, t;
    bitvec_t sel;
    int32 *count;
    float32 *gmean;
    ptmr_t tm;
    
    assert ((rows >= vqrows) && (maxiter >= 0) && (epsilon > 0.0));
    
    sel = bitvec_alloc (rows);
    
    ptmr_init (&tm);
    ptmr_start (&tm);
    
    /* Pick a random initial set of centroids */
#   ifdef SPEC_CPU
      spec_srand(seed);
#   else
#      ifndef WIN32 /* RAH */
          srandom (seed);
          seed ^= random();
#      else  /* RAH */
            srand ((unsigned) time(NULL)); /* RAH */
#      endif
#   endif  /* SPEC_CPU */

    for (i = 0; i < vqrows; i++) {
	/* Find r = a random, previously unselected row from the input */

#       ifdef SPEC_CPU
            r = ((int) (spec_rand() * (int32)0x7fffffff)) % rows;
#       else
#          ifndef WIN32			/* RAH */
	      r = (random() & (int32)0x7fffffff) % rows;
#          else  /* RAH */
	     r = (rand() & (int32)0x7fffffff) % rows; /* RAH */
#          endif /* RAH */
#      endif /* SPEC_CPU */

	while (bitvec_is_set (sel, r)) {	/* BUG: possible infinite loop!! */
	    if (++r >= rows)
		r = 0;
	}
	bitvec_set (sel, r);
	
	memcpy ((void *)(mean[i]), (void *)(data[r]), cols * sizeof(float32));
	/* BUG: What if two randomly selected rows are identical in content?? */
    }
    bitvec_free (sel);
    
    count = (int32 *) ckd_calloc (vqrows, sizeof(int32));
    
    /* In k-means, unmapped means in any iteration are a problem.  Replace them with gmean */
    gmean = (float32 *) ckd_calloc (cols, sizeof(float32));
    vector_mean (gmean, mean, vqrows, cols);

    for (it = 0;; it++) {		/* Iterations of k-means algorithm */
	/* Find the current data->mean mappings (labels) */
	sqerr = 0.0;
	for (i = 0; i < rows; i++) {
	    map[i] = vector_vqlabel (data[i], mean, vqrows, cols, &t);
	    sqerr += t;
	}
	ptmr_stop(&tm);
	
	if (it == 0)
	    E_INFO("Iter %4d: %.1fs CPU; sqerr= %e\n", it, tm.t_cpu, sqerr);
	else
	    E_INFO("Iter %4d: %.1fs CPU; sqerr= %e; delta= %e\n",
		   it, tm.t_cpu, sqerr, (prev_sqerr-sqerr)/prev_sqerr);
	
	/* Check if exit condition satisfied */
	if ((sqerr == 0.0) || (it >= maxiter-1) ||
	    ((it > 0) && ( ((prev_sqerr - sqerr) / prev_sqerr) < epsilon )) )
	    break;
	prev_sqerr = sqerr;
	
	ptmr_start(&tm);
	
	/* Update (reestimate) means */
	for (i = 0; i < vqrows; i++) {
	    for (j = 0; j < cols; j++)
		mean[i][j] = 0.0;
	    count[i] = 0;
	}
	for (i = 0; i < rows; i++) {
	    vector_accum (mean[map[i]], data[i], cols);
	    count[map[i]]++;
	}
	for (i = 0; i < vqrows; i++) {
	    if (count[i] > 1) {
		t = 1.0 / (float64)(count[i]);
		for (j = 0; j < cols; j++)
		  /*		  mean[i][j] *= t; */ /* RAH, compiler was complaining about this,  */
		  mean[i][j] = (float32) ((float64) mean[i][j] * (float64) t); /*  */
	    } else if (count[i] == 0) {
		E_ERROR("Iter %d: mean[%d] unmapped\n", it, i);
		memcpy (mean[i], gmean, cols * sizeof(float32));
	    }
	}
    }
    
    ckd_free (count);
    ckd_free (gmean);
    
    return sqerr;
}


float64 vector_pdf_entropy (float32 *p, int32 len)
{
    float64 sum;
    int32 i;
    
    sum = 0.0;
    for (i = 0; i < len; i++) {
	if (p[i] > 0.0)
	    sum -= p[i] * log(p[i]);
    }
    sum /= log(2.0);
    
    return sum;
}


float64 vector_pdf_cross_entropy (float32 *p1, float32 *p2, int32 len)
{
    float64 sum;
    int32 i;
    
    sum = 0.0;
    for (i = 0; i < len; i++) {
	if (p2[i] > 0.0)
	    sum -= p1[i] * log(p2[i]);
    }
    sum /= log(2.0);
    
    return sum;
}


void vector_gautbl_alloc (vector_gautbl_t *gautbl, int32 n_gau, int32 veclen)
{
    gautbl->n_gau = n_gau;
    gautbl->veclen = veclen;
    gautbl->mean = (float32 **) ckd_calloc_2d (n_gau, veclen, sizeof(float32));
    gautbl->var = (float32 **) ckd_calloc_2d (n_gau, veclen, sizeof(float32));
    gautbl->lrd = (float32 *) ckd_calloc (n_gau, sizeof(float32));
    gautbl->distfloor = logs3_to_log (S3_LOGPROB_ZERO);
}


void vector_gautbl_free (vector_gautbl_t *gautbl)
{
    ckd_free_2d ((void **) gautbl->mean);
    ckd_free_2d ((void **) gautbl->var);
    ckd_free ((void *) gautbl->lrd);
}


void vector_gautbl_var_floor (vector_gautbl_t *gautbl, float64 floor)
{
    int32 g;
    
    for (g = 0; g < gautbl->n_gau; g++)
	vector_floor (gautbl->var[g], gautbl->veclen, floor);
}


void vector_gautbl_maha_precomp (vector_gautbl_t *gautbl)
{
    int32 g;
    
    for (g = 0; g < gautbl->n_gau; g++)
	gautbl->lrd[g] = (float32) vector_maha_precomp (gautbl->var[g], gautbl->veclen);
}


void vector_gautbl_eval_logs3 (vector_gautbl_t *gautbl,
			       int32 offset,
			       int32 count,
			       float32 *x,
			       int32 *score)
{
    int32 i, r;
    float64 f;
    int32 end, veclen;
    float32 *m1, *m2, *v1, *v2;
    float64 dval1, dval2, diff1, diff2;
    
    f = log_to_logs3_factor();
    
    /* Interleave evaluation of two vectors at a time for speed on pipelined machines */
    end = offset + count;
    veclen = gautbl->veclen;
    
    for (r = offset; r < end-1; r += 2) {
	m1 = gautbl->mean[r];
	m2 = gautbl->mean[r+1];
	v1 = gautbl->var[r];
	v2 = gautbl->var[r+1];
	dval1 = gautbl->lrd[r];
	dval2 = gautbl->lrd[r+1];

	for (i = 0; i < veclen; i++) {
	    diff1 = x[i] - m1[i];
	    dval1 -= diff1 * diff1 * v1[i];
	    diff2 = x[i] - m2[i];
	    dval2 -= diff2 * diff2 * v2[i];
	}
	
	if (dval1 < gautbl->distfloor)
	    dval1 = gautbl->distfloor;
	if (dval2 < gautbl->distfloor)
	    dval2 = gautbl->distfloor;

	score[r] = (int32)(f * dval1);
	score[r+1] = (int32)(f * dval2);
    }
    
    if (r < end) {
	m1 = gautbl->mean[r];
	v1 = gautbl->var[r];
	dval1 = gautbl->lrd[r];
	
	for (i = 0; i < veclen; i++) {
	    diff1 = x[i] - m1[i];
	    dval1 -= diff1 * diff1 * v1[i];
	}
	
	if (dval1 < gautbl->distfloor)
	    dval1 = gautbl->distfloor;

	score[r] = (int32)(f * dval1);
    }
}

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
97	### 1	3031.43069306931	7.652000e+03	1.616000e+03	78.8
77	### 2	3817.9792746114	7.882000e+03	1.544000e+03	80.4
94	### 3	3315.96216216216	7.002000e+03	1.480000e+03	78.8
155	### 4	2151.24645892351	7.416000e+03	1.412000e+03	80.9

                /* ====================================================================
 * Copyright (c) 1999-2001 Carnegie Mellon University.  All rights
 * reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer. 
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 * This work was supported in part by funding from the Defense Advanced 
 * Research Projects Agency and the National Science Foundation of the 
 * United States of America, and the CMU Sphinx Speech Consortium.
 *
 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * ====================================================================
 *
 */
/*
 * mdef.c -- HMM model definition: base (CI) phones and triphones
 *
 * **********************************************
 * CMU ARPA Speech Project
 *
 * Copyright (c) 1999 Carnegie Mellon University.
 * ALL RIGHTS RESERVED.
 * **********************************************
 * 
 * HISTORY
 * 
 * 19.Apr-2001  Ricky Houghton, added code for free allocated memory
 *
 * 14-Oct-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon
 * 		Added mdef_sseq2sen_active().
 * 
 * 06-May-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon
 * 		In mdef_phone_id(), added backing off to silence phone context from filler
 * 		context if original triphone not found.
 * 
 * 30-Apr-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon
 * 		Added senone-sequence id (ssid) to phone_t and appropriate functions to
 * 		maintain it.  Instead, moved state sequence info to mdef_t.
 * 
 * 13-Jul-96	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Added mdef_phone_str().
 * 
 * 01-Jan-96	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Allowed mdef_phone_id_nearest to return base phone id if either
 * 		left or right context (or both) is undefined.
 * 
 * 01-Jan-96	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Created.
 */


/*
 * Major assumptions:
 *   All phones have same #states, same topology.
 *   Every phone has exactly one non-emitting, final state--the last one.
 *   CI phones must appear first in model definition file.
 */


#include "mdef.h"


#define MODEL_DEF_VERSION	"0.3"


void mdef_dump (FILE *fp, mdef_t *m)
{
    int32 i, j;
    int32 ssid;
    char buf[1024];
    
    fprintf (fp, "%d ciphone\n", m->n_ciphone);
    fprintf (fp, "%d phone\n", m->n_phone);
    fprintf (fp, "%d emitstate\n", m->n_emit_state);
    fprintf (fp, "%d cisen\n", m->n_ci_sen);
    fprintf (fp, "%d sen\n", m->n_sen);
    fprintf (fp, "%d tmat\n", m->n_tmat);
    
    for (i = 0; i < m->n_phone; i++) {
	mdef_phone_str (m, i, buf);
	ssid = m->phone[i].ssid;
	
	fprintf (fp, "%3d %5d", m->phone[i].tmat, ssid);
	for (j = 0; j < m->n_emit_state; j++)
	    fprintf (fp, " %5d", m->sseq[ssid][j]);
	fprintf (fp, "\t");
	for (j = 0; j < m->n_emit_state; j++)
	    fprintf (fp, " %3d", m->cd2cisen[m->sseq[ssid][j]]);
	fprintf (fp, "\t%s\n", buf);
    }
    
    fflush (fp);
}


#if 0
int32 mdef_hmm_cmp (mdef_t *m, s3pid_t p1, s3pid_t p2)
{
    int32 i;
    
    if (m->phone[p1].tmat != m->phone[p2].tmat)
	return -1;
    
    for (i = 0; i < m->n_emit_state; i++)
	if (m->phone[p1].state[i] != m->phone[p2].state[i])
	    return -1;
    
    return 0;
}
#endif


static void ciphone_add (mdef_t *m, char *ci, s3pid_t p)
{
    assert (p < m->n_ciphone);

    m->ciphone[p].name = (char *) ckd_salloc (ci); /* freed in mdef_free */
    if (hash_enter (m->ciphone_ht, m->ciphone[p].name, p) != p)
	E_FATAL("hash_enter(%s) failed; duplicate CIphone?\n", m->ciphone[p].name);
}


static ph_lc_t *find_ph_lc (ph_lc_t *lclist, s3cipid_t lc)
{
    ph_lc_t *lcptr;

    for (lcptr = lclist; lcptr && (lcptr->lc != lc); lcptr = lcptr->next);
    return lcptr;
}


static ph_rc_t *find_ph_rc (ph_rc_t *rclist, s3cipid_t rc)
{
    ph_rc_t *rcptr;

    for (rcptr = rclist; rcptr && (rcptr->rc != rc); rcptr = rcptr->next);
    return rcptr;
}


static void triphone_add (mdef_t *m,
			  s3cipid_t ci, s3cipid_t lc, s3cipid_t rc, word_posn_t wpos,
			  s3pid_t p)
{
    ph_lc_t *lcptr;
    ph_rc_t *rcptr;
    
    assert (p < m->n_phone);

    /* Fill in phone[p] information (state and tmat mappings added later) */
    m->phone[p].ci = ci;
    m->phone[p].lc = lc;
    m->phone[p].rc = rc;
    m->phone[p].wpos = wpos;

    /* Create <ci,lc,rc,wpos> -> p mapping if not a CI phone */
    if (p >= m->n_ciphone) {
	if ((lcptr = find_ph_lc (m->wpos_ci_lclist[wpos][(int)ci], lc)) 
	    == NULL)
	{
	lcptr = (ph_lc_t *) ckd_calloc (1, sizeof(ph_lc_t)); /* freed at mdef_free, I believe */
	    lcptr->lc = lc;
	    lcptr->next = m->wpos_ci_lclist[wpos][(int)ci];
	    m->wpos_ci_lclist[wpos][(int)ci] = lcptr; /* This is what needs to be freed */
	}
	if ((rcptr = find_ph_rc (lcptr->rclist, rc)) != NULL) {
	    char buf[4096];
	    
	    mdef_phone_str (m, rcptr->pid, buf);
	    E_FATAL("Duplicate triphone: %s\n", buf);
	}
	
      rcptr = (ph_rc_t *) ckd_calloc (1, sizeof(ph_rc_t)); /* freed in mdef_free, I believe */
	rcptr->rc = rc;
	rcptr->pid = p;
	rcptr->next = lcptr->rclist;
	lcptr->rclist = rcptr;
    }
}


s3cipid_t mdef_ciphone_id (mdef_t *m, char *ci)
{
    int32 id;
    
    assert (m);
    assert (ci);
    
    if (hash_lookup (m->ciphone_ht, ci, &id) < 0)
	return (BAD_S3CIPID);
    return ((s3cipid_t) id);
}


const char *mdef_ciphone_str (mdef_t *m, s3cipid_t id)
{
    assert (m);
    assert ((id >= 0) && (id < m->n_ciphone));
    
    return (m->ciphone[(int)id].name);
}


int32 mdef_phone_str (mdef_t *m, s3pid_t pid, char *buf)
{
    char *wpos_name;
    
    assert (m);
    assert ((pid >= 0) && (pid < m->n_phone));
    wpos_name = WPOS_NAME;
    
    buf[0] = '\0';
    if (pid < m->n_ciphone)
	sprintf (buf, "%s", mdef_ciphone_str (m, (s3cipid_t) pid));
    else {
	sprintf (buf, "(%s,%s,%s,%c)",
		 mdef_ciphone_str(m, m->phone[pid].ci),
		 mdef_ciphone_str(m, m->phone[pid].lc),
		 mdef_ciphone_str(m, m->phone[pid].rc),
		 wpos_name[m->phone[pid].wpos]);
    }
    return 0;
}


s3pid_t mdef_phone_id (mdef_t *m, 
		       s3cipid_t ci, s3cipid_t lc, s3cipid_t rc, word_posn_t wpos)
{
    ph_lc_t *lcptr;
    ph_rc_t *rcptr;
    s3cipid_t newl, newr;
    
    assert (m);
    assert ((ci >= 0) && (ci < m->n_ciphone));
    assert ((lc >= 0) && (lc < m->n_ciphone));
    assert ((rc >= 0) && (rc < m->n_ciphone));
    assert ((wpos >= 0) && (wpos < N_WORD_POSN));

    if (((lcptr = find_ph_lc (m->wpos_ci_lclist[wpos][(int)ci], lc)) == NULL) ||
	((rcptr = find_ph_rc (lcptr->rclist, rc)) == NULL)) {
	/* Not found; backoff to silence context if non-silence filler context */
	if (NOT_S3CIPID(m->sil))
	    return BAD_S3PID;
	
	newl = m->ciphone[(int)lc].filler ? m->sil : lc;
	newr = m->ciphone[(int)rc].filler ? m->sil : rc;
	if ((newl == lc) && (newr == rc))
	    return BAD_S3PID;
	
	return (mdef_phone_id (m, ci, newl, newr, wpos));
    }
    
    return (rcptr->pid);
}


s3pid_t mdef_phone_id_nearest (mdef_t *m, 
			       s3cipid_t b, s3cipid_t l, s3cipid_t r, word_posn_t pos)
{
    word_posn_t tmppos;
    s3pid_t p;
    s3cipid_t newl, newr;
    char *wpos_name;
    
    assert (m);
    assert ((b >= 0) && (b < m->n_ciphone));
    assert ((pos >= 0) && (pos < N_WORD_POSN));

    if ((NOT_S3CIPID(l)) || (NOT_S3CIPID(r)))
	return ((s3pid_t) b);
    
    assert ((l >= 0) && (l < m->n_ciphone));
    assert ((r >= 0) && (r < m->n_ciphone));
    
    p = mdef_phone_id (m, b, l, r, pos);
    if (IS_S3PID(p))
	return p;
    
    /* Exact triphone not found; backoff to other word positions */
    for (tmppos = 0; tmppos < N_WORD_POSN; tmppos++) {
	if (tmppos != pos) {
	    p = mdef_phone_id (m, b, l, r, tmppos);
	    if (IS_S3PID(p))
		return p;
	}
    }
    
    /* Nothing yet; backoff to silence phone if non-silence filler context */
    if (IS_S3CIPID(m->sil)) {
	newl = m->ciphone[(int)l].filler ? m->sil : l;
	newr = m->ciphone[(int)r].filler ? m->sil : r;
	if ((newl != l) || (newr != r)) {
	    p = mdef_phone_id (m, b, newl, newr, pos);
	    if (IS_S3PID(p))
		return p;
	    
	    for (tmppos = 0; tmppos < N_WORD_POSN; tmppos++) {
		if (tmppos != pos) {
		    p = mdef_phone_id (m, b, newl, newr, tmppos);
		    if (IS_S3PID(p))
			return p;
		}
	    }
	}
    }
    
    /* Nothing yet; backoff to base phone */
    if ((m->n_phone > m->n_ciphone) && (! m->ciphone[(int)b].filler)) {
	wpos_name = WPOS_NAME;
#if 0
	E_WARN("Triphone(%s,%s,%s,%c) not found; backing off to CIphone\n",
	       mdef_ciphone_str(m, b),
	       mdef_ciphone_str(m, l),
	       mdef_ciphone_str(m, r),
	       wpos_name[pos]);
#endif
    }
    return ((s3pid_t) b);
}


int32 mdef_phone_components (mdef_t *m,
			     s3pid_t p,
			     s3cipid_t *b,
			     s3cipid_t *l,
			     s3cipid_t *r,
			     word_posn_t *pos)
{
    assert (m);
    assert ((p >= 0) && (p < m->n_phone));

    *b = m->phone[p].ci;
    *l = m->phone[p].lc;
    *r = m->phone[p].rc;
    *pos = m->phone[p].wpos;

    return 0;
}


int32 mdef_is_ciphone (mdef_t *m, s3pid_t p)
{
    assert (m);
    assert ((p >= 0) && (p < m->n_phone));
    
    return ((p < m->n_ciphone) ? 1 : 0);
}


/* Parse tmat and state->senone mappings for phone p and fill in structure */
static void parse_tmat_senmap (mdef_t *m, char *line, int32 off, s3pid_t p)
{
    int32 wlen, n, s;
    char word[1024], *lp;

    lp = line + off;
    
    /* Read transition matrix id */
    if ((sscanf (lp, "%d%n", &n, &wlen) != 1) || (n < 0))
        E_FATAL("Missing or bad transition matrix id: %s\n", line);
    m->phone[p].tmat = n;
    if (m->n_tmat <= n)
	E_FATAL("tmat-id(%d) > #tmat in header(%d): %s\n", n, m->n_tmat, line);
    lp += wlen;
    
    /* Read senone mappings for each emitting state */
    for (n = 0; n < m->n_emit_state; n++) {
        if ((sscanf (lp, "%d%n", &s, &wlen) != 1) || (s < 0))
	    E_FATAL("Missing or bad state[%d]->senone mapping: %s\n", n, line);
	if ((p < m->n_ciphone) && (m->n_ci_sen <= s))
	    E_FATAL("CI-senone-id(%d) > #CI-senones(%d): %s\n", s, m->n_ci_sen, line);
	if (m->n_sen <= s)
	    E_FATAL("Senone-id(%d) > #senones(%d): %s\n", s, m->n_sen, line);

	m->sseq[p][n] = s;
	lp += wlen;
    }

    /* Check for the last non-emitting state N */
    if ((sscanf (lp, "%s%n", word, &wlen) != 1) || (strcmp (word, "N") != 0))
        E_FATAL("Missing non-emitting state spec: %s\n", line);
    lp += wlen;

    /* Check for end of line */
    if (sscanf (lp, "%s%n", word, &wlen) == 1)
        E_FATAL("Non-empty beyond non-emitting final state: %s\n", line);
}


static void parse_base_line (mdef_t *m, char *line, s3pid_t p)
{
    int32 wlen, n;
    char word[1024], *lp;
    s3cipid_t ci;

    lp = line;
    
    /* Read base phone name */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing base phone name: %s\n", line);
    lp += wlen;
    
    /* Make sure it's not a duplicate */
    ci = mdef_ciphone_id (m, word);
    if (IS_S3CIPID(ci))
        E_FATAL("Duplicate base phone: %s\n", line);

    /* Add ciphone to ciphone table with id p */
    ciphone_add (m, word, p);
    ci = (s3cipid_t) p;

    /* Read and skip "-" for lc, rc, wpos */
    for (n = 0; n < 3; n++) {
	if ((sscanf (lp, "%s%n", word, &wlen) != 1) || (strcmp (word, "-") != 0))
	    E_FATAL("Bad context info for base phone: %s\n", line);
	lp += wlen;
    }
    
    /* Read filler attribute, if present */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing filler atribute field: %s\n", line);
    lp += wlen;
    if (strcmp (word, "filler") == 0)
        m->ciphone[(int)ci].filler = 1;
    else if (strcmp (word, "n/a") == 0)
        m->ciphone[(int)ci].filler = 0;
    else
        E_FATAL("Bad filler attribute field: %s\n", line);

    triphone_add (m, ci, BAD_S3CIPID, BAD_S3CIPID, WORD_POSN_UNDEFINED, p);

    /* Parse remainder of line: transition matrix and state->senone mappings */
    parse_tmat_senmap (m, line, lp-line, p);
}


static void parse_tri_line (mdef_t *m, char *line, s3pid_t p)
{
    int32 wlen;
    char word[1024], *lp;
    s3cipid_t ci, lc, rc;
    word_posn_t wpos = WORD_POSN_BEGIN;

    lp = line;
    
    /* Read base phone name */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing base phone name: %s\n", line);
    lp += wlen;

    ci = mdef_ciphone_id (m, word);
    if (NOT_S3CIPID(ci))
        E_FATAL("Unknown base phone: %s\n", line);

    /* Read lc */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing left context: %s\n", line);
    lp += wlen;
    lc = mdef_ciphone_id (m, word);
    if (NOT_S3CIPID(lc))
        E_FATAL("Unknown left context: %s\n", line);

    /* Read rc */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing right context: %s\n", line);
    lp += wlen;
    rc = mdef_ciphone_id (m, word);
    if (NOT_S3CIPID(rc))
        E_FATAL("Unknown right  context: %s\n", line);
    
    /* Read tripone word-position within word */
    if ((sscanf (lp, "%s%n", word, &wlen) != 1) || (word[1] != '\0'))
        E_FATAL("Missing or bad word-position spec: %s\n", line);
    lp += wlen;
    switch (word[0]) {
    case 'b': wpos = WORD_POSN_BEGIN; break;
    case 'e': wpos = WORD_POSN_END; break;
    case 's': wpos = WORD_POSN_SINGLE; break;
    case 'i': wpos = WORD_POSN_INTERNAL; break;
    default: E_FATAL("Bad word-position spec: %s\n", line);
    }

    /* Read filler attribute, if present.  Must match base phone attribute */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing filler attribute field: %s\n", line);
    lp += wlen;
    if (((strcmp (word, "filler") == 0) && (m->ciphone[(int)ci].filler)) ||
	((strcmp (word, "n/a") == 0) && (! m->ciphone[(int)ci].filler))) {
	/* Everything is fine */
    } else
        E_FATAL("Bad filler attribute field: %s\n", line);
    
    triphone_add (m, ci, lc, rc, wpos, p);

    /* Parse remainder of line: transition matrix and state->senone mappings */
    parse_tmat_senmap (m, line, lp-line, p);
}


static void sseq_compress (mdef_t *m)
{
    hash_table_t *h;
    s3senid_t **sseq;
    int32 n_sseq;
    int32 p, j, k;
    glist_t g;
    gnode_t *gn;
    hash_entry_t *he;
    
    k = m->n_emit_state * sizeof(s3senid_t);
    
    h = hash_new (m->n_phone, HASH_CASE_YES);
    n_sseq = 0;
    
    /* Identify unique senone-sequence IDs.  BUG: tmat-id not being considered!! */
    for (p = 0; p < m->n_phone; p++) {
	/* Add senone sequence to hash table */
	if ((j = hash_enter_bkey (h, (char *)(m->sseq[p]), k, n_sseq)) == n_sseq)
	    n_sseq++;
	
	m->phone[p].ssid = j;
    }
    
    /* Generate compacted sseq table */
    sseq = (s3senid_t **) ckd_calloc_2d (n_sseq, m->n_emit_state, sizeof(s3senid_t));/* freed in mdef_free() */
    
    g = hash_tolist (h, &j);
    assert (j == n_sseq);
    
    for (gn = g; gn; gn = gnode_next(gn)) {
	he = (hash_entry_t *) gnode_ptr (gn);
	j = hash_entry_val(he);
	memcpy (sseq[j], hash_entry_key(he), k);
    }
    glist_free (g);
    
    /* Free the old, temporary senone sequence table, replace with compacted one */
    ckd_free_2d ((void **) m->sseq);
    m->sseq = sseq;
    m->n_sseq = n_sseq;
    
    hash_free (h);
}


static int32 noncomment_line(char *line, int32 size, FILE *fp)
{
    while (fgets (line, size, fp) != NULL) {
        if (line[0] != '#')
	    return 0;
    }
    return -1;
}


/*
 * Initialize phones (ci and triphones) and state->senone mappings from .mdef file.
 */
mdef_t *mdef_init (char *mdeffile)
{
    FILE *fp;
    int32 n_ci, n_tri, n_map, n;
    char tag[1024], buf[1024];
    s3senid_t **senmap;
    s3pid_t p;
    int32 s, ci, cd;
    mdef_t *m;
    int32 *cdsen_start, *cdsen_end;

    if (! mdeffile)
	E_FATAL("No mdef-file\n");

    E_INFO("Reading model definition: %s\n", mdeffile);

    m = (mdef_t *) ckd_calloc (1, sizeof(mdef_t)); /* freed in mdef_free */
    
    if ((fp = fopen(mdeffile, "r")) == NULL)
        E_FATAL_SYSTEM("fopen(%s,r) failed\n", mdeffile);

    if (noncomment_line(buf, sizeof(buf), fp) < 0)
        E_FATAL("Empty file: %s\n", mdeffile);

    if (strncmp(buf, MODEL_DEF_VERSION, strlen(MODEL_DEF_VERSION)) != 0)
        E_FATAL("Version error: Expecing %s, but read %s\n", MODEL_DEF_VERSION, buf);

    /* Read #base phones, #triphones, #senone mappings defined in header */
    n_ci = -1;
    n_tri = -1;
    n_map = -1;
    m->n_ci_sen = -1;
    m->n_sen = -1;
    m->n_tmat = -1;
    do {
	if (noncomment_line(buf, sizeof(buf), fp) < 0)
	    E_FATAL("Incomplete header\n");

	if ((sscanf(buf, "%d %s", &n, tag) != 2) || (n < 0))
	    E_FATAL("Error in header: %s\n", buf);

	if (strcmp(tag, "n_base") == 0)
	    n_ci = n;
	else if (strcmp(tag, "n_tri") == 0)
	    n_tri = n;
	else if (strcmp(tag, "n_state_map") == 0)
	    n_map = n;
	else if (strcmp(tag, "n_tied_ci_state") == 0)
	    m->n_ci_sen = n;
	else if (strcmp(tag, "n_tied_state") == 0)
	    m->n_sen = n;
	else if (strcmp(tag, "n_tied_tmat") == 0)
	    m->n_tmat = n;
	else
	    E_FATAL("Unknown header line: %s\n", buf);
    } while ((n_ci < 0) || (n_tri < 0) || (n_map < 0) ||
	     (m->n_ci_sen < 0) || (m->n_sen < 0) || (m->n_tmat < 0));

    if ((n_ci == 0) || (m->n_ci_sen == 0) || (m->n_tmat == 0) || (m->n_ci_sen > m->n_sen))
        E_FATAL("%s: Error in header\n", mdeffile);
    
    /* Check typesize limits */
    if (n_ci >= MAX_S3CIPID)
	E_FATAL("%s: #CI phones (%d) exceeds limit (%d)\n", mdeffile, n_ci, MAX_S3CIPID);
    if (n_ci + n_tri >= MAX_S3PID)
	E_FATAL("%s: #Phones (%d) exceeds limit (%d)\n", mdeffile, n_ci+n_tri, MAX_S3PID);
    if (m->n_sen >= MAX_S3SENID)
	E_FATAL("%s: #senones (%d) exceeds limit (%d)\n", mdeffile, m->n_sen, MAX_S3SENID);
    if (m->n_tmat >= MAX_S3TMATID)
	E_FATAL("%s: #tmats (%d) exceeds limit (%d)\n", mdeffile, m->n_tmat, MAX_S3TMATID);
    
    m->n_emit_state = (n_map / (n_ci+n_tri)) - 1;
    if ((m->n_emit_state+1) * (n_ci+n_tri) != n_map)
        E_FATAL("Header error: n_state_map not a multiple of n_ci*n_tri\n");

    /* Initialize ciphone info */
    m->n_ciphone = n_ci;
    m->ciphone_ht = hash_new (n_ci, 1);	/* With case-insensitive string names */ /* freed in mdef_free */
    m->ciphone = (ciphone_t *) ckd_calloc (n_ci, sizeof(ciphone_t)); /* freed in mdef_free */
    /* RAH, let's null the pointers so that we can reliably deallocate them */
    /*    for (i=0;i<m->n_ciphone;i++) { */ /*  */
    /*m->ciphone[i].name = NULL; */ /*  */
    /*} */ /*  */

    /* Initialize phones info (ciphones + triphones) */
    m->n_phone = n_ci + n_tri;
    m->phone = (phone_t *) ckd_calloc (m->n_phone, sizeof(phone_t)); /* freed in mdef_free */

    /* Allocate space for state->senone map for each phone */
    senmap = (s3senid_t **) ckd_calloc_2d (m->n_phone, m->n_emit_state, sizeof(s3senid_t));/* freed in mdef_free */
    m->sseq = senmap;	/* TEMPORARY; until it is compressed into just the unique ones */
    
    /* Allocate initial space for <ci,lc,rc,wpos> -> pid mapping */
    m->wpos_ci_lclist = (ph_lc_t ***) ckd_calloc_2d (N_WORD_POSN, m->n_ciphone, sizeof(ph_lc_t *)); /* freed in mdef_free */

    /*
     * Read base phones and triphones.  They'll simply be assigned a running sequence
     * number as their "phone-id".  If the phone-id < n_ci, it's a ciphone.
     */

    /* Read base phones */
    for (p = 0; p < n_ci; p++) {
        if (noncomment_line(buf, sizeof(buf), fp) < 0)
	    E_FATAL("Premature EOF reading CIphone %d\n", p);
        parse_base_line (m, buf, p);
    }
    m->sil = mdef_ciphone_id (m, S3_SILENCE_CIPHONE);
    
    /* Read triphones, if any */
    for (; p < m->n_phone; p++) {
        if (noncomment_line(buf, sizeof(buf), fp) < 0)
	    E_FATAL("Premature EOF reading phone %d\n", p);
        parse_tri_line (m, buf, p);
    }

    if (noncomment_line(buf, sizeof(buf), fp) >= 0)
	E_ERROR("Non-empty file beyond expected #phones (%d)\n", m->n_phone);

    /* Build CD senones to CI senones map */
    if (m->n_ciphone * m->n_emit_state != m->n_ci_sen)
	E_FATAL("#CI-senones(%d) != #CI-phone(%d) x #emitting-states(%d)\n",
		m->n_ci_sen, m->n_ciphone, m->n_emit_state);
    m->cd2cisen = (s3senid_t *) ckd_calloc (m->n_sen, sizeof(s3senid_t)); /* freed in mdef_free */

    m->sen2cimap = (s3cipid_t *) ckd_calloc (m->n_sen, sizeof(s3cipid_t)); /* freed in mdef_free */

    for (s = 0; s < m->n_sen; s++)
	m->sen2cimap[s] = BAD_S3CIPID;
    for (s = 0; s < m->n_ci_sen; s++) {		/* CI senones */
	m->cd2cisen[s] = (s3senid_t) s;
	m->sen2cimap[s] = s / m->n_emit_state;
    }
    for (p = n_ci; p < m->n_phone; p++) {	/* CD senones */
	for (s = 0; s < m->n_emit_state; s++) {
	    cd = m->sseq[p][s];
	    ci = m->sseq[(int)m->phone[p].ci][s];
	    m->cd2cisen[cd] = (s3senid_t) ci;
	    m->sen2cimap[cd] = m->phone[p].ci;
	}
    }
    
    /*
     * Count #senones (CI+CD) for each CI phone.
     * HACK!!  For handling holes in senone-CIphone mappings.  Does not work if holes
     * are present at the beginning or end of senones for a given CIphone.
     */
    cdsen_start = (int32 *) ckd_calloc (m->n_ciphone, sizeof(int32)); /* freed locally */

    cdsen_end = (int32 *) ckd_calloc (m->n_ciphone, sizeof(int32)); /* freed locally */

    for (s = m->n_ci_sen; s < m->n_sen; s++) {
	if (NOT_S3CIPID(m->sen2cimap[s]))
	    continue;
	
	if (! cdsen_start[(int)m->sen2cimap[s]])
	    cdsen_start[(int)m->sen2cimap[s]] = s;
	cdsen_end[(int)m->sen2cimap[s]] = s;
    }

    /* Fill up holes */
    for (s = m->n_ci_sen; s < m->n_sen; s++) {
	if (IS_S3CIPID(m->sen2cimap[s]))
	    continue;

	/* Check if properly inside the observed ranges above */
	for (p = 0; p < m->n_ciphone; p++) {
	    if ((s > cdsen_start[p]) && (s < cdsen_end[p]))
		break;
	}
	if (p >= m->n_ciphone)
	    E_FATAL("Unreferenced senone %d; cannot determine parent CIphone\n", s);
	m->sen2cimap[s] = p;
    }

    /* Build #CD-senones for each CIphone */
    m->ciphone2n_cd_sen = (int32 *) ckd_calloc (m->n_ciphone, sizeof(int32));/* freed mdef_free */
    n = 0;
    for (p = 0; p < m->n_ciphone; p++) {
	if (cdsen_start[p] > 0) {
	    m->ciphone2n_cd_sen[p] = cdsen_end[p] - cdsen_start[p] + 1;
	    n += m->ciphone2n_cd_sen[p];
	}
    }
    n += m->n_ci_sen;
    assert (n == m->n_sen);

    ckd_free (cdsen_start);
    ckd_free (cdsen_end);

    sseq_compress (m);
    
    E_INFO("%d CI-phone, %d CD-phone, %d emitstate/phone, %d CI-sen, %d Sen, %d Sen-Seq\n",
	   m->n_ciphone, m->n_phone - m->n_ciphone, m->n_emit_state,
	   m->n_ci_sen, m->n_sen, m->n_sseq);
    
    fclose (fp);
    
    return m;
}


void mdef_sseq2sen_active (mdef_t *mdef, int32 *sseq, int32 *sen)
{
    int32 ss, i;
    s3senid_t *sp;
    
    for (ss = 0; ss < mdef_n_sseq(mdef); ss++) {
	if (sseq[ss]) {
	    sp = mdef->sseq[ss];
	    for (i = 0; i < mdef_n_emit_state(mdef); i++)
		sen[sp[i]] = 1;
	}
    }
}

/* RAH 4.23.01, Need to step down the ->next list to see if there are
   any more things to free
 */



/* RAH 4.19.01, Attempt to free memory that was allocated within this module
   I have not verified that all the memory has been freed. I've taken only a 
   reasonable effort for now.
   RAH 4.24.01 - verified that all memory is released.
 */
void mdef_free_recursive_lc (ph_lc_t *lc)
{
  if (lc == NULL) return;

  if (lc->rclist) 
    mdef_free_recursive_rc (lc->rclist);

  if (lc->next) 
    mdef_free_recursive_lc (lc->next);

  ckd_free ((void *) lc);
}

void mdef_free_recursive_rc (ph_rc_t *rc)
{
  if (rc == NULL) return;

    if (rc->next) 
      mdef_free_recursive_rc (rc->next);

    ckd_free ((void *) rc);
}


/* RAH, Free memory that was allocated in mdef_init 
   Rational purify shows that no leaks exist
 */
   
void mdef_free (mdef_t *m)
{
  int i,j;

  if (m) { 
    if (m->ciphone2n_cd_sen)
      ckd_free    ((void *)m->ciphone2n_cd_sen);
    if (m->sen2cimap)
      ckd_free    ((void *)m->sen2cimap);
    if (m->cd2cisen)
      ckd_free    ((void *)m->cd2cisen);

    /* RAH, go down the ->next list and delete all the pieces */
    for (i=0;i<N_WORD_POSN;i++)
      for (j=0;j<m->n_ciphone;j++) 
	if (m->wpos_ci_lclist[i][j]) {
	  mdef_free_recursive_lc (m->wpos_ci_lclist[i][j]->next);
	  mdef_free_recursive_rc (m->wpos_ci_lclist[i][j]->rclist);
	}
    
    for (i=0;i<N_WORD_POSN;i++)
      for (j=0;j<m->n_ciphone;j++) 
	if (m->wpos_ci_lclist[i][j])  
	  ckd_free ((void *) m->wpos_ci_lclist[i][j]);

    
    if (m->wpos_ci_lclist)
      ckd_free_2d ((void *)m->wpos_ci_lclist);
    if (m->sseq) 
      ckd_free_2d ((void *)m->sseq);
    /* Free phone context */
    if (m->phone) 
      ckd_free    ((void *)m->phone);    
    if (m->ciphone_ht)
      hash_free (m->ciphone_ht);

    for (i=0;i<m->n_ciphone;i++) {
      if (m->ciphone[i].name) 
	ckd_free    ((void *)m->ciphone[i].name);
    }
    if (m->ciphone) 
      ckd_free    ((void *)m->ciphone);
    
    ckd_free    ((void *)m);
  }
}

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
33522	### 1	93095.8904721724	2.240000e+02	7.200000e+01	67.8
12119	### 2	89153.026256781	7.820000e+02	1.480000e+02	81.0

                /* ====================================================================
 * Copyright (c) 1999-2001 Carnegie Mellon University.  All rights
 * reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer. 
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 * This work was supported in part by funding from the Defense Advanced 
 * Research Projects Agency and the National Science Foundation of the 
 * United States of America, and the CMU Sphinx Speech Consortium.
 *
 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * ====================================================================
 *
 */
/*
 * subvq.c
 * 
 * **********************************************
 * CMU ARPA Speech Project
 *
 * Copyright (c) 1999 Carnegie Mellon University.
 * ALL RIGHTS RESERVED.
 * **********************************************
 * 
 * HISTORY
 * 
 * 20.Apr.2001  RAH (rhoughton@mediasite.com, ricky.houghton@cs.cmu.edu)
 *              Updated subvq_free () to free allocated memory
 * 
 * 17-Dec-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Added handling of a single sub-vector in subvq_mgau_shortlist().
 * 
 * 15-Dec-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Changes due to moving subvq_t.{frm_sen_eval,frm_gau_eval} to cont_mgau.h.
 * 
 * 14-Dec-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Added subvq_t.{frm_sen_eval,frm_gau_eval}.  Changed subvq_frame_eval to
 * 		return the normalization factor.
 * 
 * 06-Dec-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Added subvq_subvec_eval_logs3().
 * 
 * 14-Oct-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Changed ci_active flags input to sen_active in subvq_frame_eval().
 * 
 * 21-Jul-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Bugfix in subvq_init() that used g for allocating working areas, even
 *		though g wasn't defined.
 * 
 * 20-Jul-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Added subvq_gautbl_eval_logs3() and used it in subvq_frame_eval().
 * 
 * 12-Mar-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 * 		Started.
 */


#include "subvq.h"
/* #include "cmd_ln_args.h"	*/ /* RAH, added so we can allow for -vqeval parameter */
#include "s3types.h"

/* RAH, 5.8.01, VQ_EVAL determines how many vectors are used to
 * compute the shortlist, for now this value is only relevant when n_sv =3.
 * Setting it to 1 means that only the CEP values are estimated, 2 means that 
 * CEP and delta values are estimated, 3 means all three are estimated.
 * Note, we must adjust the beam widths as we muck around with these.
 */
static int VQ_EVAL = 3;		


/*
 * Precompute variances/(covariance-matrix-determinants) to simplify Mahalanobis distance
 * calculation.  Also, calculate 1/(det) for the original codebooks, based on the VQ vars.
 * (See comment in libmisc/vector.h.)
 */
static void subvq_maha_precomp (subvq_t *svq, float64 floor)
{
    int32 s;
    float32 *lrd;
    vector_gautbl_t *gautbl;
    
    E_INFO("Precomputing Mahalanobis distance invariants\n");
    
    lrd = (float32 *) ckd_calloc (svq->n_sv * svq->vqsize, sizeof(float32));
    
    for (s = 0; s < svq->n_sv; s++) {
	gautbl = &(svq->gautbl[s]);

	vector_gautbl_var_floor (gautbl, floor);
	
	gautbl->lrd = lrd;
	lrd += svq->vqsize;
	vector_gautbl_maha_precomp (gautbl);
    }
}


static void subvq_map_compact (subvq_t *vq, mgau_model_t *g)
{
    int32 r, c, c2, s;
    
    if (g) {
	if ((g->n_mgau != vq->origsize.r) || (g->max_comp != vq->origsize.c))
	    E_FATAL("Model size conflict: %d x %d (SubVQ) vs %d x %d (Original)\n",
		    vq->origsize.r, vq->origsize.c, g->n_mgau, g->max_comp);
    }
    
    /*
     * Compress map entries by removing invalid ones.  NOTE: There's not much of a consistency
     * check to ensure that the entries remaining do map correctly on to the valid entries in
     * the original (parent) mixture Gaussian model g.  The only check we have is to verify
     * the NUMBER of valid entries (components) in each mixture.
     */
    for (r = 0; r < vq->origsize.r; r++) {
	for (c = 0, c2 = 0; c < vq->origsize.c; c++) {
	    if (vq->map[r][c][0] < 0) {
		/* All ought to be < 0 */
		for (s = 1; s < vq->n_sv; s++) {
		    if (vq->map[r][c][s] >= 0)
			E_FATAL("Partially undefined map[%d][%d]\n", r, c);
		}
	    } else {
		if (c2 != c) {
		    for (s = 0; s < vq->n_sv; s++) {
			if (vq->map[r][c][s] < 0)
			    E_FATAL("Partially undefined map[%d][%d]\n", r, c);
			vq->map[r][c2][s] = vq->map[r][c][s];
		    }
		}
		c2++;
	    }
	}
	
	if (g && (c2 != mgau_n_comp (g, r)))
	    E_FATAL("Mixture %d: #Valid components conflict: %d (SubVQ) vs %d (Original)\n",
		    r, c2, mgau_n_comp(g,r));
	
	/* Invalidate the remaining map entries, for good measure */
	for (; c2 < vq->origsize.c; c2++) {
	    for (s = 0; s < vq->n_sv; s++)
		vq->map[r][c2][s] = -1;
	}
    }
}


/*
 * At this point, a map entries is an index within the subvector; i.e., map[r][c][s] points to
 * a subvq codeword within subvector s.  In evaluating a complete density using it subvq
 * component scores, these maps are used, with 2 lookups for each sub-vector.  To reduce the
 * lookups, linearize the map entries by viewing the subvq scores as a linear array rather than
 * a 2-D array.
 */
static void subvq_map_linearize (subvq_t *vq)
{
    int32 r, c, s;
    
    for (r = 0; r < vq->origsize.r; r++) {
	for (c = 0; (c < vq->origsize.c) && (vq->map[r][c][0] >= 0); c++) {
	    for (s = 0; s < vq->n_sv; s++)
		vq->map[r][c][s] = (s * vq->vqsize) + vq->map[r][c][s];
	}
    }
}


subvq_t *subvq_init (char *file, float64 varfloor, int32 max_sv, mgau_model_t *g)
{
    FILE *fp;
    char line[16384];
    int32 n_sv;		/* #Subvectors in file, as opposed to that used */
    int32 s, k, l, n, r, c;
    char *strp;
    subvq_t *vq;
    
    VQ_EVAL = cmd_ln_int32 ("-vqeval");	/* RAH, Currently only works when n_sv = 3, values computed but ignored in other cases */

    E_INFO("Loading Mixture Gaussian sub-VQ file '%s' (vq_eval: %d)\n", file,VQ_EVAL);
    
    vq = (subvq_t *) ckd_calloc (1, sizeof(subvq_t));
    
    fp = myfopen(file, "r");
    
    /* Read until "Sub-vectors" */
    for (;;) {
	if (fgets (line, sizeof(line), fp) == NULL)
	    E_FATAL("Failed to read VQParam header\n");
	if (sscanf (line, "VQParam %d %d -> %d %d",
		    &(vq->origsize.r), &(vq->origsize.c), &(vq->n_sv), &(vq->vqsize)) == 4)
	    break;
    }
    
    if (g) {
	if ((g->n_mgau != vq->origsize.r) || (g->max_comp != vq->origsize.c))
	    E_FATAL("Model size conflict: %d x %d (SubVQ) vs %d x %d (Original)\n",
		    vq->origsize.r, vq->origsize.c, g->n_mgau, g->max_comp);
    }
    
    if (max_sv < 0)
      max_sv = vq->n_sv;
    if (max_sv < vq->n_sv)
      E_INFO("Using %d subvectors out of %d\n", max_sv, vq->n_sv);
    else if (max_sv > vq->n_sv) {
      E_WARN("#Subvectors specified(%d) > available(%d); using latter\n", max_sv, vq->n_sv);
      max_sv = vq->n_sv;
    }
    
    n_sv = vq->n_sv;
    vq->n_sv = max_sv;
    if (vq->n_sv < VQ_EVAL)	/* RAH, 5.9.01, sanity check to make sure VQ_EVAL isn't higher than the n_sv */
      VQ_EVAL = vq->n_sv;
    vq->featdim = (int32 **) ckd_calloc (vq->n_sv, sizeof(int32 *));
    vq->gautbl = (vector_gautbl_t *) ckd_calloc (vq->n_sv, sizeof(vector_gautbl_t));
    vq->map = (int32 ***) ckd_calloc_3d (vq->origsize.r, vq->origsize.c, vq->n_sv,
					 sizeof(int32));
    
    /* Read subvector sizes and feature dimension maps */
    for (s = 0; s < n_sv; s++) {
	if ((fgets (line, sizeof(line), fp) == NULL) ||
	    (sscanf (line, "Subvector %d length %d%n", &k, &l, &n) != 2) || (k != s))
	    E_FATAL("Error reading length(subvector %d)\n", s);
	
	if (s < vq->n_sv) {
	    vq->gautbl[s].veclen = l;
	    vq->featdim[s] = (int32 *) ckd_calloc (vq->gautbl[s].veclen, sizeof(int32));
	
	    for (strp = line+n, c = 0; c < vq->gautbl[s].veclen; c++) {
		if (sscanf (strp, "%d%n", &(vq->featdim[s][c]), &n) != 1)
		    E_FATAL("Error reading subvector(%d).featdim(%d)\n", s, c);
		strp += n;
	    }
	    
	    vector_gautbl_alloc (&(vq->gautbl[s]), vq->vqsize, vq->gautbl[s].veclen);
	}
    }
    
    /* Echo info for sanity check */
    E_INFO("Original #codebooks(states)/codewords: %d x %d\n", vq->origsize.r, vq->origsize.c);
    E_INFO("Subvectors: %d, VQsize: %d\n", vq->n_sv, vq->vqsize);
    for (s = 0; s < vq->n_sv; s++) {
	E_INFO("SV %d feature dims(%d): ", s, vq->gautbl[s].veclen);
	for (c = 0; c < vq->gautbl[s].veclen; c++)
	    fprintf (stdout, " %2d", vq->featdim[s][c]);
	fprintf (stdout, "\n");
    }
    //SPEC fflush (stderr);
    
    /* Read VQ codebooks and maps for each subvector */
    for (s = 0; s < n_sv; s++) {
	E_INFO("Reading subvq %d%s\n", s, (s < vq->n_sv) ? "" : " (not used)");
	
	E_INFO("Reading codebook\n");
	if ((fgets (line, sizeof(line), fp) == NULL) ||
	    (sscanf (line, "Codebook %d", &k) != 1) || (k != s))
	    E_FATAL("Error reading codebook header\n", s);
	
	for (r = 0; r < vq->vqsize; r++) {
	    if (fgets (line, sizeof(line), fp) == NULL)
		E_FATAL("Error reading row(%d)\n", r);
	    
	    if (s >= vq->n_sv)
		continue;
	    
	    for (strp = line, c = 0; c < vq->gautbl[s].veclen; c++) {
		if (sscanf (strp, "%f %f%n",
			    &(vq->gautbl[s].mean[r][c]), &(vq->gautbl[s].var[r][c]), &k) != 2)
		    E_FATAL("Error reading row(%d) col(%d)\n", r, c);
		strp += k;
	    }
	}
	
	E_INFO("Reading map\n");
	if ((fgets (line, sizeof(line), fp) == NULL) ||
	    (sscanf (line, "Map %d", &k) != 1) || (k != s))
	    E_FATAL("Error reading map header\n", s);
	
	for (r = 0; r < vq->origsize.r; r++) {
	    if (fgets (line, sizeof(line), fp) == NULL)
		E_FATAL("Error reading row(%d)\n", r);
	    
	    if (s >= vq->n_sv)
		continue;
	    
	    for (strp = line, c = 0; c < vq->origsize.c; c++) {
		if (sscanf (strp, "%d%n", &(vq->map[r][c][s]), &k) != 1)
		    E_FATAL("Error reading row(%d) col(%d)\n", r, c);
		strp += k;
	    }
	}
	
	fflush (stdout);
    }
    
    if ((fscanf (fp, "%s", line) != 1) || (strcmp (line, "End") != 0))
	E_FATAL("Error reading 'End' token\n");
    
    fclose (fp);
    
    subvq_maha_precomp (vq, varfloor);
    subvq_map_compact (vq, g);
    subvq_map_linearize (vq);
    
    n = 0;
    for (s = 0; s < n_sv; s++) {
	if (vq->gautbl[s].veclen > n)
	    n = vq->gautbl[s].veclen;
    }
    assert (n > 0);
    vq->subvec = (float32 *) ckd_calloc (n, sizeof(float32));
    vq->vqdist = (int32 **) ckd_calloc_2d (vq->n_sv, vq->vqsize, sizeof(int32));
    vq->gauscore = (int32 *) ckd_calloc (vq->origsize.c, sizeof(int32));
    vq->mgau_sl = (int32 *) ckd_calloc (vq->origsize.c + 1, sizeof(int32));
    
    return vq;
}




/*
 * Based on previously computed subvq scores (Mahalanobis distances), determine the active
 * components in the given mixture (using the vq->map).
 * Return value: #Candidates in the returned shortlist.
 */
int32 subvq_mgau_shortlist (subvq_t *vq,
				   int32 m,	/* In: Mixture index */
				   int32 n,	/* In: #Components in specified mixture */
				   int32 beam)	/* In: Threshold to select active components */
{
    int32 *gauscore;
    int32 *map;
    int32 i, v, bv, th, nc;
    int32 *sl;
    int32 *vqdist;
    int32 sv_id;
    
    vqdist = vq->vqdist[0];	/* Since map is linearized for efficiency, must also
				   look at vqdist[][] as vqdist[] */
    gauscore = vq->gauscore;
    sl = vq->mgau_sl;
    /* Special case when vq->n_sv == 3; for speed */
    map = vq->map[m][0];
    bv = MAX_NEG_INT32;
    
    switch (vq->n_sv) {
    case 3:
	for (i = 0; i < n; i++) {
	  if (VQ_EVAL == 1) {
	    v = (int32) vqdist[*map];/* If we are not weighting the cep values, we need to adjust the subvqbeam */
	    map += 3;
	  } else {
	    /* RAH, we are ignoring the delta-delta, scoring the delta twice, strangely this works better than weighting the scores  */
	    /* I believe it has to do with the beam widths */
	    if (VQ_EVAL == 2) {
	    v = vqdist[*(map++)];
	      v += 2 * vqdist[*map]; /* RAH Count delta twice, we can keep the same subvqbeam as vq_eval = 3 if we double the delta*/
	      map += 2;
	    } else {
	      v = vqdist[*(map++)];/* Standard way */
	      v += vqdist[*(map++)]; /*  */
	      v += vqdist[*(map++)]; /*  */
	    }
	  }

	  gauscore[i] = v;
	    
	  if (bv < v)
	    bv = v;
	}
	break;
    case 2:
	for (i = 0; i < n; i++) {
	    v = vqdist[*(map++)];
	    v += vqdist[*(map++)];
	    gauscore[i] = v;
	   
	    if (bv < v)
		bv = v;
	}
	break;
    case 1:
      for (i = 0; i < n; i++) {
	v = vqdist[*(map++)];
	gauscore[i] = v;
	
	if (bv < v)
	  bv = v;
      }
      break;
    default:
      for (i = 0; i < n; i++) {
	v=0;
	for(sv_id =0 ; sv_id<vq->n_sv; sv_id++){
	  v+=vqdist[*(map++)];
	}
	gauscore[i] = v;
	
	if (bv < v)
	  bv = v;
      }
    }
    
    th = bv + beam;
    nc = 0;
    for (i = 0; i < n; i++) {
      if (gauscore[i] >= th)
	sl[nc++] = i;
    }
    sl[nc] = -1;
    
    return nc;
}


void subvq_subvec_eval_logs3 (subvq_t *vq, float32 *feat, int32 s)
{
    int32 i;
    int32 *featdim;
    
    /* Extract subvector from feat */
    featdim = vq->featdim[s];
    for (i = 0; i < vq->gautbl[s].veclen; i++)
	vq->subvec[i] = feat[featdim[i]];
    
    /* Evaluate distances between extracted subvector and corresponding codebook */
    vector_gautbl_eval_logs3(&(vq->gautbl[s]), 0, vq->vqsize, vq->subvec, vq->vqdist[s]);
}


void subvq_gautbl_eval_logs3 (subvq_t *vq, float32 *feat)
{
    int32 s, i;
    int32 *featdim;
    
    for (s = 0; s < vq->n_sv; s++) {
	/* Extract subvector from feat */
	featdim = vq->featdim[s];
	for (i = 0; i < vq->gautbl[s].veclen; i++)
	    vq->subvec[i] = feat[featdim[i]];
	
	/* Evaluate distances between extracted subvector and corresponding codebook */
    /* RAH, only evaluate the first VQ_EVAL set of features */
    if (s < VQ_EVAL) 
	vector_gautbl_eval_logs3(&(vq->gautbl[s]), 0, vq->vqsize, vq->subvec, vq->vqdist[s]);
    }
}


int32 subvq_frame_eval (subvq_t *vq, mgau_model_t *g, int32 beam, float32 *feat,
			int32 *sen_active, int32 *senscr)
{
  int32 s;
  int32 best, ns, ng;
  
  best = MAX_NEG_INT32;
  ns = 0;
  ng = 0;
  if (! vq) {
    /* No subvq model, use the original (SLOW!!) */
    for (s = 0; s < g->n_mgau; s++) {
      if ((! sen_active) || sen_active[s]) {
	senscr[s] = mgau_eval (g, s, NULL, feat);
	if (best < senscr[s])
	  best = senscr[s];
	ns++;
	ng += mgau_n_comp (g, s);
      } else
	senscr[s] = S3_LOGPROB_ZERO;
    }
  } else {
    /* Evaluate subvq model for given feature vector */
    subvq_gautbl_eval_logs3 (vq, feat);
    
    /* Find mixture component shortlists using subvq scores, and evaluate senones */
    for (s = 0; s < g->n_mgau; s++) {
      if ((! sen_active) || sen_active[s]) {
	ng += subvq_mgau_shortlist (vq, s, mgau_n_comp(g,s), beam);
	
	senscr[s] = mgau_eval (g, s, vq->mgau_sl, feat);
	if (best < senscr[s])
	  best = senscr[s];
	
	ns++;
      } else
	senscr[s] = S3_LOGPROB_ZERO;
    }
  }
    
    /* Normalize senone scores */
    for (s = 0; s < g->n_mgau; s++)
	senscr[s] -= best;
    
    g->frm_sen_eval = ns;
    g->frm_gau_eval = ng;
    
    return best;
}

int32 subvq_mgau_eval (mgau_model_t *g, subvq_t *vq, int32 m, int32 n, int32 *active)
{
    mgau_t *mgau;
    int32 *map;
    int32 i,j, v, sv_id;
    int32 c;
    int32 *vqdist;
    int32 score;
	int32 last_active;

    float64 f;
    f = log_to_logs3_factor();

    vqdist = vq->vqdist[0];	
    score = S3_LOGPROB_ZERO;
    mgau = &(g->mgau[m]);
    map = vq->map[m][0];

    if(!active){
      for (i = 0; i < n; i++) {
	v=0;
	for(sv_id =0 ; sv_id<vq->n_sv; sv_id++){
	  v+=vqdist[*(map++)];
	}
	score = logs3_add (score, v + mgau->mixw[i]);
	/*v=vqdist[*(map++)];
	  v+=vqdist[*(map++)];
	  v+=vqdist[*(map++)];*/
      }
    }else{

      last_active=0;
      for (i = 0; active[i] >=0; i++) {
	c=active[i];
      }
      for (i = 0; active[i] >=0; i++) {

	/*	E_INFO("Value of c %d\n",c);
		E_INFO("Value of last_active %d\n",last_active);*/
	c=active[i];
	map+=(c-last_active)*vq->n_sv;
	v=0;
	for(sv_id =0 ; sv_id<vq->n_sv; sv_id++){
	  v+=vqdist[*(map++)];
	}
	/*v=vqdist[*(map++)];
	v+=vqdist[*(map++)];
	v+=vqdist[*(map++)];*/

	last_active=c+1;

	score = logs3_add (score, v + mgau->mixw[i]);
      }
    }

    if(score == S3_LOGPROB_ZERO){
      E_INFO("Warning!! Score is S3_LOGPROB_ZERO\n");
    }

    return score;

}


/* RAH, free memory allocated by subvq_init() */
void subvq_free (subvq_t *s)
{
  int i;

  if (s) {
    
    for (i=0;i<s->n_sv;i++) {
      //      vector_gautbl_free (&(s->gautbl[i]));
      if (s->featdim[i]) ckd_free ((void *) s->featdim[i]);
    }


    if (s->featdim) 
      ckd_free ((void *) s->featdim);

    /* Free gaussian table */
    if (s->gautbl) 
      ckd_free ((void *)s->gautbl);      


    /* Free map */
    if (s->map)
      ckd_free_3d ((void ***) s->map);

    if (s->subvec) 
      ckd_free ((void *) s->subvec);

    if (s->vqdist) 
      ckd_free_2d ((void **) s->vqdist);

    if (s->gauscore) 
      ckd_free ((void *) s->gauscore);

    if (s->mgau_sl) 
      ckd_free ((void *) s->mgau_sl);

	
    ckd_free ((void *)s);


  }
}

​x
 
/* ====================================================================
 * Copyright (c) 1999-2001 Carnegie Mellon University.  All rights
 * reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer. 
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 * This work was supported in part by funding from the Defense Advanced 
 * Research Projects Agency and the National Science Foundation of the 
 * United States of America, and the CMU Sphinx Speech Consortium.
 *
 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * ====================================================================
 *
 */
/*
 * subvq.c
 * 
 * **********************************************
 * CMU ARPA Speech Project
 *
 * Copyright (c) 1999 Carnegie Mellon University.
 * ALL RIGHTS RESERVED.
 * **********************************************
 * 
 * HISTORY
 * 
 * 20.Apr.2001  RAH (rhoughton@mediasite.com, ricky.houghton@cs.cmu.edu)
 *              Updated subvq_free () to free allocated memory
 * 
 * 17-Dec-1999  M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 *      Added handling of a single sub-vector in subvq_mgau_shortlist().
 * 
 * 15-Dec-1999  M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 *      Changes due to moving subvq_t.{frm_sen_eval,frm_gau_eval} to cont_mgau.h.
 * 
 * 14-Dec-1999  M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 *      Added subvq_t.{frm_sen_eval,frm_gau_eval}.  Changed subvq_frame_eval to
 *      return the normalization factor.
 * 
 * 06-Dec-1999  M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 *      Added subvq_subvec_eval_logs3().
 * 
 * 14-Oct-1999  M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 *      Changed ci_active flags input to sen_active in subvq_frame_eval().
 * 
 * 21-Jul-1999  M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 *      Bugfix in subvq_init() that used g for allocating working areas, even
 *      though g wasn't defined.
 * 
 * 20-Jul-1999  M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 *      Added subvq_gautbl_eval_logs3() and used it in subvq_frame_eval().
 * 
 * 12-Mar-1999  M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 *      Started.
 */
​
​
#include "subvq.h"
/* #include "cmd_ln_args.h" */ /* RAH, added so we can allow for -vqeval parameter */
#include "s3types.h"
​
/* RAH, 5.8.01, VQ_EVAL determines how many vectors are used to
 * compute the shortlist, for now this value is only relevant when n_sv =3.
 * Setting it to 1 means that only the CEP values are estimated, 2 means that 
 * CEP and delta values are estimated, 3 means all three are estimated.
 * Note, we must adjust the beam widths as we muck around with these.
 */
static int VQ_EVAL = 3;     
​
​
/*
 * Precompute variances/(covariance-matrix-determinants) to simplify Mahalanobis distance
 * calculation.  Also, calculate 1/(det) for the original codebooks, based on the VQ vars.
 * (See comment in libmisc/vector.h.)
 */
static void subvq_maha_precomp (subvq_t *svq, float64 floor)
{
    int32 s;
    float32 *lrd;
    vector_gautbl_t *gautbl;
    
    E_INFO("Precomputing Mahalanobis distance invariants\n");
    
    lrd = (float32 *) ckd_calloc (svq->n_sv * svq->vqsize, sizeof(float32));
    
    for (s = 0; s < svq->n_sv; s++) {
    gautbl = &(svq->gautbl[s]);
​
    vector_gautbl_var_floor (gautbl, floor);
    
    gautbl->lrd = lrd;
    lrd += svq->vqsize;
    vector_gautbl_maha_precomp (gautbl);
    }
}
​
​
static void subvq_map_compact (subvq_t *vq, mgau_model_t *g)
{
    int32 r, c, c2, s;
    
    if (g) {
    if ((g->n_mgau != vq->origsize.r) || (g->max_comp != vq->origsize.c))
        E_FATAL("Model size conflict: %d x %d (SubVQ) vs %d x %d (Original)\n",
            vq->origsize.r, vq->origsize.c, g->n_mgau, g->max_comp);
    }
    
    /*
     * Compress map entries by removing invalid ones.  NOTE: There's not much of a consistency
     * check to ensure that the entries remaining do map correctly on to the valid entries in
     * the original (parent) mixture Gaussian model g.  The only check we have is to verify
     * the NUMBER of valid entries (components) in each mixture.
     */
    for (r = 0; r < vq->origsize.r; r++) {
    for (c = 0, c2 = 0; c < vq->origsize.c; c++) {
        if (vq->map[r][c][0] < 0) {
        /* All ought to be < 0 */
        for (s = 1; s < vq->n_sv; s++) {
            if (vq->map[r][c][s] >= 0)
            E_FATAL("Partially undefined map[%d][%d]\n", r, c);
        }
        } else {
        if (c2 != c) {
            for (s = 0; s < vq->n_sv; s++) {
            if (vq->map[r][c][s] < 0)
                E_FATAL("Partially undefined map[%d][%d]\n", r, c);
            vq->map[r][c2][s] = vq->map[r][c][s];
            }
        }
        c2++;
        }
    }
    
    if (g && (c2 != mgau_n_comp (g, r)))
        E_FATAL("Mixture %d: #Valid components conflict: %d (SubVQ) vs %d (Original)\n",
            r, c2, mgau_n_comp(g,r));
    
    /* Invalidate the remaining map entries, for good measure */
    for (; c2 < vq->origsize.c; c2++) {
        for (s = 0; s < vq->n_sv; s++)
        vq->map[r][c2][s] = -1;
    }
    }
}
​
​
/*
 * At this point, a map entries is an index within the subvector; i.e., map[r][c][s] points to
 * a subvq codeword within subvector s.  In evaluating a complete density using it subvq
 * component scores, these maps are used, with 2 lookups for each sub-vector.  To reduce the
 * lookups, linearize the map entries by viewing the subvq scores as a linear array rather than
 * a 2-D array.
 */
static void subvq_map_linearize (subvq_t *vq)
{
    int32 r, c, s;
    
    for (r = 0; r < vq->origsize.r; r++) {
    for (c = 0; (c < vq->origsize.c) && (vq->map[r][c][0] >= 0); c++) {
        for (s = 0; s < vq->n_sv; s++)
        vq->map[r][c][s] = (s * vq->vqsize) + vq->map[r][c][s];
    }
    }
}
​
​
subvq_t *subvq_init (char *file, float64 varfloor, int32 max_sv, mgau_model_t *g)
{
    FILE *fp;
    char line[16384];
    int32 n_sv;     /* #Subvectors in file, as opposed to that used */
    int32 s, k, l, n, r, c;
    char *strp;
    subvq_t *vq;
    
    VQ_EVAL = cmd_ln_int32 ("-vqeval"); /* RAH, Currently only works when n_sv = 3, values computed but ignored in other cases */
​
    E_INFO("Loading Mixture Gaussian sub-VQ file '%s' (vq_eval: %d)\n", file,VQ_EVAL);
    
    vq = (subvq_t *) ckd_calloc (1, sizeof(subvq_t));
    
    fp = myfopen(file, "r");
    
    /* Read until "Sub-vectors" */
    for (;;) {
    if (fgets (line, sizeof(line), fp) == NULL)
        E_FATAL("Failed to read VQParam header\n");
    if (sscanf (line, "VQParam %d %d -> %d %d",
            &(vq->origsize.r), &(vq->origsize.c), &(vq->n_sv), &(vq->vqsize)) == 4)
        break;
    }
    
    if (g) {
    if ((g->n_mgau != vq->origsize.r) || (g->max_comp != vq->origsize.c))
        E_FATAL("Model size conflict: %d x %d (SubVQ) vs %d x %d (Original)\n",
            vq->origsize.r, vq->origsize.c, g->n_mgau, g->max_comp);
    }
    
    if (max_sv < 0)
      max_sv = vq->n_sv;
    if (max_sv < vq->n_sv)
      E_INFO("Using %d subvectors out of %d\n", max_sv, vq->n_sv);
    else if (max_sv > vq->n_sv) {
      E_WARN("#Subvectors specified(%d) > available(%d); using latter\n", max_sv, vq->n_sv);
      max_sv = vq->n_sv;
    }
    
    n_sv = vq->n_sv;
    vq->n_sv = max_sv;
    if (vq->n_sv < VQ_EVAL) /* RAH, 5.9.01, sanity check to make sure VQ_EVAL isn't higher than the n_sv */
      VQ_EVAL = vq->n_sv;
    vq->featdim = (int32 **) ckd_calloc (vq->n_sv, sizeof(int32 *));
    vq->gautbl = (vector_gautbl_t *) ckd_calloc (vq->n_sv, sizeof(vector_gautbl_t));
    vq->map = (int32 ***) ckd_calloc_3d (vq->origsize.r, vq->origsize.c, vq->n_sv,
                     sizeof(int32));
    
    /* Read subvector sizes and feature dimension maps */
    for (s = 0; s < n_sv; s++) {
    if ((fgets (line, sizeof(line), fp) == NULL) ||
        (sscanf (line, "Subvector %d length %d%n", &k, &l, &n) != 2) || (k != s))
        E_FATAL("Error reading length(subvector %d)\n", s);
    
    if (s < vq->n_sv) {
        vq->gautbl[s].veclen = l;
        vq->featdim[s] = (int32 *) ckd_calloc (vq->gautbl[s].veclen, sizeof(int32));
    
        for (strp = line+n, c = 0; c < vq->gautbl[s].veclen; c++) {
        if (sscanf (strp, "%d%n", &(vq->featdim[s][c]), &n) != 1)
            E_FATAL("Error reading subvector(%d).featdim(%d)\n", s, c);
        strp += n;
        }
        
        vector_gautbl_alloc (&(vq->gautbl[s]), vq->vqsize, vq->gautbl[s].veclen);
    }
    }
    
    /* Echo info for sanity check */
    E_INFO("Original #codebooks(states)/codewords: %d x %d\n", vq->origsize.r, vq->origsize.c);
    E_INFO("Subvectors: %d, VQsize: %d\n", vq->n_sv, vq->vqsize);
    for (s = 0; s < vq->n_sv; s++) {
    E_INFO("SV %d feature dims(%d): ", s, vq->gautbl[s].veclen);
    for (c = 0; c < vq->gautbl[s].veclen; c++)
        fprintf (stdout, " %2d", vq->featdim[s][c]);
    fprintf (stdout, "\n");
    }
    //SPEC fflush (stderr);
    
    /* Read VQ codebooks and maps for each subvector */
    for (s = 0; s < n_sv; s++) {
    E_INFO("Reading subvq %d%s\n", s, (s < vq->n_sv) ? "" : " (not used)");
    
    E_INFO("Reading codebook\n");
    if ((fgets (line, sizeof(line), fp) == NULL) ||
        (sscanf (line, "Codebook %d", &k) != 1) || (k != s))
        E_FATAL("Error reading codebook header\n", s);
    
    for (r = 0; r < vq->vqsize; r++) {
        if (fgets (line, sizeof(line), fp) == NULL)
        E_FATAL("Error reading row(%d)\n", r);
        
        if (s >= vq->n_sv)
        continue;
        
        for (strp = line, c = 0; c < vq->gautbl[s].veclen; c++) {
        if (sscanf (strp, "%f %f%n",
                &(vq->gautbl[s].mean[r][c]), &(vq->gautbl[s].var[r][c]), &k) != 2)
            E_FATAL("Error reading row(%d) col(%d)\n", r, c);
        strp += k;
        }
    }
    
    E_INFO("Reading map\n");
    if ((fgets (line, sizeof(line), fp) == NULL) ||
        (sscanf (line, "Map %d", &k) != 1) || (k != s))
        E_FATAL("Error reading map header\n", s);
    
    for (r = 0; r < vq->origsize.r; r++) {
        if (fgets (line, sizeof(line), fp) == NULL)
        E_FATAL("Error reading row(%d)\n", r);
        
        if (s >= vq->n_sv)
        continue;
        
        for (strp = line, c = 0; c < vq->origsize.c; c++) {
        if (sscanf (strp, "%d%n", &(vq->map[r][c][s]), &k) != 1)
            E_FATAL("Error reading row(%d) col(%d)\n", r, c);
        strp += k;
        }
    }
    
    fflush (stdout);
    }
    
    if ((fscanf (fp, "%s", line) != 1) || (strcmp (line, "End") != 0))
    E_FATAL("Error reading 'End' token\n");
    
    fclose (fp);
    
    subvq_maha_precomp (vq, varfloor);
    subvq_map_compact (vq, g);
    subvq_map_linearize (vq);
    
    n = 0;
    for (s = 0; s < n_sv; s++) {
    if (vq->gautbl[s].veclen > n)
        n = vq->gautbl[s].veclen;
    }
    assert (n > 0);
    vq->subvec = (float32 *) ckd_calloc (n, sizeof(float32));
    vq->vqdist = (int32 **) ckd_calloc_2d (vq->n_sv, vq->vqsize, sizeof(int32));
    vq->gauscore = (int32 *) ckd_calloc (vq->origsize.c, sizeof(int32));
    vq->mgau_sl = (int32 *) ckd_calloc (vq->origsize.c + 1, sizeof(int32));
    
    return vq;
}
​
​
​
​
/*
 * Based on previously computed subvq scores (Mahalanobis distances), determine the active
 * components in the given mixture (using the vq->map).
 * Return value: #Candidates in the returned shortlist.
 */
int32 subvq_mgau_shortlist (subvq_t *vq,
                   int32 m, /* In: Mixture index */
                   int32 n, /* In: #Components in specified mixture */
                   int32 beam)  /* In: Threshold to select active components */
{
    int32 *gauscore;
    int32 *map;
    int32 i, v, bv, th, nc;
    int32 *sl;
    int32 *vqdist;
    int32 sv_id;
    
    vqdist = vq->vqdist[0]; /* Since map is linearized for efficiency, must also
                   look at vqdist[][] as vqdist[] */
    gauscore = vq->gauscore;
    sl = vq->mgau_sl;
    /* Special case when vq->n_sv == 3; for speed */
    map = vq->map[m][0];
    bv = MAX_NEG_INT32;
    
    switch (vq->n_sv) {
    case 3:
    for (i = 0; i < n; i++) {
      if (VQ_EVAL == 1) {
        v = (int32) vqdist[*map];/* If we are not weighting the cep values, we need to adjust the subvqbeam */
        map += 3;
      } else {
        /* RAH, we are ignoring the delta-delta, scoring the delta twice, strangely this works better than weighting the scores  */
        /* I believe it has to do with the beam widths */
        if (VQ_EVAL == 2) {
        v = vqdist[*(map++)];
          v += 2 * vqdist[*map]; /* RAH Count delta twice, we can keep the same subvqbeam as vq_eval = 3 if we double the delta*/
          map += 2;
        } else {
          v = vqdist[*(map++)];/* Standard way */
          v += vqdist[*(map++)]; /*  */
          v += vqdist[*(map++)]; /*  */
        }
      }
​
      gauscore[i] = v;
        
      if (bv < v)
        bv = v;
    }
    break;
    case 2:
    for (i = 0; i < n; i++) {
        v = vqdist[*(map++)];
        v += vqdist[*(map++)];
        gauscore[i] = v;
       
        if (bv < v)
        bv = v;
    }
    break;
    case 1:
      for (i = 0; i < n; i++) {
    v = vqdist[*(map++)];
    gauscore[i] = v;
    
    if (bv < v)
      bv = v;
      }
      break;
    default:
      for (i = 0; i < n; i++) {
    v=0;
    for(sv_id =0 ; sv_id<vq->n_sv; sv_id++){
      v+=vqdist[*(map++)];
    }
    gauscore[i] = v;
    
    if (bv < v)
      bv = v;
      }
    }
    
    th = bv + beam;
    nc = 0;
    for (i = 0; i < n; i++) {
      if (gauscore[i] >= th)
    sl[nc++] = i;
    }
    sl[nc] = -1;
    
    return nc;
}
​
​
void subvq_subvec_eval_logs3 (subvq_t *vq, float32 *feat, int32 s)
{
    int32 i;
    int32 *featdim;
    
    /* Extract subvector from feat */
    featdim = vq->featdim[s];
    for (i = 0; i < vq->gautbl[s].veclen; i++)
    vq->subvec[i] = feat[featdim[i]];
    
    /* Evaluate distances between extracted subvector and corresponding codebook */
    vector_gautbl_eval_logs3(&(vq->gautbl[s]), 0, vq->vqsize, vq->subvec, vq->vqdist[s]);
}
​
​
void subvq_gautbl_eval_logs3 (subvq_t *vq, float32 *feat)
{
    int32 s, i;
    int32 *featdim;
    
    for (s = 0; s < vq->n_sv; s++) {
    /* Extract subvector from feat */
    featdim = vq->featdim[s];
    for (i = 0; i < vq->gautbl[s].veclen; i++)
        vq->subvec[i] = feat[featdim[i]];
    
    /* Evaluate distances between extracted subvector and corresponding codebook */
    /* RAH, only evaluate the first VQ_EVAL set of features */
    if (s < VQ_EVAL) 
    vector_gautbl_eval_logs3(&(vq->gautbl[s]), 0, vq->vqsize, vq->subvec, vq->vqdist[s]);
    }
}
​
​
int32 subvq_frame_eval (subvq_t *vq, mgau_model_t *g, int32 beam, float32 *feat,
            int32 *sen_active, int32 *senscr)
{
  int32 s;
  int32 best, ns, ng;
  
  best = MAX_NEG_INT32;
  ns = 0;
  ng = 0;
  if (! vq) {
    /* No subvq model, use the original (SLOW!!) */
    for (s = 0; s < g->n_mgau; s++) {
      if ((! sen_active) || sen_active[s]) {
    senscr[s] = mgau_eval (g, s, NULL, feat);
    if (best < senscr[s])
      best = senscr[s];
    ns++;
    ng += mgau_n_comp (g, s);
      } else
    senscr[s] = S3_LOGPROB_ZERO;
    }
  } else {
    /* Evaluate subvq model for given feature vector */
    subvq_gautbl_eval_logs3 (vq, feat);
    
    /* Find mixture component shortlists using subvq scores, and evaluate senones */
    for (s = 0; s < g->n_mgau; s++) {
      if ((! sen_active) || sen_active[s]) {
    ng += subvq_mgau_shortlist (vq, s, mgau_n_comp(g,s), beam);
    
    senscr[s] = mgau_eval (g, s, vq->mgau_sl, feat);
    if (best < senscr[s])
      best = senscr[s];
    
    ns++;
      } else
    senscr[s] = S3_LOGPROB_ZERO;
    }
  }
    
    /* Normalize senone scores */
    for (s = 0; s < g->n_mgau; s++)
    senscr[s] -= best;
    
    g->frm_sen_eval = ns;
    g->frm_gau_eval = ng;
    
    return best;
}
​
int32 subvq_mgau_eval (mgau_model_t *g, subvq_t *vq, int32 m, int32 n, int32 *active)
{
    mgau_t *mgau;
    int32 *map;
    int32 i,j, v, sv_id;
    int32 c;
    int32 *vqdist;
    int32 score;
    int32 last_active;
​
    float64 f;
    f = log_to_logs3_factor();
​
    vqdist = vq->vqdist[0]; 
    score = S3_LOGPROB_ZERO;
    mgau = &(g->mgau[m]);
    map = vq->map[m][0];
​
    if(!active){
      for (i = 0; i < n; i++) {
    v=0;
    for(sv_id =0 ; sv_id<vq->n_sv; sv_id++){
      v+=vqdist[*(map++)];
    }
    score = logs3_add (score, v + mgau->mixw[i]);
    /*v=vqdist[*(map++)];
      v+=vqdist[*(map++)];
      v+=vqdist[*(map++)];*/
      }
    }else{
​
      last_active=0;
      for (i = 0; active[i] >=0; i++) {
    c=active[i];
      }
      for (i = 0; active[i] >=0; i++) {
​
    /*  E_INFO("Value of c %d\n",c);
        E_INFO("Value of last_active %d\n",last_active);*/
    c=active[i];
    map+=(c-last_active)*vq->n_sv;
    v=0;
    for(sv_id =0 ; sv_id<vq->n_sv; sv_id++){
      v+=vqdist[*(map++)];
    }
    /*v=vqdist[*(map++)];
    v+=vqdist[*(map++)];
    v+=vqdist[*(map++)];*/
​
    last_active=c+1;
​
    score = logs3_add (score, v + mgau->mixw[i]);
      }
    }
​
    if(score == S3_LOGPROB_ZERO){
      E_INFO("Warning!! Score is S3_LOGPROB_ZERO\n");
    }
​
    return score;
​
}
​
​
/* RAH, free memory allocated by subvq_init() */
void subvq_free (subvq_t *s)
{
  int i;
​
  if (s) {
    
    for (i=0;i<s->n_sv;i++) {
      //      vector_gautbl_free (&(s->gautbl[i]));
      if (s->featdim[i]) ckd_free ((void *) s->featdim[i]);
    }
​
​
    if (s->featdim) 
      ckd_free ((void *) s->featdim);
​
    /* Free gaussian table */
    if (s->gautbl) 
      ckd_free ((void *)s->gautbl);      
​
​
    /* Free map */
    if (s->map)
      ckd_free_3d ((void ***) s->map);
​
    if (s->subvec) 
      ckd_free ((void *) s->subvec);
​
    if (s->vqdist) 
      ckd_free_2d ((void **) s->vqdist);
​
    if (s->gauscore) 
      ckd_free ((void *) s->gauscore);
​
    if (s->mgau_sl) 
      ckd_free ((void *) s->mgau_sl);
​
    
    ckd_free ((void *)s);
​
​
  }
}

One invocation only

invivosubvq_subvq_init_270

Real time: 2.926047e+08 -- Predicted time: 0.000000e+00

Callcount: 1

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	0.000000e+00	2.927055e+08	100.0

                /* ====================================================================
 * Copyright (c) 1999-2001 Carnegie Mellon University.  All rights
 * reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer. 
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 * This work was supported in part by funding from the Defense Advanced 
 * Research Projects Agency and the National Science Foundation of the 
 * United States of America, and the CMU Sphinx Speech Consortium.
 *
 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * ====================================================================
 *
 */
/*
 * mdef.c -- HMM model definition: base (CI) phones and triphones
 *
 * **********************************************
 * CMU ARPA Speech Project
 *
 * Copyright (c) 1999 Carnegie Mellon University.
 * ALL RIGHTS RESERVED.
 * **********************************************
 * 
 * HISTORY
 * 
 * 19.Apr-2001  Ricky Houghton, added code for free allocated memory
 *
 * 14-Oct-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon
 * 		Added mdef_sseq2sen_active().
 * 
 * 06-May-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon
 * 		In mdef_phone_id(), added backing off to silence phone context from filler
 * 		context if original triphone not found.
 * 
 * 30-Apr-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon
 * 		Added senone-sequence id (ssid) to phone_t and appropriate functions to
 * 		maintain it.  Instead, moved state sequence info to mdef_t.
 * 
 * 13-Jul-96	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Added mdef_phone_str().
 * 
 * 01-Jan-96	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Allowed mdef_phone_id_nearest to return base phone id if either
 * 		left or right context (or both) is undefined.
 * 
 * 01-Jan-96	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Created.
 */


/*
 * Major assumptions:
 *   All phones have same #states, same topology.
 *   Every phone has exactly one non-emitting, final state--the last one.
 *   CI phones must appear first in model definition file.
 */


#include "mdef.h"


#define MODEL_DEF_VERSION	"0.3"


void mdef_dump (FILE *fp, mdef_t *m)
{
    int32 i, j;
    int32 ssid;
    char buf[1024];
    
    fprintf (fp, "%d ciphone\n", m->n_ciphone);
    fprintf (fp, "%d phone\n", m->n_phone);
    fprintf (fp, "%d emitstate\n", m->n_emit_state);
    fprintf (fp, "%d cisen\n", m->n_ci_sen);
    fprintf (fp, "%d sen\n", m->n_sen);
    fprintf (fp, "%d tmat\n", m->n_tmat);
    
    for (i = 0; i < m->n_phone; i++) {
	mdef_phone_str (m, i, buf);
	ssid = m->phone[i].ssid;
	
	fprintf (fp, "%3d %5d", m->phone[i].tmat, ssid);
	for (j = 0; j < m->n_emit_state; j++)
	    fprintf (fp, " %5d", m->sseq[ssid][j]);
	fprintf (fp, "\t");
	for (j = 0; j < m->n_emit_state; j++)
	    fprintf (fp, " %3d", m->cd2cisen[m->sseq[ssid][j]]);
	fprintf (fp, "\t%s\n", buf);
    }
    
    fflush (fp);
}


#if 0
int32 mdef_hmm_cmp (mdef_t *m, s3pid_t p1, s3pid_t p2)
{
    int32 i;
    
    if (m->phone[p1].tmat != m->phone[p2].tmat)
	return -1;
    
    for (i = 0; i < m->n_emit_state; i++)
	if (m->phone[p1].state[i] != m->phone[p2].state[i])
	    return -1;
    
    return 0;
}
#endif


static void ciphone_add (mdef_t *m, char *ci, s3pid_t p)
{
    assert (p < m->n_ciphone);

    m->ciphone[p].name = (char *) ckd_salloc (ci); /* freed in mdef_free */
    if (hash_enter (m->ciphone_ht, m->ciphone[p].name, p) != p)
	E_FATAL("hash_enter(%s) failed; duplicate CIphone?\n", m->ciphone[p].name);
}


static ph_lc_t *find_ph_lc (ph_lc_t *lclist, s3cipid_t lc)
{
    ph_lc_t *lcptr;

    for (lcptr = lclist; lcptr && (lcptr->lc != lc); lcptr = lcptr->next);
    return lcptr;
}


static ph_rc_t *find_ph_rc (ph_rc_t *rclist, s3cipid_t rc)
{
    ph_rc_t *rcptr;

    for (rcptr = rclist; rcptr && (rcptr->rc != rc); rcptr = rcptr->next);
    return rcptr;
}


static void triphone_add (mdef_t *m,
			  s3cipid_t ci, s3cipid_t lc, s3cipid_t rc, word_posn_t wpos,
			  s3pid_t p)
{
    ph_lc_t *lcptr;
    ph_rc_t *rcptr;
    
    assert (p < m->n_phone);

    /* Fill in phone[p] information (state and tmat mappings added later) */
    m->phone[p].ci = ci;
    m->phone[p].lc = lc;
    m->phone[p].rc = rc;
    m->phone[p].wpos = wpos;

    /* Create <ci,lc,rc,wpos> -> p mapping if not a CI phone */
    if (p >= m->n_ciphone) {
	if ((lcptr = find_ph_lc (m->wpos_ci_lclist[wpos][(int)ci], lc)) 
	    == NULL)
	{
	lcptr = (ph_lc_t *) ckd_calloc (1, sizeof(ph_lc_t)); /* freed at mdef_free, I believe */
	    lcptr->lc = lc;
	    lcptr->next = m->wpos_ci_lclist[wpos][(int)ci];
	    m->wpos_ci_lclist[wpos][(int)ci] = lcptr; /* This is what needs to be freed */
	}
	if ((rcptr = find_ph_rc (lcptr->rclist, rc)) != NULL) {
	    char buf[4096];
	    
	    mdef_phone_str (m, rcptr->pid, buf);
	    E_FATAL("Duplicate triphone: %s\n", buf);
	}
	
      rcptr = (ph_rc_t *) ckd_calloc (1, sizeof(ph_rc_t)); /* freed in mdef_free, I believe */
	rcptr->rc = rc;
	rcptr->pid = p;
	rcptr->next = lcptr->rclist;
	lcptr->rclist = rcptr;
    }
}


s3cipid_t mdef_ciphone_id (mdef_t *m, char *ci)
{
    int32 id;
    
    assert (m);
    assert (ci);
    
    if (hash_lookup (m->ciphone_ht, ci, &id) < 0)
	return (BAD_S3CIPID);
    return ((s3cipid_t) id);
}


const char *mdef_ciphone_str (mdef_t *m, s3cipid_t id)
{
    assert (m);
    assert ((id >= 0) && (id < m->n_ciphone));
    
    return (m->ciphone[(int)id].name);
}


int32 mdef_phone_str (mdef_t *m, s3pid_t pid, char *buf)
{
    char *wpos_name;
    
    assert (m);
    assert ((pid >= 0) && (pid < m->n_phone));
    wpos_name = WPOS_NAME;
    
    buf[0] = '\0';
    if (pid < m->n_ciphone)
	sprintf (buf, "%s", mdef_ciphone_str (m, (s3cipid_t) pid));
    else {
	sprintf (buf, "(%s,%s,%s,%c)",
		 mdef_ciphone_str(m, m->phone[pid].ci),
		 mdef_ciphone_str(m, m->phone[pid].lc),
		 mdef_ciphone_str(m, m->phone[pid].rc),
		 wpos_name[m->phone[pid].wpos]);
    }
    return 0;
}


s3pid_t mdef_phone_id (mdef_t *m, 
		       s3cipid_t ci, s3cipid_t lc, s3cipid_t rc, word_posn_t wpos)
{
    ph_lc_t *lcptr;
    ph_rc_t *rcptr;
    s3cipid_t newl, newr;
    
    assert (m);
    assert ((ci >= 0) && (ci < m->n_ciphone));
    assert ((lc >= 0) && (lc < m->n_ciphone));
    assert ((rc >= 0) && (rc < m->n_ciphone));
    assert ((wpos >= 0) && (wpos < N_WORD_POSN));

    if (((lcptr = find_ph_lc (m->wpos_ci_lclist[wpos][(int)ci], lc)) == NULL) ||
	((rcptr = find_ph_rc (lcptr->rclist, rc)) == NULL)) {
	/* Not found; backoff to silence context if non-silence filler context */
	if (NOT_S3CIPID(m->sil))
	    return BAD_S3PID;
	
	newl = m->ciphone[(int)lc].filler ? m->sil : lc;
	newr = m->ciphone[(int)rc].filler ? m->sil : rc;
	if ((newl == lc) && (newr == rc))
	    return BAD_S3PID;
	
	return (mdef_phone_id (m, ci, newl, newr, wpos));
    }
    
    return (rcptr->pid);
}


s3pid_t mdef_phone_id_nearest (mdef_t *m, 
			       s3cipid_t b, s3cipid_t l, s3cipid_t r, word_posn_t pos)
{
    word_posn_t tmppos;
    s3pid_t p;
    s3cipid_t newl, newr;
    char *wpos_name;
    
    assert (m);
    assert ((b >= 0) && (b < m->n_ciphone));
    assert ((pos >= 0) && (pos < N_WORD_POSN));

    if ((NOT_S3CIPID(l)) || (NOT_S3CIPID(r)))
	return ((s3pid_t) b);
    
    assert ((l >= 0) && (l < m->n_ciphone));
    assert ((r >= 0) && (r < m->n_ciphone));
    
    p = mdef_phone_id (m, b, l, r, pos);
    if (IS_S3PID(p))
	return p;
    
    /* Exact triphone not found; backoff to other word positions */
    for (tmppos = 0; tmppos < N_WORD_POSN; tmppos++) {
	if (tmppos != pos) {
	    p = mdef_phone_id (m, b, l, r, tmppos);
	    if (IS_S3PID(p))
		return p;
	}
    }
    
    /* Nothing yet; backoff to silence phone if non-silence filler context */
    if (IS_S3CIPID(m->sil)) {
	newl = m->ciphone[(int)l].filler ? m->sil : l;
	newr = m->ciphone[(int)r].filler ? m->sil : r;
	if ((newl != l) || (newr != r)) {
	    p = mdef_phone_id (m, b, newl, newr, pos);
	    if (IS_S3PID(p))
		return p;
	    
	    for (tmppos = 0; tmppos < N_WORD_POSN; tmppos++) {
		if (tmppos != pos) {
		    p = mdef_phone_id (m, b, newl, newr, tmppos);
		    if (IS_S3PID(p))
			return p;
		}
	    }
	}
    }
    
    /* Nothing yet; backoff to base phone */
    if ((m->n_phone > m->n_ciphone) && (! m->ciphone[(int)b].filler)) {
	wpos_name = WPOS_NAME;
#if 0
	E_WARN("Triphone(%s,%s,%s,%c) not found; backing off to CIphone\n",
	       mdef_ciphone_str(m, b),
	       mdef_ciphone_str(m, l),
	       mdef_ciphone_str(m, r),
	       wpos_name[pos]);
#endif
    }
    return ((s3pid_t) b);
}


int32 mdef_phone_components (mdef_t *m,
			     s3pid_t p,
			     s3cipid_t *b,
			     s3cipid_t *l,
			     s3cipid_t *r,
			     word_posn_t *pos)
{
    assert (m);
    assert ((p >= 0) && (p < m->n_phone));

    *b = m->phone[p].ci;
    *l = m->phone[p].lc;
    *r = m->phone[p].rc;
    *pos = m->phone[p].wpos;

    return 0;
}


int32 mdef_is_ciphone (mdef_t *m, s3pid_t p)
{
    assert (m);
    assert ((p >= 0) && (p < m->n_phone));
    
    return ((p < m->n_ciphone) ? 1 : 0);
}


/* Parse tmat and state->senone mappings for phone p and fill in structure */
static void parse_tmat_senmap (mdef_t *m, char *line, int32 off, s3pid_t p)
{
    int32 wlen, n, s;
    char word[1024], *lp;

    lp = line + off;
    
    /* Read transition matrix id */
    if ((sscanf (lp, "%d%n", &n, &wlen) != 1) || (n < 0))
        E_FATAL("Missing or bad transition matrix id: %s\n", line);
    m->phone[p].tmat = n;
    if (m->n_tmat <= n)
	E_FATAL("tmat-id(%d) > #tmat in header(%d): %s\n", n, m->n_tmat, line);
    lp += wlen;
    
    /* Read senone mappings for each emitting state */
    for (n = 0; n < m->n_emit_state; n++) {
        if ((sscanf (lp, "%d%n", &s, &wlen) != 1) || (s < 0))
	    E_FATAL("Missing or bad state[%d]->senone mapping: %s\n", n, line);
	if ((p < m->n_ciphone) && (m->n_ci_sen <= s))
	    E_FATAL("CI-senone-id(%d) > #CI-senones(%d): %s\n", s, m->n_ci_sen, line);
	if (m->n_sen <= s)
	    E_FATAL("Senone-id(%d) > #senones(%d): %s\n", s, m->n_sen, line);

	m->sseq[p][n] = s;
	lp += wlen;
    }

    /* Check for the last non-emitting state N */
    if ((sscanf (lp, "%s%n", word, &wlen) != 1) || (strcmp (word, "N") != 0))
        E_FATAL("Missing non-emitting state spec: %s\n", line);
    lp += wlen;

    /* Check for end of line */
    if (sscanf (lp, "%s%n", word, &wlen) == 1)
        E_FATAL("Non-empty beyond non-emitting final state: %s\n", line);
}


static void parse_base_line (mdef_t *m, char *line, s3pid_t p)
{
    int32 wlen, n;
    char word[1024], *lp;
    s3cipid_t ci;

    lp = line;
    
    /* Read base phone name */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing base phone name: %s\n", line);
    lp += wlen;
    
    /* Make sure it's not a duplicate */
    ci = mdef_ciphone_id (m, word);
    if (IS_S3CIPID(ci))
        E_FATAL("Duplicate base phone: %s\n", line);

    /* Add ciphone to ciphone table with id p */
    ciphone_add (m, word, p);
    ci = (s3cipid_t) p;

    /* Read and skip "-" for lc, rc, wpos */
    for (n = 0; n < 3; n++) {
	if ((sscanf (lp, "%s%n", word, &wlen) != 1) || (strcmp (word, "-") != 0))
	    E_FATAL("Bad context info for base phone: %s\n", line);
	lp += wlen;
    }
    
    /* Read filler attribute, if present */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing filler atribute field: %s\n", line);
    lp += wlen;
    if (strcmp (word, "filler") == 0)
        m->ciphone[(int)ci].filler = 1;
    else if (strcmp (word, "n/a") == 0)
        m->ciphone[(int)ci].filler = 0;
    else
        E_FATAL("Bad filler attribute field: %s\n", line);

    triphone_add (m, ci, BAD_S3CIPID, BAD_S3CIPID, WORD_POSN_UNDEFINED, p);

    /* Parse remainder of line: transition matrix and state->senone mappings */
    parse_tmat_senmap (m, line, lp-line, p);
}


static void parse_tri_line (mdef_t *m, char *line, s3pid_t p)
{
    int32 wlen;
    char word[1024], *lp;
    s3cipid_t ci, lc, rc;
    word_posn_t wpos = WORD_POSN_BEGIN;

    lp = line;
    
    /* Read base phone name */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing base phone name: %s\n", line);
    lp += wlen;

    ci = mdef_ciphone_id (m, word);
    if (NOT_S3CIPID(ci))
        E_FATAL("Unknown base phone: %s\n", line);

    /* Read lc */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing left context: %s\n", line);
    lp += wlen;
    lc = mdef_ciphone_id (m, word);
    if (NOT_S3CIPID(lc))
        E_FATAL("Unknown left context: %s\n", line);

    /* Read rc */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing right context: %s\n", line);
    lp += wlen;
    rc = mdef_ciphone_id (m, word);
    if (NOT_S3CIPID(rc))
        E_FATAL("Unknown right  context: %s\n", line);
    
    /* Read tripone word-position within word */
    if ((sscanf (lp, "%s%n", word, &wlen) != 1) || (word[1] != '\0'))
        E_FATAL("Missing or bad word-position spec: %s\n", line);
    lp += wlen;
    switch (word[0]) {
    case 'b': wpos = WORD_POSN_BEGIN; break;
    case 'e': wpos = WORD_POSN_END; break;
    case 's': wpos = WORD_POSN_SINGLE; break;
    case 'i': wpos = WORD_POSN_INTERNAL; break;
    default: E_FATAL("Bad word-position spec: %s\n", line);
    }

    /* Read filler attribute, if present.  Must match base phone attribute */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing filler attribute field: %s\n", line);
    lp += wlen;
    if (((strcmp (word, "filler") == 0) && (m->ciphone[(int)ci].filler)) ||
	((strcmp (word, "n/a") == 0) && (! m->ciphone[(int)ci].filler))) {
	/* Everything is fine */
    } else
        E_FATAL("Bad filler attribute field: %s\n", line);
    
    triphone_add (m, ci, lc, rc, wpos, p);

    /* Parse remainder of line: transition matrix and state->senone mappings */
    parse_tmat_senmap (m, line, lp-line, p);
}


static void sseq_compress (mdef_t *m)
{
    hash_table_t *h;
    s3senid_t **sseq;
    int32 n_sseq;
    int32 p, j, k;
    glist_t g;
    gnode_t *gn;
    hash_entry_t *he;
    
    k = m->n_emit_state * sizeof(s3senid_t);
    
    h = hash_new (m->n_phone, HASH_CASE_YES);
    n_sseq = 0;
    
    /* Identify unique senone-sequence IDs.  BUG: tmat-id not being considered!! */
    for (p = 0; p < m->n_phone; p++) {
	/* Add senone sequence to hash table */
	if ((j = hash_enter_bkey (h, (char *)(m->sseq[p]), k, n_sseq)) == n_sseq)
	    n_sseq++;
	
	m->phone[p].ssid = j;
    }
    
    /* Generate compacted sseq table */
    sseq = (s3senid_t **) ckd_calloc_2d (n_sseq, m->n_emit_state, sizeof(s3senid_t));/* freed in mdef_free() */
    
    g = hash_tolist (h, &j);
    assert (j == n_sseq);
    
    for (gn = g; gn; gn = gnode_next(gn)) {
	he = (hash_entry_t *) gnode_ptr (gn);
	j = hash_entry_val(he);
	memcpy (sseq[j], hash_entry_key(he), k);
    }
    glist_free (g);
    
    /* Free the old, temporary senone sequence table, replace with compacted one */
    ckd_free_2d ((void **) m->sseq);
    m->sseq = sseq;
    m->n_sseq = n_sseq;
    
    hash_free (h);
}


static int32 noncomment_line(char *line, int32 size, FILE *fp)
{
    while (fgets (line, size, fp) != NULL) {
        if (line[0] != '#')
	    return 0;
    }
    return -1;
}


/*
 * Initialize phones (ci and triphones) and state->senone mappings from .mdef file.
 */
mdef_t *mdef_init (char *mdeffile)
{
    FILE *fp;
    int32 n_ci, n_tri, n_map, n;
    char tag[1024], buf[1024];
    s3senid_t **senmap;
    s3pid_t p;
    int32 s, ci, cd;
    mdef_t *m;
    int32 *cdsen_start, *cdsen_end;

    if (! mdeffile)
	E_FATAL("No mdef-file\n");

    E_INFO("Reading model definition: %s\n", mdeffile);

    m = (mdef_t *) ckd_calloc (1, sizeof(mdef_t)); /* freed in mdef_free */
    
    if ((fp = fopen(mdeffile, "r")) == NULL)
        E_FATAL_SYSTEM("fopen(%s,r) failed\n", mdeffile);

    if (noncomment_line(buf, sizeof(buf), fp) < 0)
        E_FATAL("Empty file: %s\n", mdeffile);

    if (strncmp(buf, MODEL_DEF_VERSION, strlen(MODEL_DEF_VERSION)) != 0)
        E_FATAL("Version error: Expecing %s, but read %s\n", MODEL_DEF_VERSION, buf);

    /* Read #base phones, #triphones, #senone mappings defined in header */
    n_ci = -1;
    n_tri = -1;
    n_map = -1;
    m->n_ci_sen = -1;
    m->n_sen = -1;
    m->n_tmat = -1;
    do {
	if (noncomment_line(buf, sizeof(buf), fp) < 0)
	    E_FATAL("Incomplete header\n");

	if ((sscanf(buf, "%d %s", &n, tag) != 2) || (n < 0))
	    E_FATAL("Error in header: %s\n", buf);

	if (strcmp(tag, "n_base") == 0)
	    n_ci = n;
	else if (strcmp(tag, "n_tri") == 0)
	    n_tri = n;
	else if (strcmp(tag, "n_state_map") == 0)
	    n_map = n;
	else if (strcmp(tag, "n_tied_ci_state") == 0)
	    m->n_ci_sen = n;
	else if (strcmp(tag, "n_tied_state") == 0)
	    m->n_sen = n;
	else if (strcmp(tag, "n_tied_tmat") == 0)
	    m->n_tmat = n;
	else
	    E_FATAL("Unknown header line: %s\n", buf);
    } while ((n_ci < 0) || (n_tri < 0) || (n_map < 0) ||
	     (m->n_ci_sen < 0) || (m->n_sen < 0) || (m->n_tmat < 0));

    if ((n_ci == 0) || (m->n_ci_sen == 0) || (m->n_tmat == 0) || (m->n_ci_sen > m->n_sen))
        E_FATAL("%s: Error in header\n", mdeffile);
    
    /* Check typesize limits */
    if (n_ci >= MAX_S3CIPID)
	E_FATAL("%s: #CI phones (%d) exceeds limit (%d)\n", mdeffile, n_ci, MAX_S3CIPID);
    if (n_ci + n_tri >= MAX_S3PID)
	E_FATAL("%s: #Phones (%d) exceeds limit (%d)\n", mdeffile, n_ci+n_tri, MAX_S3PID);
    if (m->n_sen >= MAX_S3SENID)
	E_FATAL("%s: #senones (%d) exceeds limit (%d)\n", mdeffile, m->n_sen, MAX_S3SENID);
    if (m->n_tmat >= MAX_S3TMATID)
	E_FATAL("%s: #tmats (%d) exceeds limit (%d)\n", mdeffile, m->n_tmat, MAX_S3TMATID);
    
    m->n_emit_state = (n_map / (n_ci+n_tri)) - 1;
    if ((m->n_emit_state+1) * (n_ci+n_tri) != n_map)
        E_FATAL("Header error: n_state_map not a multiple of n_ci*n_tri\n");

    /* Initialize ciphone info */
    m->n_ciphone = n_ci;
    m->ciphone_ht = hash_new (n_ci, 1);	/* With case-insensitive string names */ /* freed in mdef_free */
    m->ciphone = (ciphone_t *) ckd_calloc (n_ci, sizeof(ciphone_t)); /* freed in mdef_free */
    /* RAH, let's null the pointers so that we can reliably deallocate them */
    /*    for (i=0;i<m->n_ciphone;i++) { */ /*  */
    /*m->ciphone[i].name = NULL; */ /*  */
    /*} */ /*  */

    /* Initialize phones info (ciphones + triphones) */
    m->n_phone = n_ci + n_tri;
    m->phone = (phone_t *) ckd_calloc (m->n_phone, sizeof(phone_t)); /* freed in mdef_free */

    /* Allocate space for state->senone map for each phone */
    senmap = (s3senid_t **) ckd_calloc_2d (m->n_phone, m->n_emit_state, sizeof(s3senid_t));/* freed in mdef_free */
    m->sseq = senmap;	/* TEMPORARY; until it is compressed into just the unique ones */
    
    /* Allocate initial space for <ci,lc,rc,wpos> -> pid mapping */
    m->wpos_ci_lclist = (ph_lc_t ***) ckd_calloc_2d (N_WORD_POSN, m->n_ciphone, sizeof(ph_lc_t *)); /* freed in mdef_free */

    /*
     * Read base phones and triphones.  They'll simply be assigned a running sequence
     * number as their "phone-id".  If the phone-id < n_ci, it's a ciphone.
     */

    /* Read base phones */
    for (p = 0; p < n_ci; p++) {
        if (noncomment_line(buf, sizeof(buf), fp) < 0)
	    E_FATAL("Premature EOF reading CIphone %d\n", p);
        parse_base_line (m, buf, p);
    }
    m->sil = mdef_ciphone_id (m, S3_SILENCE_CIPHONE);
    
    /* Read triphones, if any */
    for (; p < m->n_phone; p++) {
        if (noncomment_line(buf, sizeof(buf), fp) < 0)
	    E_FATAL("Premature EOF reading phone %d\n", p);
        parse_tri_line (m, buf, p);
    }

    if (noncomment_line(buf, sizeof(buf), fp) >= 0)
	E_ERROR("Non-empty file beyond expected #phones (%d)\n", m->n_phone);

    /* Build CD senones to CI senones map */
    if (m->n_ciphone * m->n_emit_state != m->n_ci_sen)
	E_FATAL("#CI-senones(%d) != #CI-phone(%d) x #emitting-states(%d)\n",
		m->n_ci_sen, m->n_ciphone, m->n_emit_state);
    m->cd2cisen = (s3senid_t *) ckd_calloc (m->n_sen, sizeof(s3senid_t)); /* freed in mdef_free */

    m->sen2cimap = (s3cipid_t *) ckd_calloc (m->n_sen, sizeof(s3cipid_t)); /* freed in mdef_free */

    for (s = 0; s < m->n_sen; s++)
	m->sen2cimap[s] = BAD_S3CIPID;
    for (s = 0; s < m->n_ci_sen; s++) {		/* CI senones */
	m->cd2cisen[s] = (s3senid_t) s;
	m->sen2cimap[s] = s / m->n_emit_state;
    }
    for (p = n_ci; p < m->n_phone; p++) {	/* CD senones */
	for (s = 0; s < m->n_emit_state; s++) {
	    cd = m->sseq[p][s];
	    ci = m->sseq[(int)m->phone[p].ci][s];
	    m->cd2cisen[cd] = (s3senid_t) ci;
	    m->sen2cimap[cd] = m->phone[p].ci;
	}
    }
    
    /*
     * Count #senones (CI+CD) for each CI phone.
     * HACK!!  For handling holes in senone-CIphone mappings.  Does not work if holes
     * are present at the beginning or end of senones for a given CIphone.
     */
    cdsen_start = (int32 *) ckd_calloc (m->n_ciphone, sizeof(int32)); /* freed locally */

    cdsen_end = (int32 *) ckd_calloc (m->n_ciphone, sizeof(int32)); /* freed locally */

    for (s = m->n_ci_sen; s < m->n_sen; s++) {
	if (NOT_S3CIPID(m->sen2cimap[s]))
	    continue;
	
	if (! cdsen_start[(int)m->sen2cimap[s]])
	    cdsen_start[(int)m->sen2cimap[s]] = s;
	cdsen_end[(int)m->sen2cimap[s]] = s;
    }

    /* Fill up holes */
    for (s = m->n_ci_sen; s < m->n_sen; s++) {
	if (IS_S3CIPID(m->sen2cimap[s]))
	    continue;

	/* Check if properly inside the observed ranges above */
	for (p = 0; p < m->n_ciphone; p++) {
	    if ((s > cdsen_start[p]) && (s < cdsen_end[p]))
		break;
	}
	if (p >= m->n_ciphone)
	    E_FATAL("Unreferenced senone %d; cannot determine parent CIphone\n", s);
	m->sen2cimap[s] = p;
    }

    /* Build #CD-senones for each CIphone */
    m->ciphone2n_cd_sen = (int32 *) ckd_calloc (m->n_ciphone, sizeof(int32));/* freed mdef_free */
    n = 0;
    for (p = 0; p < m->n_ciphone; p++) {
	if (cdsen_start[p] > 0) {
	    m->ciphone2n_cd_sen[p] = cdsen_end[p] - cdsen_start[p] + 1;
	    n += m->ciphone2n_cd_sen[p];
	}
    }
    n += m->n_ci_sen;
    assert (n == m->n_sen);

    ckd_free (cdsen_start);
    ckd_free (cdsen_end);

    sseq_compress (m);
    
    E_INFO("%d CI-phone, %d CD-phone, %d emitstate/phone, %d CI-sen, %d Sen, %d Sen-Seq\n",
	   m->n_ciphone, m->n_phone - m->n_ciphone, m->n_emit_state,
	   m->n_ci_sen, m->n_sen, m->n_sseq);
    
    fclose (fp);
    
    return m;
}


void mdef_sseq2sen_active (mdef_t *mdef, int32 *sseq, int32 *sen)
{
    int32 ss, i;
    s3senid_t *sp;
    
    for (ss = 0; ss < mdef_n_sseq(mdef); ss++) {
	if (sseq[ss]) {
	    sp = mdef->sseq[ss];
	    for (i = 0; i < mdef_n_emit_state(mdef); i++)
		sen[sp[i]] = 1;
	}
    }
}

/* RAH 4.23.01, Need to step down the ->next list to see if there are
   any more things to free
 */



/* RAH 4.19.01, Attempt to free memory that was allocated within this module
   I have not verified that all the memory has been freed. I've taken only a 
   reasonable effort for now.
   RAH 4.24.01 - verified that all memory is released.
 */
void mdef_free_recursive_lc (ph_lc_t *lc)
{
  if (lc == NULL) return;

  if (lc->rclist) 
    mdef_free_recursive_rc (lc->rclist);

  if (lc->next) 
    mdef_free_recursive_lc (lc->next);

  ckd_free ((void *) lc);
}

void mdef_free_recursive_rc (ph_rc_t *rc)
{
  if (rc == NULL) return;

    if (rc->next) 
      mdef_free_recursive_rc (rc->next);

    ckd_free ((void *) rc);
}


/* RAH, Free memory that was allocated in mdef_init 
   Rational purify shows that no leaks exist
 */
   
void mdef_free (mdef_t *m)
{
  int i,j;

  if (m) { 
    if (m->ciphone2n_cd_sen)
      ckd_free    ((void *)m->ciphone2n_cd_sen);
    if (m->sen2cimap)
      ckd_free    ((void *)m->sen2cimap);
    if (m->cd2cisen)
      ckd_free    ((void *)m->cd2cisen);

    /* RAH, go down the ->next list and delete all the pieces */
    for (i=0;i<N_WORD_POSN;i++)
      for (j=0;j<m->n_ciphone;j++) 
	if (m->wpos_ci_lclist[i][j]) {
	  mdef_free_recursive_lc (m->wpos_ci_lclist[i][j]->next);
	  mdef_free_recursive_rc (m->wpos_ci_lclist[i][j]->rclist);
	}
    
    for (i=0;i<N_WORD_POSN;i++)
      for (j=0;j<m->n_ciphone;j++) 
	if (m->wpos_ci_lclist[i][j])  
	  ckd_free ((void *) m->wpos_ci_lclist[i][j]);

    
    if (m->wpos_ci_lclist)
      ckd_free_2d ((void *)m->wpos_ci_lclist);
    if (m->sseq) 
      ckd_free_2d ((void *)m->sseq);
    /* Free phone context */
    if (m->phone) 
      ckd_free    ((void *)m->phone);    
    if (m->ciphone_ht)
      hash_free (m->ciphone_ht);

    for (i=0;i<m->n_ciphone;i++) {
      if (m->ciphone[i].name) 
	ckd_free    ((void *)m->ciphone[i].name);
    }
    if (m->ciphone) 
      ckd_free    ((void *)m->ciphone);
    
    ckd_free    ((void *)m);
  }
}

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	3.784801e+07	3.627406e+07	4.1

                /* ====================================================================
 * Copyright (c) 1999-2001 Carnegie Mellon University.  All rights
 * reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer. 
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 * This work was supported in part by funding from the Defense Advanced 
 * Research Projects Agency and the National Science Foundation of the 
 * United States of America, and the CMU Sphinx Speech Consortium.
 *
 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * ====================================================================
 *
 */
/*
 * mdef.c -- HMM model definition: base (CI) phones and triphones
 *
 * **********************************************
 * CMU ARPA Speech Project
 *
 * Copyright (c) 1999 Carnegie Mellon University.
 * ALL RIGHTS RESERVED.
 * **********************************************
 * 
 * HISTORY
 * 
 * 19.Apr-2001  Ricky Houghton, added code for free allocated memory
 *
 * 14-Oct-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon
 * 		Added mdef_sseq2sen_active().
 * 
 * 06-May-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon
 * 		In mdef_phone_id(), added backing off to silence phone context from filler
 * 		context if original triphone not found.
 * 
 * 30-Apr-1999	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon
 * 		Added senone-sequence id (ssid) to phone_t and appropriate functions to
 * 		maintain it.  Instead, moved state sequence info to mdef_t.
 * 
 * 13-Jul-96	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Added mdef_phone_str().
 * 
 * 01-Jan-96	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Allowed mdef_phone_id_nearest to return base phone id if either
 * 		left or right context (or both) is undefined.
 * 
 * 01-Jan-96	M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University.
 * 		Created.
 */


/*
 * Major assumptions:
 *   All phones have same #states, same topology.
 *   Every phone has exactly one non-emitting, final state--the last one.
 *   CI phones must appear first in model definition file.
 */


#include "mdef.h"


#define MODEL_DEF_VERSION	"0.3"


void mdef_dump (FILE *fp, mdef_t *m)
{
    int32 i, j;
    int32 ssid;
    char buf[1024];
    
    fprintf (fp, "%d ciphone\n", m->n_ciphone);
    fprintf (fp, "%d phone\n", m->n_phone);
    fprintf (fp, "%d emitstate\n", m->n_emit_state);
    fprintf (fp, "%d cisen\n", m->n_ci_sen);
    fprintf (fp, "%d sen\n", m->n_sen);
    fprintf (fp, "%d tmat\n", m->n_tmat);
    
    for (i = 0; i < m->n_phone; i++) {
	mdef_phone_str (m, i, buf);
	ssid = m->phone[i].ssid;
	
	fprintf (fp, "%3d %5d", m->phone[i].tmat, ssid);
	for (j = 0; j < m->n_emit_state; j++)
	    fprintf (fp, " %5d", m->sseq[ssid][j]);
	fprintf (fp, "\t");
	for (j = 0; j < m->n_emit_state; j++)
	    fprintf (fp, " %3d", m->cd2cisen[m->sseq[ssid][j]]);
	fprintf (fp, "\t%s\n", buf);
    }
    
    fflush (fp);
}


#if 0
int32 mdef_hmm_cmp (mdef_t *m, s3pid_t p1, s3pid_t p2)
{
    int32 i;
    
    if (m->phone[p1].tmat != m->phone[p2].tmat)
	return -1;
    
    for (i = 0; i < m->n_emit_state; i++)
	if (m->phone[p1].state[i] != m->phone[p2].state[i])
	    return -1;
    
    return 0;
}
#endif


static void ciphone_add (mdef_t *m, char *ci, s3pid_t p)
{
    assert (p < m->n_ciphone);

    m->ciphone[p].name = (char *) ckd_salloc (ci); /* freed in mdef_free */
    if (hash_enter (m->ciphone_ht, m->ciphone[p].name, p) != p)
	E_FATAL("hash_enter(%s) failed; duplicate CIphone?\n", m->ciphone[p].name);
}


static ph_lc_t *find_ph_lc (ph_lc_t *lclist, s3cipid_t lc)
{
    ph_lc_t *lcptr;

    for (lcptr = lclist; lcptr && (lcptr->lc != lc); lcptr = lcptr->next);
    return lcptr;
}


static ph_rc_t *find_ph_rc (ph_rc_t *rclist, s3cipid_t rc)
{
    ph_rc_t *rcptr;

    for (rcptr = rclist; rcptr && (rcptr->rc != rc); rcptr = rcptr->next);
    return rcptr;
}


static void triphone_add (mdef_t *m,
			  s3cipid_t ci, s3cipid_t lc, s3cipid_t rc, word_posn_t wpos,
			  s3pid_t p)
{
    ph_lc_t *lcptr;
    ph_rc_t *rcptr;
    
    assert (p < m->n_phone);

    /* Fill in phone[p] information (state and tmat mappings added later) */
    m->phone[p].ci = ci;
    m->phone[p].lc = lc;
    m->phone[p].rc = rc;
    m->phone[p].wpos = wpos;

    /* Create <ci,lc,rc,wpos> -> p mapping if not a CI phone */
    if (p >= m->n_ciphone) {
	if ((lcptr = find_ph_lc (m->wpos_ci_lclist[wpos][(int)ci], lc)) 
	    == NULL)
	{
	lcptr = (ph_lc_t *) ckd_calloc (1, sizeof(ph_lc_t)); /* freed at mdef_free, I believe */
	    lcptr->lc = lc;
	    lcptr->next = m->wpos_ci_lclist[wpos][(int)ci];
	    m->wpos_ci_lclist[wpos][(int)ci] = lcptr; /* This is what needs to be freed */
	}
	if ((rcptr = find_ph_rc (lcptr->rclist, rc)) != NULL) {
	    char buf[4096];
	    
	    mdef_phone_str (m, rcptr->pid, buf);
	    E_FATAL("Duplicate triphone: %s\n", buf);
	}
	
      rcptr = (ph_rc_t *) ckd_calloc (1, sizeof(ph_rc_t)); /* freed in mdef_free, I believe */
	rcptr->rc = rc;
	rcptr->pid = p;
	rcptr->next = lcptr->rclist;
	lcptr->rclist = rcptr;
    }
}


s3cipid_t mdef_ciphone_id (mdef_t *m, char *ci)
{
    int32 id;
    
    assert (m);
    assert (ci);
    
    if (hash_lookup (m->ciphone_ht, ci, &id) < 0)
	return (BAD_S3CIPID);
    return ((s3cipid_t) id);
}


const char *mdef_ciphone_str (mdef_t *m, s3cipid_t id)
{
    assert (m);
    assert ((id >= 0) && (id < m->n_ciphone));
    
    return (m->ciphone[(int)id].name);
}


int32 mdef_phone_str (mdef_t *m, s3pid_t pid, char *buf)
{
    char *wpos_name;
    
    assert (m);
    assert ((pid >= 0) && (pid < m->n_phone));
    wpos_name = WPOS_NAME;
    
    buf[0] = '\0';
    if (pid < m->n_ciphone)
	sprintf (buf, "%s", mdef_ciphone_str (m, (s3cipid_t) pid));
    else {
	sprintf (buf, "(%s,%s,%s,%c)",
		 mdef_ciphone_str(m, m->phone[pid].ci),
		 mdef_ciphone_str(m, m->phone[pid].lc),
		 mdef_ciphone_str(m, m->phone[pid].rc),
		 wpos_name[m->phone[pid].wpos]);
    }
    return 0;
}


s3pid_t mdef_phone_id (mdef_t *m, 
		       s3cipid_t ci, s3cipid_t lc, s3cipid_t rc, word_posn_t wpos)
{
    ph_lc_t *lcptr;
    ph_rc_t *rcptr;
    s3cipid_t newl, newr;
    
    assert (m);
    assert ((ci >= 0) && (ci < m->n_ciphone));
    assert ((lc >= 0) && (lc < m->n_ciphone));
    assert ((rc >= 0) && (rc < m->n_ciphone));
    assert ((wpos >= 0) && (wpos < N_WORD_POSN));

    if (((lcptr = find_ph_lc (m->wpos_ci_lclist[wpos][(int)ci], lc)) == NULL) ||
	((rcptr = find_ph_rc (lcptr->rclist, rc)) == NULL)) {
	/* Not found; backoff to silence context if non-silence filler context */
	if (NOT_S3CIPID(m->sil))
	    return BAD_S3PID;
	
	newl = m->ciphone[(int)lc].filler ? m->sil : lc;
	newr = m->ciphone[(int)rc].filler ? m->sil : rc;
	if ((newl == lc) && (newr == rc))
	    return BAD_S3PID;
	
	return (mdef_phone_id (m, ci, newl, newr, wpos));
    }
    
    return (rcptr->pid);
}


s3pid_t mdef_phone_id_nearest (mdef_t *m, 
			       s3cipid_t b, s3cipid_t l, s3cipid_t r, word_posn_t pos)
{
    word_posn_t tmppos;
    s3pid_t p;
    s3cipid_t newl, newr;
    char *wpos_name;
    
    assert (m);
    assert ((b >= 0) && (b < m->n_ciphone));
    assert ((pos >= 0) && (pos < N_WORD_POSN));

    if ((NOT_S3CIPID(l)) || (NOT_S3CIPID(r)))
	return ((s3pid_t) b);
    
    assert ((l >= 0) && (l < m->n_ciphone));
    assert ((r >= 0) && (r < m->n_ciphone));
    
    p = mdef_phone_id (m, b, l, r, pos);
    if (IS_S3PID(p))
	return p;
    
    /* Exact triphone not found; backoff to other word positions */
    for (tmppos = 0; tmppos < N_WORD_POSN; tmppos++) {
	if (tmppos != pos) {
	    p = mdef_phone_id (m, b, l, r, tmppos);
	    if (IS_S3PID(p))
		return p;
	}
    }
    
    /* Nothing yet; backoff to silence phone if non-silence filler context */
    if (IS_S3CIPID(m->sil)) {
	newl = m->ciphone[(int)l].filler ? m->sil : l;
	newr = m->ciphone[(int)r].filler ? m->sil : r;
	if ((newl != l) || (newr != r)) {
	    p = mdef_phone_id (m, b, newl, newr, pos);
	    if (IS_S3PID(p))
		return p;
	    
	    for (tmppos = 0; tmppos < N_WORD_POSN; tmppos++) {
		if (tmppos != pos) {
		    p = mdef_phone_id (m, b, newl, newr, tmppos);
		    if (IS_S3PID(p))
			return p;
		}
	    }
	}
    }
    
    /* Nothing yet; backoff to base phone */
    if ((m->n_phone > m->n_ciphone) && (! m->ciphone[(int)b].filler)) {
	wpos_name = WPOS_NAME;
#if 0
	E_WARN("Triphone(%s,%s,%s,%c) not found; backing off to CIphone\n",
	       mdef_ciphone_str(m, b),
	       mdef_ciphone_str(m, l),
	       mdef_ciphone_str(m, r),
	       wpos_name[pos]);
#endif
    }
    return ((s3pid_t) b);
}


int32 mdef_phone_components (mdef_t *m,
			     s3pid_t p,
			     s3cipid_t *b,
			     s3cipid_t *l,
			     s3cipid_t *r,
			     word_posn_t *pos)
{
    assert (m);
    assert ((p >= 0) && (p < m->n_phone));

    *b = m->phone[p].ci;
    *l = m->phone[p].lc;
    *r = m->phone[p].rc;
    *pos = m->phone[p].wpos;

    return 0;
}


int32 mdef_is_ciphone (mdef_t *m, s3pid_t p)
{
    assert (m);
    assert ((p >= 0) && (p < m->n_phone));
    
    return ((p < m->n_ciphone) ? 1 : 0);
}


/* Parse tmat and state->senone mappings for phone p and fill in structure */
static void parse_tmat_senmap (mdef_t *m, char *line, int32 off, s3pid_t p)
{
    int32 wlen, n, s;
    char word[1024], *lp;

    lp = line + off;
    
    /* Read transition matrix id */
    if ((sscanf (lp, "%d%n", &n, &wlen) != 1) || (n < 0))
        E_FATAL("Missing or bad transition matrix id: %s\n", line);
    m->phone[p].tmat = n;
    if (m->n_tmat <= n)
	E_FATAL("tmat-id(%d) > #tmat in header(%d): %s\n", n, m->n_tmat, line);
    lp += wlen;
    
    /* Read senone mappings for each emitting state */
    for (n = 0; n < m->n_emit_state; n++) {
        if ((sscanf (lp, "%d%n", &s, &wlen) != 1) || (s < 0))
	    E_FATAL("Missing or bad state[%d]->senone mapping: %s\n", n, line);
	if ((p < m->n_ciphone) && (m->n_ci_sen <= s))
	    E_FATAL("CI-senone-id(%d) > #CI-senones(%d): %s\n", s, m->n_ci_sen, line);
	if (m->n_sen <= s)
	    E_FATAL("Senone-id(%d) > #senones(%d): %s\n", s, m->n_sen, line);

	m->sseq[p][n] = s;
	lp += wlen;
    }

    /* Check for the last non-emitting state N */
    if ((sscanf (lp, "%s%n", word, &wlen) != 1) || (strcmp (word, "N") != 0))
        E_FATAL("Missing non-emitting state spec: %s\n", line);
    lp += wlen;

    /* Check for end of line */
    if (sscanf (lp, "%s%n", word, &wlen) == 1)
        E_FATAL("Non-empty beyond non-emitting final state: %s\n", line);
}


static void parse_base_line (mdef_t *m, char *line, s3pid_t p)
{
    int32 wlen, n;
    char word[1024], *lp;
    s3cipid_t ci;

    lp = line;
    
    /* Read base phone name */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing base phone name: %s\n", line);
    lp += wlen;
    
    /* Make sure it's not a duplicate */
    ci = mdef_ciphone_id (m, word);
    if (IS_S3CIPID(ci))
        E_FATAL("Duplicate base phone: %s\n", line);

    /* Add ciphone to ciphone table with id p */
    ciphone_add (m, word, p);
    ci = (s3cipid_t) p;

    /* Read and skip "-" for lc, rc, wpos */
    for (n = 0; n < 3; n++) {
	if ((sscanf (lp, "%s%n", word, &wlen) != 1) || (strcmp (word, "-") != 0))
	    E_FATAL("Bad context info for base phone: %s\n", line);
	lp += wlen;
    }
    
    /* Read filler attribute, if present */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing filler atribute field: %s\n", line);
    lp += wlen;
    if (strcmp (word, "filler") == 0)
        m->ciphone[(int)ci].filler = 1;
    else if (strcmp (word, "n/a") == 0)
        m->ciphone[(int)ci].filler = 0;
    else
        E_FATAL("Bad filler attribute field: %s\n", line);

    triphone_add (m, ci, BAD_S3CIPID, BAD_S3CIPID, WORD_POSN_UNDEFINED, p);

    /* Parse remainder of line: transition matrix and state->senone mappings */
    parse_tmat_senmap (m, line, lp-line, p);
}


static void parse_tri_line (mdef_t *m, char *line, s3pid_t p)
{
    int32 wlen;
    char word[1024], *lp;
    s3cipid_t ci, lc, rc;
    word_posn_t wpos = WORD_POSN_BEGIN;

    lp = line;
    
    /* Read base phone name */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing base phone name: %s\n", line);
    lp += wlen;

    ci = mdef_ciphone_id (m, word);
    if (NOT_S3CIPID(ci))
        E_FATAL("Unknown base phone: %s\n", line);

    /* Read lc */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing left context: %s\n", line);
    lp += wlen;
    lc = mdef_ciphone_id (m, word);
    if (NOT_S3CIPID(lc))
        E_FATAL("Unknown left context: %s\n", line);

    /* Read rc */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing right context: %s\n", line);
    lp += wlen;
    rc = mdef_ciphone_id (m, word);
    if (NOT_S3CIPID(rc))
        E_FATAL("Unknown right  context: %s\n", line);
    
    /* Read tripone word-position within word */
    if ((sscanf (lp, "%s%n", word, &wlen) != 1) || (word[1] != '\0'))
        E_FATAL("Missing or bad word-position spec: %s\n", line);
    lp += wlen;
    switch (word[0]) {
    case 'b': wpos = WORD_POSN_BEGIN; break;
    case 'e': wpos = WORD_POSN_END; break;
    case 's': wpos = WORD_POSN_SINGLE; break;
    case 'i': wpos = WORD_POSN_INTERNAL; break;
    default: E_FATAL("Bad word-position spec: %s\n", line);
    }

    /* Read filler attribute, if present.  Must match base phone attribute */
    if (sscanf (lp, "%s%n", word, &wlen) != 1)
	E_FATAL("Missing filler attribute field: %s\n", line);
    lp += wlen;
    if (((strcmp (word, "filler") == 0) && (m->ciphone[(int)ci].filler)) ||
	((strcmp (word, "n/a") == 0) && (! m->ciphone[(int)ci].filler))) {
	/* Everything is fine */
    } else
        E_FATAL("Bad filler attribute field: %s\n", line);
    
    triphone_add (m, ci, lc, rc, wpos, p);

    /* Parse remainder of line: transition matrix and state->senone mappings */
    parse_tmat_senmap (m, line, lp-line, p);
}


static void sseq_compress (mdef_t *m)
{
    hash_table_t *h;
    s3senid_t **sseq;
    int32 n_sseq;
    int32 p, j, k;
    glist_t g;
    gnode_t *gn;
    hash_entry_t *he;
    
    k = m->n_emit_state * sizeof(s3senid_t);
    
    h = hash_new (m->n_phone, HASH_CASE_YES);
    n_sseq = 0;
    
    /* Identify unique senone-sequence IDs.  BUG: tmat-id not being considered!! */
    for (p = 0; p < m->n_phone; p++) {
	/* Add senone sequence to hash table */
	if ((j = hash_enter_bkey (h, (char *)(m->sseq[p]), k, n_sseq)) == n_sseq)
	    n_sseq++;
	
	m->phone[p].ssid = j;
    }
    
    /* Generate compacted sseq table */
    sseq = (s3senid_t **) ckd_calloc_2d (n_sseq, m->n_emit_state, sizeof(s3senid_t));/* freed in mdef_free() */
    
    g = hash_tolist (h, &j);
    assert (j == n_sseq);
    
    for (gn = g; gn; gn = gnode_next(gn)) {
	he = (hash_entry_t *) gnode_ptr (gn);
	j = hash_entry_val(he);
	memcpy (sseq[j], hash_entry_key(he), k);
    }
    glist_free (g);
    
    /* Free the old, temporary senone sequence table, replace with compacted one */
    ckd_free_2d ((void **) m->sseq);
    m->sseq = sseq;
    m->n_sseq = n_sseq;
    
    hash_free (h);
}


static int32 noncomment_line(char *line, int32 size, FILE *fp)
{
    while (fgets (line, size, fp) != NULL) {
        if (line[0] != '#')
	    return 0;
    }
    return -1;
}


/*
 * Initialize phones (ci and triphones) and state->senone mappings from .mdef file.
 */
mdef_t *mdef_init (char *mdeffile)
{
    FILE *fp;
    int32 n_ci, n_tri, n_map, n;
    char tag[1024], buf[1024];
    s3senid_t **senmap;
    s3pid_t p;
    int32 s, ci, cd;
    mdef_t *m;
    int32 *cdsen_start, *cdsen_end;

    if (! mdeffile)
	E_FATAL("No mdef-file\n");

    E_INFO("Reading model definition: %s\n", mdeffile);

    m = (mdef_t *) ckd_calloc (1, sizeof(mdef_t)); /* freed in mdef_free */
    
    if ((fp = fopen(mdeffile, "r")) == NULL)
        E_FATAL_SYSTEM("fopen(%s,r) failed\n", mdeffile);

    if (noncomment_line(buf, sizeof(buf), fp) < 0)
        E_FATAL("Empty file: %s\n", mdeffile);

    if (strncmp(buf, MODEL_DEF_VERSION, strlen(MODEL_DEF_VERSION)) != 0)
        E_FATAL("Version error: Expecing %s, but read %s\n", MODEL_DEF_VERSION, buf);

    /* Read #base phones, #triphones, #senone mappings defined in header */
    n_ci = -1;
    n_tri = -1;
    n_map = -1;
    m->n_ci_sen = -1;
    m->n_sen = -1;
    m->n_tmat = -1;
    do {
	if (noncomment_line(buf, sizeof(buf), fp) < 0)
	    E_FATAL("Incomplete header\n");

	if ((sscanf(buf, "%d %s", &n, tag) != 2) || (n < 0))
	    E_FATAL("Error in header: %s\n", buf);

	if (strcmp(tag, "n_base") == 0)
	    n_ci = n;
	else if (strcmp(tag, "n_tri") == 0)
	    n_tri = n;
	else if (strcmp(tag, "n_state_map") == 0)
	    n_map = n;
	else if (strcmp(tag, "n_tied_ci_state") == 0)
	    m->n_ci_sen = n;
	else if (strcmp(tag, "n_tied_state") == 0)
	    m->n_sen = n;
	else if (strcmp(tag, "n_tied_tmat") == 0)
	    m->n_tmat = n;
	else
	    E_FATAL("Unknown header line: %s\n", buf);
    } while ((n_ci < 0) || (n_tri < 0) || (n_map < 0) ||
	     (m->n_ci_sen < 0) || (m->n_sen < 0) || (m->n_tmat < 0));

    if ((n_ci == 0) || (m->n_ci_sen == 0) || (m->n_tmat == 0) || (m->n_ci_sen > m->n_sen))
        E_FATAL("%s: Error in header\n", mdeffile);
    
    /* Check typesize limits */
    if (n_ci >= MAX_S3CIPID)
	E_FATAL("%s: #CI phones (%d) exceeds limit (%d)\n", mdeffile, n_ci, MAX_S3CIPID);
    if (n_ci + n_tri >= MAX_S3PID)
	E_FATAL("%s: #Phones (%d) exceeds limit (%d)\n", mdeffile, n_ci+n_tri, MAX_S3PID);
    if (m->n_sen >= MAX_S3SENID)
	E_FATAL("%s: #senones (%d) exceeds limit (%d)\n", mdeffile, m->n_sen, MAX_S3SENID);
    if (m->n_tmat >= MAX_S3TMATID)
	E_FATAL("%s: #tmats (%d) exceeds limit (%d)\n", mdeffile, m->n_tmat, MAX_S3TMATID);
    
    m->n_emit_state = (n_map / (n_ci+n_tri)) - 1;
    if ((m->n_emit_state+1) * (n_ci+n_tri) != n_map)
        E_FATAL("Header error: n_state_map not a multiple of n_ci*n_tri\n");

    /* Initialize ciphone info */
    m->n_ciphone = n_ci;
    m->ciphone_ht = hash_new (n_ci, 1);	/* With case-insensitive string names */ /* freed in mdef_free */
    m->ciphone = (ciphone_t *) ckd_calloc (n_ci, sizeof(ciphone_t)); /* freed in mdef_free */
    /* RAH, let's null the pointers so that we can reliably deallocate them */
    /*    for (i=0;i<m->n_ciphone;i++) { */ /*  */
    /*m->ciphone[i].name = NULL; */ /*  */
    /*} */ /*  */

    /* Initialize phones info (ciphones + triphones) */
    m->n_phone = n_ci + n_tri;
    m->phone = (phone_t *) ckd_calloc (m->n_phone, sizeof(phone_t)); /* freed in mdef_free */

    /* Allocate space for state->senone map for each phone */
    senmap = (s3senid_t **) ckd_calloc_2d (m->n_phone, m->n_emit_state, sizeof(s3senid_t));/* freed in mdef_free */
    m->sseq = senmap;	/* TEMPORARY; until it is compressed into just the unique ones */
    
    /* Allocate initial space for <ci,lc,rc,wpos> -> pid mapping */
    m->wpos_ci_lclist = (ph_lc_t ***) ckd_calloc_2d (N_WORD_POSN, m->n_ciphone, sizeof(ph_lc_t *)); /* freed in mdef_free */

    /*
     * Read base phones and triphones.  They'll simply be assigned a running sequence
     * number as their "phone-id".  If the phone-id < n_ci, it's a ciphone.
     */

    /* Read base phones */
    for (p = 0; p < n_ci; p++) {
        if (noncomment_line(buf, sizeof(buf), fp) < 0)
	    E_FATAL("Premature EOF reading CIphone %d\n", p);
        parse_base_line (m, buf, p);
    }
    m->sil = mdef_ciphone_id (m, S3_SILENCE_CIPHONE);
    
    /* Read triphones, if any */
    for (; p < m->n_phone; p++) {
        if (noncomment_line(buf, sizeof(buf), fp) < 0)
	    E_FATAL("Premature EOF reading phone %d\n", p);
        parse_tri_line (m, buf, p);
    }

    if (noncomment_line(buf, sizeof(buf), fp) >= 0)
	E_ERROR("Non-empty file beyond expected #phones (%d)\n", m->n_phone);

    /* Build CD senones to CI senones map */
    if (m->n_ciphone * m->n_emit_state != m->n_ci_sen)
	E_FATAL("#CI-senones(%d) != #CI-phone(%d) x #emitting-states(%d)\n",
		m->n_ci_sen, m->n_ciphone, m->n_emit_state);
    m->cd2cisen = (s3senid_t *) ckd_calloc (m->n_sen, sizeof(s3senid_t)); /* freed in mdef_free */

    m->sen2cimap = (s3cipid_t *) ckd_calloc (m->n_sen, sizeof(s3cipid_t)); /* freed in mdef_free */

    for (s = 0; s < m->n_sen; s++)
	m->sen2cimap[s] = BAD_S3CIPID;
    for (s = 0; s < m->n_ci_sen; s++) {		/* CI senones */
	m->cd2cisen[s] = (s3senid_t) s;
	m->sen2cimap[s] = s / m->n_emit_state;
    }
    for (p = n_ci; p < m->n_phone; p++) {	/* CD senones */
	for (s = 0; s < m->n_emit_state; s++) {
	    cd = m->sseq[p][s];
	    ci = m->sseq[(int)m->phone[p].ci][s];
	    m->cd2cisen[cd] = (s3senid_t) ci;
	    m->sen2cimap[cd] = m->phone[p].ci;
	}
    }
    
    /*
     * Count #senones (CI+CD) for each CI phone.
     * HACK!!  For handling holes in senone-CIphone mappings.  Does not work if holes
     * are present at the beginning or end of senones for a given CIphone.
     */
    cdsen_start = (int32 *) ckd_calloc (m->n_ciphone, sizeof(int32)); /* freed locally */

    cdsen_end = (int32 *) ckd_calloc (m->n_ciphone, sizeof(int32)); /* freed locally */

    for (s = m->n_ci_sen; s < m->n_sen; s++) {
	if (NOT_S3CIPID(m->sen2cimap[s]))
	    continue;
	
	if (! cdsen_start[(int)m->sen2cimap[s]])
	    cdsen_start[(int)m->sen2cimap[s]] = s;
	cdsen_end[(int)m->sen2cimap[s]] = s;
    }

    /* Fill up holes */
    for (s = m->n_ci_sen; s < m->n_sen; s++) {
	if (IS_S3CIPID(m->sen2cimap[s]))
	    continue;

	/* Check if properly inside the observed ranges above */
	for (p = 0; p < m->n_ciphone; p++) {
	    if ((s > cdsen_start[p]) && (s < cdsen_end[p]))
		break;
	}
	if (p >= m->n_ciphone)
	    E_FATAL("Unreferenced senone %d; cannot determine parent CIphone\n", s);
	m->sen2cimap[s] = p;
    }

    /* Build #CD-senones for each CIphone */
    m->ciphone2n_cd_sen = (int32 *) ckd_calloc (m->n_ciphone, sizeof(int32));/* freed mdef_free */
    n = 0;
    for (p = 0; p < m->n_ciphone; p++) {
	if (cdsen_start[p] > 0) {
	    m->ciphone2n_cd_sen[p] = cdsen_end[p] - cdsen_start[p] + 1;
	    n += m->ciphone2n_cd_sen[p];
	}
    }
    n += m->n_ci_sen;
    assert (n == m->n_sen);

    ckd_free (cdsen_start);
    ckd_free (cdsen_end);

    sseq_compress (m);
    
    E_INFO("%d CI-phone, %d CD-phone, %d emitstate/phone, %d CI-sen, %d Sen, %d Sen-Seq\n",
	   m->n_ciphone, m->n_phone - m->n_ciphone, m->n_emit_state,
	   m->n_ci_sen, m->n_sen, m->n_sseq);
    
    fclose (fp);
    
    return m;
}


void mdef_sseq2sen_active (mdef_t *mdef, int32 *sseq, int32 *sen)
{
    int32 ss, i;
    s3senid_t *sp;
    
    for (ss = 0; ss < mdef_n_sseq(mdef); ss++) {
	if (sseq[ss]) {
	    sp = mdef->sseq[ss];
	    for (i = 0; i < mdef_n_emit_state(mdef); i++)
		sen[sp[i]] = 1;
	}
    }
}

/* RAH 4.23.01, Need to step down the ->next list to see if there are
   any more things to free
 */



/* RAH 4.19.01, Attempt to free memory that was allocated within this module
   I have not verified that all the memory has been freed. I've taken only a 
   reasonable effort for now.
   RAH 4.24.01 - verified that all memory is released.
 */
void mdef_free_recursive_lc (ph_lc_t *lc)
{
  if (lc == NULL) return;

  if (lc->rclist) 
    mdef_free_recursive_rc (lc->rclist);

  if (lc->next) 
    mdef_free_recursive_lc (lc->next);

  ckd_free ((void *) lc);
}

void mdef_free_recursive_rc (ph_rc_t *rc)
{
  if (rc == NULL) return;

    if (rc->next) 
      mdef_free_recursive_rc (rc->next);

    ckd_free ((void *) rc);
}


/* RAH, Free memory that was allocated in mdef_init 
   Rational purify shows that no leaks exist
 */
   
void mdef_free (mdef_t *m)
{
  int i,j;

  if (m) { 
    if (m->ciphone2n_cd_sen)
      ckd_free    ((void *)m->ciphone2n_cd_sen);
    if (m->sen2cimap)
      ckd_free    ((void *)m->sen2cimap);
    if (m->cd2cisen)
      ckd_free    ((void *)m->cd2cisen);

    /* RAH, go down the ->next list and delete all the pieces */
    for (i=0;i<N_WORD_POSN;i++)
      for (j=0;j<m->n_ciphone;j++) 
	if (m->wpos_ci_lclist[i][j]) {
	  mdef_free_recursive_lc (m->wpos_ci_lclist[i][j]->next);
	  mdef_free_recursive_rc (m->wpos_ci_lclist[i][j]->rclist);
	}
    
    for (i=0;i<N_WORD_POSN;i++)
      for (j=0;j<m->n_ciphone;j++) 
	if (m->wpos_ci_lclist[i][j])  
	  ckd_free ((void *) m->wpos_ci_lclist[i][j]);

    
    if (m->wpos_ci_lclist)
      ckd_free_2d ((void *)m->wpos_ci_lclist);
    if (m->sseq) 
      ckd_free_2d ((void *)m->sseq);
    /* Free phone context */
    if (m->phone) 
      ckd_free    ((void *)m->phone);    
    if (m->ciphone_ht)
      hash_free (m->ciphone_ht);

    for (i=0;i<m->n_ciphone;i++) {
      if (m->ciphone[i].name) 
	ckd_free    ((void *)m->ciphone[i].name);
    }
    if (m->ciphone) 
      ckd_free    ((void *)m->ciphone);
    
    ckd_free    ((void *)m);
  }
}

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	0.000000e+00	4.469710e+08	100.0

482-sphinx3

Exec Time : 1.834394e+09 cycles

Matching : 14.9%

invivodict2pid_dict2pid_build_297

Real time: 2.867914e+07 -- Predicted time: 2.916783e+07

Callcount: 1

invivocont_mgau_mgau_precomp_476

Real time: 2.397192e+08 -- Predicted time: 2.452782e+08

Callcount: 1

invivomdef_parse_tmat_senmap_388

Real time: 1.335019e+08 -- Predicted time: 1.009084e+09

Callcount: 133548

invivosubvq_subvq_maha_precomp_107

Real time: 1.848929e+07 -- Predicted time: 1.856825e+07

Callcount: 1

invivovector_vector_maha_precomp_256

Real time: 1.851739e+07 -- Predicted time: 9.246183e+07

Callcount: 12288

invivomdef_find_ph_rc_157

Real time: 1.933746e+07 -- Predicted time: 9.057115e+07

Callcount: 170512

invivosubvq_subvq_init_270

Real time: 2.926047e+08 -- Predicted time: 0.000000e+00

Callcount: 1

invivomdef_sseq_compress_537

Real time: 3.581604e+07 -- Predicted time: 3.784801e+07

Callcount: 1

invivomdef_mdef_init_689

Real time: 4.408045e+08 -- Predicted time: 0.000000e+00

Callcount: 1

482-sphinx3

Exec Time : 1.834394e+09 cycles

Matching : 14.9%

__invivo__dict2pid_dict2pid_build_297

Real time: 2.867914e+07 -- Predicted time: 2.916783e+07

Callcount: 1

__invivo__cont_mgau_mgau_precomp_476

Real time: 2.397192e+08 -- Predicted time: 2.452782e+08

Callcount: 1

__invivo__mdef_parse_tmat_senmap_388

Real time: 1.335019e+08 -- Predicted time: 1.009084e+09

Callcount: 133548

__invivo__subvq_subvq_maha_precomp_107

Real time: 1.848929e+07 -- Predicted time: 1.856825e+07

Callcount: 1

__invivo__vector_vector_maha_precomp_256

Real time: 1.851739e+07 -- Predicted time: 9.246183e+07

Callcount: 12288

__invivo__mdef_find_ph_rc_157

Real time: 1.933746e+07 -- Predicted time: 9.057115e+07

Callcount: 170512

__invivo__subvq_subvq_init_270

Real time: 2.926047e+08 -- Predicted time: 0.000000e+00

Callcount: 1

__invivo__mdef_sseq_compress_537

Real time: 3.581604e+07 -- Predicted time: 3.784801e+07

Callcount: 1

__invivo__mdef_mdef_init_689

Real time: 4.408045e+08 -- Predicted time: 0.000000e+00

Callcount: 1

invivodict2pid_dict2pid_build_297

invivocont_mgau_mgau_precomp_476

invivomdef_parse_tmat_senmap_388

invivosubvq_subvq_maha_precomp_107

invivovector_vector_maha_precomp_256

invivomdef_find_ph_rc_157

invivosubvq_subvq_init_270

invivomdef_sseq_compress_537

invivomdef_mdef_init_689