4 * Created by Mark Levy on 12/02/2006.
5 * Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
7 This program is free software; you can redistribute it and/or
8 modify it under the terms of the GNU General Public License as
9 published by the Free Software Foundation; either version 2 of the
10 License, or (at your option) any later version. See the file
11 COPYING included with this distribution for more information.
19 #include <time.h> /* to seed random number generator */
21 #include <clapack.h> /* LAPACK for matrix inversion */
23 #include "maths/nan-inf.h"
30 // Using ATLAS C interface to LAPACK
31 #define dgetrf_(m, n, a, lda, ipiv, info) \
32 clapack_dgetrf(CblasColMajor, *m, *n, a, *lda, ipiv)
33 #define dgetri_(n, a, lda, ipiv, work, lwork, info) \
34 clapack_dgetri(CblasColMajor, *n, a, *lda, ipiv)
38 #include <vecLib/cblas.h>
40 #include <cblas.h> /* BLAS for matrix multiplication */
45 model_t* hmm_init(double** x, int T, int L, int N)
51 model = (model_t*) malloc(sizeof(model_t));
54 model->p0 = (double*) malloc(N*sizeof(double));
55 model->a = (double**) malloc(N*sizeof(double*));
56 model->mu = (double**) malloc(N*sizeof(double*));
57 for (i = 0; i < N; i++)
59 model->a[i] = (double*) malloc(N*sizeof(double));
60 model->mu[i] = (double*) malloc(L*sizeof(double));
62 model->cov = (double**) malloc(L*sizeof(double*));
63 for (i = 0; i < L; i++)
64 model->cov[i] = (double*) malloc(L*sizeof(double));
67 double* global_mean = (double*) malloc(L*sizeof(double));
69 /* find global mean */
70 for (d = 0; d < L; d++)
73 for (t = 0; t < T; t++)
74 global_mean[d] += x[t][d];
78 /* calculate global diagonal covariance */
79 for (d = 0; d < L; d++)
81 for (e = 0; e < L; e++)
83 for (t = 0; t < T; t++)
84 model->cov[d][d] += (x[t][d] - global_mean[d]) * (x[t][d] - global_mean[d]);
85 model->cov[d][d] /= T-1;
88 /* set all means close to global mean */
89 for (i = 0; i < N; i++)
91 for (d = 0; d < L; d++)
93 /* add some random noise related to covariance */
94 /* ideally the random number would be Gaussian(0,1), as a hack we make it uniform on [-0.25,0.25] */
95 model->mu[i][d] = global_mean[d] + (0.5 * rand() / (double) RAND_MAX - 0.25) * sqrt(model->cov[d][d]);
99 /* random intial and transition probs */
101 for (i = 0; i < N; i++)
103 model->p0[i] = 1 + rand() / (double) RAND_MAX;
106 for (j = 0; j < N; j++)
108 model->a[i][j] = 1 + rand() / (double) RAND_MAX;
109 ss += model->a[i][j];
111 for (j = 0; j < N; j++)
113 model->a[i][j] /= ss;
116 for (i = 0; i < N; i++)
124 void hmm_close(model_t* model)
128 for (i = 0; i < model->N; i++)
135 for (i = 0; i < model->L; i++)
141 void hmm_train(double** x, int T, model_t* model)
144 double loglik; /* overall log-likelihood at each iteration */
148 double* p0 = model->p0;
149 double** a = model->a;
150 double** mu = model->mu;
151 double** cov = model->cov;
153 /* allocate memory */
154 double** gamma = (double**) malloc(T*sizeof(double*));
155 double*** xi = (double***) malloc(T*sizeof(double**));
156 for (t = 0; t < T; t++)
158 gamma[t] = (double*) malloc(N*sizeof(double));
159 xi[t] = (double**) malloc(N*sizeof(double*));
160 for (i = 0; i < N; i++)
161 xi[t][i] = (double*) malloc(N*sizeof(double));
164 /* temporary memory */
165 double* gauss_y = (double*) malloc(L*sizeof(double));
166 double* gauss_z = (double*) malloc(L*sizeof(double));
168 /* obs probs P(j|{x}) */
169 double** b = (double**) malloc(T*sizeof(double*));
170 for (t = 0; t < T; t++)
171 b[t] = (double*) malloc(N*sizeof(double));
173 /* inverse covariance and its determinant */
174 double** icov = (double**) malloc(L*sizeof(double*));
175 for (i = 0; i < L; i++)
176 icov[i] = (double*) malloc(L*sizeof(double));
179 double thresh = 0.0001;
182 double loglik1, loglik2;
185 while (iter < niter && !foundnan && !(iter > 1 && (loglik - loglik1) < thresh * (loglik1 - loglik2)))
189 fprintf(stderr, "calculating obsprobs...\n");
192 /* precalculate obs probs */
193 invert(cov, L, icov, &detcov);
195 for (t = 0; t < T; t++)
198 for (i = 0; i < N; i++)
200 b[t][i] = exp(loggauss(x[t], L, mu[i], icov, detcov, gauss_y, gauss_z));
208 printf("all the b[t][i] were zero for t = %d, correcting...\n", t);
209 for (i = 0; i < N; i++)
217 fprintf(stderr, "forwards-backwards...\n");
220 forward_backwards(xi, gamma, &loglik, &loglik1, &loglik2, iter, N, T, p0, a, b);
222 fprintf(stderr, "iteration %d: loglik = %f\n", iter, loglik);
223 fprintf(stderr, "re-estimation...\n");
231 baum_welch(p0, a, mu, cov, N, T, L, x, xi, gamma);
235 for (i = 0; i < model->N; i++)
237 for (j = 0; j < model->N; j++)
238 printf("%f ", model->a[i][j]);
246 /* deallocate memory */
247 for (t = 0; t < T; t++)
251 for (i = 0; i < N; i++)
259 for (i = 0; i < L; i++)
267 void mlss_reestimate(double* p0, double** a, double** mu, double** cov, int N, int T, int L, int* q, double** x)
269 /* fit a single Gaussian to observations in each state */
271 /* calculate the mean observation in each state */
273 /* calculate the overall covariance */
275 /* count transitions */
277 /* estimate initial probs from transitions (???) */
280 void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, int L, double** x, double*** xi, double** gamma)
284 double* sum_gamma = (double*) malloc(N*sizeof(double));
286 /* temporary memory */
287 double* u = (double*) malloc(L*L*sizeof(double));
288 double* yy = (double*) malloc(T*L*sizeof(double));
289 double* yy2 = (double*) malloc(T*L*sizeof(double));
291 /* re-estimate transition probs */
292 for (i = 0; i < N; i++)
295 for (t = 0; t < T-1; t++)
296 sum_gamma[i] += gamma[t][i];
299 for (i = 0; i < N; i++)
301 if (sum_gamma[i] == 0)
303 /* fprintf(stderr, "sum_gamma[%d] was zero...\n", i); */
306 for (j = 0; j < N; j++)
309 if (sum_gamma[i] == 0.) continue;
310 for (t = 0; t < T-1; t++)
311 a[i][j] += xi[t][i][j];
313 a[i][j] /= sum_gamma[i];
316 for (j = 0; j < N; j++)
323 /* NB: now we need to sum gamma over all t */
324 for (i = 0; i < N; i++)
325 sum_gamma[i] += gamma[T-1][i];
327 /* re-estimate initial probs */
328 for (i = 0; i < N; i++)
331 /* re-estimate covariance */
333 double sum_sum_gamma = 0;
334 for (i = 0; i < N; i++)
335 sum_sum_gamma += sum_gamma[i];
338 for (d = 0; d < L; d++)
340 for (e = d; e < L; e++)
343 for (t = 0; t < T; t++)
344 for (j = 0; j < N; j++)
345 cov[d][e] += gamma[t][j] * (x[t][d] - mu[j][d]) * (x[t][e] - mu[j][e]);
347 cov[d][e] /= sum_sum_gamma;
349 if (ISNAN(cov[d][e]))
351 printf("cov[%d][%d] was nan\n", d, e);
352 for (j = 0; j < N; j++)
353 for (i = 0; i < L; i++)
355 printf("mu[%d][%d] was nan\n", j, i);
356 for (t = 0; t < T; t++)
357 for (j = 0; j < N; j++)
358 if (ISNAN(gamma[t][j]))
359 printf("gamma[%d][%d] was nan\n", t, j);
364 for (d = 0; d < L; d++)
365 for (e = 0; e < d; e++)
366 cov[d][e] = cov[e][d];
370 for (d = 0; d < L; d++)
371 for (e = 0; e < L; e++)
374 for (j = 0; j < N; j++)
376 for (d = 0; d < L; d++)
377 for (t = 0; t < T; t++)
379 yy[d*T+t] = x[t][d] - mu[j][d];
380 yy2[d*T+t] = gamma[t][j] * (x[t][d] - mu[j][d]);
383 cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, L, L, T, 1.0, yy, T, yy2, T, 0, u, L);
385 for (e = 0; e < L; e++)
386 for (d = 0; d < L; d++)
387 cov[d][e] += u[e*L+d];
390 for (d = 0; d < L; d++)
391 for (e = 0; e < L; e++)
392 cov[d][e] /= T; /* sum_sum_gamma; */
394 //printf("sum_sum_gamma = %f\n", sum_sum_gamma); /* fine, = T IS THIS ALWAYS TRUE with pooled cov?? */
396 /* re-estimate means */
397 for (j = 0; j < N; j++)
399 for (d = 0; d < L; d++)
402 for (t = 0; t < T; t++)
403 mu[j][d] += gamma[t][j] * x[t][d];
404 mu[j][d] /= sum_gamma[j];
408 /* deallocate memory */
415 void forward_backwards(double*** xi, double** gamma, double* loglik, double* loglik1, double* loglik2, int iter, int N, int T, double* p0, double** a, double** b)
417 /* forwards-backwards with scaling */
420 double** alpha = (double**) malloc(T*sizeof(double*));
421 double** beta = (double**) malloc(T*sizeof(double*));
422 for (t = 0; t < T; t++)
424 alpha[t] = (double*) malloc(N*sizeof(double));
425 beta[t] = (double*) malloc(N*sizeof(double));
428 /* scaling coefficients */
429 double* c = (double*) malloc(T*sizeof(double));
431 /* calculate forward probs and scale coefficients */
433 for (i = 0; i < N; i++)
435 alpha[0][i] = p0[i] * b[0][i];
438 //printf("p0[%d] = %f, b[0][%d] = %f\n", i, p0[i], i, b[0][i]);
441 for (i = 0; i < N; i++)
445 //printf("alpha[0][%d] = %f\n", i, alpha[0][i]); /* OK agrees with Matlab */
449 *loglik = -log(c[0]);
453 for (t = 1; t < T; t++)
456 for (j = 0; j < N; j++)
459 for (i = 0; i < N; i++)
460 alpha[t][j] += alpha[t-1][i] * a[i][j];
461 alpha[t][j] *= b[t][j];
469 printf("c[%d] = 0, going to blow up so exiting\n", t);
470 for (i = 0; i < N; i++)
472 fprintf(stderr, "b[%d][%d] was zero\n", t, i);
473 fprintf(stderr, "x[t] was \n");
474 for (i = 0; i < L; i++)
475 fprintf(stderr, "%f ", x[t][i]);
476 fprintf(stderr, "\n\n");
482 for (j = 0; j < N; j++)
485 //printf("c[%d] = %e\n", t, c[t]);
487 *loglik -= log(c[t]);
490 /* calculate backwards probs using same coefficients */
491 for (i = 0; i < N; i++)
496 for (i = 0; i < N; i++)
502 for (i = 0; i < N; i++)
505 for (j = 0; j < N; j++)
506 beta[t-1][i] += a[i][j] * b[t][j] * beta[t][j];
514 for (t = 0; t < T; t++)
516 for (i = 0; i < N; i++)
517 printf("%4.4e\t\t", alpha[t][i]);
520 printf("\n\n");printf("beta:\n");
521 for (t = 0; t < T; t++)
523 for (i = 0; i < N; i++)
524 printf("%4.4e\t\t", beta[t][i]);
530 /* calculate posterior probs */
532 for (t = 0; t < T; t++)
535 for (i = 0; i < N; i++)
537 gamma[t][i] = alpha[t][i] * beta[t][i];
540 for (i = 0; i < N; i++)
544 //printf("gamma[%d][%d] = %f\n", t, i, gamma[t][i]);
548 for (t = 0; t < T-1; t++)
551 for (i = 0; i < N; i++)
553 for (j = 0; j < N; j++)
555 xi[t][i][j] = alpha[t][i] * a[i][j] * b[t+1][j] * beta[t+1][j];
559 for (i = 0; i < N; i++)
560 for (j = 0; j < N; j++)
566 // gamma[t][i] = \sum_j{xi[t][i][j]}
568 for (j = 0; j < N; j++)
570 printf("gamma[3][1] = %f, sum_j(xi[3][1][j]) = %f\n", gamma[3][1], tot);
573 for (t = 0; t < T; t++)
583 void viterbi_decode(double** x, int T, model_t* model, int* q)
591 double* p0 = model->p0;
592 double** a = model->a;
593 double** mu = model->mu;
594 double** cov = model->cov;
596 /* inverse covariance and its determinant */
597 double** icov = (double**) malloc(L*sizeof(double*));
598 for (i = 0; i < L; i++)
599 icov[i] = (double*) malloc(L*sizeof(double));
602 double** logb = (double**) malloc(T*sizeof(double*));
603 double** phi = (double**) malloc(T*sizeof(double*));
604 int** psi = (int**) malloc(T*sizeof(int*));
605 for (t = 0; t < T; t++)
607 logb[t] = (double*) malloc(N*sizeof(double));
608 phi[t] = (double*) malloc(N*sizeof(double));
609 psi[t] = (int*) malloc(N*sizeof(int));
612 /* temporary memory */
613 double* gauss_y = (double*) malloc(L*sizeof(double));
614 double* gauss_z = (double*) malloc(L*sizeof(double));
616 /* calculate observation logprobs */
617 invert(cov, L, icov, &detcov);
618 for (t = 0; t < T; t++)
619 for (i = 0; i < N; i++)
620 logb[t][i] = loggauss(x[t], L, mu[i], icov, detcov, gauss_y, gauss_z);
623 for (i = 0; i < N; i++)
625 phi[0][i] = log(p0[i]) + logb[0][i];
629 for (t = 1; t < T; t++)
631 for (j = 0; j < N; j++)
637 for (i = 0; i < N; i++)
639 if (phi[t-1][i] + log(a[i][j]) > max || !havemax)
641 max = phi[t-1][i] + log(a[i][j]);
642 phi[t][j] = max + logb[t][j];
650 /* find maximising state at time T-1 */
653 for (i = 1; i < N; i++)
655 if (phi[T-1][i] > max)
667 q[t] = psi[t+1][q[t+1]];
671 /* de-allocate memory */
672 for (i = 0; i < L; i++)
675 for (t = 0; t < T; t++)
689 /* invert matrix and calculate determinant using LAPACK */
690 void invert(double** cov, int L, double** icov, double* detcov)
692 /* copy square matrix into a vector in column-major order */
693 double* a = (double*) malloc(L*L*sizeof(double));
696 for (i=0; i < L; i++)
697 a[j*L+i] = cov[i][j];
700 int* ipiv = (int *) malloc(L*L*sizeof(int));
703 /* LU decomposition */
704 ret = dgetrf_(&M, &M, a, &M, ipiv, &ret); /* ret should be zero, negative if cov is singular */
707 fprintf(stderr, "Covariance matrix was singular, couldn't invert\n");
711 /* find determinant */
713 for(i = 0; i < L; i++)
715 // TODO: get this to work!!! If detcov < 0 then cov is bad anyway...
718 for (i = 0; i < L; i++)
727 /* allocate required working storage */
731 dgetri_(&M, a, &M, ipiv, &lwbest, &lwork, &ret);
732 lwork = (int) lwbest;
733 double* work = (double*) malloc(lwork*sizeof(double));
737 dgetri_(&M, a, &M, ipiv, work, &lwork, &ret);
740 for (i=0; i < L; i++)
741 icov[i][j] = a[j*L+i];
749 /* probability of multivariate Gaussian given mean, inverse and determinant of covariance */
750 double gauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z)
754 for (i = 0; i < L; i++)
756 for (i = 0; i < L; i++)
759 //for (j = 0; j < L; j++)
760 // z[i] += icov[i][j] * y[j];
761 z[i] = cblas_ddot(L, &icov[i][0], 1, y, 1);
763 s = cblas_ddot(L, z, 1, y, 1);
764 //for (i = 0; i < L; i++)
767 return exp(-s/2.0) / (pow(2*PI, L/2.0) * sqrt(detcov));
770 /* log probability of multivariate Gaussian given mean, inverse and determinant of covariance */
771 double loggauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z)
776 for (i = 0; i < L; i++)
778 for (i = 0; i < L; i++)
781 //for (j = 0; j < L; j++)
782 // z[i] += icov[i][j] * y[j];
783 z[i] = cblas_ddot(L, &icov[i][0], 1, y, 1);
785 s = cblas_ddot(L, z, 1, y, 1);
786 //for (i = 0; i < L; i++)
789 ret = -0.5 * (s + L * log(2*PI) + log(detcov));
794 printf("loggauss returning infinity\n");
796 printf("loggauss returning -infinity\n");
798 printf("loggauss returning nan\n");
804 void hmm_print(model_t* model)
808 for (i = 0; i < model->N; i++)
809 printf("%f ", model->p0[i]);
812 for (i = 0; i < model->N; i++)
814 for (j = 0; j < model->N; j++)
815 printf("%f ", model->a[i][j]);
820 for (i = 0; i < model->N; i++)
822 for (j = 0; j < model->L; j++)
823 printf("%f ", model->mu[i][j]);
828 for (i = 0; i < model->L; i++)
830 for (j = 0; j < model->L; j++)
831 printf("%f ", model->cov[i][j]);