diff options
Diffstat (limited to 'libs/qm-dsp/hmm/hmm.c')
-rw-r--r-- | libs/qm-dsp/hmm/hmm.c | 837 |
1 files changed, 0 insertions, 837 deletions
diff --git a/libs/qm-dsp/hmm/hmm.c b/libs/qm-dsp/hmm/hmm.c deleted file mode 100644 index fbe202613d..0000000000 --- a/libs/qm-dsp/hmm/hmm.c +++ /dev/null @@ -1,837 +0,0 @@ -/* - * hmm.c - * - * Created by Mark Levy on 12/02/2006. - * Copyright 2006 Centre for Digital Music, Queen Mary, University of London. - - This program is free software; you can redistribute it and/or - modify it under the terms of the GNU General Public License as - published by the Free Software Foundation; either version 2 of the - License, or (at your option) any later version. See the file - COPYING included with this distribution for more information. - * - */ - -#include <stdio.h> -#include <math.h> -#include <stdlib.h> -#include <float.h> -#include <time.h> /* to seed random number generator */ - -#include <clapack.h> /* LAPACK for matrix inversion */ - -#include "maths/nan-inf.h" - -#ifdef ATLAS_ORDER -#define HAVE_ATLAS 1 -#endif - -#ifdef HAVE_ATLAS -// Using ATLAS C interface to LAPACK -#define dgetrf_(m, n, a, lda, ipiv, info) \ - clapack_dgetrf(CblasColMajor, *m, *n, a, *lda, ipiv) -#define dgetri_(n, a, lda, ipiv, work, lwork, info) \ - clapack_dgetri(CblasColMajor, *n, a, *lda, ipiv) -#endif - -#ifdef _MAC_OS_X -#include <vecLib/cblas.h> -#else -#include <cblas.h> /* BLAS for matrix multiplication */ -#endif - -#include "hmm.h" - -model_t* hmm_init(double** x, int T, int L, int N) -{ - int i, j, d, e, t; - double s, ss; - - model_t* model; - model = (model_t*) malloc(sizeof(model_t)); - model->N = N; - model->L = L; - model->p0 = (double*) malloc(N*sizeof(double)); - model->a = (double**) malloc(N*sizeof(double*)); - model->mu = (double**) malloc(N*sizeof(double*)); - for (i = 0; i < N; i++) - { - model->a[i] = (double*) malloc(N*sizeof(double)); - model->mu[i] = (double*) malloc(L*sizeof(double)); - } - model->cov = (double**) malloc(L*sizeof(double*)); - for (i = 0; i < L; i++) - model->cov[i] = (double*) malloc(L*sizeof(double)); - - srand(time(0)); - double* global_mean = (double*) malloc(L*sizeof(double)); - - /* find global mean */ - for (d = 0; d < L; d++) - { - global_mean[d] = 0; - for (t = 0; t < T; t++) - global_mean[d] += x[t][d]; - global_mean[d] /= T; - } - - /* calculate global diagonal covariance */ - for (d = 0; d < L; d++) - { - for (e = 0; e < L; e++) - model->cov[d][e] = 0; - for (t = 0; t < T; t++) - model->cov[d][d] += (x[t][d] - global_mean[d]) * (x[t][d] - global_mean[d]); - model->cov[d][d] /= T-1; - } - - /* set all means close to global mean */ - for (i = 0; i < N; i++) - { - for (d = 0; d < L; d++) - { - /* add some random noise related to covariance */ - /* ideally the random number would be Gaussian(0,1), as a hack we make it uniform on [-0.25,0.25] */ - model->mu[i][d] = global_mean[d] + (0.5 * rand() / (double) RAND_MAX - 0.25) * sqrt(model->cov[d][d]); - } - } - - /* random intial and transition probs */ - s = 0; - for (i = 0; i < N; i++) - { - model->p0[i] = 1 + rand() / (double) RAND_MAX; - s += model->p0[i]; - ss = 0; - for (j = 0; j < N; j++) - { - model->a[i][j] = 1 + rand() / (double) RAND_MAX; - ss += model->a[i][j]; - } - for (j = 0; j < N; j++) - { - model->a[i][j] /= ss; - } - } - for (i = 0; i < N; i++) - model->p0[i] /= s; - - free(global_mean); - - return model; -} - -void hmm_close(model_t* model) -{ - int i; - - for (i = 0; i < model->N; i++) - { - free(model->a[i]); - free(model->mu[i]); - } - free(model->a); - free(model->mu); - for (i = 0; i < model->L; i++) - free(model->cov[i]); - free(model->cov); - free(model); -} - -void hmm_train(double** x, int T, model_t* model) -{ - int i, t; - double loglik; /* overall log-likelihood at each iteration */ - - int N = model->N; - int L = model->L; - double* p0 = model->p0; - double** a = model->a; - double** mu = model->mu; - double** cov = model->cov; - - /* allocate memory */ - double** gamma = (double**) malloc(T*sizeof(double*)); - double*** xi = (double***) malloc(T*sizeof(double**)); - for (t = 0; t < T; t++) - { - gamma[t] = (double*) malloc(N*sizeof(double)); - xi[t] = (double**) malloc(N*sizeof(double*)); - for (i = 0; i < N; i++) - xi[t][i] = (double*) malloc(N*sizeof(double)); - } - - /* temporary memory */ - double* gauss_y = (double*) malloc(L*sizeof(double)); - double* gauss_z = (double*) malloc(L*sizeof(double)); - - /* obs probs P(j|{x}) */ - double** b = (double**) malloc(T*sizeof(double*)); - for (t = 0; t < T; t++) - b[t] = (double*) malloc(N*sizeof(double)); - - /* inverse covariance and its determinant */ - double** icov = (double**) malloc(L*sizeof(double*)); - for (i = 0; i < L; i++) - icov[i] = (double*) malloc(L*sizeof(double)); - double detcov; - - double thresh = 0.0001; - int niter = 50; - int iter = 0; - double loglik1, loglik2; - int foundnan = 0; - - while (iter < niter && !foundnan && !(iter > 1 && (loglik - loglik1) < thresh * (loglik1 - loglik2))) - { - ++iter; -/* - fprintf(stderr, "calculating obsprobs...\n"); - fflush(stderr); -*/ - /* precalculate obs probs */ - invert(cov, L, icov, &detcov); - - for (t = 0; t < T; t++) - { - //int allzero = 1; - for (i = 0; i < N; i++) - { - b[t][i] = exp(loggauss(x[t], L, mu[i], icov, detcov, gauss_y, gauss_z)); - - //if (b[t][i] != 0) - // allzero = 0; - } - /* - if (allzero) - { - printf("all the b[t][i] were zero for t = %d, correcting...\n", t); - for (i = 0; i < N; i++) - { - b[t][i] = 0.00001; - } - } - */ - } -/* - fprintf(stderr, "forwards-backwards...\n"); - fflush(stderr); -*/ - forward_backwards(xi, gamma, &loglik, &loglik1, &loglik2, iter, N, T, p0, a, b); -/* - fprintf(stderr, "iteration %d: loglik = %f\n", iter, loglik); - fprintf(stderr, "re-estimation...\n"); - fflush(stderr); -*/ - if (ISNAN(loglik)) { - foundnan = 1; - continue; - } - - baum_welch(p0, a, mu, cov, N, T, L, x, xi, gamma); - - /* - printf("a:\n"); - for (i = 0; i < model->N; i++) - { - for (j = 0; j < model->N; j++) - printf("%f ", model->a[i][j]); - printf("\n"); - } - printf("\n\n"); - */ - //hmm_print(model); - } - - /* deallocate memory */ - for (t = 0; t < T; t++) - { - free(gamma[t]); - free(b[t]); - for (i = 0; i < N; i++) - free(xi[t][i]); - free(xi[t]); - } - free(gamma); - free(xi); - free(b); - - for (i = 0; i < L; i++) - free(icov[i]); - free(icov); - - free(gauss_y); - free(gauss_z); -} - -void mlss_reestimate(double* p0, double** a, double** mu, double** cov, int N, int T, int L, int* q, double** x) -{ - /* fit a single Gaussian to observations in each state */ - - /* calculate the mean observation in each state */ - - /* calculate the overall covariance */ - - /* count transitions */ - - /* estimate initial probs from transitions (???) */ -} - -void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, int L, double** x, double*** xi, double** gamma) -{ - int i, j, t; - - double* sum_gamma = (double*) malloc(N*sizeof(double)); - - /* temporary memory */ - double* u = (double*) malloc(L*L*sizeof(double)); - double* yy = (double*) malloc(T*L*sizeof(double)); - double* yy2 = (double*) malloc(T*L*sizeof(double)); - - /* re-estimate transition probs */ - for (i = 0; i < N; i++) - { - sum_gamma[i] = 0; - for (t = 0; t < T-1; t++) - sum_gamma[i] += gamma[t][i]; - } - - for (i = 0; i < N; i++) - { - if (sum_gamma[i] == 0) - { -/* fprintf(stderr, "sum_gamma[%d] was zero...\n", i); */ - } - //double s = 0; - for (j = 0; j < N; j++) - { - a[i][j] = 0; - if (sum_gamma[i] == 0.) continue; - for (t = 0; t < T-1; t++) - a[i][j] += xi[t][i][j]; - //s += a[i][j]; - a[i][j] /= sum_gamma[i]; - } - /* - for (j = 0; j < N; j++) - { - a[i][j] /= s; - } - */ - } - - /* NB: now we need to sum gamma over all t */ - for (i = 0; i < N; i++) - sum_gamma[i] += gamma[T-1][i]; - - /* re-estimate initial probs */ - for (i = 0; i < N; i++) - p0[i] = gamma[0][i]; - - /* re-estimate covariance */ - int d, e; - double sum_sum_gamma = 0; - for (i = 0; i < N; i++) - sum_sum_gamma += sum_gamma[i]; - - /* - for (d = 0; d < L; d++) - { - for (e = d; e < L; e++) - { - cov[d][e] = 0; - for (t = 0; t < T; t++) - for (j = 0; j < N; j++) - cov[d][e] += gamma[t][j] * (x[t][d] - mu[j][d]) * (x[t][e] - mu[j][e]); - - cov[d][e] /= sum_sum_gamma; - - if (ISNAN(cov[d][e])) - { - printf("cov[%d][%d] was nan\n", d, e); - for (j = 0; j < N; j++) - for (i = 0; i < L; i++) - if (ISNAN(mu[j][i])) - printf("mu[%d][%d] was nan\n", j, i); - for (t = 0; t < T; t++) - for (j = 0; j < N; j++) - if (ISNAN(gamma[t][j])) - printf("gamma[%d][%d] was nan\n", t, j); - exit(-1); - } - } - } - for (d = 0; d < L; d++) - for (e = 0; e < d; e++) - cov[d][e] = cov[e][d]; - */ - - /* using BLAS */ - for (d = 0; d < L; d++) - for (e = 0; e < L; e++) - cov[d][e] = 0; - - for (j = 0; j < N; j++) - { - for (d = 0; d < L; d++) - for (t = 0; t < T; t++) - { - yy[d*T+t] = x[t][d] - mu[j][d]; - yy2[d*T+t] = gamma[t][j] * (x[t][d] - mu[j][d]); - } - - cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, L, L, T, 1.0, yy, T, yy2, T, 0, u, L); - - for (e = 0; e < L; e++) - for (d = 0; d < L; d++) - cov[d][e] += u[e*L+d]; - } - - for (d = 0; d < L; d++) - for (e = 0; e < L; e++) - cov[d][e] /= T; /* sum_sum_gamma; */ - - //printf("sum_sum_gamma = %f\n", sum_sum_gamma); /* fine, = T IS THIS ALWAYS TRUE with pooled cov?? */ - - /* re-estimate means */ - for (j = 0; j < N; j++) - { - for (d = 0; d < L; d++) - { - mu[j][d] = 0; - for (t = 0; t < T; t++) - mu[j][d] += gamma[t][j] * x[t][d]; - mu[j][d] /= sum_gamma[j]; - } - } - - /* deallocate memory */ - free(sum_gamma); - free(yy); - free(yy2); - free(u); -} - -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) -{ - /* forwards-backwards with scaling */ - int i, j, t; - - double** alpha = (double**) malloc(T*sizeof(double*)); - double** beta = (double**) malloc(T*sizeof(double*)); - for (t = 0; t < T; t++) - { - alpha[t] = (double*) malloc(N*sizeof(double)); - beta[t] = (double*) malloc(N*sizeof(double)); - } - - /* scaling coefficients */ - double* c = (double*) malloc(T*sizeof(double)); - - /* calculate forward probs and scale coefficients */ - c[0] = 0; - for (i = 0; i < N; i++) - { - alpha[0][i] = p0[i] * b[0][i]; - c[0] += alpha[0][i]; - - //printf("p0[%d] = %f, b[0][%d] = %f\n", i, p0[i], i, b[0][i]); - } - c[0] = 1 / c[0]; - for (i = 0; i < N; i++) - { - alpha[0][i] *= c[0]; - - //printf("alpha[0][%d] = %f\n", i, alpha[0][i]); /* OK agrees with Matlab */ - } - - *loglik1 = *loglik; - *loglik = -log(c[0]); - if (iter == 2) - *loglik2 = *loglik; - - for (t = 1; t < T; t++) - { - c[t] = 0; - for (j = 0; j < N; j++) - { - alpha[t][j] = 0; - for (i = 0; i < N; i++) - alpha[t][j] += alpha[t-1][i] * a[i][j]; - alpha[t][j] *= b[t][j]; - - c[t] += alpha[t][j]; - } - - /* - if (c[t] == 0) - { - printf("c[%d] = 0, going to blow up so exiting\n", t); - for (i = 0; i < N; i++) - if (b[t][i] == 0) - fprintf(stderr, "b[%d][%d] was zero\n", t, i); - fprintf(stderr, "x[t] was \n"); - for (i = 0; i < L; i++) - fprintf(stderr, "%f ", x[t][i]); - fprintf(stderr, "\n\n"); - exit(-1); - } - */ - - c[t] = 1 / c[t]; - for (j = 0; j < N; j++) - alpha[t][j] *= c[t]; - - //printf("c[%d] = %e\n", t, c[t]); - - *loglik -= log(c[t]); - } - - /* calculate backwards probs using same coefficients */ - for (i = 0; i < N; i++) - beta[T-1][i] = 1; - t = T - 1; - while (1) - { - for (i = 0; i < N; i++) - beta[t][i] *= c[t]; - - if (t == 0) - break; - - for (i = 0; i < N; i++) - { - beta[t-1][i] = 0; - for (j = 0; j < N; j++) - beta[t-1][i] += a[i][j] * b[t][j] * beta[t][j]; - } - - t--; - } - - /* - printf("alpha:\n"); - for (t = 0; t < T; t++) - { - for (i = 0; i < N; i++) - printf("%4.4e\t\t", alpha[t][i]); - printf("\n"); - } - printf("\n\n");printf("beta:\n"); - for (t = 0; t < T; t++) - { - for (i = 0; i < N; i++) - printf("%4.4e\t\t", beta[t][i]); - printf("\n"); - } - printf("\n\n"); - */ - - /* calculate posterior probs */ - double tot; - for (t = 0; t < T; t++) - { - tot = 0; - for (i = 0; i < N; i++) - { - gamma[t][i] = alpha[t][i] * beta[t][i]; - tot += gamma[t][i]; - } - for (i = 0; i < N; i++) - { - gamma[t][i] /= tot; - - //printf("gamma[%d][%d] = %f\n", t, i, gamma[t][i]); - } - } - - for (t = 0; t < T-1; t++) - { - tot = 0; - for (i = 0; i < N; i++) - { - for (j = 0; j < N; j++) - { - xi[t][i][j] = alpha[t][i] * a[i][j] * b[t+1][j] * beta[t+1][j]; - tot += xi[t][i][j]; - } - } - for (i = 0; i < N; i++) - for (j = 0; j < N; j++) - xi[t][i][j] /= tot; - } - - /* - // CHECK - fine - // gamma[t][i] = \sum_j{xi[t][i][j]} - tot = 0; - for (j = 0; j < N; j++) - tot += xi[3][1][j]; - printf("gamma[3][1] = %f, sum_j(xi[3][1][j]) = %f\n", gamma[3][1], tot); - */ - - for (t = 0; t < T; t++) - { - free(alpha[t]); - free(beta[t]); - } - free(alpha); - free(beta); - free(c); -} - -void viterbi_decode(double** x, int T, model_t* model, int* q) -{ - int i, j, t; - double max; - int havemax; - - int N = model->N; - int L = model->L; - double* p0 = model->p0; - double** a = model->a; - double** mu = model->mu; - double** cov = model->cov; - - /* inverse covariance and its determinant */ - double** icov = (double**) malloc(L*sizeof(double*)); - for (i = 0; i < L; i++) - icov[i] = (double*) malloc(L*sizeof(double)); - double detcov; - - double** logb = (double**) malloc(T*sizeof(double*)); - double** phi = (double**) malloc(T*sizeof(double*)); - int** psi = (int**) malloc(T*sizeof(int*)); - for (t = 0; t < T; t++) - { - logb[t] = (double*) malloc(N*sizeof(double)); - phi[t] = (double*) malloc(N*sizeof(double)); - psi[t] = (int*) malloc(N*sizeof(int)); - } - - /* temporary memory */ - double* gauss_y = (double*) malloc(L*sizeof(double)); - double* gauss_z = (double*) malloc(L*sizeof(double)); - - /* calculate observation logprobs */ - invert(cov, L, icov, &detcov); - for (t = 0; t < T; t++) - for (i = 0; i < N; i++) - logb[t][i] = loggauss(x[t], L, mu[i], icov, detcov, gauss_y, gauss_z); - - /* initialise */ - for (i = 0; i < N; i++) - { - phi[0][i] = log(p0[i]) + logb[0][i]; - psi[0][i] = 0; - } - - for (t = 1; t < T; t++) - { - for (j = 0; j < N; j++) - { - max = -1000000; - havemax = 0; - - psi[t][j] = 0; - for (i = 0; i < N; i++) - { - if (phi[t-1][i] + log(a[i][j]) > max || !havemax) - { - max = phi[t-1][i] + log(a[i][j]); - phi[t][j] = max + logb[t][j]; - psi[t][j] = i; - havemax = 1; - } - } - } - } - - /* find maximising state at time T-1 */ - max = phi[T-1][0]; - q[T-1] = 0; - for (i = 1; i < N; i++) - { - if (phi[T-1][i] > max) - { - max = phi[T-1][i]; - q[T-1] = i; - } - } - - - /* track back */ - t = T - 2; - while (t >= 0) - { - q[t] = psi[t+1][q[t+1]]; - t--; - } - - /* de-allocate memory */ - for (i = 0; i < L; i++) - free(icov[i]); - free(icov); - for (t = 0; t < T; t++) - { - free(logb[t]); - free(phi[t]); - free(psi[t]); - } - free(logb); - free(phi); - free(psi); - - free(gauss_y); - free(gauss_z); -} - -/* invert matrix and calculate determinant using LAPACK */ -void invert(double** cov, int L, double** icov, double* detcov) -{ - /* copy square matrix into a vector in column-major order */ - double* a = (double*) malloc(L*L*sizeof(double)); - int i, j; - for(j=0; j < L; j++) - for (i=0; i < L; i++) - a[j*L+i] = cov[i][j]; - - int M = (int) L; - int* ipiv = (int *) malloc(L*L*sizeof(int)); - int ret; - - /* LU decomposition */ - ret = dgetrf_(&M, &M, a, &M, ipiv, &ret); /* ret should be zero, negative if cov is singular */ - if (ret < 0) - { - fprintf(stderr, "Covariance matrix was singular, couldn't invert\n"); - exit(-1); - } - - /* find determinant */ - double det = 1; - for(i = 0; i < L; i++) - det *= a[i*L+i]; - // TODO: get this to work!!! If detcov < 0 then cov is bad anyway... - /* - int sign = 1; - for (i = 0; i < L; i++) - if (ipiv[i] != i) - sign = -sign; - det *= sign; - */ - if (det < 0) - det = -det; - *detcov = det; - - /* allocate required working storage */ -#ifndef HAVE_ATLAS - int lwork = -1; - double lwbest = 0.0; - dgetri_(&M, a, &M, ipiv, &lwbest, &lwork, &ret); - lwork = (int) lwbest; - double* work = (double*) malloc(lwork*sizeof(double)); -#endif - - /* find inverse */ - dgetri_(&M, a, &M, ipiv, work, &lwork, &ret); - - for(j=0; j < L; j++) - for (i=0; i < L; i++) - icov[i][j] = a[j*L+i]; - -#ifndef HAVE_ATLAS - free(work); -#endif - free(a); -} - -/* probability of multivariate Gaussian given mean, inverse and determinant of covariance */ -double gauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z) -{ - int i, j; - double s = 0; - for (i = 0; i < L; i++) - y[i] = x[i] - mu[i]; - for (i = 0; i < L; i++) - { - //z[i] = 0; - //for (j = 0; j < L; j++) - // z[i] += icov[i][j] * y[j]; - z[i] = cblas_ddot(L, &icov[i][0], 1, y, 1); - } - s = cblas_ddot(L, z, 1, y, 1); - //for (i = 0; i < L; i++) - // s += z[i] * y[i]; - - return exp(-s/2.0) / (pow(2*PI, L/2.0) * sqrt(detcov)); -} - -/* log probability of multivariate Gaussian given mean, inverse and determinant of covariance */ -double loggauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z) -{ - int i, j; - double s = 0; - double ret; - for (i = 0; i < L; i++) - y[i] = x[i] - mu[i]; - for (i = 0; i < L; i++) - { - //z[i] = 0; - //for (j = 0; j < L; j++) - // z[i] += icov[i][j] * y[j]; - z[i] = cblas_ddot(L, &icov[i][0], 1, y, 1); - } - s = cblas_ddot(L, z, 1, y, 1); - //for (i = 0; i < L; i++) - // s += z[i] * y[i]; - - ret = -0.5 * (s + L * log(2*PI) + log(detcov)); - - /* - // TEST - if (ISINF(ret) > 0) - printf("loggauss returning infinity\n"); - if (ISINF(ret) < 0) - printf("loggauss returning -infinity\n"); - if (ISNAN(ret)) - printf("loggauss returning nan\n"); - */ - - return ret; -} - -void hmm_print(model_t* model) -{ - int i, j; - printf("p0:\n"); - for (i = 0; i < model->N; i++) - printf("%f ", model->p0[i]); - printf("\n\n"); - printf("a:\n"); - for (i = 0; i < model->N; i++) - { - for (j = 0; j < model->N; j++) - printf("%f ", model->a[i][j]); - printf("\n"); - } - printf("\n\n"); - printf("mu:\n"); - for (i = 0; i < model->N; i++) - { - for (j = 0; j < model->L; j++) - printf("%f ", model->mu[i][j]); - printf("\n"); - } - printf("\n\n"); - printf("cov:\n"); - for (i = 0; i < model->L; i++) - { - for (j = 0; j < model->L; j++) - printf("%f ", model->cov[i][j]); - printf("\n"); - } - printf("\n\n"); -} - - |