apaam/mlpack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

HMM codes

Browse Source
This commit is contained in:
tqlong
2008-01-17 22:34:47 +00:00
parent 48e66040d6
commit fa472f35ba
13 changed files with 706 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/u/tqlong/hmm/build.py
+8
View File
@@ -26,6 +26,14 @@ binrule(
deplibs = ["fastlib:fastlib_int"] # depends on example in this folder
)
binrule(
name = "hmm", # the executable name
sources = ["hmm.cc","support.cc","discreteHMM.cc","gaussianHMM.cc","mixgaussHMM.cc","mixtureDST.cc"],
headers = ["support.h", "discreteHMM.h","gaussianHMM.h","mixgaussHMM.h","mixtureDST.h"],
deplibs = ["fastlib:fastlib_int"] # depends on example in this folder
)
# to build:
# 1. make sure have environment variables set up:
# $ source /full/path/to/fastlib/script/fl-env /full/path/to/fastlib
fastlib/u/tqlong/hmm/discreteHMM.cc
+70
View File
@@ -181,6 +181,10 @@ double hmm_decodeD(const Vector& seq, const Matrix& trans, const Matrix& emis, M
double hmm_viterbiD_init(const Vector& seq, const Matrix& trans, const Matrix& emis, Vector* states) {
int L = seq.length();
return hmm_viterbiD_init(L, seq, trans, emis, states);
}
double hmm_viterbiD_init(int L, const Vector& seq, const Matrix& trans, const Matrix& emis, Vector* states) {
int M = trans.n_rows();
int N = emis.n_cols();
DEBUG_ASSERT_MSG((M==trans.n_cols() && M==emis.n_rows()),"hmm_viterbiD: sizes do not match");
@@ -312,3 +316,69 @@ void hmm_trainD(const ArrayList<Vector>& seqs, Matrix* guessTR, Matrix* guessEM,
oldlog = loglik;
}
}
void hmm_train_viterbiD(const ArrayList<Vector>& seqs, Matrix* guessTR, Matrix* guessEM, int max_iter, double tol) {
int L = -1;
int M = guessTR->n_rows();
int N = guessEM->n_cols();
DEBUG_ASSERT_MSG((M==guessTR->n_cols() && M==guessEM->n_rows()),"hmm_trainD: sizes do not match");
for (int i = 0; i < seqs.size(); i++)
if (seqs[i].length() > L) L = seqs[i].length();
Matrix &gTR = *guessTR, &gEM = *guessEM;
Matrix TR, EM; // guess transition and emission matrix
TR.Init(M, M);
EM.Init(M, N);
double loglik = 0, oldlog;
for (int iter = 0; iter < max_iter; iter++) {
oldlog = loglik;
loglik = 0;
TR.SetAll(1e-4);
EM.SetAll(1e-4);
for (int idx = 0; idx < seqs.size(); idx++) {
Vector states;
L = seqs[idx].length();
loglik += hmm_viterbiD_init(L, seqs[idx], gTR, gEM, &states);
for (int t = 0; t < L-1; t++) {
int i = (int) states[t];
int j = (int) states[t+1];
TR.ref(i, j) ++;
}
for (int t = 0; t < L; t++) {
int e = (int) seqs[idx][t];
int i = (int) states[t];
EM.ref(i, e) ++;
}
}
double s;
print_matrix(TR, "TR");
for (int i = 0; i < M; i++) {
s = 0;
for (int j = 0; j < M; j++) s += TR.get(i, j);
if (s == 0) {
for (int j = 0; j < M; j++) gTR.ref(i, j) = 0;
gTR.ref(i, i) = 1;
}
else {
for (int j = 0; j < M; j++) gTR.ref(i, j) = TR.get(i, j) / s;
}
s = 0;
for (int j = 0; j < N; j++) s += EM.get(i, j);
for (int j = 0; j < N; j++) gEM.ref(i, j) = EM.get(i, j) / s;
}
printf("Iter = %d Loglik = %8.4f\n", iter, loglik);
if (fabs(oldlog - loglik) < tol) {
printf("\nConverged after %d iterations\n", iter);
break;
}
oldlog = loglik;
}
}
fastlib/u/tqlong/hmm/discreteHMM.h
+5 -2
View File
@@ -44,11 +44,14 @@ double hmm_decodeD(const Vector& seq, const Matrix& trans, const Matrix& emis, M
RETURN: log probability of the most probable sequence
*/
double hmm_viterbiD_init(const Vector& seq, const Matrix& trans, const Matrix& emis, Vector* states);
double hmm_viterbiD_init(int L, const Vector& seq, const Matrix& trans, const Matrix& emis, Vector* states);
/** Baum-Welch estimation of transition and emission probabilities
*/
void hmm_trainD(const ArrayList<Vector>& seqs, Matrix* guessTR, Matrix* guessEM, int max_iter = 500, double tol = 1e-3);
/** Viterbi estimation of transition and emission probabilities
*/
void hmm_train_viterbiD(const ArrayList<Vector>& seqs, Matrix* guessTR, Matrix* guessEM, int max_iter = 500, double tol = 1e-3);
#endif
fastlib/u/tqlong/hmm/gaussianHMM.cc
+223
View File
@@ -2,6 +2,52 @@
#include "support.h"
#include "gaussianHMM.h"
success_t load_profileG(const char* profile, Matrix* trans, ArrayList<Vector>* means, ArrayList<Matrix>* covs) {
ArrayList<Matrix> matlst;
if (!PASSED(load_matrix_list(profile, &matlst))) {
NONFATAL("Couldn't open '%s' for reading.", profile);
return SUCCESS_FAIL;
}
DEBUG_ASSERT(matlst.size() > 0);
trans->Copy(matlst[0]);
means->Init();
covs->Init();
int M = trans->n_rows(); // num of states
DEBUG_ASSERT(matlst.size() == 2*M+1);
int N = matlst[1].n_rows(); // dimension
for (int i = 1; i < 2*M+1; i+=2) {
DEBUG_ASSERT(matlst[i].n_rows()==N && matlst[i].n_cols()==1);
DEBUG_ASSERT(matlst[i+1].n_rows()==N && matlst[i+1].n_cols()==N);
Vector m;
matlst[i].MakeColumnVector(0, &m);
means->AddBackItem(m);
covs->AddBackItem(matlst[i+1]);
}
return SUCCESS_PASS;
}
success_t save_profileG(const char* profile, const Matrix& trans, const ArrayList<Vector>& means, const ArrayList<Matrix>& covs) {
TextWriter w_pro;
if (!PASSED(w_pro.Open(profile))) {
NONFATAL("Couldn't open '%s' for writing.", profile);
return SUCCESS_FAIL;
}
int M = trans.n_rows(); // num of states
DEBUG_ASSERT(means.size() == M && covs.size() == M);
int N = means[0].length(); // dimension
print_matrix(w_pro, trans, "% transmission", "%E,");
for (int i = 0; i < M; i++) {
DEBUG_ASSERT(means[i].length() == N);
DEBUG_ASSERT(covs[i].n_rows()==N && covs[i].n_cols()==N);
char s[100];
sprintf(s, "%% mean - state %d", i);
print_vector(w_pro, means[i], s, "%E,");
sprintf(s, "%% covariance - state%d", i);
print_matrix(w_pro, covs[i], s, "%E,");
}
return SUCCESS_PASS;
}
void hmm_generateG_init(int L, const Matrix& trans, const ArrayList<Vector>& means, const ArrayList<Matrix>& covs, Matrix* seq, Vector* states){
DEBUG_ASSERT_MSG((trans.n_rows()==trans.n_cols() && trans.n_rows()==means.size() && trans.n_rows()==covs.size()), "hmm_generateG_init: matrices sizes do not match");
Matrix trsum;
@@ -200,6 +246,10 @@ double hmm_decodeG(const Matrix& trans, const Matrix& emis_prob, Matrix* pstates
double hmm_viterbiG_init(const Matrix& trans, const Matrix& emis_prob, Vector* states) {
int L = emis_prob.n_cols();
return hmm_viterbiG_init(L, trans, emis_prob, states);
}
double hmm_viterbiG_init(int L, const Matrix& trans, const Matrix& emis_prob, Vector* states) {
int M = trans.n_rows();
DEBUG_ASSERT_MSG((M==trans.n_cols() && M==emis_prob.n_rows()),"hmm_viterbiG: sizes do not match");
@@ -255,6 +305,179 @@ double hmm_viterbiG_init(const Matrix& trans, const Matrix& emis_prob, Vector* s
return bestVal;
}
void hmm_cal_emis_prob(const Matrix& seq, const ArrayList<Vector>& means, const ArrayList<Matrix>& inv_covs, const Vector& det, Matrix* emis_prob) {
int L = seq.n_cols();
int M = means.size();
for (int t = 0; t < L; t++) {
Vector e;
seq.MakeColumnVector(t, &e);
for (int i = 0; i < M; i++)
emis_prob->ref(i, t) = NORMAL_DENSITY(e, means[i], inv_covs[i], det[i]);
}
}
void init_gauss_param(int M, const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO) {
int N = seqs[0].n_rows();
Matrix& gTR = *guessTR;
ArrayList<Vector>& gME = *guessME;
ArrayList<Matrix>& gCO = *guessCO;
ArrayList<int> labels;
Vector sumState;
kmeans(seqs, M, &labels, &gME, 1000, 1e-5);
//for (int i = 0; i < labels.size(); i++) printf("%8d", labels[i]);
//printf("---1---\n");
gTR.Init(M, M); gTR.SetZero();
sumState.Init(M); sumState.SetZero();
gCO.Init();
for (int i = 0; i < M; i++) {
Matrix m;
m.Init(N, N); m.SetZero();
gCO.AddBackItem(m);
}
//printf("---2---\n");
int t = 0;
for (int p=0; p < seqs.size(); p++) {
for (int q=0; q < seqs[p].n_cols(); q++,t++) {
if (q == seqs[p].n_cols() -1) continue;
int i = labels[t];
int j = labels[t+1];
gTR.ref(i, j)++;
sumState[i]++;
Vector data_j_Vec, sub_Vec;
Matrix tmp_cov;
seqs[p].MakeColumnVector(q, &data_j_Vec);
la::SubInit(gME[i], data_j_Vec, &sub_Vec);
tmp_cov.AliasColVector(sub_Vec);
//printf("t = %d x = %8.3f\n", t, sub_Vec[0]);
la::MulExpert(1, false, tmp_cov, true, tmp_cov, 1, &gCO[i]);
}
}
for (int i = 0; i < M; i++)
if (sumState[i] == 0) {
for (int j = 0; j < M; j++) gTR.ref(i, j) = 0;
gTR.ref(i, i) = 1;
gME[i].SetZero();
gCO[i].SetZero();
for (int j = 0; j < N; j++) gCO[i].ref(j, j) = 1;
}
else {
for (int j = 0; j < M; j++) gTR.ref(i, j) /= sumState[i];
la::Scale(1.0/sumState[i], &gCO[i]);
for (int j = 0; j < N; j++) gCO[i].ref(j, j) += 1e-3; // make sure the diagonal elements are not too small
}
}
void hmm_train_viterbiG(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO, int max_iter, double tol) {
Matrix &gTR = *guessTR;
ArrayList<Vector>& gME = *guessME;
ArrayList<Matrix>& gCO = *guessCO;
int L = -1;
int M = gTR.n_rows();
int N = gME[0].length();
DEBUG_ASSERT_MSG((M==gTR.n_cols() && M==gME.size() && M == gCO.size()),"hmm_trainD: sizes do not match");
for (int i = 0; i < seqs.size(); i++)
if (seqs[i].n_cols() > L) L = seqs[i].n_cols();
Matrix TR; // accumulating transition
ArrayList<Vector> ME; // accumulating mean
ArrayList<Matrix> CO; // accumulating covariance
ArrayList<Matrix> INV_CO; // inverse matrix of the covariance
Vector DET; // the determinant * constant of the Normal PDF formula
TR.Init(M, M);
ME.Copy(gME);
CO.Copy(gCO);
INV_CO.Copy(CO);
DET.Init(M);
Matrix emis_prob;
Vector sumState; // the denominator for each state
emis_prob.Init(M, L);
sumState.Init(M);
double loglik = 0, oldlog;
for (int iter = 0; iter < max_iter; iter++) {
oldlog = loglik;
loglik = 0;
// set the accumulating values to zeros and compute the inverse matrices and determinant constants
TR.SetZero();
for (int i = 0; i < M; i++) {
ME[i].SetZero();
CO[i].SetZero();
la::InverseOverwrite(gCO[i], &INV_CO[i]);
DET[i] = pow(2.0*math::PI, -N/2.0) * pow(la::Determinant(gCO[i]), -0.5);
}
sumState.SetZero();
// for each sequence, we will use forward-backward procedure and then accumulate
for (int idx = 0; idx < seqs.size(); idx++) {
L = seqs[idx].n_cols();
Vector states;
hmm_cal_emis_prob(seqs[idx], gME, INV_CO, DET, &emis_prob); // first calculate the emission probabilities of the sequence
loglik += hmm_viterbiG_init(L, gTR, emis_prob, &states); // get the most probable state sequence
// accumulate expected transition & mean & covariance
for (int t = 0; t < L-1; t++) {
int i = (int) states[t];
int j = (int) states[t+1];
TR.ref(i, j) ++;
}
for (int t = 0; t < L; t++) {
Vector e;
seqs[idx].MakeColumnVector(t, &e);
int i = (int) states[t];
sumState[i] ++;
la::AddTo(e, &ME[i]);
Vector d;
la::SubInit(e, gME[i], &d);
Matrix D;
D.AliasColVector(d);
la::MulExpert(1, false, D, true, D, 1.0, &CO[i]);
}
// end accumulate
}
// after accumulate all sequences: re-estimate transition & mean & covariance for the next iteration
for (int i = 0; i < M; i++) {
double s = 0;
for (int j = 0; j < M; j++) s += TR.get(i, j);
if (s == 0) {
for (int j = 0; j < M; j++) gTR.ref(i, j) = 0;
gTR.ref(i, i) = 1;
}
else {
for (int j = 0; j < M; j++) gTR.ref(i, j) = TR.get(i, j) / s;
}
if (sumState[i] != 0) {
la::ScaleOverwrite(1.0/sumState[i], ME[i], &gME[i]);
la::ScaleOverwrite(1.0/sumState[i], CO[i], &gCO[i]);
}
}
// end re-estimate
printf("Iter = %d Loglik = %8.4f\n", iter, loglik);
if (fabs(oldlog - loglik) < tol) {
printf("\nConverged after %d iterations\n", iter);
break;
}
oldlog = loglik;
}
}
void hmm_trainG(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO, int max_iter, double tol) {
Matrix &gTR = *guessTR;
ArrayList<Vector>& gME = *guessME;
fastlib/u/tqlong/hmm/gaussianHMM.h
+8 -4
View File
@@ -3,6 +3,8 @@
#include "fastlib/fastlib.h"
success_t load_profileG(const char* profile, Matrix* trans, ArrayList<Vector>* means, ArrayList<Matrix>* covs);
success_t save_profileG(const char* profile, const Matrix& trans, const ArrayList<Vector>& means, const ArrayList<Matrix>& covs);
/**
Generating a sequence and states using transition and emission probabilities.
L: sequence length
@@ -35,6 +37,7 @@ void hmm_estimateG_init(int numStates, const Matrix& seq, const Vector& states,
*/
double hmm_decodeG(const Matrix& trans, const Matrix& emis_prob, Matrix* pstates, Matrix* fs, Matrix* bs, Vector* scales);
double hmm_decodeG(int L, const Matrix& trans, const Matrix& emis_prob, Matrix* pstates, Matrix* fs, Matrix* bs, Vector* scales);
void hmm_cal_emis_prob(const Matrix& seq, const ArrayList<Vector>& means, const ArrayList<Matrix>& inv_covs, const Vector& det, Matrix* emis_prob);
/** Calculate the most probable states for a sequence
Viterbi algorithm
@@ -45,13 +48,14 @@ double hmm_decodeG(int L, const Matrix& trans, const Matrix& emis_prob, Matrix*
RETURN: log probability of the most probable sequence
*/
double hmm_viterbiG_init(const Matrix& trans, const Matrix& emis_prob, Vector* states);
//double hmm_viterbiD_init(const Vector& seq, const Matrix& trans, const Matrix& emis, Vector* states);
/** Baum-Welch estimation of transition and emission distribution (Gaussian)
double hmm_viterbiG_init(int L, const Matrix& trans, const Matrix& emis_prob, Vector* states);
/** Baum-Welch and Viterbi estimation of transition and emission distribution (Gaussian)
*/
void init_gauss_param(int M, const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO);
void hmm_trainG(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO, int max_iter, double tol);
void hmm_train_viterbiG(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO, int max_iter, double tol);
#endif
fastlib/u/tqlong/hmm/main.cc
+24 -2
View File
@@ -829,8 +829,30 @@ int main(int argc, char* argv[]) {
//test_hmm_estimateM();
//test_hmm_decodeM();
//test_hmm_viterbiM();
test_hmm_trainM();
//test_hmm_trainM();
/*
{
ArrayList<Matrix> matlst;
load_matrix_list("seq.out", &matlst);
printf("NO SEQ = %d\n", matlst.size());
for (int i = 0; i < matlst.size(); i++) {
char s[100];
sprintf(s, "SEQ #%d:\n", i);
print_matrix(matlst[i],s);
}
}
*/
ArrayList<Vector> veclst;
load_vector_list("seq.out", &veclst);
printf("NO SEQ = %d\n", veclst.size());
for (int i = 0; i < veclst.size(); i++) {
char s[100];
sprintf(s, "SEQ #%d:\n", i);
print_vector(veclst[i],s);
}
/*
MixtureGauss mix;
mix.InitFromFile("meansG.txt", "covsG.txt");
fastlib/u/tqlong/hmm/mixgaussHMM.cc
+157
View File
@@ -3,6 +3,53 @@
#include "mixgaussHMM.h"
#include "gaussianHMM.h"
success_t load_profileM(const char* profile, Matrix* trans, ArrayList<MixtureGauss>* mixs) {
ArrayList<Matrix> matlst;
if (!PASSED(load_matrix_list(profile, &matlst))) {
NONFATAL("Couldn't open '%s' for reading.", profile);
return SUCCESS_FAIL;
}
DEBUG_ASSERT(matlst.size() >= 4); // at least 1 trans, 1 prior, 1 mean, 1 cov
trans->Copy(matlst[0]);
mixs->Init();
int M = trans->n_rows(); // num of states
int N = matlst[2].n_rows(); // dimension
int p = 1;
for (int i = 0; i < M; i++) {
int K = matlst[p].n_rows(); // num of clusters
//printf("load p=%d K=%d\n", p, K);
DEBUG_ASSERT(matlst.size() > p+2*K);
MixtureGauss mix;
mix.InitFromProfile(matlst, p, N);
mixs->AddBackItem(mix);
p += 2*K+1;
}
return SUCCESS_PASS;
}
success_t save_profileM(const char* profile, const Matrix& trans, const ArrayList<MixtureGauss>& mixs) {
TextWriter w_pro;
if (!PASSED(w_pro.Open(profile))) {
NONFATAL("Couldn't open '%s' for writing.", profile);
return SUCCESS_FAIL;
}
int M = trans.n_rows(); // num of states
print_matrix(w_pro, trans, "% transmission", "%E,");
for (int i = 0; i < M; i++) {
int K = mixs[i].n_clusters(); // num of clusters
char s[100];
sprintf(s, "%% prior - state %d", i);
print_vector(w_pro, mixs[i].get_prior(), s, "%E,");
for (int k=0; k < K; k++) {
sprintf(s, "%% mean %d - state %d", k, i);
print_vector(w_pro, mixs[i].get_mean(k), s, "%E,");
sprintf(s, "%% covariance %d - state %d", k, i);
print_matrix(w_pro, mixs[i].get_cov(k), s, "%E,");
}
}
return SUCCESS_PASS;
}
void hmm_generateM_init(int L, const Matrix& trans, const ArrayList<MixtureGauss>& mixs, Matrix* seq, Vector* states){
DEBUG_ASSERT_MSG((trans.n_rows()==trans.n_cols() && trans.n_rows()==mixs.size()), "hmm_generateM_init: matrices sizes do not match");
Matrix trsum;
@@ -252,6 +299,17 @@ double hmm_viterbiG_init(const Matrix& trans, const Matrix& emis_prob, Vector* s
}
*/
void hmm_cal_emis_probM(const Matrix& seq, const ArrayList<MixtureGauss>& mixs, Matrix* emis_prob) {
int M = mixs.size();
int L = seq.n_cols();
for (int t = 0; t < L; t++) {
Vector e;
seq.MakeColumnVector(t, &e);
for (int i = 0; i < M; i++)
emis_prob->ref(i, t) = mixs[i].getPDF(e);
}
}
void hmm_trainM(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<MixtureGauss>* guessMG, int max_iter, double tol) {
Matrix &gTR = *guessTR;
ArrayList<MixtureGauss>& gMG = *guessMG;
@@ -361,3 +419,102 @@ void hmm_trainM(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Mixtur
}
}
void hmm_train_viterbiM(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<MixtureGauss>* guessMG, int max_iter, double tol) {
Matrix &gTR = *guessTR;
ArrayList<MixtureGauss>& gMG = *guessMG;
int L = -1;
int M = gTR.n_rows();
DEBUG_ASSERT_MSG((M==gTR.n_cols() && M==gMG.size()),"hmm_trainM: sizes do not match");
for (int i = 0; i < seqs.size(); i++)
if (seqs[i].n_cols() > L) L = seqs[i].n_cols();
Matrix TR; // guess transition and emission matrix
TR.Init(M, M);
Matrix emis_prob; // to hold hmm_decodeG results
ArrayList<Matrix> emis_prob_cluster;
emis_prob.Init(M, L);
emis_prob_cluster.Init();
for (int i = 0; i < M; i++) {
Matrix m;
int K = gMG[i].n_clusters();
m.Init(K, L);
emis_prob_cluster.AddBackItem(m);
}
double loglik = 0, oldlog;
for (int iter = 0; iter < max_iter; iter++) {
oldlog = loglik;
loglik = 0;
// set the accumulating values to zeros and compute the inverse matrices and determinant constants
TR.SetZero();
for (int i = 0; i < M; i++)
gMG[i].start_accumulate();
// for each sequence, we will use viterbi procedure to find the most probable state sequence and then accumulate
for (int idx = 0; idx < seqs.size(); idx++) {
Vector states;
// first calculate the emission probabilities of the sequence
L = seqs[idx].n_cols();
for (int t = 0; t < L; t++) {
Vector e;
seqs[idx].MakeColumnVector(t, &e);
for (int i = 0; i < M; i++) {
double s = 0;
int K = gMG[i].n_clusters();
for (int j = 0; j < K; j++) {
emis_prob_cluster[i].ref(j, t) = gMG[i].getPDF(j, e);
s += emis_prob_cluster[i].ref(j, t);
}
emis_prob.ref(i, t) = s;
}
}
loglik += hmm_viterbiG_init(L, gTR, emis_prob, &states); // viterbi procedure
// accumulate expected transition & gaussian mixture parameters
for (int t = 0; t < L-1; t++) {
int i = (int) states[t];
int j = (int) states[t+1];
TR.ref(i, j)++;
}
for (int t = 0; t < L; t++) {
Vector e;
int i = (int) states[t];
seqs[idx].MakeColumnVector(t, &e);
int K = gMG[i].n_clusters();
for (int j = 0; j < K; j++)
gMG[i].accumulate(emis_prob_cluster[i].get(j, t) / emis_prob.get(i, t), j, e);
}
// end accumulate
}
// after accumulate all sequences: re-estimate transition & mean & covariance for the next iteration
for (int i = 0; i < M; i++) {
double s = 0;
for (int j = 0; j < M; j++) s += TR.get(i, j);
if (s == 0) {
for (int j = 0; j < M; j++) gTR.ref(i, j) = 0;
gTR.ref(i, i) = 1;
}
else {
for (int j = 0; j < M; j++) gTR.ref(i, j) = TR.get(i, j) / s;
}
gMG[i].end_accumulate();
}
// end re-estimate
printf("Iter = %d Loglik = %8.4f\n", iter, loglik);
if (fabs(oldlog - loglik) < tol) {
printf("\nConverged after %d iterations\n", iter);
break;
}
oldlog = loglik;
}
}
fastlib/u/tqlong/hmm/mixgaussHMM.h
+5 -3
View File
@@ -4,6 +4,9 @@
#include "fastlib/fastlib.h"
#include "mixtureDST.h"
success_t load_profileM(const char* profile, Matrix* trans, ArrayList<MixtureGauss>* mixs);
success_t save_profileM(const char* profile, const Matrix& trans, const ArrayList<MixtureGauss>& mixs);
/**
Generating a sequence and states using transition and emission probabilities.
L: sequence length
@@ -48,10 +51,9 @@ void hmm_estimateM_init(int numStates, int NumClusters, const Matrix& seq, const
//double hmm_viterbiD_init(const Vector& seq, const Matrix& trans, const Matrix& emis, Vector* states);
/** Baum-Welch estimation of transition and emission distribution (Gaussian)
*/
void hmm_cal_emis_probM(const Matrix& seq, const ArrayList<MixtureGauss>& mixs, Matrix* emis_prob);
void hmm_trainM(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<MixtureGauss>* guessMG, int max_iter=500, double tol=1e-3);
void hmm_train_viterbiM(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<MixtureGauss>* guessMG, int max_iter=500, double tol=1e-3);
#endif
fastlib/u/tqlong/hmm/mixtureDST.cc
+29
View File
@@ -113,6 +113,35 @@ void MixtureGauss::InitFromFile(const char* mean_fn, const char* covs_fn, const
}
}
void MixtureGauss::InitFromProfile(const ArrayList<Matrix>& matlst, int start, int N) {
DEBUG_ASSERT(matlst[start].n_cols()==1);
Vector tmp;
matlst[start].MakeColumnVector(0, &tmp);
prior.Copy(tmp);
means.Init();
covs.Init();
int K = prior.length();
for (int i = start+1; i < start+2*K+1; i+=2) {
DEBUG_ASSERT(matlst[i].n_rows()==N && matlst[i].n_cols()==1);
DEBUG_ASSERT(matlst[i+1].n_rows()==N && matlst[i+1].n_cols()==N);
Vector m;
matlst[i].MakeColumnVector(0, &m);
means.AddBackItem(m);
covs.AddBackItem(matlst[i+1]);
}
ACC_means.Copy(means);
ACC_covs.Copy(covs);
ACC_prior.Init(K);
inv_covs.Copy(covs);
det_covs.Init(covs.size());
for (int i = 0; i < K; i++) {
double det = la::Determinant(covs[i]);
la::InverseOverwrite(covs[i], &inv_covs[i]);
det_covs[i] = pow(2.0*math::PI, -N/2.0) * pow(det, -0.5);
}
}
void MixtureGauss::print_mixture(const char* s) const {
int K = means.size();
printf("%s - Mixture (%d)\n", s, K);
fastlib/u/tqlong/hmm/mixtureDST.h
+4
View File
@@ -16,12 +16,16 @@ class MixtureGauss {
double total;
public:
void InitFromFile(const char* mean_fn, const char* covs_fn = NULL, const char* prior_fn = NULL);
void InitFromProfile(const ArrayList<Matrix>& matlst, int start, int N);
void Init(int K, int N);
void Init(int K, const Matrix& data, const ArrayList<int>& labels);
void print_mixture(const char* s) const;
void generate(Vector* v) const;
double MixtureGauss::getPDF(const Vector& v) const;
double MixtureGauss::getPDF(int cluster, const Vector& v) const;
const Vector& get_prior() const { return prior; }
const Vector& get_mean(int k) const { return means[k]; }
const Matrix& get_cov(int k) const { return covs[k]; }
int n_clusters() const { return means.size(); }
int v_length() const { return means[0].length(); }
fastlib/u/tqlong/hmm/seq.out
+4 -1
View File
File diff suppressed because one or more lines are too long
fastlib/u/tqlong/hmm/support.cc
+164
View File
@@ -10,6 +10,15 @@ void print_matrix(const Matrix& a, const char* msg) {
}
}
void print_matrix(TextWriter& writer, const Matrix& a, const char* msg, const char* format) {
writer.Printf("%s - Matrix (%d x %d) = \n", msg, a.n_rows(), a.n_cols());
for (int j = 0; j < a.n_cols(); j++) {
for (int i = 0; i < a.n_rows(); i++)
writer.Printf(format, a.get(i, j));
writer.Printf("\n");
}
}
void print_vector(const Vector& a, const char* msg) {
printf("%s - Vector (%d) = \n", msg, a.length());
for (int i = 0; i < a.length(); i++)
@@ -17,6 +26,13 @@ void print_vector(const Vector& a, const char* msg) {
printf("\n");
}
void print_vector(TextWriter& writer, const Vector& a, const char* msg, const char* format) {
writer.Printf("%s - Vector (%d) = \n", msg, a.length());
for (int i = 0; i < a.length(); i++)
writer.Printf(format, a[i]);
writer.Printf("\n");
}
double RAND_NORMAL_01() {
double r = 2, u, v;
while (r > 1) {
@@ -258,3 +274,151 @@ void mat2arrlstmat(int N, Matrix& a, ArrayList<Matrix> * seqs) {
s_.AddBackItem(b);
}
}
bool skip_blank(TextLineReader& reader) {
for (;;){
if (!reader.MoreLines()) return false;
char* pos = reader.Peek().begin();
while (*pos == ' ' || *pos == ',' || *pos == '\t')
pos++;
if (*pos == '\0' || *pos == '%') reader.Gobble();
else break;
}
return true;
}
success_t read_matrix(TextLineReader& reader, Matrix* matrix) {
if (!skip_blank(reader)) { // EOF ?
matrix->Init(0,0);
return SUCCESS_FAIL;
}
else {
int n_rows = 0;
int n_cols = 0;
bool is_done;
{// How many columns ?
ArrayList<String> num_str;
num_str.Init();
reader.Peek().Split(", \t", &num_str);
n_cols = num_str.size();
}
ArrayList<double> num_double;
num_double.Init();
for(;;) { // read each rows
n_rows++;
double* point = num_double.AddBack(n_cols);
ArrayList<String> num_str;
num_str.Init();
reader.Peek().Split(", \t", &num_str);
DEBUG_ASSERT(num_str.size() == n_cols);
for (int i = 0; i < n_cols; i++)
*(point+i) = strtod(num_str[i], NULL);
is_done = false;
reader.Gobble();
for (;;){
if (!reader.MoreLines()) {
is_done = true;
break;
}
char* pos = reader.Peek().begin();
while (*pos == ' ' || *pos == '\t')
pos++;
if (*pos == '\0') reader.Gobble();
else if (*pos == '%') {
is_done = true;
break;
}
else break;
}
if (is_done) {
num_double.Trim();
matrix->Own(num_double.ReleasePointer(), n_cols, n_rows);
return SUCCESS_PASS;
}
}
}
}
success_t read_vector(TextLineReader& reader, Vector* vec) {
if (!skip_blank(reader)) { // EOF ?
vec->Init(0);
return SUCCESS_FAIL;
}
else {
ArrayList<double> num_double;
num_double.Init();
for(;;) { // read each rows
bool is_done = false;
ArrayList<String> num_str;
num_str.Init();
reader.Peek().Split(", \t", &num_str);
double* point = num_double.AddBack(num_str.size());
for (int i = 0; i < num_str.size(); i++)
*(point+i) = strtod(num_str[i], NULL);
reader.Gobble();
for (;;){
if (!reader.MoreLines()) {
is_done = true;
break;
}
char* pos = reader.Peek().begin();
while (*pos == ' ' || *pos == '\t')
pos++;
if (*pos == '\0') reader.Gobble();
else if (*pos == '%') {
is_done = true;
break;
}
else break;
}
if (is_done) {
num_double.Trim();
int length = num_double.size();
vec->Own(num_double.ReleasePointer(), length);
return SUCCESS_PASS;
}
}
}
}
success_t load_matrix_list(const char* filename, ArrayList<Matrix> *matlst) {
TextLineReader reader;
matlst->Init();
if (!PASSED(reader.Open(filename))) return SUCCESS_FAIL;
do {
Matrix tmp;
if (read_matrix(reader, &tmp) == SUCCESS_PASS) {
matlst->AddBackItem(tmp);
}
else break;
} while (1);
return SUCCESS_PASS;
}
success_t load_vector_list(const char* filename, ArrayList<Vector> *veclst) {
TextLineReader reader;
veclst->Init();
if (!PASSED(reader.Open(filename))) return SUCCESS_FAIL;
do {
Vector vec;
if (read_vector(reader, &vec) == SUCCESS_PASS) {
veclst->AddBackItem(vec);
}
else break;
} while (1);
return SUCCESS_PASS;
}
fastlib/u/tqlong/hmm/support.h
+5
View File
@@ -14,8 +14,10 @@ double MyMulExpert(const Vector& x, const Matrix& A, const Vector& y);
double NORMAL_DENSITY(const Vector& x, const Vector& mean, const Matrix& inv_cov, double det_cov);
void print_matrix(const Matrix& a, const char* msg);
void print_matrix(TextWriter& writer, const Matrix& a, const char* msg, const char* format = "%f,");
void print_vector(const Vector& a, const char* msg);
void print_vector(TextWriter& writer, const Vector& a, const char* msg, const char* format = "%f,");
bool kmeans(const ArrayList<Matrix>& data, int num_clusters,
ArrayList<int> *labels_, ArrayList<Vector> *cetroids_,
@@ -29,5 +31,8 @@ void mat2arrlst(Matrix& a, ArrayList<Vector> * seqs);
void mat2arrlstmat(int N, Matrix& a, ArrayList<Matrix> * seqs);
success_t load_matrix_list(const char* filename, ArrayList<Matrix> *matlst);
success_t load_vector_list(const char* filename, ArrayList<Vector> *veclst);
#endif