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-22 21:47:49 +00:00
parent 5cc07d1a02
commit bb1bc193db
20 changed files with 1363 additions and 288 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/u/tqlong/hmm/#mixtureDST.cc#
-214
View File
@@ -1,214 +0,0 @@
#include "fastlib/fastlib.h"
#include "support.h"
#include "mixtureDST.h"
void MixtureGauss::Init(int K, int N) {
means.Init();
for (int i = 0; i < K; i++) {
Vector v;
RAND_NORMAL_01_INIT(N, &v);
means.AddBackItem(v);
}
covs.Init();
for (int i = 0; i < means.size(); i++) {
Matrix m;
m.Init(N, N); m.SetZero();
for (int j = 0; j < N; j++) m.ref(j, j) = 1.0;
covs.AddBackItem(m);
}
prior.Init(means.size());
for (int i = 0; i < prior.length(); i++) prior[i] = 1.0/K;
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::Init(int K, const Matrix& data, const ArrayList<int>& labels) {
means.Init();
int N = data.n_rows();
for (int i = 0; i < K; i++) {
Vector v;
v.Init(N);
means.AddBackItem(v);
}
covs.Init();
for (int i = 0; i < means.size(); i++) {
Matrix m;
m.Init(N, N);
covs.AddBackItem(m);
}
prior.Init(means.size());
ACC_means.Copy(means);
ACC_covs.Copy(covs);
ACC_prior.Init(K);
inv_covs.Copy(covs);
det_covs.Init(covs.size());
start_accumulate();
//printf("cols = %d rows = %d\n", data.n_cols(), data.n_rows());
for (int i = 0; i < data.n_cols(); i++) {
Vector v;
data.MakeColumnVector(i, &v);
//printf("%d\n", i);
accumulate_cluster(labels[i], v);
}
end_accumulate_cluster();
}
void MixtureGauss::InitFromFile(const char* mean_fn, const char* covs_fn, const char* prior_fn) {
Matrix meansmat;
data::Load(mean_fn, &meansmat);
mat2arrlst(meansmat, &means);
int N = means[0].length();
int K = means.size();
if (covs_fn != NULL) {
Matrix covsmat;
data::Load(covs_fn, &covsmat);
mat2arrlstmat(N, covsmat, &covs);
DEBUG_ASSERT_MSG(K==covs.size(), "InitFromFile: sizes do not match !");
}
else {
covs.Init();
for (int i = 0; i < means.size(); i++) {
Matrix m;
m.Init(N, N); m.SetZero();
for (int j = 0; j < N; j++) m.ref(j, j) = 1.0;
covs.AddBackItem(m);
}
}
if (prior_fn != NULL) {
Matrix priormat;
data::Load(prior_fn, &priormat);
DEBUG_ASSERT_MSG(K==priormat.n_cols(), "InitFromFile: sizes do not match !!");
prior.Init(K);
for (int i = 0; i < K; i++) prior[i] = priormat.get(0, i);
}
else {
prior.Init(means.size());
for (int i = 0; i < prior.length(); i++) prior[i] = 1.0/K;
}
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 Mixture::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);
int K = prior.length();
for (int i = start+1; i < start+2*K+1; i+=2) {
DEBUG_ASSERT(matlst[i].n_rows()==N && matlst[starti].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]);
}
Matrix meansmat;
data::Load(mean_fn, &meansmat);
mat2arrlst(meansmat, &means);
int N = means[0].length();
int K = means.size();
if (covs_fn != NULL) {
Matrix covsmat;
data::Load(covs_fn, &covsmat);
mat2arrlstmat(N, covsmat, &covs);
DEBUG_ASSERT_MSG(K==covs.size(), "InitFromFile: sizes do not match !");
}
else {
covs.Init();
for (int i = 0; i < means.size(); i++) {
Matrix m;
m.Init(N, N); m.SetZero();
for (int j = 0; j < N; j++) m.ref(j, j) = 1.0;
covs.AddBackItem(m);
}
}
if (prior_fn != NULL) {
Matrix priormat;
data::Load(prior_fn, &priormat);
DEBUG_ASSERT_MSG(K==priormat.n_cols(), "InitFromFile: sizes do not match !!");
prior.Init(K);
for (int i = 0; i < K; i++) prior[i] = priormat.get(0, i);
}
else {
prior.Init(means.size());
for (int i = 0; i < prior.length(); i++) prior[i] = 1.0/K;
}
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);
print_vector(prior, " PRIOR");
for (int i = 0; i < K; i++) {
printf(" CLUSTER %d:\n", i);
print_vector(means[i], " MEANS");
print_matrix(covs[i], " COVS");
}
}
void MixtureGauss::generate(Vector* v) const {
int K = means.size();
double r = RAND_UNIFORM_01;
int cluster = K-1;
double s = 0;
for (int i = 0; i < K; i++) {
s += prior[i];
if (s >= r) {
cluster = i;
break;
}
}
RAND_NORMAL_INIT(means[cluster], covs[cluster], v);
}
double MixtureGauss::getPDF(const Vector& v) const {
int K = means.size();
double s = 0;
for (int i = 0; i < K; i++)
s += getPDF(i, v);
return s;
}
double MixtureGauss::getPDF(int cluster, const Vector& v) const {
return prior[cluster]*NORMAL_DENSITY(v, means[cluster], inv_covs[cluster], det_covs[cluster]);
}
fastlib/u/tqlong/hmm/build.py
+26 -3
View File
@@ -1,10 +1,33 @@
binrule(
librule(
name = "hmm", # the executable name
sources = ["hmm.cc","support.cc","discreteHMM.cc","gaussianHMM.cc","mixgaussHMM.cc","mixtureDST.cc"],
sources = ["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
deplibs = ["fastlib:fastlib_int"] # depends on faslib core library
)
binrule(
name = "generate", # the executable name
sources = ["generate.cc"],
deplibs = [":hmm"] # depends on hmm library in this folder
)
binrule(
name = "loglik", # the executable name
sources = ["loglik.cc"],
deplibs = [":hmm"] # depends on hmm library in this folder
)
binrule(
name = "viterbi", # the executable name
sources = ["viterbi.cc"],
deplibs = [":hmm"] # depends on hmm library in this folder
)
binrule(
name = "train", # the executable name
sources = ["train.cc"],
deplibs = [":hmm"] # depends on hmm library in this folder
)
# to build:
# 1. make sure have environment variables set up:
fastlib/u/tqlong/hmm/discreteHMM.cc
+138 -2
View File
@@ -2,6 +2,143 @@
#include "support.h"
#include "discreteHMM.h"
void DiscreteHMM::setModel(const Matrix& transmission, const Matrix& emission) {
DEBUG_ASSERT(transmission.n_rows() == transmission.n_cols());
DEBUG_ASSERT(transmission.n_rows() == emission.n_rows());
transmission_.Destruct();
emission_.Destruct();
transmission_.Copy(transmission);
emission_.Copy(emission);
}
void DiscreteHMM::Init(const Matrix& transmission, const Matrix& emission) {
transmission_.Copy(transmission);
emission_.Copy(emission);
DEBUG_ASSERT(transmission.n_rows() == transmission.n_cols());
DEBUG_ASSERT(transmission.n_rows() == emission.n_rows());
}
void DiscreteHMM::InitFromFile(const char* profile) {
ArrayList<Matrix> list_mat;
load_matrix_list(profile, &list_mat);
if (list_mat.size() < 2)
FATAL("Number of matrices in the file should be at least 2.");
else if (list_mat.size() > 2)
NONFATAL("Number of matrices in the file should be 2.");
transmission_.Own(&list_mat[0]);
emission_.Own(&list_mat[1]);
DEBUG_ASSERT(transmission_.n_rows() == transmission_.n_cols());
DEBUG_ASSERT(transmission_.n_rows() == emission_.n_rows());
}
void DiscreteHMM::InitFromData(const ArrayList<Vector>& list_data_seq, int numstate) {
int numsymbol = 0;
int maxseq = 0;
for (int i = 0; i < list_data_seq.size(); i++)
if (list_data_seq[i].length() > list_data_seq[maxseq].length()) maxseq = i;
for (int i = 0; i < list_data_seq[maxseq].length(); i++)
if (list_data_seq[maxseq][i] > numsymbol) numsymbol = (int) list_data_seq[maxseq][i];
numsymbol++;
Vector states;
int L = list_data_seq[maxseq].length();
states.Init(L);
for (int i = 0; i < L; i++) states[i] = rand() % numstate;
hmm_estimateD_init(numsymbol, numstate, list_data_seq[maxseq], states, &transmission_, &emission_);
}
void DiscreteHMM::LoadProfile(const char* profile) {
transmission_.Destruct();
emission_.Destruct();
InitFromFile(profile);
}
void DiscreteHMM::SaveProfile(const char* profile) const {
TextWriter w_pro;
if (!PASSED(w_pro.Open(profile))) {
NONFATAL("Couldn't open '%s' for writing.", profile);
return;
}
print_matrix(w_pro, transmission_, "%% transmision", "%f,");
print_matrix(w_pro, emission_, "%% emission", "%f,");
}
void DiscreteHMM::GenerateSequence(int length, Vector* data_seq, Vector* state_seq) const {
hmm_generateD_init(length, transmission_, emission_, data_seq, state_seq);
}
void DiscreteHMM::EstimateModel(const Vector& data_seq, const Vector& state_seq) {
transmission_.Destruct();
emission_.Destruct();
hmm_estimateD_init(data_seq, state_seq, &transmission_, &emission_);
}
void DiscreteHMM::EstimateModel(int numstate, int numsymbol, const Vector& data_seq, const Vector& state_seq) {
transmission_.Destruct();
emission_.Destruct();
hmm_estimateD_init(numsymbol, numstate, data_seq, state_seq, &transmission_, &emission_);
}
void DiscreteHMM::DecodeOverwrite(const Vector& data_seq, Matrix* state_prob_mat, Matrix* forward_prob_mat, Matrix* backward_prob_mat, Vector* scale_vec) const {
hmm_decodeD(data_seq, transmission_, emission_, state_prob_mat, forward_prob_mat, backward_prob_mat, scale_vec);
}
void DiscreteHMM::DecodeInit(const Vector& data_seq, Matrix* state_prob_mat, Matrix* forward_prob_mat, Matrix* backward_prob_mat, Vector* scale_vec) const {
int M = transmission_.n_rows();
int L = data_seq.length();
state_prob_mat->Init(M, L);
forward_prob_mat->Init(M, L);
backward_prob_mat->Init(M, L);
scale_vec->Init(L);
hmm_decodeD(data_seq, transmission_, emission_, state_prob_mat, forward_prob_mat, backward_prob_mat, scale_vec);
}
void forward_procedure(const Vector& seq, const Matrix& trans, const Matrix& emis, Vector *scales, Matrix* fs);
double DiscreteHMM::ComputeLogLikelihood(const Vector& data_seq) const {
int L = data_seq.length();
int M = transmission_.n_rows();
Matrix fs(M, L);
Vector sc;
sc.Init(L);
forward_procedure(data_seq, transmission_, emission_, &sc, &fs);
double loglik = 0;
for (int t = 0; t < L; t++)
loglik += log(sc[t]);
return loglik;
}
void DiscreteHMM::ComputeLogLikelihood(const ArrayList<Vector>& list_data_seq, ArrayList<double>* list_likelihood) const {
int L = 0;
for (int i = 0; i < list_data_seq.size(); i++)
if (list_data_seq[i].length() > L) L = list_data_seq[i].length();
int M = transmission_.n_rows();
Matrix fs(M, L);
Vector sc;
sc.Init(L);
list_likelihood->Init();
for (int i = 0; i < list_data_seq.size(); i++) {
forward_procedure(list_data_seq[i], transmission_, emission_, &sc, &fs);
int L = list_data_seq[i].length();
double loglik = 0;
for (int t = 0; t < L; t++)
loglik += log(sc[t]);
list_likelihood->AddBackItem(loglik);
}
}
void DiscreteHMM::ComputeViterbiStateSequence(const Vector& data_seq, Vector* state_seq) const {
hmm_viterbiD_init(data_seq, transmission_, emission_, state_seq);
}
void DiscreteHMM::TrainBaumWelch(const ArrayList<Vector>& list_data_seq, int max_iteration, double tolerance) {
hmm_trainD(list_data_seq, &transmission_, &emission_, max_iteration, tolerance);
}
void DiscreteHMM::TrainViterbi(const ArrayList<Vector>& list_data_seq, int max_iteration, double tolerance) {
hmm_train_viterbiD(list_data_seq, &transmission_, &emission_, max_iteration, tolerance);
}
void hmm_generateD_init(int L, const Matrix& trans, const Matrix& emis, Vector* seq, Vector* states) {
DEBUG_ASSERT_MSG((trans.n_rows()==trans.n_cols() && trans.n_rows()==emis.n_rows()), "hmm_generateD_init: matrices sizes do not match");
Matrix trsum, esum;
@@ -110,6 +247,7 @@ void forward_procedure(const Vector& seq, const Matrix& trans, const Matrix& emi
Vector& s_ = *scales;
fs_.SetZero();
s_.SetZero();
// NOTE: start state is 0
// time t = 0
int e = (int) seq[0];
@@ -161,9 +299,7 @@ double hmm_decodeD(const Vector& seq, const Matrix& trans, const Matrix& emis, M
M==trans.n_cols() && M==emis.n_rows()),"hmm_decodeD: sizes do not match");
Matrix& ps_ = *pstates;
Vector& s_ = *scales;
s_.SetZero();
forward_procedure(seq, trans, emis, &s_, fs);
backward_procedure(seq, trans, emis, s_, bs);
fastlib/u/tqlong/hmm/discreteHMM.h
+92 -1
View File
@@ -1,8 +1,100 @@
/**
* @file discreteHMM.h
*
* This file contains functions of a discrete Hidden Markov Models (in C) and a
* wrapper class (in C++) for these functions. It implements log-likelihood
* computation, viterbi algorithm for the most probable sequence, Baum-Welch
* algorithm and Viterbi-like algorithm for parameter estimation. It can also
* generate sequences from a Hidden Markov Model.
*/
#ifndef FASTLIB_DISCRETE_HMM_H
#define FASTLIB_DISCRETE_HMM_H
#include "fastlib/fastlib.h"
/**
* A wrapper class for HMM functionals in discrete case
*
* This class maintains transition probabilities and emission probabilities
* matrices and performs basic HMM functionals, more details below.
*
*/
class DiscreteHMM {
/////////// Member variables /////////////////////////////////////
private:
/** Transmission probabilities matrix between states */
Matrix transmission_;
/** Emission probabilities in each state */
Matrix emission_;
public:
/** Basic getters */
const Matrix& getTransmission() const { return transmission_; }
const Matrix& getEmission() const { return emission_; }
/** Setters used when already initialized */
void setModel(const Matrix& transmission, const Matrix& emission);
/** Initializes from computed transmission and emission matrices */
void Init(const Matrix& transmission, const Matrix& emission);
/** Initializes by loading from a file */
void InitFromFile(const char* profile);
/** Initializes randomly using data as a guide */
void InitFromData(const ArrayList<Vector>& list_data_seq, int numstate);
/** Load from file, used when already initialized */
void LoadProfile(const char* profile);
/** Save matrices to file */
void SaveProfile(const char* profile) const;
/** Generate a random data sequence of a given length */
void GenerateSequence(int length, Vector* data_seq, Vector* state_seq) const;
/**
* Estimate the matrices by a data sequence and a state sequence
* Must be already initialized
*/
void EstimateModel(const Vector& data_seq, const Vector& state_seq);
void EstimateModel(int numstate, int numsymbol, const Vector& data_seq, const Vector& state_seq);
/**
* Decode a sequence into probabilities of each state at each time step
* using scaled forward-backward algorithm.
* Also return forward, backward probabilities and scale factors
*/
void DecodeOverwrite(const Vector& data_seq, Matrix* state_prob_mat, Matrix* forward_prob_mat, Matrix* backward_prob_mat, Vector* scale_vec) const;
/** A decode version that initialized the out matrices */
void DecodeInit(const Vector& data_seq, Matrix* state_prob_mat, Matrix* forward_prob_mat, Matrix* backward_prob_mat, Vector* scale_vec) const;
/** Compute the log-likelihood of a sequence */
double ComputeLogLikelihood(const Vector& data_seq) const;
/** Compute the log-likelihood of a list of sequences */
void ComputeLogLikelihood(const ArrayList<Vector>& list_data_seq, ArrayList<double>* list_likelihood) const;
/** Compute the most probable sequence (Viterbi) */
void ComputeViterbiStateSequence(const Vector& data_seq, Vector* state_seq) const;
/**
* Train the model with a list of sequences, must be already initialized
* using Baum-Welch EM algorithm
*/
void TrainBaumWelch(const ArrayList<Vector>& list_data_seq, int max_iteration = 500, double tolerance = 1e-3);
/**
* Train the model with a list of sequences, must be already initialized
* using Viterbi algorithm to determine the state sequence of each sequence
*/
void TrainViterbi(const ArrayList<Vector>& list_data_seq, int max_iteration = 500, double tolerance = 1e-3);
};
/**
Generating a sequence and states using transition and emission probabilities.
L: sequence length
@@ -13,7 +105,6 @@ states: uninitialized vector, will have length L
*/
void hmm_generateD_init(int L, const Matrix& trans, const Matrix& emis, Vector* seq, Vector * states);
/** Estimate transition and emission probabilities from sequence and states
*/
void hmm_estimateD_init(const Vector& seq, const Vector& states, Matrix* trans, Matrix* emis);
fastlib/u/tqlong/hmm/gaussianHMM.cc
+107 -13
View File
@@ -2,12 +2,104 @@
#include "support.h"
#include "gaussianHMM.h"
void GaussianHMM::InitFromFile(const char* profile) {
if (!PASSED(load_profileG(profile, &transmission_, &list_mean_vec_, &list_covariance_mat_)))
FATAL("Couldn't open '%s' for reading.", profile);
list_inverse_cov_mat_.Copy(list_covariance_mat_);
gauss_const_vec_.Init(list_covariance_mat_.size());
CalculateInverse();
}
void GaussianHMM::InitFromData(const ArrayList<Matrix>& list_data_seq, int numstate) {
init_gauss_param(numstate, list_data_seq, &transmission_, &list_mean_vec_, &list_covariance_mat_);
list_inverse_cov_mat_.Copy(list_covariance_mat_);
gauss_const_vec_.Init(list_covariance_mat_.size());
CalculateInverse();
}
void GaussianHMM::LoadProfile(const char* profile) {
transmission_.Destruct();
list_mean_vec_.Destruct();
list_covariance_mat_.Destruct();
list_inverse_cov_mat_.Destruct();
gauss_const_vec_.Destruct();
InitFromFile(profile);
}
void GaussianHMM::SaveProfile(const char* profile) {
save_profileG(profile, transmission_, list_mean_vec_, list_covariance_mat_);
}
void GaussianHMM::CalculateInverse() {
int M = transmission_.n_rows();
int N = list_mean_vec_[0].length();
for (int i = 0; i < M; i++) {
la::InverseOverwrite(list_covariance_mat_[i], &list_inverse_cov_mat_[i]);
gauss_const_vec_[i] = pow(2.0*math::PI, -N/2.0) * pow(la::Determinant(list_covariance_mat_[i]), -0.5);
}
}
void GaussianHMM::GenerateSequence(int L, Matrix* data_seq, Vector* state_seq) {
hmm_generateG_init(L, transmission_, list_mean_vec_, list_covariance_mat_, data_seq, state_seq);
}
void forward_procedure(int L, const Matrix& trans, const Matrix& emis_prob, Vector *scales, Matrix* fs);
double GaussianHMM::ComputeLogLikelihood(const Matrix& data_seq) {
int L = data_seq.n_cols();
int M = transmission_.n_rows();
Matrix fs(M, L), emis_prob(M, L);
Vector sc;
sc.Init(L);
hmm_cal_emis_prob(data_seq, list_mean_vec_, list_inverse_cov_mat_, gauss_const_vec_, &emis_prob);
forward_procedure(L, transmission_, emis_prob, &sc, &fs);
double loglik = 0;
for (int t = 0; t < L; t++)
loglik += log(sc[t]);
return loglik;
}
void GaussianHMM::ComputeLogLikelihood(const ArrayList<Matrix>& list_data_seq, ArrayList<double>* list_likelihood) {
int L = 0;
for (int i = 0; i < list_data_seq.size(); i++)
if (list_data_seq[i].n_cols() > L) L = list_data_seq[i].n_cols();
int M = transmission_.n_rows();
Matrix fs(M, L), emis_prob(M, L);
Vector sc;
sc.Init(L);
list_likelihood->Init();
for (int i = 0; i < list_data_seq.size(); i++) {
int L = list_data_seq[i].n_cols();
hmm_cal_emis_prob(list_data_seq[i], list_mean_vec_, list_inverse_cov_mat_, gauss_const_vec_, &emis_prob);
forward_procedure(L, transmission_, emis_prob, &sc, &fs);
double loglik = 0;
for (int t = 0; t < L; t++)
loglik += log(sc[t]);
list_likelihood->AddBackItem(loglik);
}
}
void GaussianHMM::ComputeViterbiStateSequence(const Matrix& data_seq, Vector* state_seq) {
int M = transmission_.n_rows();
int L = data_seq.n_cols();
Matrix emis_prob(M, L);
hmm_cal_emis_prob(data_seq, list_mean_vec_, list_inverse_cov_mat_, gauss_const_vec_, &emis_prob);
hmm_viterbiG_init(transmission_, emis_prob, state_seq);
}
void GaussianHMM::TrainBaumWelch(const ArrayList<Matrix>& list_data_seq, int max_iteration, double tolerance) {
hmm_trainG(list_data_seq, &transmission_, &list_mean_vec_, &list_covariance_mat_, max_iteration, tolerance);
}
void GaussianHMM::TrainViterbi(const ArrayList<Matrix>& list_data_seq, int max_iteration, double tolerance) {
hmm_train_viterbiG(list_data_seq, &transmission_, &list_mean_vec_, &list_covariance_mat_, max_iteration, tolerance);
}
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);
if (!PASSED(load_matrix_list(profile, &matlst)))
return SUCCESS_FAIL;
}
DEBUG_ASSERT(matlst.size() > 0);
trans->Copy(matlst[0]);
means->Init();
@@ -35,15 +127,15 @@ success_t save_profileG(const char* profile, const Matrix& trans, const ArrayLis
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,");
print_matrix(w_pro, trans, "% transmission", "%f,");
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,");
print_vector(w_pro, means[i], s, "%f,");
sprintf(s, "%% covariance - state%d", i);
print_matrix(w_pro, covs[i], s, "%E,");
print_matrix(w_pro, covs[i], s, "%f,");
}
return SUCCESS_PASS;
}
@@ -175,6 +267,7 @@ void forward_procedure(int L, const Matrix& trans, const Matrix& emis_prob, Vect
Vector& s_ = *scales;
fs_.SetZero();
s_.SetZero();
// NOTE: start state is 0
// time t = 0
for (int i = 0; i < M; i++) {
@@ -221,9 +314,7 @@ double hmm_decodeG(int L, const Matrix& trans, const Matrix& emis_prob, Matrix*
M==trans.n_cols() && M==emis_prob.n_rows()),"hmm_decodeG: sizes do not match");
Matrix& ps_ = *pstates;
Vector& s_ = *scales;
s_.SetZero();
forward_procedure(L, trans, emis_prob, &s_, fs);
backward_procedure(L, trans, emis_prob, s_, bs);
@@ -359,7 +450,8 @@ void init_gauss_param(int M, const ArrayList<Matrix>& seqs, Matrix* guessTR, Arr
la::MulExpert(1, false, tmp_cov, true, tmp_cov, 1, &gCO[i]);
}
}
//printf("---3---\n");
for (int i = 0; i < M; i++)
if (sumState[i] == 0) {
for (int j = 0; j < M; j++) gTR.ref(i, j) = 0;
@@ -373,6 +465,7 @@ void init_gauss_param(int M, const ArrayList<Matrix>& seqs, Matrix* guessTR, Arr
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
}
//printf("---4---\n");
}
void hmm_train_viterbiG(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO, int max_iter, double tol) {
@@ -554,10 +647,8 @@ void hmm_trainG(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Vector
sumState[i] += ps.get(i, t);
la::AddExpert(ps.get(i, t), e, &ME[i]);
Vector d;
la::SubInit(e, gME[i], &d);
Matrix D;
D.AliasColVector(d);
D.AliasColVector(e);
la::MulExpert(ps.get(i, t), false, D, true, D, 1.0, &CO[i]);
}
}
@@ -578,7 +669,10 @@ void hmm_trainG(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Vector
if (sumState[i] != 0) {
la::ScaleOverwrite(1.0/sumState[i], ME[i], &gME[i]);
la::ScaleOverwrite(1.0/sumState[i], CO[i], &gCO[i]);
Matrix D;
D.AliasColVector(gME[i]);
la::MulExpert(-1.0, false, D, true, D, 1.0/sumState[i], &CO[i]);
gCO[i].CopyValues(CO[i]);
}
}
// end re-estimate
fastlib/u/tqlong/hmm/gaussianHMM.h
+31
View File
@@ -3,6 +3,37 @@
#include "fastlib/fastlib.h"
/**
* A wrapper class for HMM functionals in single Gaussian case
*/
class GaussianHMM {
private:
Matrix transmission_;
ArrayList<Vector> list_mean_vec_;
ArrayList<Matrix> list_covariance_mat_;
ArrayList<Matrix> list_inverse_cov_mat_;
Vector gauss_const_vec_;
void CalculateInverse();
public:
void InitFromFile(const char* profile);
void InitFromData(const ArrayList<Matrix>& list_data_seq, int numstate);
void LoadProfile(const char* profile);
void SaveProfile(const char* profile);
void GenerateSequence(int L, Matrix* data_seq, Vector* state_seq);
double ComputeLogLikelihood(const Matrix& data_seq);
void ComputeLogLikelihood(const ArrayList<Matrix>& list_data_seq, ArrayList<double>* list_likelihood);
void ComputeViterbiStateSequence(const Matrix& data_seq, Vector* state_seq);
void TrainBaumWelch(const ArrayList<Matrix>& list_data_seq, int max_iteration = 500, double tolerance = 1e-3);
void TrainViterbi(const ArrayList<Matrix>& list_data_seq, int max_iteration = 500, double tolerance = 1e-3);
};
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);
/**
fastlib/u/tqlong/hmm/generate.cc
+188
View File
@@ -0,0 +1,188 @@
#include "fastlib/fastlib.h"
#include "support.h"
#include "discreteHMM.h"
#include "gaussianHMM.h"
#include "mixgaussHMM.h"
#include "mixtureDST.h"
success_t generate_discrete();
success_t generate_gaussian();
success_t generate_mixture();
void usage();
int main(int argc, char* argv[]) {
fx_init(argc, argv);
success_t s = SUCCESS_PASS;
if (fx_param_exists(NULL,"type")) {
const char* type = fx_param_str_req(NULL, "type");
if (strcmp(type, "discrete")==0)
s = generate_discrete();
else if (strcmp(type, "gaussian")==0)
s = generate_gaussian();
else if (strcmp(type, "mixture")==0)
s = generate_mixture();
else {
printf("Unrecognized type: must be: discrete | gaussian | mixture !!!\n");
return SUCCESS_PASS;
}
}
else {
printf("Unrecognized type: must be: discrete | gaussian | mixture !!!\n");
s = SUCCESS_FAIL;
}
if (!PASSED(s)) usage();
fx_done();
}
void usage() {
printf("\nUsage:\n");
printf(" generate --type=={discrete|gaussian|mixture} OPTIONS\n");
printf("[OPTIONS]\n");
printf(" --profile=file : file contains HMM profile\n");
printf(" --length=NUM : sequence length\n");
printf(" --lenmax=NUM : maximum sequence length, default = length\n");
printf(" --numseq=NUM : number of sequence\n");
printf(" --seqfile=file : output file for generated sequences\n");
printf(" --statefile=file : output file for generated state sequences\n");
}
success_t generate_mixture() {
if (!fx_param_exists(NULL, "profile")) {
printf("--profile must be defined.\n");
return SUCCESS_FAIL;
}
const char* profile = fx_param_str_req(NULL, "profile");
const int seqlen = fx_param_int(NULL, "length", 10);
const int seqlmax = fx_param_int(NULL, "lenmax", seqlen);
const int numseq = fx_param_int(NULL, "numseq", 10);
const char* seqout = fx_param_str(NULL, "seqfile", "seq.mix.out");
const char* stateout = fx_param_str(NULL, "statefile", "state.mix.out");
DEBUG_ASSERT_MSG(seqlen <= seqlmax, "LENMAX must bigger than LENGTH");
DEBUG_ASSERT_MSG(numseq > 0, "NUMSEQ must be positive");
double step = (double) (seqlmax-seqlen) / numseq;
MixtureofGaussianHMM hmm;
hmm.InitFromFile(profile);
TextWriter w_seq, w_state;
if (!PASSED(w_seq.Open(seqout))) {
NONFATAL("Couldn't open '%s' for writing.", seqout);
return SUCCESS_FAIL;
}
if (!PASSED(w_state.Open(stateout))) {
NONFATAL("Couldn't open '%s' for writing.", stateout);
return SUCCESS_FAIL;
}
double L = seqlen;
for (int i = 0; i < numseq; i++, L+=step) {
Matrix seq;
Vector states;
char s[100];
hmm.GenerateSequence((int)L, &seq, &states);
sprintf(s, "%% sequence %d", i);
print_matrix(w_seq, seq, s, "%E,");
sprintf(s, "%% state sequence %d", i);
print_vector(w_state, states, s, "%.0f,");
}
//printf("---END---");
return SUCCESS_PASS;
}
success_t generate_gaussian() {
if (!fx_param_exists(NULL, "profile")) {
printf("--profile must be defined.\n");
return SUCCESS_FAIL;
}
const char* profile = fx_param_str_req(NULL, "profile");
const int seqlen = fx_param_int(NULL, "length", 10);
const int seqlmax = fx_param_int(NULL, "lenmax", seqlen);
const int numseq = fx_param_int(NULL, "numseq", 10);
const char* seqout = fx_param_str(NULL, "seqfile", "seq.gauss.out");
const char* stateout = fx_param_str(NULL, "statefile", "state.gauss.out");
DEBUG_ASSERT_MSG(seqlen <= seqlmax, "LENMAX must bigger than LENGTH");
DEBUG_ASSERT_MSG(numseq > 0, "NUMSEQ must be positive");
double step = (double) (seqlmax-seqlen) / numseq;
GaussianHMM hmm;
hmm.InitFromFile(profile);
TextWriter w_seq, w_state;
if (!PASSED(w_seq.Open(seqout))) {
NONFATAL("Couldn't open '%s' for writing.", seqout);
return SUCCESS_FAIL;
}
if (!PASSED(w_state.Open(stateout))) {
NONFATAL("Couldn't open '%s' for writing.", stateout);
return SUCCESS_FAIL;
}
double L = seqlen;
for (int i = 0; i < numseq; i++, L+=step) {
Matrix seq;
Vector states;
char s[100];
hmm.GenerateSequence((int)L, &seq, &states);
sprintf(s, "%% sequence %d", i);
print_matrix(w_seq, seq, s, "%E,");
sprintf(s, "%% state sequence %d", i);
print_vector(w_state, states, s, "%.0f,");
}
return SUCCESS_PASS;
}
success_t generate_discrete() {
if (!fx_param_exists(NULL, "profile")) {
printf("--profile must be defined.\n");
return SUCCESS_FAIL;
}
const char* profile = fx_param_str_req(NULL, "profile");
const int seqlen = fx_param_int(NULL, "length", 10);
const int seqlmax = fx_param_int(NULL, "lenmax", seqlen);
const int numseq = fx_param_int(NULL, "numseq", 10);
const char* seqout = fx_param_str(NULL, "seqfile", "seq.out");
const char* stateout = fx_param_str(NULL, "statefile", "state.out");
DEBUG_ASSERT_MSG(seqlen <= seqlmax, "LENMAX must bigger than LENGTH");
DEBUG_ASSERT_MSG(numseq > 0, "NUMSEQ must be positive");
double step = (double) (seqlmax-seqlen) / numseq;
DiscreteHMM hmm;
hmm.InitFromFile(profile);
TextWriter w_seq, w_state;
if (!PASSED(w_seq.Open(seqout))) {
NONFATAL("Couldn't open '%s' for writing.", seqout);
return SUCCESS_FAIL;
}
if (!PASSED(w_state.Open(stateout))) {
NONFATAL("Couldn't open '%s' for writing.", stateout);
return SUCCESS_FAIL;
}
double L = seqlen;
for (int i = 0; i < numseq; i++, L+=step) {
Vector seq, states;
char s[100];
hmm.GenerateSequence((int)L, &seq, &states);
sprintf(s, "%% sequence %d", i);
print_vector(w_seq, seq, s, "%.0f,");
sprintf(s, "%% state sequence %d", i);
print_vector(w_state, states, s, "%.0f,");
}
return SUCCESS_PASS;
}
fastlib/u/tqlong/hmm/log.out
+10 -4
View File
@@ -1,4 +1,10 @@
% Loglik sequence 0
-8188.898628
% Loglik sequence 1
-8206.808088
-1.209996E+01
-1.570194E+01
-1.793114E+01
-1.236031E+01
-1.290610E+01
-1.787109E+01
-1.623056E+01
-1.853745E+01
-1.786704E+01
-1.615386E+01
fastlib/u/tqlong/hmm/loglik.cc
+136
View File
@@ -0,0 +1,136 @@
#include "fastlib/fastlib.h"
#include "support.h"
#include "discreteHMM.h"
#include "gaussianHMM.h"
#include "mixgaussHMM.h"
#include "mixtureDST.h"
success_t loglik_discrete();
success_t loglik_gaussian();
success_t loglik_mixture();
void usage();
int main(int argc, char* argv[]) {
fx_init(argc, argv);
success_t s = SUCCESS_PASS;
if (fx_param_exists(NULL,"type")) {
const char* type = fx_param_str_req(NULL, "type");
if (strcmp(type, "discrete")==0)
s = loglik_discrete();
else if (strcmp(type, "gaussian")==0)
s = loglik_gaussian();
else if (strcmp(type, "mixture")==0)
s = loglik_mixture();
else {
printf("Unrecognized type: must be: discrete | gaussian | mixture !!!\n");
s = SUCCESS_FAIL;
}
}
else {
printf("Unrecognized type: must be: discrete | gaussian | mixture !!!\n");
s = SUCCESS_FAIL;
}
if (!PASSED(s)) usage();
fx_done();
}
void usage() {
printf("\n"
"Usage:\n"
" loglik --type=={discrete|gaussian|mixture} OPTIONS\n"
"[OPTIONS]\n"
" --profile==file : file contains HMM profile\n"
" --seqfile==file : file contains input sequences\n"
" --logfile==file : output file for log-likelihood of the sequences\n"
);
}
success_t loglik_mixture() {
if (!fx_param_exists(NULL, "profile")) {
printf("--profile must be defined.\n");
return SUCCESS_FAIL;
}
const char* profile = fx_param_str_req(NULL, "profile");
const char* seqin = fx_param_str(NULL, "seqfile", "seq.mix.out");
const char* logout = fx_param_str(NULL, "logfile", "log.mix.out");
MixtureofGaussianHMM hmm;
hmm.InitFromFile(profile);
ArrayList<Matrix> seqs;
load_matrix_list(seqin, &seqs);
TextWriter w_log;
if (!PASSED(w_log.Open(logout))) {
NONFATAL("Couldn't open '%s' for writing.", logout);
return SUCCESS_FAIL;
}
ArrayList<double> list_loglik;
hmm.ComputeLogLikelihood(seqs, &list_loglik);
for (int i = 0; i < seqs.size(); i++)
w_log.Printf("%f\n", list_loglik[i]);
return SUCCESS_PASS;
}
success_t loglik_gaussian() {
if (!fx_param_exists(NULL, "profile")) {
printf("--profile must be defined.\n");
return SUCCESS_FAIL;
}
const char* profile = fx_param_str_req(NULL, "profile");
const char* seqin = fx_param_str(NULL, "seqfile", "seq.gauss.out");
const char* logout = fx_param_str(NULL, "logfile", "log.gauss.out");
GaussianHMM hmm;
hmm.InitFromFile(profile);
ArrayList<Matrix> seqs;
load_matrix_list(seqin, &seqs);
TextWriter w_log;
if (!PASSED(w_log.Open(logout))) {
NONFATAL("Couldn't open '%s' for writing.", logout);
return SUCCESS_FAIL;
}
ArrayList<double> list_loglik;
hmm.ComputeLogLikelihood(seqs, &list_loglik);
for (int i = 0; i < seqs.size(); i++)
w_log.Printf("%f\n", list_loglik[i]);
return SUCCESS_PASS;
}
success_t loglik_discrete() {
if (!fx_param_exists(NULL, "profile")) {
printf("--profile must be defined.\n");
return SUCCESS_FAIL;
}
const char* profile = fx_param_str_req(NULL, "profile");
const char* seqin = fx_param_str(NULL, "seqfile", "seq.out");
const char* logout = fx_param_str(NULL, "logfile", "log.out");
DiscreteHMM hmm;
hmm.InitFromFile(profile);
ArrayList<Vector> seqs;
load_vector_list(seqin, &seqs);
TextWriter w_log;
if (!PASSED(w_log.Open(logout))) {
NONFATAL("Couldn't open '%s' for writing.", logout);
return SUCCESS_FAIL;
}
ArrayList<double> list_loglik;
hmm.ComputeLogLikelihood(seqs, &list_loglik);
for (int i = 0; i < seqs.size(); i++)
w_log.Printf("%f\n", list_loglik[i]);
return SUCCESS_PASS;
}
fastlib/u/tqlong/hmm/mixgaussHMM.cc
+71
View File
@@ -3,6 +3,77 @@
#include "mixgaussHMM.h"
#include "gaussianHMM.h"
void MixtureofGaussianHMM::InitFromFile(const char* profile) {
if (!PASSED(load_profileM(profile, &transmission_, &list_mixture_gauss_)))
FATAL("Couldn't open '%s' for reading.", profile);
}
void MixtureofGaussianHMM::LoadProfile(const char* profile) {
transmission_.Destruct();
list_mixture_gauss_.Destruct();
InitFromFile(profile);
}
void MixtureofGaussianHMM::SaveProfile(const char* profile) {
save_profileM(profile, transmission_, list_mixture_gauss_);
}
void MixtureofGaussianHMM::GenerateSequence(int L, Matrix* data_seq, Vector* state_seq) {
hmm_generateM_init(L, transmission_, list_mixture_gauss_, data_seq, state_seq);
}
void forward_procedure(int L, const Matrix& trans, const Matrix& emis_prob, Vector *scales, Matrix* fs);
double MixtureofGaussianHMM::ComputeLogLikelihood(const Matrix& data_seq) {
int L = data_seq.n_cols();
int M = transmission_.n_rows();
Matrix fs(M, L), emis_prob(M, L);
Vector sc;
sc.Init(L);
hmm_cal_emis_probM(data_seq, list_mixture_gauss_, &emis_prob);
forward_procedure(L, transmission_, emis_prob, &sc, &fs);
double loglik = 0;
for (int t = 0; t < L; t++)
loglik += log(sc[t]);
return loglik;
}
void MixtureofGaussianHMM::ComputeLogLikelihood(const ArrayList<Matrix>& list_data_seq, ArrayList<double>* list_likelihood) {
int L = 0;
for (int i = 0; i < list_data_seq.size(); i++)
if (list_data_seq[i].n_cols() > L) L = list_data_seq[i].n_cols();
int M = transmission_.n_rows();
Matrix fs(M, L), emis_prob(M, L);
Vector sc;
sc.Init(L);
list_likelihood->Init();
for (int i = 0; i < list_data_seq.size(); i++) {
int L = list_data_seq[i].n_cols();
hmm_cal_emis_probM(list_data_seq[i], list_mixture_gauss_, &emis_prob);
forward_procedure(L, transmission_, emis_prob, &sc, &fs);
double loglik = 0;
for (int t = 0; t < L; t++)
loglik += log(sc[t]);
list_likelihood->AddBackItem(loglik);
}
}
void MixtureofGaussianHMM::ComputeViterbiStateSequence(const Matrix& data_seq, Vector* state_seq) {
int M = transmission_.n_rows();
int L = data_seq.n_cols();
Matrix emis_prob(M, L);
hmm_cal_emis_probM(data_seq, list_mixture_gauss_, &emis_prob);
hmm_viterbiG_init(transmission_, emis_prob, state_seq);
}
void MixtureofGaussianHMM::TrainBaumWelch(const ArrayList<Matrix>& list_data_seq, int max_iteration, double tolerance) {
hmm_trainM(list_data_seq, &transmission_, &list_mixture_gauss_, max_iteration, tolerance);
}
void MixtureofGaussianHMM::TrainViterbi(const ArrayList<Matrix>& list_data_seq, int max_iteration, double tolerance) {
hmm_train_viterbiM(list_data_seq, &transmission_, &list_mixture_gauss_, max_iteration, tolerance);
}
success_t load_profileM(const char* profile, Matrix* trans, ArrayList<MixtureGauss>* mixs) {
ArrayList<Matrix> matlst;
if (!PASSED(load_matrix_list(profile, &matlst))) {
fastlib/u/tqlong/hmm/mixgaussHMM.h
+25
View File
@@ -4,6 +4,31 @@
#include "fastlib/fastlib.h"
#include "mixtureDST.h"
class MixtureofGaussianHMM {
private:
Matrix transmission_;
ArrayList<MixtureGauss> list_mixture_gauss_;
public:
void InitFromFile(const char* profile);
void Init() {
transmission_.Init(0, 0);
list_mixture_gauss_.Init();
}
void LoadProfile(const char* profile);
void SaveProfile(const char* profile);
void GenerateSequence(int L, Matrix* data_seq, Vector* state_seq);
double ComputeLogLikelihood(const Matrix& data_seq);
void ComputeLogLikelihood(const ArrayList<Matrix>& list_data_seq, ArrayList<double>* list_likelihood);
void ComputeViterbiStateSequence(const Matrix& data_seq, Vector* state_seq);
void TrainBaumWelch(const ArrayList<Matrix>& list_data_seq, int max_iteration = 500, double tolerance = 1e-3);
void TrainViterbi(const ArrayList<Matrix>& list_data_seq, int max_iteration = 500, double tolerance = 1e-3);
};
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);
fastlib/u/tqlong/hmm/mixtureDST.h
+8 -7
View File
@@ -60,13 +60,11 @@ class MixtureGauss {
total++;
}
void accumulate(double p, int i, const Vector& v) {
ACC_prior[i] += p;
la::AddExpert(p, v, &ACC_means[i]);
Vector d;
la::SubInit(v, means[i], &d);
Matrix D;
D.AliasColVector(d);
la::MulExpert(p, false, D, true, D, 1.0, &ACC_covs[i]);
Matrix V;
V.AliasColVector(v);
la::MulExpert(p, false, V, true, V, 1.0, &ACC_covs[i]);
ACC_prior[i] += p;
total += p;
}
void end_accumulate_cluster() {
@@ -88,7 +86,10 @@ class MixtureGauss {
for (int i = 0; i < means.size(); i++) {
if (ACC_prior[i] != 0) {
la::ScaleOverwrite(1.0/ACC_prior[i], ACC_means[i], &means[i]);
la::ScaleOverwrite(1.0/ACC_prior[i], ACC_covs[i], &covs[i]);
Matrix M;
M.AliasColVector(means[i]);
la::MulExpert(-1.0, false, M, true, M, 1.0/ACC_prior[i], &ACC_covs[i]);
covs[i].CopyValues(ACC_covs[i]);
prior[i] = ACC_prior[i]/total;
double det = la::Determinant(covs[i]);
fastlib/u/tqlong/hmm/pro.dis.out
+8 -8
View File
@@ -1,10 +1,10 @@
%% transmision - Matrix (2 x 2) =
0.951217,0.000000,
0.048783,1.000000,
0.913559,0.000041,
0.086441,0.999959,
%% emission - Matrix (2 x 6) =
0.097562,0.124008,
0.146342,0.122000,
0.195122,0.122201,
0.219511,0.121398,
0.146342,0.117984,
0.195122,0.392409,
0.176788,0.032479,
0.117033,0.002811,
0.207122,0.000096,
0.093345,0.052184,
0.169223,0.048246,
0.236488,0.864183,
fastlib/u/tqlong/hmm/pro.gauss.out
+8 -8
View File
@@ -1,13 +1,13 @@
% transmission - Matrix (2 x 2) =
9.008915E-01,5.015528E-02,
9.910846E-02,9.498447E-01,
0.900891,0.050155,
0.099109,0.949845,
% mean - state 0 - Vector (2) =
-1.558390E-02,-5.167563E-03,
-0.015584,-0.005167,
% covariance - state0 - Matrix (2 x 2) =
1.030753E+00,8.553333E-02,
8.553333E-02,9.869338E-01,
1.030753,0.085533,
0.085533,0.986934,
% mean - state 1 - Vector (2) =
1.971526E+00,1.991312E+00,
1.971527,1.991312,
% covariance - state1 - Matrix (2 x 2) =
1.027874E+00,1.014292E-01,
1.014292E-01,1.015718E+00,
1.027874,0.101429,
0.101429,1.015717,
fastlib/u/tqlong/hmm/pro.mix.out
+16 -16
View File
@@ -1,27 +1,27 @@
% transmission - Matrix (2 x 2) =
9.013158E-01,4.244604E-02,
9.868421E-02,9.575540E-01,
9.011601E-01,4.249988E-02,
9.883990E-02,9.575001E-01,
% prior - state 0 - Vector (2) =
5.186639E-01,4.813361E-01,
5.192414E-01,4.807586E-01,
% mean 0 - state 0 - Vector (2) =
-1.188965E-02,-3.753936E-03,
-8.598510E-03,-2.119347E-03,
% covariance 0 - state 0 - Matrix (2 x 2) =
9.505489E-01,5.312182E-02,
5.312182E-02,1.016244E+00,
9.602378E-01,5.660176E-02,
5.660176E-02,1.017587E+00,
% mean 1 - state 0 - Vector (2) =
2.486977E-03,4.987819E+00,
-3.370962E-04,4.986145E+00,
% covariance 1 - state 0 - Matrix (2 x 2) =
9.356824E-01,-2.942855E-02,
-2.942855E-02,1.143288E+00,
9.294454E-01,-3.376060E-02,
-3.376060E-02,1.142541E+00,
% prior - state 1 - Vector (2) =
5.009455E-01,4.990545E-01,
5.007006E-01,4.992994E-01,
% mean 0 - state 1 - Vector (2) =
5.005855E+00,-5.875820E-03,
5.006180E+00,-6.596133E-03,
% covariance 0 - state 1 - Matrix (2 x 2) =
9.461827E-01,4.419038E-02,
4.419038E-02,1.027436E+00,
9.477496E-01,4.588493E-02,
4.588493E-02,1.026848E+00,
% mean 1 - state 1 - Vector (2) =
4.978663E+00,5.000389E+00,
4.977038E+00,5.001125E+00,
% covariance 1 - state 1 - Matrix (2 x 2) =
9.377226E-01,-5.077517E-02,
-5.077517E-02,9.670612E-01,
9.423989E-01,-5.292953E-02,
-5.292953E-02,9.673613E-01,
fastlib/u/tqlong/hmm/seq.out
+20 -4
View File
File diff suppressed because one or more lines are too long
fastlib/u/tqlong/hmm/state.out
+20 -4
View File
File diff suppressed because one or more lines are too long
fastlib/u/tqlong/hmm/state_viterbi.out
+20 -4
View File
File diff suppressed because one or more lines are too long
fastlib/u/tqlong/hmm/train.cc
+295
View File
@@ -0,0 +1,295 @@
#include "fastlib/fastlib.h"
#include "support.h"
#include "discreteHMM.h"
#include "gaussianHMM.h"
#include "mixgaussHMM.h"
#include "mixtureDST.h"
success_t train_baumwelch();
success_t train_viterbi();
void usage();
int main(int argc, char* argv[]) {
fx_init(argc, argv);
success_t s = SUCCESS_PASS;
if (fx_param_exists(NULL,"type")) {
const char* algorithm = fx_param_str(NULL, "algorithm", "baumwelch");
if (strcmp(algorithm,"baumwelch")==0)
s = train_baumwelch();
else if (strcmp(algorithm,"viterbi")==0)
s = train_viterbi();
else {
printf("Unrecognized algorithm: must be baumwelch or viterbi !!!\n");
s = SUCCESS_FAIL;
}
}
else {
printf("Unrecognized type: must be: discrete | gaussian | mixture !!!\n");
s = SUCCESS_FAIL;
}
if (!PASSED(s)) usage();
fx_done();
}
success_t train_baumwelch_discrete();
success_t train_baumwelch_gaussian();
success_t train_baumwelch_mixture();
success_t train_baumwelch() {
const char* type = fx_param_str_req(NULL, "type");
if (strcmp(type, "discrete")==0)
return train_baumwelch_discrete();
else if (strcmp(type, "gaussian")==0)
return train_baumwelch_gaussian();
else if (strcmp(type, "mixture")==0)
return train_baumwelch_mixture();
else {
printf("Unrecognized type: must be: discrete | gaussian | mixture !!!\n");
return SUCCESS_FAIL;
}
}
success_t train_viterbi_discrete();
success_t train_viterbi_gaussian();
success_t train_viterbi_mixture();
success_t train_viterbi() {
const char* type = fx_param_str_req(NULL, "type");
if (strcmp(type, "discrete")==0)
return train_viterbi_discrete();
else if (strcmp(type, "gaussian")==0)
return train_viterbi_gaussian();
else if (strcmp(type, "mixture")==0)
return train_viterbi_mixture();
else {
printf("Unrecognized type: must be: discrete | gaussian | mixture !!!\n");
return SUCCESS_FAIL;
}
}
void usage() {
printf("\nUsage:\n"
" train --type=={discrete|gaussian|mixture} OPTION\n"
"[OPTIONS]\n"
" --algorithm={baumwelch|viterbi} : algorithm used for training, default Baum-Welch\n"
" --seqfile=file : file contains input sequences\n"
" --guess=file : file contains guess HMM profile\n"
" --numstate=NUM : if no guess profile is specified, at least specify the number of state\n"
" --profile=file : output file for estimated HMM profile\n"
" --maxiter=NUM : maximum number of iteration, default=500\n"
" --tolerance=NUM : error tolerance on log-likelihood, default=1e-3\n"
);
}
success_t train_baumwelch_mixture() {
if (!fx_param_exists(NULL, "seqfile")) {
printf("--seqfile must be defined.\n");
return SUCCESS_FAIL;
}
MixtureofGaussianHMM hmm;
ArrayList<Matrix> seqs;
const char* seqin = fx_param_str_req(NULL, "seqfile");
const char* proout = fx_param_str(NULL, "profile", "pro.mix.out");
load_matrix_list(seqin, &seqs);
if (fx_param_exists(NULL, "guess")) { // guessed parameters in a file
const char* guess = fx_param_str_req(NULL, "guess");
printf("Load parameters from file %s\n", guess);
hmm.InitFromFile(guess);
}
else {
hmm.Init();
printf("Automatic initialization not supported !!!");
return SUCCESS_FAIL;
}
int maxiter = fx_param_int(NULL, "maxiter", 500);
double tol = fx_param_double(NULL, "tolerance", 1e-3);
hmm.TrainBaumWelch(seqs, maxiter, tol);
hmm.SaveProfile(proout);
return SUCCESS_PASS;
}
success_t train_baumwelch_gaussian() {
if (!fx_param_exists(NULL, "seqfile")) {
printf("--seqfile must be defined.\n");
return SUCCESS_FAIL;
}
GaussianHMM hmm;
ArrayList<Matrix> seqs;
const char* seqin = fx_param_str_req(NULL, "seqfile");
const char* proout = fx_param_str(NULL, "profile", "pro.gauss.out");
load_matrix_list(seqin, &seqs);
if (fx_param_exists(NULL, "guess")) { // guessed parameters in a file
const char* guess = fx_param_str_req(NULL, "guess");
printf("Load parameters from file %s\n", guess);
hmm.InitFromFile(guess);
}
else { // otherwise initialized using information from the data
int numstate = fx_param_int_req(NULL, "numstate");
printf("Generate HMM parameters: NUMSTATE = %d\n", numstate);
hmm.InitFromData(seqs, numstate);
printf("Done.\n");
}
int maxiter = fx_param_int(NULL, "maxiter", 500);
double tol = fx_param_double(NULL, "tolerance", 1e-3);
printf("Training ...\n");
hmm.TrainBaumWelch(seqs, maxiter, tol);
printf("Done.\n");
hmm.SaveProfile(proout);
return SUCCESS_PASS;
}
success_t train_baumwelch_discrete() {
if (!fx_param_exists(NULL, "seqfile")) {
printf("--seqfile must be defined.\n");
return SUCCESS_FAIL;
}
const char* seqin = fx_param_str_req(NULL, "seqfile");
const char* proout = fx_param_str(NULL, "profile", "pro.dis.out");
ArrayList<Vector> seqs;
load_vector_list(seqin, &seqs);
DiscreteHMM hmm;
if (fx_param_exists(NULL, "guess")) { // guessed parameters in a file
const char* guess = fx_param_str_req(NULL, "guess");
printf("Load HMM parameters from file %s\n", guess);
hmm.InitFromFile(guess);
}
else { // otherwise randomly initialized using information from the data
int numstate = fx_param_int_req(NULL, "numstate");
printf("Randomly generate parameters: NUMSTATE = %d\n", numstate);
hmm.InitFromData(seqs, numstate);
}
int maxiter = fx_param_int(NULL, "maxiter", 500);
double tol = fx_param_double(NULL, "tolerance", 1e-3);
hmm.TrainBaumWelch(seqs, maxiter, tol);
hmm.SaveProfile(proout);
return SUCCESS_PASS;
}
success_t train_viterbi_mixture() {
if (!fx_param_exists(NULL, "seqfile")) {
printf("--seqfile must be defined.\n");
return SUCCESS_FAIL;
}
MixtureofGaussianHMM hmm;
ArrayList<Matrix> seqs;
const char* seqin = fx_param_str_req(NULL, "seqfile");
const char* proout = fx_param_str(NULL, "profile", "pro.mix.out");
load_matrix_list(seqin, &seqs);
if (fx_param_exists(NULL, "guess")) { // guessed parameters in a file
const char* guess = fx_param_str_req(NULL, "guess");
printf("Load parameters from file %s\n", guess);
hmm.InitFromFile(guess);
}
else {
hmm.Init();
printf("Automatic initialization not supported !!!");
return SUCCESS_FAIL;
}
int maxiter = fx_param_int(NULL, "maxiter", 500);
double tol = fx_param_double(NULL, "tolerance", 1e-3);
hmm.TrainViterbi(seqs, maxiter, tol);
hmm.SaveProfile(proout);
return SUCCESS_PASS;
}
success_t train_viterbi_gaussian() {
if (!fx_param_exists(NULL, "seqfile")) {
printf("--seqfile must be defined.\n");
return SUCCESS_FAIL;
}
GaussianHMM hmm;
ArrayList<Matrix> seqs;
const char* seqin = fx_param_str_req(NULL, "seqfile");
const char* proout = fx_param_str(NULL, "profile", "pro.gauss.viterbi.out");
load_matrix_list(seqin, &seqs);
if (fx_param_exists(NULL, "guess")) { // guessed parameters in a file
const char* guess = fx_param_str_req(NULL, "guess");
printf("Load parameters from file %s\n", guess);
hmm.InitFromFile(guess);
}
else { // otherwise initialized using information from the data
int numstate = fx_param_int_req(NULL, "numstate");
printf("Generate parameters: NUMSTATE = %d\n", numstate);
hmm.InitFromData(seqs, numstate);
}
int maxiter = fx_param_int(NULL, "maxiter", 500);
double tol = fx_param_double(NULL, "tolerance", 1e-3);
hmm.TrainViterbi(seqs, maxiter, tol);
hmm.SaveProfile(proout);
return SUCCESS_PASS;
}
success_t train_viterbi_discrete() {
if (!fx_param_exists(NULL, "seqfile")) {
printf("--seqfile must be defined.\n");
return SUCCESS_FAIL;
}
DiscreteHMM hmm;
ArrayList<Vector> seqs;
const char* seqin = fx_param_str_req(NULL, "seqfile");
const char* proout = fx_param_str(NULL, "profile", "pro.dis.viterbi.out");
load_vector_list(seqin, &seqs);
if (fx_param_exists(NULL, "guess")) { // guessed parameters in a file
ArrayList<Matrix> matlst;
const char* guess = fx_param_str_req(NULL, "guess");
printf("Load parameters from file %s\n", guess);
hmm.InitFromFile(guess);
}
else { // otherwise randomly initialized using information from the data
int numstate = fx_param_int_req(NULL, "numstate");
printf("Generate parameters with NUMSTATE = %d\n", numstate);
hmm.InitFromData(seqs, numstate);
}
int maxiter = fx_param_int(NULL, "maxiter", 500);
double tol = fx_param_double(NULL, "tolerance", 1e-3);
hmm.TrainViterbi(seqs, maxiter, tol);
hmm.SaveProfile(proout);
return SUCCESS_PASS;
}
fastlib/u/tqlong/hmm/viterbi.cc
+144
View File
@@ -0,0 +1,144 @@
#include "fastlib/fastlib.h"
#include "support.h"
#include "discreteHMM.h"
#include "gaussianHMM.h"
#include "mixgaussHMM.h"
#include "mixtureDST.h"
success_t viterbi_discrete();
success_t viterbi_gaussian();
success_t viterbi_mixture();
void usage();
int main(int argc, char* argv[]) {
fx_init(argc, argv);
success_t s = SUCCESS_PASS;
if (fx_param_exists(NULL,"type")) {
const char* type = fx_param_str_req(NULL, "type");
if (strcmp(type, "discrete")==0)
s = viterbi_discrete();
else if (strcmp(type, "gaussian")==0)
s = viterbi_gaussian();
else if (strcmp(type, "mixture")==0)
s = viterbi_mixture();
else {
printf("Unrecognized type: must be: discrete | gaussian | mixture !!!\n");
s = SUCCESS_FAIL;
}
}
else {
printf("Unrecognized type: must be: discrete | gaussian | mixture !!!\n");
s = SUCCESS_FAIL;
}
if (!PASSED(s)) usage();
fx_done();
}
void usage() {
printf("\n"
"Usage:\n"
" viterbi --type=={discrete|gaussian|mixture} OPTIONS\n"
"[OPTIONS]\n"
" --profile=file : file contains HMM profile\n"
" --seqfile=file : file contains input sequences\n"
" --statefile=file : output file for state sequences\n"
);
}
success_t viterbi_mixture() {
if (!fx_param_exists(NULL, "profile")) {
printf("--profile must be defined.\n");
return SUCCESS_FAIL;
}
const char* profile = fx_param_str_req(NULL, "profile");
const char* seqin = fx_param_str(NULL, "seqfile", "seq.mix.out");
const char* stateout = fx_param_str(NULL, "statefile", "state.viterbi.mix.out");
MixtureofGaussianHMM hmm;
hmm.InitFromFile(profile);
ArrayList<Matrix> seqs;
load_matrix_list(seqin, &seqs);
TextWriter w_state;
if (!PASSED(w_state.Open(stateout))) {
NONFATAL("Couldn't open '%s' for writing.", stateout);
return SUCCESS_FAIL;
}
for (int i = 0; i < seqs.size(); i++) {
Vector states;
char s[100];
hmm.ComputeViterbiStateSequence(seqs[i], &states);
sprintf(s, "%% viterbi state sequence %d", i);
print_vector(w_state, states, s, "%.0f,");
}
return SUCCESS_PASS;
}
success_t viterbi_gaussian() {
if (!fx_param_exists(NULL, "profile")) {
printf("--profile must be defined.\n");
return SUCCESS_FAIL;
}
const char* profile = fx_param_str_req(NULL, "profile");
const char* seqin = fx_param_str(NULL, "seqfile", "seq.gauss.out");
const char* stateout = fx_param_str(NULL, "statefile", "state.viterbi.gauss.out");
GaussianHMM hmm;
hmm.InitFromFile(profile);
ArrayList<Matrix> seqs;
load_matrix_list(seqin, &seqs);
TextWriter w_state;
if (!PASSED(w_state.Open(stateout))) {
NONFATAL("Couldn't open '%s' for writing.", stateout);
return SUCCESS_FAIL;
}
for (int i = 0; i < seqs.size(); i++) {
Vector states;
char s[100];
hmm.ComputeViterbiStateSequence(seqs[i], &states);
sprintf(s, "%% viterbi state sequence %d", i);
print_vector(w_state, states, s, "%.0f,");
}
return SUCCESS_PASS;
}
success_t viterbi_discrete() {
if (!fx_param_exists(NULL, "profile")) {
printf("--profile must be defined.\n");
return SUCCESS_FAIL;
}
const char* profile = fx_param_str_req(NULL, "profile");
const char* seqin = fx_param_str(NULL, "seqfile", "seq.out");
const char* stateout = fx_param_str(NULL, "statefile", "state.viterbi.out");
DiscreteHMM hmm;
hmm.InitFromFile(profile);
ArrayList<Vector> seqs;
load_vector_list(seqin, &seqs);
TextWriter w_state;
if (!PASSED(w_state.Open(stateout))) {
NONFATAL("Couldn't open '%s' for writing.", stateout);
return SUCCESS_FAIL;
}
for (int i = 0; i < seqs.size(); i++) {
Vector states;
char s[100];
hmm.ComputeViterbiStateSequence(seqs[i], &states);
sprintf(s, "%% viterbi state sequence %d", i);
print_vector(w_state, states, s, "%.0f,");
}
return SUCCESS_PASS;
}