1 /* -*- c-basic-offset: 4; indent-tabs-mode: nil -*- */
2 /* ====================================================================
3  * Copyright (c) 1999-2010 Carnegie Mellon University.  All rights
4  * reserved.
5  *
6  * Redistribution and use in source and binary forms, with or without
7  * modification, are permitted provided that the following conditions
8  * are met:
9  *
10  * 1. Redistributions of source code must retain the above copyright
11  *    notice, this list of conditions and the following disclaimer.
12  *
13  * 2. Redistributions in binary form must reproduce the above copyright
14  *    notice, this list of conditions and the following disclaimer in
15  *    the documentation and/or other materials provided with the
16  *    distribution.
17  *
18  * This work was supported in part by funding from the Defense Advanced
19  * Research Projects Agency and the National Science Foundation of the
20  * United States of America, and the CMU Sphinx Speech Consortium.
21  *
22  * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND
23  * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
24  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
25  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
26  * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
27  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
28  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
29  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
30  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
31  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
32  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33  *
34  * ====================================================================
35  *
36  */
37 
38 /* System headers */
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <string.h>
42 #include <assert.h>
43 #include <limits.h>
44 #include <math.h>
45 #if defined(__ADSPBLACKFIN__)
46 #elif !defined(_WIN32_WCE)
47 #include <sys/types.h>
48 #endif
49 
50 /* SphinxBase headers */
51 #include <sphinx_config.h>
52 #include <sphinxbase/cmd_ln.h>
53 #include <sphinxbase/fixpoint.h>
54 #include <sphinxbase/ckd_alloc.h>
55 #include <sphinxbase/bio.h>
56 #include <sphinxbase/err.h>
57 #include <sphinxbase/prim_type.h>
58 
59 /* Local headers */
60 #include "tied_mgau_common.h"
61 #include "ptm_mgau.h"
62 
63 static ps_mgaufuncs_t ptm_mgau_funcs = {
64     "ptm",
65     ptm_mgau_frame_eval,      /* frame_eval */
66     ptm_mgau_mllr_transform,  /* transform */
67     ptm_mgau_free             /* free */
68 };
69 
70 #define COMPUTE_GMM_MAP(_idx)                           \
71     diff[_idx] = obs[_idx] - mean[_idx];                \
72     sqdiff[_idx] = MFCCMUL(diff[_idx], diff[_idx]);     \
73     compl[_idx] = MFCCMUL(sqdiff[_idx], var[_idx]);
74 #define COMPUTE_GMM_REDUCE(_idx)                \
75     d = GMMSUB(d, compl[_idx]);
76 
77 static void
insertion_sort_topn(ptm_topn_t * topn,int i,int32 d)78 insertion_sort_topn(ptm_topn_t *topn, int i, int32 d)
79 {
80     ptm_topn_t vtmp;
81     int j;
82 
83     topn[i].score = d;
84     if (i == 0)
85         return;
86     vtmp = topn[i];
87     for (j = i - 1; j >= 0 && d > topn[j].score; j--) {
88         topn[j + 1] = topn[j];
89     }
90     topn[j + 1] = vtmp;
91 }
92 
93 static int
eval_topn(ptm_mgau_t * s,int cb,int feat,mfcc_t * z)94 eval_topn(ptm_mgau_t *s, int cb, int feat, mfcc_t *z)
95 {
96     ptm_topn_t *topn;
97     int i, ceplen;
98 
99     topn = s->f->topn[cb][feat];
100     ceplen = s->g->featlen[feat];
101 
102     for (i = 0; i < s->max_topn; i++) {
103         mfcc_t *mean, diff[4], sqdiff[4], compl[4]; /* diff, diff^2, component likelihood */
104         mfcc_t *var, d;
105         mfcc_t *obs;
106         int32 cw, j;
107 
108         cw = topn[i].cw;
109         mean = s->g->mean[cb][feat][0] + cw * ceplen;
110         var = s->g->var[cb][feat][0] + cw * ceplen;
111         d = s->g->det[cb][feat][cw];
112         obs = z;
113         for (j = 0; j < ceplen % 4; ++j) {
114             diff[0] = *obs++ - *mean++;
115             sqdiff[0] = MFCCMUL(diff[0], diff[0]);
116             compl[0] = MFCCMUL(sqdiff[0], *var);
117             d = GMMSUB(d, compl[0]);
118             ++var;
119         }
120         /* We could vectorize this but it's unlikely to make much
121          * difference as the outer loop here isn't very big. */
122         for (;j < ceplen; j += 4) {
123             COMPUTE_GMM_MAP(0);
124             COMPUTE_GMM_MAP(1);
125             COMPUTE_GMM_MAP(2);
126             COMPUTE_GMM_MAP(3);
127             COMPUTE_GMM_REDUCE(0);
128             COMPUTE_GMM_REDUCE(1);
129             COMPUTE_GMM_REDUCE(2);
130             COMPUTE_GMM_REDUCE(3);
131             var += 4;
132             obs += 4;
133             mean += 4;
134         }
135         insertion_sort_topn(topn, i, (int32)d);
136     }
137 
138     return topn[0].score;
139 }
140 
141 /* This looks bad, but it actually isn't.  Less than 1% of eval_cb's
142  * time is spent doing this. */
143 static void
insertion_sort_cb(ptm_topn_t ** cur,ptm_topn_t * worst,ptm_topn_t * best,int cw,int32 intd)144 insertion_sort_cb(ptm_topn_t **cur, ptm_topn_t *worst, ptm_topn_t *best,
145                   int cw, int32 intd)
146 {
147     for (*cur = worst - 1; *cur >= best && intd >= (*cur)->score; --*cur)
148         memcpy(*cur + 1, *cur, sizeof(**cur));
149     ++*cur;
150     (*cur)->cw = cw;
151     (*cur)->score = intd;
152 }
153 
154 static int
eval_cb(ptm_mgau_t * s,int cb,int feat,mfcc_t * z)155 eval_cb(ptm_mgau_t *s, int cb, int feat, mfcc_t *z)
156 {
157     ptm_topn_t *worst, *best, *topn;
158     mfcc_t *mean;
159     mfcc_t *var, *det, *detP, *detE;
160     int32 i, ceplen;
161 
162     best = topn = s->f->topn[cb][feat];
163     worst = topn + (s->max_topn - 1);
164     mean = s->g->mean[cb][feat][0];
165     var = s->g->var[cb][feat][0];
166     det = s->g->det[cb][feat];
167     detE = det + s->g->n_density;
168     ceplen = s->g->featlen[feat];
169 
170     for (detP = det; detP < detE; ++detP) {
171         mfcc_t diff[4], sqdiff[4], compl[4]; /* diff, diff^2, component likelihood */
172         mfcc_t d, thresh;
173         mfcc_t *obs;
174         ptm_topn_t *cur;
175         int32 cw, j;
176 
177         d = *detP;
178         thresh = (mfcc_t) worst->score; /* Avoid int-to-float conversions */
179         obs = z;
180         cw = detP - det;
181 
182         /* Unroll the loop starting with the first dimension(s).  In
183          * theory this might be a bit faster if this Gaussian gets
184          * "knocked out" by C0. In practice not. */
185         for (j = 0; (j < ceplen % 4) && (d >= thresh); ++j) {
186             diff[0] = *obs++ - *mean++;
187             sqdiff[0] = MFCCMUL(diff[0], diff[0]);
188             compl[0] = MFCCMUL(sqdiff[0], *var++);
189             d = GMMSUB(d, compl[0]);
190         }
191         /* Now do 4 dimensions at a time.  You'd think that GCC would
192          * vectorize this?  Apparently not.  And it's right, because
193          * that won't make this any faster, at least on x86-64. */
194         for (; j < ceplen && d >= thresh; j += 4) {
195             COMPUTE_GMM_MAP(0);
196             COMPUTE_GMM_MAP(1);
197             COMPUTE_GMM_MAP(2);
198             COMPUTE_GMM_MAP(3);
199             COMPUTE_GMM_REDUCE(0);
200             COMPUTE_GMM_REDUCE(1);
201             COMPUTE_GMM_REDUCE(2);
202             COMPUTE_GMM_REDUCE(3);
203             var += 4;
204             obs += 4;
205             mean += 4;
206         }
207         if (j < ceplen) {
208             /* terminated early, so not in topn */
209             mean += (ceplen - j);
210             var += (ceplen - j);
211             continue;
212         }
213         if (d < thresh)
214             continue;
215         for (i = 0; i < s->max_topn; i++) {
216             /* already there, so don't need to insert */
217             if (topn[i].cw == cw)
218                 break;
219         }
220         if (i < s->max_topn)
221             continue;       /* already there.  Don't insert */
222         insertion_sort_cb(&cur, worst, best, cw, (int32)d);
223     }
224 
225     return best->score;
226 }
227 
228 /**
229  * Compute top-N densities for active codebooks (and prune)
230  */
231 static int
ptm_mgau_codebook_eval(ptm_mgau_t * s,mfcc_t ** z,int frame)232 ptm_mgau_codebook_eval(ptm_mgau_t *s, mfcc_t **z, int frame)
233 {
234     int i, j;
235 
236     /* First evaluate top-N from previous frame. */
237     for (i = 0; i < s->g->n_mgau; ++i)
238         for (j = 0; j < s->g->n_feat; ++j)
239             eval_topn(s, i, j, z[j]);
240 
241     /* If frame downsampling is in effect, possibly do nothing else. */
242     if (frame % s->ds_ratio)
243         return 0;
244 
245     /* Evaluate remaining codebooks. */
246     for (i = 0; i < s->g->n_mgau; ++i) {
247         if (bitvec_is_clear(s->f->mgau_active, i))
248             continue;
249         for (j = 0; j < s->g->n_feat; ++j) {
250             eval_cb(s, i, j, z[j]);
251         }
252     }
253 
254     /* Normalize densities to produce "posterior probabilities",
255      * i.e. things with a reasonable dynamic range, then scale and
256      * clamp them to the acceptable range.  This is actually done
257      * solely to ensure that we can use fast_logmath_add().  Note that
258      * unless we share the same normalizer across all codebooks for
259      * each feature stream we get defective scores (that's why these
260      * loops are inside out - doing it per-feature should give us
261      * greater precision). */
262     for (j = 0; j < s->g->n_feat; ++j) {
263         int32 norm = 0x7fffffff;
264         for (i = 0; i < s->g->n_mgau; ++i) {
265             if (bitvec_is_clear(s->f->mgau_active, i))
266                 continue;
267             if (norm > s->f->topn[i][j][0].score >> SENSCR_SHIFT)
268                 norm = s->f->topn[i][j][0].score >> SENSCR_SHIFT;
269         }
270         assert(norm != 0x7fffffff);
271         for (i = 0; i < s->g->n_mgau; ++i) {
272             int32 k;
273             if (bitvec_is_clear(s->f->mgau_active, i))
274                 continue;
275             for (k = 0; k < s->max_topn; ++k) {
276                 s->f->topn[i][j][k].score >>= SENSCR_SHIFT;
277                 s->f->topn[i][j][k].score -= norm;
278                 s->f->topn[i][j][k].score = -s->f->topn[i][j][k].score;
279                 if (s->f->topn[i][j][k].score > MAX_NEG_ASCR)
280                     s->f->topn[i][j][k].score = MAX_NEG_ASCR;
281             }
282         }
283     }
284 
285     return 0;
286 }
287 
288 static int
ptm_mgau_calc_cb_active(ptm_mgau_t * s,uint8 * senone_active,int32 n_senone_active,int compallsen)289 ptm_mgau_calc_cb_active(ptm_mgau_t *s, uint8 *senone_active,
290                         int32 n_senone_active, int compallsen)
291 {
292     int i, lastsen;
293 
294     if (compallsen) {
295         bitvec_set_all(s->f->mgau_active, s->g->n_mgau);
296         return 0;
297     }
298     bitvec_clear_all(s->f->mgau_active, s->g->n_mgau);
299     for (lastsen = i = 0; i < n_senone_active; ++i) {
300         int sen = senone_active[i] + lastsen;
301         int cb = s->sen2cb[sen];
302         bitvec_set(s->f->mgau_active, cb);
303         lastsen = sen;
304     }
305     E_DEBUG(1, ("Active codebooks:"));
306     for (i = 0; i < s->g->n_mgau; ++i) {
307         if (bitvec_is_clear(s->f->mgau_active, i))
308             continue;
309         E_DEBUGCONT(1, (" %d", i));
310     }
311     E_DEBUGCONT(1, ("\n"));
312     return 0;
313 }
314 
315 /**
316  * Compute senone scores from top-N densities for active codebooks.
317  */
318 static int
ptm_mgau_senone_eval(ptm_mgau_t * s,int16 * senone_scores,uint8 * senone_active,int32 n_senone_active,int compall)319 ptm_mgau_senone_eval(ptm_mgau_t *s, int16 *senone_scores,
320                      uint8 *senone_active, int32 n_senone_active,
321                      int compall)
322 {
323     int i, lastsen, bestscore;
324 
325     memset(senone_scores, 0, s->n_sen * sizeof(*senone_scores));
326     /* FIXME: This is the non-cache-efficient way to do this.  We want
327      * to evaluate one codeword at a time but this requires us to have
328      * a reverse codebook to senone mapping, which we don't have
329      * (yet), since different codebooks have different top-N
330      * codewords. */
331     if (compall)
332         n_senone_active = s->n_sen;
333     bestscore = 0x7fffffff;
334     for (lastsen = i = 0; i < n_senone_active; ++i) {
335         int sen, f, cb;
336         int ascore;
337 
338         if (compall)
339             sen = i;
340         else
341             sen = senone_active[i] + lastsen;
342         lastsen = sen;
343         cb = s->sen2cb[sen];
344 
345         if (bitvec_is_clear(s->f->mgau_active, cb)) {
346             int j;
347             /* Because senone_active is deltas we can't really "knock
348              * out" senones from pruned codebooks, and in any case,
349              * it wouldn't make any difference to the search code,
350              * which doesn't expect senone_active to change. */
351             for (f = 0; f < s->g->n_feat; ++f) {
352                 for (j = 0; j < s->max_topn; ++j) {
353                     s->f->topn[cb][f][j].score = MAX_NEG_ASCR;
354                 }
355             }
356         }
357         /* For each feature, log-sum codeword scores + mixw to get
358          * feature density, then sum (multiply) to get ascore */
359         ascore = 0;
360         for (f = 0; f < s->g->n_feat; ++f) {
361             ptm_topn_t *topn;
362             int j, fden = 0;
363             topn = s->f->topn[cb][f];
364             for (j = 0; j < s->max_topn; ++j) {
365                 int mixw;
366                 /* Find mixture weight for this codeword. */
367                 if (s->mixw_cb) {
368                     int dcw = s->mixw[f][topn[j].cw][sen/2];
369                     dcw = (dcw & 1) ? dcw >> 4 : dcw & 0x0f;
370                     mixw = s->mixw_cb[dcw];
371                 }
372                 else {
373                     mixw = s->mixw[f][topn[j].cw][sen];
374                 }
375                 if (j == 0)
376                     fden = mixw + topn[j].score;
377                 else
378                     fden = fast_logmath_add(s->lmath_8b, fden,
379                                        mixw + topn[j].score);
380                 E_DEBUG(3, ("fden[%d][%d] l+= %d + %d = %d\n",
381                             sen, f, mixw, topn[j].score, fden));
382             }
383             ascore += fden;
384         }
385         if (ascore < bestscore) bestscore = ascore;
386         senone_scores[sen] = ascore;
387     }
388     /* Normalize the scores again (finishing the job we started above
389      * in ptm_mgau_codebook_eval...) */
390     for (i = 0; i < s->n_sen; ++i) {
391         senone_scores[i] -= bestscore;
392     }
393 
394     return 0;
395 }
396 
397 /**
398  * Compute senone scores for the active senones.
399  */
400 int32
ptm_mgau_frame_eval(ps_mgau_t * ps,int16 * senone_scores,uint8 * senone_active,int32 n_senone_active,mfcc_t ** featbuf,int32 frame,int32 compallsen)401 ptm_mgau_frame_eval(ps_mgau_t *ps,
402                     int16 *senone_scores,
403                     uint8 *senone_active,
404                     int32 n_senone_active,
405                     mfcc_t ** featbuf, int32 frame,
406                     int32 compallsen)
407 {
408     ptm_mgau_t *s = (ptm_mgau_t *)ps;
409     int fast_eval_idx;
410 
411     /* Find the appropriate frame in the rotating history buffer
412      * corresponding to the requested input frame.  No bounds checking
413      * is done here, which just means you'll get semi-random crap if
414      * you request a frame in the future or one that's too far in the
415      * past.  Since the history buffer is just used for fast match
416      * that might not be fatal. */
417     fast_eval_idx = frame % s->n_fast_hist;
418     s->f = s->hist + fast_eval_idx;
419     /* Compute the top-N codewords for every codebook, unless this
420      * is a past frame, in which case we already have them (we
421      * hope!) */
422     if (frame >= ps_mgau_base(ps)->frame_idx) {
423         ptm_fast_eval_t *lastf;
424         /* Get the previous frame's top-N information (on the
425          * first frame of the input this is just all WORST_DIST,
426          * no harm in that) */
427         if (fast_eval_idx == 0)
428             lastf = s->hist + s->n_fast_hist - 1;
429         else
430             lastf = s->hist + fast_eval_idx - 1;
431         /* Copy in initial top-N info */
432         memcpy(s->f->topn[0][0], lastf->topn[0][0],
433                s->g->n_mgau * s->g->n_feat * s->max_topn * sizeof(ptm_topn_t));
434         /* Generate initial active codebook list (this might not be
435          * necessary) */
436         ptm_mgau_calc_cb_active(s, senone_active, n_senone_active, compallsen);
437         /* Now evaluate top-N, prune, and evaluate remaining codebooks. */
438         ptm_mgau_codebook_eval(s, featbuf, frame);
439     }
440     /* Evaluate intersection of active senones and active codebooks. */
441     ptm_mgau_senone_eval(s, senone_scores, senone_active,
442                          n_senone_active, compallsen);
443 
444     return 0;
445 }
446 
447 static int32
read_sendump(ptm_mgau_t * s,bin_mdef_t * mdef,char const * file)448 read_sendump(ptm_mgau_t *s, bin_mdef_t *mdef, char const *file)
449 {
450     FILE *fp;
451     char line[1000];
452     int32 i, n, r, c;
453     int32 do_swap, do_mmap;
454     size_t offset;
455     int n_clust = 0;
456     int n_feat = s->g->n_feat;
457     int n_density = s->g->n_density;
458     int n_sen = bin_mdef_n_sen(mdef);
459     int n_bits = 8;
460 
461     s->n_sen = n_sen; /* FIXME: Should have been done earlier */
462     do_mmap = cmd_ln_boolean_r(s->config, "-mmap");
463 
464     if ((fp = fopen(file, "rb")) == NULL)
465         return -1;
466 
467     E_INFO("Loading senones from dump file %s\n", file);
468     /* Read title size, title */
469     if (fread(&n, sizeof(int32), 1, fp) != 1) {
470         E_ERROR_SYSTEM("Failed to read title size from %s", file);
471         goto error_out;
472     }
473     /* This is extremely bogus */
474     do_swap = 0;
475     if (n < 1 || n > 999) {
476         SWAP_INT32(&n);
477         if (n < 1 || n > 999) {
478             E_ERROR("Title length %x in dump file %s out of range\n", n, file);
479             goto error_out;
480         }
481         do_swap = 1;
482     }
483     if (fread(line, sizeof(char), n, fp) != n) {
484         E_ERROR_SYSTEM("Cannot read title");
485         goto error_out;
486     }
487     if (line[n - 1] != '\0') {
488         E_ERROR("Bad title in dump file\n");
489         goto error_out;
490     }
491     E_INFO("%s\n", line);
492 
493     /* Read header size, header */
494     if (fread(&n, sizeof(n), 1, fp) != 1) {
495         E_ERROR_SYSTEM("Failed to read header size from %s", file);
496         goto error_out;
497     }
498     if (do_swap) SWAP_INT32(&n);
499     if (fread(line, sizeof(char), n, fp) != n) {
500         E_ERROR_SYSTEM("Cannot read header");
501         goto error_out;
502     }
503     if (line[n - 1] != '\0') {
504         E_ERROR("Bad header in dump file\n");
505         goto error_out;
506     }
507 
508     /* Read other header strings until string length = 0 */
509     for (;;) {
510         if (fread(&n, sizeof(n), 1, fp) != 1) {
511             E_ERROR_SYSTEM("Failed to read header string size from %s", file);
512             goto error_out;
513         }
514         if (do_swap) SWAP_INT32(&n);
515         if (n == 0)
516             break;
517         if (fread(line, sizeof(char), n, fp) != n) {
518             E_ERROR_SYSTEM("Cannot read header");
519             goto error_out;
520         }
521         /* Look for a cluster count, if present */
522         if (!strncmp(line, "feature_count ", strlen("feature_count "))) {
523             n_feat = atoi(line + strlen("feature_count "));
524         }
525         if (!strncmp(line, "mixture_count ", strlen("mixture_count "))) {
526             n_density = atoi(line + strlen("mixture_count "));
527         }
528         if (!strncmp(line, "model_count ", strlen("model_count "))) {
529             n_sen = atoi(line + strlen("model_count "));
530         }
531         if (!strncmp(line, "cluster_count ", strlen("cluster_count "))) {
532             n_clust = atoi(line + strlen("cluster_count "));
533         }
534         if (!strncmp(line, "cluster_bits ", strlen("cluster_bits "))) {
535             n_bits = atoi(line + strlen("cluster_bits "));
536         }
537     }
538 
539     /* Defaults for #rows, #columns in mixw array. */
540     c = n_sen;
541     r = n_density;
542     if (n_clust == 0) {
543         /* Older mixw files have them here, and they might be padded. */
544         if (fread(&r, sizeof(r), 1, fp) != 1) {
545             E_ERROR_SYSTEM("Cannot read #rows");
546             goto error_out;
547         }
548         if (do_swap) SWAP_INT32(&r);
549         if (fread(&c, sizeof(c), 1, fp) != 1) {
550             E_ERROR_SYSTEM("Cannot read #columns");
551             goto error_out;
552         }
553         if (do_swap) SWAP_INT32(&c);
554         E_INFO("Rows: %d, Columns: %d\n", r, c);
555     }
556 
557     if (n_feat != s->g->n_feat) {
558         E_ERROR("Number of feature streams mismatch: %d != %d\n",
559                 n_feat, s->g->n_feat);
560         goto error_out;
561     }
562     if (n_density != s->g->n_density) {
563         E_ERROR("Number of densities mismatch: %d != %d\n",
564                 n_density, s->g->n_density);
565         goto error_out;
566     }
567     if (n_sen != s->n_sen) {
568         E_ERROR("Number of senones mismatch: %d != %d\n",
569                 n_sen, s->n_sen);
570         goto error_out;
571     }
572 
573     if (!((n_clust == 0) || (n_clust == 15) || (n_clust == 16))) {
574         E_ERROR("Cluster count must be 0, 15, or 16\n");
575         goto error_out;
576     }
577     if (n_clust == 15)
578         ++n_clust;
579 
580     if (!((n_bits == 8) || (n_bits == 4))) {
581         E_ERROR("Cluster count must be 4 or 8\n");
582         goto error_out;
583     }
584 
585     if (do_mmap) {
586             E_INFO("Using memory-mapped I/O for senones\n");
587     }
588     offset = ftell(fp);
589 
590     /* Allocate memory for pdfs (or memory map them) */
591     if (do_mmap) {
592         s->sendump_mmap = mmio_file_read(file);
593         /* Get cluster codebook if any. */
594         if (n_clust) {
595             s->mixw_cb = ((uint8 *) mmio_file_ptr(s->sendump_mmap)) + offset;
596             offset += n_clust;
597         }
598     }
599     else {
600         /* Get cluster codebook if any. */
601         if (n_clust) {
602             s->mixw_cb = ckd_calloc(1, n_clust);
603             if (fread(s->mixw_cb, 1, n_clust, fp) != (size_t) n_clust) {
604                 E_ERROR("Failed to read %d bytes from sendump\n", n_clust);
605                 goto error_out;
606             }
607         }
608     }
609 
610     /* Set up pointers, or read, or whatever */
611     if (s->sendump_mmap) {
612         s->mixw = ckd_calloc_2d(n_feat, n_density, sizeof(*s->mixw));
613         for (n = 0; n < n_feat; n++) {
614             int step = c;
615             if (n_bits == 4)
616                 step = (step + 1) / 2;
617             for (i = 0; i < r; i++) {
618                 s->mixw[n][i] = ((uint8 *) mmio_file_ptr(s->sendump_mmap)) + offset;
619                 offset += step;
620             }
621         }
622     }
623     else {
624         s->mixw = ckd_calloc_3d(n_feat, n_density, n_sen, sizeof(***s->mixw));
625         /* Read pdf values and ids */
626         for (n = 0; n < n_feat; n++) {
627             int step = c;
628             if (n_bits == 4)
629                 step = (step + 1) / 2;
630             for (i = 0; i < r; i++) {
631                 if (fread(s->mixw[n][i], sizeof(***s->mixw), step, fp)
632                     != (size_t) step) {
633                     E_ERROR("Failed to read %d bytes from sendump\n", step);
634                     goto error_out;
635                 }
636             }
637         }
638     }
639 
640     fclose(fp);
641     return 0;
642 error_out:
643     fclose(fp);
644     return -1;
645 }
646 
647 static int32
read_mixw(ptm_mgau_t * s,char const * file_name,double SmoothMin)648 read_mixw(ptm_mgau_t * s, char const *file_name, double SmoothMin)
649 {
650     char **argname, **argval;
651     char eofchk;
652     FILE *fp;
653     int32 byteswap, chksum_present;
654     uint32 chksum;
655     float32 *pdf;
656     int32 i, f, c, n;
657     int32 n_sen;
658     int32 n_feat;
659     int32 n_comp;
660     int32 n_err;
661 
662     E_INFO("Reading mixture weights file '%s'\n", file_name);
663 
664     if ((fp = fopen(file_name, "rb")) == NULL)
665         E_FATAL_SYSTEM("Failed to open mixture file '%s' for reading", file_name);
666 
667     /* Read header, including argument-value info and 32-bit byteorder magic */
668     if (bio_readhdr(fp, &argname, &argval, &byteswap) < 0)
669         E_FATAL("Failed to read header from '%s'\n", file_name);
670 
671     /* Parse argument-value list */
672     chksum_present = 0;
673     for (i = 0; argname[i]; i++) {
674         if (strcmp(argname[i], "version") == 0) {
675             if (strcmp(argval[i], MGAU_MIXW_VERSION) != 0)
676                 E_WARN("Version mismatch(%s): %s, expecting %s\n",
677                        file_name, argval[i], MGAU_MIXW_VERSION);
678         }
679         else if (strcmp(argname[i], "chksum0") == 0) {
680             chksum_present = 1; /* Ignore the associated value */
681         }
682     }
683     bio_hdrarg_free(argname, argval);
684     argname = argval = NULL;
685 
686     chksum = 0;
687 
688     /* Read #senones, #features, #codewords, arraysize */
689     if ((bio_fread(&n_sen, sizeof(int32), 1, fp, byteswap, &chksum) != 1)
690         || (bio_fread(&n_feat, sizeof(int32), 1, fp, byteswap, &chksum) !=
691             1)
692         || (bio_fread(&n_comp, sizeof(int32), 1, fp, byteswap, &chksum) !=
693             1)
694         || (bio_fread(&n, sizeof(int32), 1, fp, byteswap, &chksum) != 1)) {
695         E_FATAL("bio_fread(%s) (arraysize) failed\n", file_name);
696     }
697     if (n_feat != s->g->n_feat)
698         E_FATAL("#Features streams(%d) != %d\n", n_feat, s->g->n_feat);
699     if (n != n_sen * n_feat * n_comp) {
700         E_FATAL
701             ("%s: #float32s(%d) doesn't match header dimensions: %d x %d x %d\n",
702              file_name, i, n_sen, n_feat, n_comp);
703     }
704 
705     /* n_sen = number of mixture weights per codeword, which is
706      * fixed at the number of senones since we have only one codebook.
707      */
708     s->n_sen = n_sen;
709 
710     /* Quantized mixture weight arrays. */
711     s->mixw = ckd_calloc_3d(s->g->n_feat, s->g->n_density,
712                             n_sen, sizeof(***s->mixw));
713 
714     /* Temporary structure to read in floats before conversion to (int32) logs3 */
715     pdf = (float32 *) ckd_calloc(n_comp, sizeof(float32));
716 
717     /* Read senone probs data, normalize, floor, convert to logs3, truncate to 8 bits */
718     n_err = 0;
719     for (i = 0; i < n_sen; i++) {
720         for (f = 0; f < n_feat; f++) {
721             if (bio_fread((void *) pdf, sizeof(float32),
722                           n_comp, fp, byteswap, &chksum) != n_comp) {
723                 E_FATAL("bio_fread(%s) (arraydata) failed\n", file_name);
724             }
725 
726             /* Normalize and floor */
727             if (vector_sum_norm(pdf, n_comp) <= 0.0)
728                 n_err++;
729             vector_floor(pdf, n_comp, SmoothMin);
730             vector_sum_norm(pdf, n_comp);
731 
732             /* Convert to LOG, quantize, and transpose */
733             for (c = 0; c < n_comp; c++) {
734                 int32 qscr;
735 
736                 qscr = -logmath_log(s->lmath_8b, pdf[c]);
737                 if ((qscr > MAX_NEG_MIXW) || (qscr < 0))
738                     qscr = MAX_NEG_MIXW;
739                 s->mixw[f][c][i] = qscr;
740             }
741         }
742     }
743     if (n_err > 0)
744         E_WARN("Weight normalization failed for %d mixture weights components\n", n_err);
745 
746     ckd_free(pdf);
747 
748     if (chksum_present)
749         bio_verify_chksum(fp, byteswap, chksum);
750 
751     if (fread(&eofchk, 1, 1, fp) == 1)
752         E_FATAL("More data than expected in %s\n", file_name);
753 
754     fclose(fp);
755 
756     E_INFO("Read %d x %d x %d mixture weights\n", n_sen, n_feat, n_comp);
757     return n_sen;
758 }
759 
760 ps_mgau_t *
ptm_mgau_init(acmod_t * acmod,bin_mdef_t * mdef)761 ptm_mgau_init(acmod_t *acmod, bin_mdef_t *mdef)
762 {
763     ptm_mgau_t *s;
764     ps_mgau_t *ps;
765     char const *sendump_path;
766     int i;
767 
768     s = ckd_calloc(1, sizeof(*s));
769     s->config = acmod->config;
770 
771     s->lmath = logmath_retain(acmod->lmath);
772     /* Log-add table. */
773     s->lmath_8b = logmath_init(logmath_get_base(acmod->lmath), SENSCR_SHIFT, TRUE);
774     if (s->lmath_8b == NULL)
775         goto error_out;
776     /* Ensure that it is only 8 bits wide so that fast_logmath_add() works. */
777     if (logmath_get_width(s->lmath_8b) != 1) {
778         E_ERROR("Log base %f is too small to represent add table in 8 bits\n",
779                 logmath_get_base(s->lmath_8b));
780         goto error_out;
781     }
782 
783     /* Read means and variances. */
784     if ((s->g = gauden_init(cmd_ln_str_r(s->config, "-mean"),
785                             cmd_ln_str_r(s->config, "-var"),
786                             cmd_ln_float32_r(s->config, "-varfloor"),
787                             s->lmath)) == NULL)
788         goto error_out;
789     /* We only support 256 codebooks or less (like 640k or 2GB, this
790      * should be enough for anyone) */
791     if (s->g->n_mgau > 256) {
792         E_INFO("Number of codebooks exceeds 256: %d\n", s->g->n_mgau);
793         goto error_out;
794     }
795     if (s->g->n_mgau != bin_mdef_n_ciphone(mdef)) {
796         E_INFO("Number of codebooks doesn't match number of ciphones, doesn't look like PTM: %d != %d\n", s->g->n_mgau, bin_mdef_n_ciphone(mdef));
797         goto error_out;
798     }
799     /* Verify n_feat and veclen, against acmod. */
800     if (s->g->n_feat != feat_dimension1(acmod->fcb)) {
801         E_ERROR("Number of streams does not match: %d != %d\n",
802                 s->g->n_feat, feat_dimension1(acmod->fcb));
803         goto error_out;
804     }
805     for (i = 0; i < s->g->n_feat; ++i) {
806         if (s->g->featlen[i] != feat_dimension2(acmod->fcb, i)) {
807             E_ERROR("Dimension of stream %d does not match: %d != %d\n",
808                     s->g->featlen[i], feat_dimension2(acmod->fcb, i));
809             goto error_out;
810         }
811     }
812     /* Read mixture weights. */
813     if ((sendump_path = cmd_ln_str_r(s->config, "-sendump"))) {
814         if (read_sendump(s, acmod->mdef, sendump_path) < 0) {
815             goto error_out;
816         }
817     }
818     else {
819         if (read_mixw(s, cmd_ln_str_r(s->config, "-mixw"),
820                       cmd_ln_float32_r(s->config, "-mixwfloor")) < 0) {
821             goto error_out;
822         }
823     }
824     s->ds_ratio = cmd_ln_int32_r(s->config, "-ds");
825     s->max_topn = cmd_ln_int32_r(s->config, "-topn");
826     E_INFO("Maximum top-N: %d\n", s->max_topn);
827 
828     /* Assume mapping of senones to their base phones, though this
829      * will become more flexible in the future. */
830     s->sen2cb = ckd_calloc(s->n_sen, sizeof(*s->sen2cb));
831     for (i = 0; i < s->n_sen; ++i)
832         s->sen2cb[i] = bin_mdef_sen2cimap(acmod->mdef, i);
833 
834     /* Allocate fast-match history buffers.  We need enough for the
835      * phoneme lookahead window, plus the current frame, plus one for
836      * good measure? (FIXME: I don't remember why) */
837     s->n_fast_hist = cmd_ln_int32_r(s->config, "-pl_window") + 2;
838     s->hist = ckd_calloc(s->n_fast_hist, sizeof(*s->hist));
839     /* s->f will be a rotating pointer into s->hist. */
840     s->f = s->hist;
841     for (i = 0; i < s->n_fast_hist; ++i) {
842         int j, k, m;
843         /* Top-N codewords for every codebook and feature. */
844         s->hist[i].topn = ckd_calloc_3d(s->g->n_mgau, s->g->n_feat,
845                                         s->max_topn, sizeof(ptm_topn_t));
846         /* Initialize them to sane (yet arbitrary) defaults. */
847         for (j = 0; j < s->g->n_mgau; ++j) {
848             for (k = 0; k < s->g->n_feat; ++k) {
849                 for (m = 0; m < s->max_topn; ++m) {
850                     s->hist[i].topn[j][k][m].cw = m;
851                     s->hist[i].topn[j][k][m].score = WORST_DIST;
852                 }
853             }
854         }
855         /* Active codebook mapping (just codebook, not features,
856            at least not yet) */
857         s->hist[i].mgau_active = bitvec_alloc(s->g->n_mgau);
858         /* Start with them all on, prune them later. */
859         bitvec_set_all(s->hist[i].mgau_active, s->g->n_mgau);
860     }
861 
862     ps = (ps_mgau_t *)s;
863     ps->vt = &ptm_mgau_funcs;
864     return ps;
865 error_out:
866     ptm_mgau_free(ps_mgau_base(s));
867     return NULL;
868 }
869 
870 int
ptm_mgau_mllr_transform(ps_mgau_t * ps,ps_mllr_t * mllr)871 ptm_mgau_mllr_transform(ps_mgau_t *ps,
872                             ps_mllr_t *mllr)
873 {
874     ptm_mgau_t *s = (ptm_mgau_t *)ps;
875     return gauden_mllr_transform(s->g, mllr, s->config);
876 }
877 
878 void
ptm_mgau_free(ps_mgau_t * ps)879 ptm_mgau_free(ps_mgau_t *ps)
880 {
881     ptm_mgau_t *s = (ptm_mgau_t *)ps;
882 
883     logmath_free(s->lmath);
884     logmath_free(s->lmath_8b);
885     if (s->sendump_mmap) {
886         ckd_free_2d(s->mixw);
887         mmio_file_unmap(s->sendump_mmap);
888     }
889     else {
890         ckd_free_3d(s->mixw);
891     }
892     ckd_free(s->sen2cb);
893     gauden_free(s->g);
894     ckd_free(s);
895 }
896