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

More stuff compiles, all untested for correctness.

Browse Source
This commit is contained in:
James Cline
2010-06-15 21:31:33 +00:00
parent 6c50cd909a
commit 9b60acd638
16 changed files with 326 additions and 348 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/branches/fastlib-stl/mlpack/hmm/discreteHMM.cc
+17 -16
View File
@@ -7,6 +7,7 @@
#include "fastlib/fastlib.h"
#include "support.h"
#include "discreteHMM.h"
#include <algorithm>
using namespace hmm_support;
@@ -27,7 +28,7 @@ void DiscreteHMM::Init(const Matrix& transmission, const Matrix& emission) {
}
void DiscreteHMM::InitFromFile(const char* profile) {
ArrayList<Matrix> list_mat;
std::vector<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.");
@@ -39,19 +40,20 @@ void DiscreteHMM::InitFromFile(const char* profile) {
DEBUG_ASSERT(transmission_.n_rows() == emission_.n_rows());
}
void DiscreteHMM::InitFromData(const ArrayList<Vector>& list_data_seq, int numstate) {
void DiscreteHMM::InitFromData(const std::vector<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];
std::vector<Vector>::const_iterator maxseq = list_data_seq.begin();
for( std::vector<Vector>::const_iterator i = list_data_seq.begin();
i < list_data_seq.end(); ++i )
maxseq = (maxseq->length() < i->length())?i:maxseq;
for (int i = 0; i < maxseq->length(); i++)
if ((*maxseq)[i] > numsymbol) numsymbol = (int) (*maxseq)[i];
numsymbol++;
Vector states;
int L = list_data_seq[maxseq].length();
int L = maxseq->length();
states.Init(L);
for (int i = 0; i < L; i++) states[i] = rand() % numstate;
DiscreteHMM::EstimateInit(numsymbol, numstate, list_data_seq[maxseq], states, &transmission_, &emission_);
DiscreteHMM::EstimateInit(numsymbol, numstate, *maxseq, states, &transmission_, &emission_);
}
void DiscreteHMM::LoadProfile(const char* profile) {
@@ -116,7 +118,7 @@ double DiscreteHMM::ComputeLogLikelihood(const Vector& data_seq) const {
return loglik;
}
void DiscreteHMM::ComputeLogLikelihood(const ArrayList<Vector>& list_data_seq, ArrayList<double>* list_likelihood) const {
void DiscreteHMM::ComputeLogLikelihood(const std::vector<Vector>& list_data_seq, std::vector<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();
@@ -124,14 +126,13 @@ void DiscreteHMM::ComputeLogLikelihood(const ArrayList<Vector>& list_data_seq, A
Matrix fs(M, L);
Vector sc;
sc.Init(L);
list_likelihood->Init();
for (int i = 0; i < list_data_seq.size(); i++) {
DiscreteHMM::ForwardProcedure(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->PushBackCopy(loglik);
list_likelihood->push_back(loglik);
}
}
@@ -139,11 +140,11 @@ void DiscreteHMM::ComputeViterbiStateSequence(const Vector& data_seq, Vector* st
DiscreteHMM::ViterbiInit(data_seq, transmission_, emission_, state_seq);
}
void DiscreteHMM::TrainBaumWelch(const ArrayList<Vector>& list_data_seq, int max_iteration, double tolerance) {
void DiscreteHMM::TrainBaumWelch(const std::vector<Vector>& list_data_seq, int max_iteration, double tolerance) {
DiscreteHMM::Train(list_data_seq, &transmission_, &emission_, max_iteration, tolerance);
}
void DiscreteHMM::TrainViterbi(const ArrayList<Vector>& list_data_seq, int max_iteration, double tolerance) {
void DiscreteHMM::TrainViterbi(const std::vector<Vector>& list_data_seq, int max_iteration, double tolerance) {
DiscreteHMM::TrainViterbi(list_data_seq, &transmission_, &emission_, max_iteration, tolerance);
}
@@ -388,7 +389,7 @@ double DiscreteHMM::ViterbiInit(int L, const Vector& seq, const Matrix& trans, c
return bestVal;
}
void DiscreteHMM::Train(const ArrayList<Vector>& seqs, Matrix* guessTR, Matrix* guessEM, int max_iter, double tol) {
void DiscreteHMM::Train(const std::vector<Vector>& seqs, Matrix* guessTR, Matrix* guessEM, int max_iter, double tol) {
int L = -1;
int M = guessTR->n_rows();
int N = guessEM->n_cols();
@@ -461,7 +462,7 @@ void DiscreteHMM::Train(const ArrayList<Vector>& seqs, Matrix* guessTR, Matrix*
}
}
void DiscreteHMM::TrainViterbi(const ArrayList<Vector>& seqs, Matrix* guessTR, Matrix* guessEM, int max_iter, double tol) {
void DiscreteHMM::TrainViterbi(const std::vector<Vector>& seqs, Matrix* guessTR, Matrix* guessEM, int max_iter, double tol) {
int L = -1;
int M = guessTR->n_rows();
int N = guessEM->n_cols();
fastlib/branches/fastlib-stl/mlpack/hmm/discreteHMM.h
+6 -6
View File
@@ -50,7 +50,7 @@ class DiscreteHMM {
void InitFromFile(const char* profile);
/** Initializes randomly using data as a guide */
void InitFromData(const ArrayList<Vector>& list_data_seq, int numstate);
void InitFromData(const std::vector<Vector>& list_data_seq, int numstate);
/** Load from file, used when already initialized */
void LoadProfile(const char* profile);
@@ -82,7 +82,7 @@ class DiscreteHMM {
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;
void ComputeLogLikelihood(const std::vector<Vector>& list_data_seq, std::vector<double>* list_likelihood) const;
/** Compute the most probable sequence (Viterbi) */
void ComputeViterbiStateSequence(const Vector& data_seq, Vector* state_seq) const;
@@ -91,13 +91,13 @@ class DiscreteHMM {
* 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, double tolerance);
void TrainBaumWelch(const std::vector<Vector>& list_data_seq, int max_iteration, double tolerance);
/**
* 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, double tolerance);
void TrainViterbi(const std::vector<Vector>& list_data_seq, int max_iteration, double tolerance);
///////// Static helper functions ///////////////////////////////////////
@@ -143,10 +143,10 @@ class DiscreteHMM {
static double ViterbiInit(int L, const Vector& seq, const Matrix& trans, const Matrix& emis, Vector* states);
/** Baum-Welch estimation of transition and emission probabilities */
static void Train(const ArrayList<Vector>& seqs, Matrix* guessTR, Matrix* guessEM, int max_iter, double tol);
static void Train(const std::vector<Vector>& seqs, Matrix* guessTR, Matrix* guessEM, int max_iter, double tol);
/** Viterbi estimation of transition and emission probabilities */
static void TrainViterbi(const ArrayList<Vector>& seqs, Matrix* guessTR, Matrix* guessEM, int max_iter, double tol);
static void TrainViterbi(const std::vector<Vector>& seqs, Matrix* guessTR, Matrix* guessEM, int max_iter, double tol);
};
#endif
fastlib/branches/fastlib-stl/mlpack/hmm/gaussianHMM.cc
+51 -61
View File
@@ -10,7 +10,7 @@
using namespace hmm_support;
void GaussianHMM::setModel(const Matrix& transmission, const ArrayList<Vector>& list_mean_vec,const ArrayList<Matrix>& list_covariance_mat) {
void GaussianHMM::setModel(const Matrix& transmission, const std::vector<Vector>& list_mean_vec,const std::vector<Matrix>& list_covariance_mat) {
DEBUG_ASSERT(transmission.n_rows() == transmission.n_cols());
DEBUG_ASSERT(transmission.n_rows() == list_mean_vec.size());
DEBUG_ASSERT(transmission.n_rows() == list_covariance_mat.size());
@@ -20,18 +20,16 @@ void GaussianHMM::setModel(const Matrix& transmission, const ArrayList<Vector>&
DEBUG_ASSERT(list_mean_vec[0].length() == list_mean_vec[i].length());
}
transmission_.Destruct();
list_mean_vec_.Renew();
list_covariance_mat_.Renew();
transmission_.Copy(transmission);
list_mean_vec_.InitCopy(list_mean_vec);
list_covariance_mat_.InitCopy(list_covariance_mat);
list_mean_vec_.assign(list_mean_vec.begin(), list_mean_vec.end());
list_covariance_mat_.assign(list_covariance_mat.begin(), list_covariance_mat.end());
CalculateInverse();
}
void GaussianHMM::Init(const Matrix& transmission, const ArrayList<Vector>& list_mean_vec,const ArrayList<Matrix>& list_covariance_mat) {
void GaussianHMM::Init(const Matrix& transmission, const std::vector<Vector>& list_mean_vec,const std::vector<Matrix>& list_covariance_mat) {
transmission_.Copy(transmission);
list_mean_vec_.InitCopy(list_mean_vec);
list_covariance_mat_.InitCopy(list_covariance_mat);
list_mean_vec_.assign(list_mean_vec.begin(), list_mean_vec.end());
list_covariance_mat_.assign(list_covariance_mat.begin(),list_covariance_mat.end());
DEBUG_ASSERT(transmission.n_rows() == transmission.n_cols());
DEBUG_ASSERT(transmission.n_rows() == list_mean_vec.size());
DEBUG_ASSERT(transmission.n_rows() == list_covariance_mat.size());
@@ -46,14 +44,14 @@ void GaussianHMM::Init(const Matrix& transmission, const ArrayList<Vector>& lis
void GaussianHMM::InitFromFile(const char* profile) {
if (!PASSED(GaussianHMM::LoadProfile(profile, &transmission_, &list_mean_vec_, &list_covariance_mat_)))
FATAL("Couldn't open '%s' for reading.", profile);
list_inverse_cov_mat_.InitCopy(list_covariance_mat_);
list_inverse_cov_mat_.assign(list_covariance_mat_.begin(),list_covariance_mat_.end());
gauss_const_vec_.Init(list_covariance_mat_.size());
CalculateInverse();
}
void GaussianHMM::InitFromData(const ArrayList<Matrix>& list_data_seq, int numstate) {
void GaussianHMM::InitFromData(const std::vector<Matrix>& list_data_seq, int numstate) {
GaussianHMM::InitGaussParameter(numstate, list_data_seq, &transmission_, &list_mean_vec_, &list_covariance_mat_);
list_inverse_cov_mat_.InitCopy(list_covariance_mat_);
list_inverse_cov_mat_.assign(list_covariance_mat_.begin(), list_covariance_mat_.end());
gauss_const_vec_.Init(list_covariance_mat_.size());
CalculateInverse();
}
@@ -65,11 +63,13 @@ void GaussianHMM::InitFromData(const Matrix& data_seq, const Vector& state_seq)
}
void GaussianHMM::LoadProfile(const char* profile) {
/*
transmission_.Destruct();
list_mean_vec_.Renew();
list_covariance_mat_.Renew();
list_inverse_cov_mat_.Renew();
gauss_const_vec_.Destruct();
*/
InitFromFile(profile);
}
@@ -92,8 +92,6 @@ void GaussianHMM::GenerateSequence(int L, Matrix* data_seq, Vector* state_seq) c
void GaussianHMM::EstimateModel(const Matrix& data_seq, const Vector& state_seq) {
transmission_.Destruct();
list_mean_vec_.Renew();
list_covariance_mat_.Renew();
GaussianHMM::EstimateInit(data_seq, state_seq, &transmission_,
&list_mean_vec_, &list_covariance_mat_);
CalculateInverse();
@@ -101,8 +99,6 @@ void GaussianHMM::EstimateModel(const Matrix& data_seq, const Vector& state_seq)
void GaussianHMM::EstimateModel(int numstate, const Matrix& data_seq, const Vector& state_seq) {
transmission_.Destruct();
list_mean_vec_.Renew();
list_covariance_mat_.Renew();
GaussianHMM::EstimateInit(numstate, data_seq, state_seq, &transmission_,
&list_mean_vec_, &list_covariance_mat_);
CalculateInverse();
@@ -148,7 +144,7 @@ double GaussianHMM::ComputeLogLikelihood(const Matrix& data_seq) const {
return loglik;
}
void GaussianHMM::ComputeLogLikelihood(const ArrayList<Matrix>& list_data_seq, ArrayList<double>* list_likelihood) const {
void GaussianHMM::ComputeLogLikelihood(const std::vector<Matrix>& list_data_seq, std::vector<double>* list_likelihood) const {
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();
@@ -156,7 +152,6 @@ void GaussianHMM::ComputeLogLikelihood(const ArrayList<Matrix>& list_data_seq, A
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();
GaussianHMM::CalculateEmissionProb(list_data_seq[i], list_mean_vec_, list_inverse_cov_mat_, gauss_const_vec_, &emis_prob);
@@ -164,7 +159,7 @@ void GaussianHMM::ComputeLogLikelihood(const ArrayList<Matrix>& list_data_seq, A
double loglik = 0;
for (int t = 0; t < L; t++)
loglik += log(sc[t]);
list_likelihood->PushBackCopy(loglik);
list_likelihood->push_back(loglik);
}
}
@@ -176,25 +171,23 @@ void GaussianHMM::ComputeViterbiStateSequence(const Matrix& data_seq, Vector* st
GaussianHMM::ViterbiInit(transmission_, emis_prob, state_seq);
}
void GaussianHMM::TrainBaumWelch(const ArrayList<Matrix>& list_data_seq, int max_iteration, double tolerance) {
void GaussianHMM::TrainBaumWelch(const std::vector<Matrix>& list_data_seq, int max_iteration, double tolerance) {
GaussianHMM::Train(list_data_seq, &transmission_, &list_mean_vec_, &list_covariance_mat_, max_iteration, tolerance);
CalculateInverse();
}
void GaussianHMM::TrainViterbi(const ArrayList<Matrix>& list_data_seq, int max_iteration, double tolerance) {
void GaussianHMM::TrainViterbi(const std::vector<Matrix>& list_data_seq, int max_iteration, double tolerance) {
GaussianHMM::TrainViterbi(list_data_seq, &transmission_, &list_mean_vec_, &list_covariance_mat_, max_iteration, tolerance);
CalculateInverse();
}
success_t GaussianHMM::LoadProfile(const char* profile, Matrix* trans, ArrayList<Vector>* means, ArrayList<Matrix>* covs) {
ArrayList<Matrix> matlst;
success_t GaussianHMM::LoadProfile(const char* profile, Matrix* trans, std::vector<Vector>* means, std::vector<Matrix>* covs) {
std::vector<Matrix> matlst;
if (!PASSED(load_matrix_list(profile, &matlst)))
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
@@ -203,13 +196,13 @@ success_t GaussianHMM::LoadProfile(const char* profile, Matrix* trans, ArrayList
DEBUG_ASSERT(matlst[i+1].n_rows()==N && matlst[i+1].n_cols()==N);
Vector m;
matlst[i].MakeColumnVector(0, &m);
means->PushBackCopy(m);
covs->PushBackCopy(matlst[i+1]);
means->push_back(m);
covs->push_back(matlst[i+1]);
}
return SUCCESS_PASS;
}
success_t GaussianHMM::SaveProfile(const char* profile, const Matrix& trans, const ArrayList<Vector>& means, const ArrayList<Matrix>& covs) {
success_t GaussianHMM::SaveProfile(const char* profile, const Matrix& trans, const std::vector<Vector>& means, const std::vector<Matrix>& covs) {
TextWriter w_pro;
if (!PASSED(w_pro.Open(profile))) {
NONFATAL("Couldn't open '%s' for writing.", profile);
@@ -231,7 +224,7 @@ success_t GaussianHMM::SaveProfile(const char* profile, const Matrix& trans, con
return SUCCESS_PASS;
}
void GaussianHMM::GenerateInit(int L, const Matrix& trans, const ArrayList<Vector>& means, const ArrayList<Matrix>& covs, Matrix* seq, Vector* states){
void GaussianHMM::GenerateInit(int L, const Matrix& trans, const std::vector<Vector>& means, const std::vector<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;
Matrix& seq_ = *seq;
@@ -271,7 +264,7 @@ void GaussianHMM::GenerateInit(int L, const Matrix& trans, const ArrayList<Vecto
}
}
void GaussianHMM::EstimateInit(const Matrix& seq, const Vector& states, Matrix* trans, ArrayList<Vector>* means, ArrayList<Matrix>* covs) {
void GaussianHMM::EstimateInit(const Matrix& seq, const Vector& states, Matrix* trans, std::vector<Vector>* means, std::vector<Matrix>* covs) {
DEBUG_ASSERT_MSG((seq.n_cols()==states.length()), "hmm_estimateG_init: sequence and states length must be the same");
int M = 0;
for (int i = 0; i < seq.n_cols(); i++)
@@ -280,7 +273,7 @@ void GaussianHMM::EstimateInit(const Matrix& seq, const Vector& states, Matrix*
GaussianHMM::EstimateInit(M, seq, states, trans, means, covs);
}
void GaussianHMM::EstimateInit(int numStates, const Matrix& seq, const Vector& states, Matrix* trans, ArrayList<Vector>* means, ArrayList<Matrix>* covs) {
void GaussianHMM::EstimateInit(int numStates, const Matrix& seq, const Vector& states, Matrix* trans, std::vector<Vector>* means, std::vector<Matrix>* covs) {
DEBUG_ASSERT_MSG((seq.n_cols()==states.length()), "hmm_estimateD_init: sequence and states length must be the same");
int N = seq.n_rows(); // emission vector length
@@ -288,24 +281,22 @@ void GaussianHMM::EstimateInit(int numStates, const Matrix& seq, const Vector& s
int L = seq.n_cols(); // sequence length
Matrix &trans_ = *trans;
ArrayList<Vector>& mean_ = *means;
ArrayList<Matrix>& cov_ = *covs;
std::vector<Vector>& mean_ = *means;
std::vector<Matrix>& cov_ = *covs;
Vector stateSum;
trans_.Init(M, M);
mean_.Init();
cov_.Init();
stateSum.Init(M);
trans_.SetZero();
for (int i = 0; i < M; i++) {
Vector m;
m.Init(N); m.SetZero();
mean_.PushBackCopy(m);
mean_.push_back(m);
Matrix c;
c.Init(N, N); c.SetZero();
cov_.PushBackCopy(c);
cov_.push_back(c);
}
stateSum.SetZero();
@@ -487,7 +478,7 @@ double GaussianHMM::ViterbiInit(int L, const Matrix& trans, const Matrix& emis_p
return bestVal;
}
void GaussianHMM::CalculateEmissionProb(const Matrix& seq, const ArrayList<Vector>& means, const ArrayList<Matrix>& inv_covs, const Vector& det, Matrix* emis_prob) {
void GaussianHMM::CalculateEmissionProb(const Matrix& seq, const std::vector<Vector>& means, const std::vector<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++) {
@@ -498,12 +489,12 @@ void GaussianHMM::CalculateEmissionProb(const Matrix& seq, const ArrayList<Vecto
}
}
void GaussianHMM::InitGaussParameter(int M, const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO) {
void GaussianHMM::InitGaussParameter(int M, const std::vector<Matrix>& seqs, Matrix* guessTR, std::vector<Vector>* guessME, std::vector<Matrix>* guessCO) {
int N = seqs[0].n_rows();
Matrix& gTR = *guessTR;
ArrayList<Vector>& gME = *guessME;
ArrayList<Matrix>& gCO = *guessCO;
ArrayList<int> labels;
std::vector<Vector>& gME = *guessME;
std::vector<Matrix>& gCO = *guessCO;
std::vector<int> labels;
Vector sumState;
kmeans(seqs, M, &labels, &gME, 1000, 1e-5);
@@ -513,11 +504,10 @@ void GaussianHMM::InitGaussParameter(int M, const ArrayList<Matrix>& seqs, Matri
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.PushBackCopy(m);
gCO.push_back(m);
}
//printf("---2---\n");
@@ -559,10 +549,10 @@ void GaussianHMM::InitGaussParameter(int M, const ArrayList<Matrix>& seqs, Matri
//printf("---4---\n");
}
void GaussianHMM::TrainViterbi(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO, int max_iter, double tol) {
void GaussianHMM::TrainViterbi(const std::vector<Matrix>& seqs, Matrix* guessTR, std::vector<Vector>* guessME, std::vector<Matrix>* guessCO, int max_iter, double tol) {
Matrix &gTR = *guessTR;
ArrayList<Vector>& gME = *guessME;
ArrayList<Matrix>& gCO = *guessCO;
std::vector<Vector>& gME = *guessME;
std::vector<Matrix>& gCO = *guessCO;
int L = -1;
int M = gTR.n_rows();
int N = gME[0].length();
@@ -572,14 +562,14 @@ void GaussianHMM::TrainViterbi(const ArrayList<Matrix>& seqs, Matrix* guessTR, A
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
std::vector<Vector> ME; // accumulating mean
std::vector<Matrix> CO; // accumulating covariance
std::vector<Matrix> INV_CO; // inverse matrix of the covariance
Vector DET; // the determinant * constant of the Normal PDF formula
TR.Init(M, M);
ME.InitCopy(gME);
CO.InitCopy(gCO);
INV_CO.InitCopy(CO);
ME.assign(gME.begin(),gME.end());
CO.assign(gCO.begin(),gCO.end());
INV_CO.assign(CO.begin(),CO.end());
DET.Init(M);
Matrix emis_prob;
@@ -662,10 +652,10 @@ void GaussianHMM::TrainViterbi(const ArrayList<Matrix>& seqs, Matrix* guessTR, A
}
void GaussianHMM::Train(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO, int max_iter, double tol) {
void GaussianHMM::Train(const std::vector<Matrix>& seqs, Matrix* guessTR, std::vector<Vector>* guessME, std::vector<Matrix>* guessCO, int max_iter, double tol) {
Matrix &gTR = *guessTR;
ArrayList<Vector>& gME = *guessME;
ArrayList<Matrix>& gCO = *guessCO;
std::vector<Vector>& gME = *guessME;
std::vector<Matrix>& gCO = *guessCO;
int L = -1;
int M = gTR.n_rows();
int N = gME[0].length();
@@ -675,14 +665,14 @@ void GaussianHMM::Train(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayLis
if (seqs[i].n_cols() > L) L = seqs[i].n_cols();
Matrix TR; // guess transition and emission matrix
ArrayList<Vector> ME; // accumulating mean
ArrayList<Matrix> CO; // accumulating covariance
ArrayList<Matrix> INV_CO; // inverse matrix of the covariance
std::vector<Vector> ME; // accumulating mean
std::vector<Matrix> CO; // accumulating covariance
std::vector<Matrix> INV_CO; // inverse matrix of the covariance
Vector DET; // the determinant * constant of the Normal PDF formula
TR.Init(M, M);
ME.InitCopy(gME);
CO.InitCopy(gCO);
INV_CO.InitCopy(CO);
ME.assign(gME.begin(),gME.end());
CO.assign(gCO.begin(),gCO.end());
INV_CO.assign(CO.begin(),CO.end());
DET.Init(M);
Matrix ps, fs, bs, emis_prob; // to hold hmm_decodeG results
fastlib/branches/fastlib-stl/mlpack/hmm/gaussianHMM.h
+30 -30
View File
@@ -26,13 +26,13 @@ class GaussianHMM {
Matrix transmission_;
/** List of mean vectors */
ArrayList<Vector> list_mean_vec_;
std::vector<Vector> list_mean_vec_;
/** List of covariance matrices */
ArrayList<Matrix> list_covariance_mat_;
std::vector<Matrix> list_covariance_mat_;
/** List of inverse of the covariances */
ArrayList<Matrix> list_inverse_cov_mat_;
std::vector<Matrix> list_inverse_cov_mat_;
/** Vector of constant in the gaussian density fomular */
Vector gauss_const_vec_;
@@ -49,22 +49,22 @@ class GaussianHMM {
public:
/** Getters */
const Matrix& transmission() const { return transmission_; }
const ArrayList<Vector>& list_mean_vec() const { return list_mean_vec_; }
const ArrayList<Matrix>& list_covariance_mat() const { return list_covariance_mat_; }
const std::vector<Vector>& list_mean_vec() const { return list_mean_vec_; }
const std::vector<Matrix>& list_covariance_mat() const { return list_covariance_mat_; }
/** Setter used when already initialized */
void setModel(const Matrix& transmission, const ArrayList<Vector>& list_mean_vec,
const ArrayList<Matrix>& list_covariance_mat);
void setModel(const Matrix& transmission, const std::vector<Vector>& list_mean_vec,
const std::vector<Matrix>& list_covariance_mat);
/** Initializes from computed transmission and Gaussian parameters */
void Init(const Matrix& transmission, const ArrayList<Vector>& list_mean_vec,
const ArrayList<Matrix>& list_covariance_mat);
void Init(const Matrix& transmission, const std::vector<Vector>& list_mean_vec,
const std::vector<Matrix>& list_covariance_mat);
/** Initializes by loading from a file */
void InitFromFile(const char* profile);
/** Initializes using K-means algorithm using data as a guide */
void InitFromData(const ArrayList<Matrix>& list_data_seq, int numstate);
void InitFromData(const std::vector<Matrix>& list_data_seq, int numstate);
/** Initializes using data and state sequence as a guide */
void InitFromData(const Matrix& data_seq, const Vector& state_seq);
@@ -103,8 +103,8 @@ class GaussianHMM {
double ComputeLogLikelihood(const Matrix& data_seq) const;
/** Compute the log-likelihood of a list of sequences */
void ComputeLogLikelihood(const ArrayList<Matrix>& list_data_seq,
ArrayList<double>* list_likelihood) const;
void ComputeLogLikelihood(const std::vector<Matrix>& list_data_seq,
std::vector<double>* list_likelihood) const;
/** Compute the most probable sequence (Viterbi) */
void ComputeViterbiStateSequence(const Matrix& data_seq, Vector* state_seq) const;
@@ -113,24 +113,24 @@ class GaussianHMM {
* Train the model with a list of sequences, must be already initialized
* using Baum-Welch EM algorithm
*/
void TrainBaumWelch(const ArrayList<Matrix>& list_data_seq,
void TrainBaumWelch(const std::vector<Matrix>& list_data_seq,
int max_iteration, double tolerance);
/**
* 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<Matrix>& list_data_seq,
void TrainViterbi(const std::vector<Matrix>& list_data_seq,
int max_iteration, double tolerance);
////////// Static helper functions ///////////////////////////////////////
static success_t LoadProfile(const char* profile, Matrix* trans,
ArrayList<Vector>* means, ArrayList<Matrix>* covs);
std::vector<Vector>* means, std::vector<Matrix>* covs);
static success_t SaveProfile(const char* profile, const Matrix& trans,
const ArrayList<Vector>& means,
const ArrayList<Matrix>& covs);
const std::vector<Vector>& means,
const std::vector<Matrix>& covs);
/**
* Generating a sequence and states using transition and emission probabilities.
* L: sequence length
@@ -140,15 +140,15 @@ class GaussianHMM {
* seq: generated sequence, uninitialized matrix, will have size N x L
* states: generated states, uninitialized vector, will have length L
*/
static void GenerateInit(int L, const Matrix& trans, const ArrayList<Vector>& means,
const ArrayList<Matrix>& covs, Matrix* seq, Vector* states);
static void GenerateInit(int L, const Matrix& trans, const std::vector<Vector>& means,
const std::vector<Matrix>& covs, Matrix* seq, Vector* states);
/** Estimate transition and emission distribution from sequence and states */
static void EstimateInit(const Matrix& seq, const Vector& states, Matrix* trans,
ArrayList<Vector>* means, ArrayList<Matrix>* covs);
std::vector<Vector>* means, std::vector<Matrix>* covs);
static void EstimateInit(int numStates, const Matrix& seq, const Vector& states,
Matrix* trans, ArrayList<Vector>* means,
ArrayList<Matrix>* covs);
Matrix* trans, std::vector<Vector>* means,
std::vector<Matrix>* covs);
/**
* Calculate posteriori probabilities of states at each steps
@@ -170,8 +170,8 @@ class GaussianHMM {
Matrix* fs, Matrix* bs, Vector* scales);
static double Decode(int L, const Matrix& trans, const Matrix& emis_prob,
Matrix* pstates, Matrix* fs, Matrix* bs, Vector* scales);
static void CalculateEmissionProb(const Matrix& seq, const ArrayList<Vector>& means,
const ArrayList<Matrix>& inv_covs, const Vector& det,
static void CalculateEmissionProb(const Matrix& seq, const std::vector<Vector>& means,
const std::vector<Matrix>& inv_covs, const Vector& det,
Matrix* emis_prob);
/**
@@ -189,15 +189,15 @@ class GaussianHMM {
* Baum-Welch and Viterbi estimation of transition and emission
* distribution (Gaussian)
*/
static void InitGaussParameter(int M, const ArrayList<Matrix>& seqs,
Matrix* guessTR, ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO);
static void InitGaussParameter(int M, const std::vector<Matrix>& seqs,
Matrix* guessTR, std::vector<Vector>* guessME, std::vector<Matrix>* guessCO);
static void Train(const ArrayList<Matrix>& seqs, Matrix* guessTR,
ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO,
static void Train(const std::vector<Matrix>& seqs, Matrix* guessTR,
std::vector<Vector>* guessME, std::vector<Matrix>* guessCO,
int max_iter, double tol);
static void TrainViterbi(const ArrayList<Matrix>& seqs, Matrix* guessTR,
ArrayList<Vector>* guessME, ArrayList<Matrix>* guessCO,
static void TrainViterbi(const std::vector<Matrix>& seqs, Matrix* guessTR,
std::vector<Vector>* guessME, std::vector<Matrix>* guessCO,
int max_iter, double tol);
};
#endif
fastlib/branches/fastlib-stl/mlpack/hmm/mixgaussHMM.cc
+28 -38
View File
@@ -12,13 +12,12 @@
using namespace hmm_support;
void MixtureofGaussianHMM::setModel(const Matrix& transmission,
const ArrayList<MixtureGauss>& list_mixture_gauss) {
const std::vector<MixtureGauss>& list_mixture_gauss) {
DEBUG_ASSERT(transmission.n_rows() == transmission.n_cols());
DEBUG_ASSERT(transmission.n_rows() == list_mixture_gauss.size());
transmission_.Destruct();
list_mixture_gauss_.Renew();
transmission_.Copy(transmission);
list_mixture_gauss_.InitCopy(list_mixture_gauss);
list_mixture_gauss_.assign(list_mixture_gauss.begin(),list_mixture_gauss.end());
}
void MixtureofGaussianHMM::InitFromFile(const char* profile) {
@@ -28,7 +27,6 @@ void MixtureofGaussianHMM::InitFromFile(const char* profile) {
void MixtureofGaussianHMM::LoadProfile(const char* profile) {
transmission_.Destruct();
list_mixture_gauss_.Renew();
InitFromFile(profile);
}
@@ -43,7 +41,6 @@ void MixtureofGaussianHMM::GenerateSequence(int L, Matrix* data_seq, Vector* sta
void MixtureofGaussianHMM::EstimateModel(int numcluster, const Matrix& data_seq,
const Vector& state_seq) {
transmission_.Destruct();
list_mixture_gauss_.Renew();
MixtureofGaussianHMM::EstimateInit(numcluster, data_seq, state_seq, &transmission_,
&list_mixture_gauss_);
}
@@ -51,7 +48,6 @@ void MixtureofGaussianHMM::EstimateModel(int numcluster, const Matrix& data_seq,
void MixtureofGaussianHMM::EstimateModel(int numstate, int numcluster,
const Matrix& data_seq, const Vector& state_seq) {
transmission_.Destruct();
list_mixture_gauss_.Renew();
MixtureofGaussianHMM::EstimateInit(numstate, numcluster, data_seq, state_seq,
&transmission_, &list_mixture_gauss_);
}
@@ -99,7 +95,7 @@ double MixtureofGaussianHMM::ComputeLogLikelihood(const Matrix& data_seq) const
return loglik;
}
void MixtureofGaussianHMM::ComputeLogLikelihood(const ArrayList<Matrix>& list_data_seq, ArrayList<double>* list_likelihood) const {
void MixtureofGaussianHMM::ComputeLogLikelihood(const std::vector<Matrix>& list_data_seq, std::vector<double>* list_likelihood) const {
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();
@@ -107,7 +103,6 @@ void MixtureofGaussianHMM::ComputeLogLikelihood(const ArrayList<Matrix>& list_da
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();
MixtureofGaussianHMM::CalculateEmissionProb(list_data_seq[i], list_mixture_gauss_, &emis_prob);
@@ -115,7 +110,7 @@ void MixtureofGaussianHMM::ComputeLogLikelihood(const ArrayList<Matrix>& list_da
double loglik = 0;
for (int t = 0; t < L; t++)
loglik += log(sc[t]);
list_likelihood->PushBackCopy(loglik);
list_likelihood->push_back(loglik);
}
}
@@ -127,23 +122,22 @@ void MixtureofGaussianHMM::ComputeViterbiStateSequence(const Matrix& data_seq, V
MixtureofGaussianHMM::ViterbiInit(transmission_, emis_prob, state_seq);
}
void MixtureofGaussianHMM::TrainBaumWelch(const ArrayList<Matrix>& list_data_seq, int max_iteration, double tolerance) {
void MixtureofGaussianHMM::TrainBaumWelch(const std::vector<Matrix>& list_data_seq, int max_iteration, double tolerance) {
MixtureofGaussianHMM::Train(list_data_seq, &transmission_, &list_mixture_gauss_, max_iteration, tolerance);
}
void MixtureofGaussianHMM::TrainViterbi(const ArrayList<Matrix>& list_data_seq, int max_iteration, double tolerance) {
void MixtureofGaussianHMM::TrainViterbi(const std::vector<Matrix>& list_data_seq, int max_iteration, double tolerance) {
MixtureofGaussianHMM::TrainViterbi(list_data_seq, &transmission_, &list_mixture_gauss_, max_iteration, tolerance);
}
success_t MixtureofGaussianHMM::LoadProfile(const char* profile, Matrix* trans, ArrayList<MixtureGauss>* mixs) {
ArrayList<Matrix> matlst;
success_t MixtureofGaussianHMM::LoadProfile(const char* profile, Matrix* trans, std::vector<MixtureGauss>* mixs) {
std::vector<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;
@@ -153,13 +147,13 @@ success_t MixtureofGaussianHMM::LoadProfile(const char* profile, Matrix* trans,
DEBUG_ASSERT(matlst.size() > p+2*K);
MixtureGauss mix;
mix.InitFromProfile(matlst, p, N);
mixs->PushBackCopy(mix);
mixs->push_back(mix);
p += 2*K+1;
}
return SUCCESS_PASS;
}
success_t MixtureofGaussianHMM::SaveProfile(const char* profile, const Matrix& trans, const ArrayList<MixtureGauss>& mixs) {
success_t MixtureofGaussianHMM::SaveProfile(const char* profile, const Matrix& trans, const std::vector<MixtureGauss>& mixs) {
TextWriter w_pro;
if (!PASSED(w_pro.Open(profile))) {
NONFATAL("Couldn't open '%s' for writing.", profile);
@@ -182,7 +176,7 @@ success_t MixtureofGaussianHMM::SaveProfile(const char* profile, const Matrix& t
return SUCCESS_PASS;
}
void MixtureofGaussianHMM::GenerateInit(int L, const Matrix& trans, const ArrayList<MixtureGauss>& mixs, Matrix* seq, Vector* states){
void MixtureofGaussianHMM::GenerateInit(int L, const Matrix& trans, const std::vector<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;
Matrix& seq_ = *seq;
@@ -222,7 +216,7 @@ void MixtureofGaussianHMM::GenerateInit(int L, const Matrix& trans, const ArrayL
}
}
void MixtureofGaussianHMM::EstimateInit(int numStates, int numClusters, const Matrix& seq, const Vector& states, Matrix* trans, ArrayList<MixtureGauss>* mixs) {
void MixtureofGaussianHMM::EstimateInit(int numStates, int numClusters, const Matrix& seq, const Vector& states, Matrix* trans, std::vector<MixtureGauss>* mixs) {
DEBUG_ASSERT_MSG((seq.n_cols()==states.length()), "hmm_estimateM_init: sequence and states length must be the same");
int N = seq.n_rows(); // emission vector length
@@ -231,11 +225,10 @@ void MixtureofGaussianHMM::EstimateInit(int numStates, int numClusters, const Ma
int K = numClusters;
Matrix &trans_ = *trans;
ArrayList<MixtureGauss>& mix_ = *mixs;
std::vector<MixtureGauss>& mix_ = *mixs;
Vector stateSum;
trans_.Init(M, M);
mix_.Init();
stateSum.Init(M);
trans_.SetZero();
@@ -253,8 +246,7 @@ void MixtureofGaussianHMM::EstimateInit(int numStates, int numClusters, const Ma
trans_.ref(i, j) /= stateSum[i];
}
ArrayList<Matrix> data;
data.Init();
std::vector<Matrix> data;
Vector n_data;
n_data.Init(M); n_data.SetZero();
for (int i = 0; i < L; i++) {
@@ -265,7 +257,7 @@ void MixtureofGaussianHMM::EstimateInit(int numStates, int numClusters, const Ma
Matrix m;
//printf("n[%d]=%8.0f\n", i, n_data[i]);
m.Init(N, (int)n_data[i]);
data.PushBackCopy(m);
data.push_back(m);
}
n_data.SetZero();
for (int i = 0; i < L; i++) {
@@ -275,18 +267,18 @@ void MixtureofGaussianHMM::EstimateInit(int numStates, int numClusters, const Ma
//printf("%d %d %8.0f\n", i, state, n_data[state]);
}
for (int i = 0; i < M; i++) {
ArrayList<int> labels;
ArrayList<Vector> means;
std::vector<int> labels;
std::vector<Vector> means;
kmeans(data[i], K, &labels, &means, 500, 1e-3);
//printf("STATE #%d %d\n", i, K);
MixtureGauss m;
m.Init(K, data[i], labels);
mix_.PushBackCopy(m);
mix_.push_back(m);
}
}
void MixtureofGaussianHMM::EstimateInit(int NumClusters, const Matrix& seq, const Vector& states, Matrix* trans, ArrayList<MixtureGauss>* mixs) {
void MixtureofGaussianHMM::EstimateInit(int NumClusters, const Matrix& seq, const Vector& states, Matrix* trans, std::vector<MixtureGauss>* mixs) {
DEBUG_ASSERT_MSG((seq.n_cols()==states.length()), "hmm_estimateG_init: sequence and states length must be the same");
int M = 0;
for (int i = 0; i < seq.n_cols(); i++)
@@ -323,7 +315,7 @@ double MixtureofGaussianHMM::ViterbiInit(int L, const Matrix& trans, const Matri
return GaussianHMM::ViterbiInit(L, trans, emis_prob, states);
}
void MixtureofGaussianHMM::CalculateEmissionProb(const Matrix& seq, const ArrayList<MixtureGauss>& mixs, Matrix* emis_prob) {
void MixtureofGaussianHMM::CalculateEmissionProb(const Matrix& seq, const std::vector<MixtureGauss>& mixs, Matrix* emis_prob) {
int M = mixs.size();
int L = seq.n_cols();
for (int t = 0; t < L; t++) {
@@ -334,9 +326,9 @@ void MixtureofGaussianHMM::CalculateEmissionProb(const Matrix& seq, const ArrayL
}
}
void MixtureofGaussianHMM::Train(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<MixtureGauss>* guessMG, int max_iter, double tol) {
void MixtureofGaussianHMM::Train(const std::vector<Matrix>& seqs, Matrix* guessTR, std::vector<MixtureGauss>* guessMG, int max_iter, double tol) {
Matrix &gTR = *guessTR;
ArrayList<MixtureGauss>& gMG = *guessMG;
std::vector<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");
@@ -348,7 +340,7 @@ void MixtureofGaussianHMM::Train(const ArrayList<Matrix>& seqs, Matrix* guessTR,
TR.Init(M, M);
Matrix ps, fs, bs, emis_prob; // to hold hmm_decodeG results
ArrayList<Matrix> emis_prob_cluster;
std::vector<Matrix> emis_prob_cluster;
Vector s; // scaling factors
Vector sumState; // the denominator for each state
@@ -358,12 +350,11 @@ void MixtureofGaussianHMM::Train(const ArrayList<Matrix>& seqs, Matrix* guessTR,
s.Init(L);
emis_prob.Init(M, L);
sumState.Init(M);
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.PushBackCopy(m);
emis_prob_cluster.push_back(m);
}
double loglik = 0, oldlog;
@@ -443,9 +434,9 @@ void MixtureofGaussianHMM::Train(const ArrayList<Matrix>& seqs, Matrix* guessTR,
}
}
void MixtureofGaussianHMM::TrainViterbi(const ArrayList<Matrix>& seqs, Matrix* guessTR, ArrayList<MixtureGauss>* guessMG, int max_iter, double tol) {
void MixtureofGaussianHMM::TrainViterbi(const std::vector<Matrix>& seqs, Matrix* guessTR, std::vector<MixtureGauss>* guessMG, int max_iter, double tol) {
Matrix &gTR = *guessTR;
ArrayList<MixtureGauss>& gMG = *guessMG;
std::vector<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");
@@ -457,15 +448,14 @@ void MixtureofGaussianHMM::TrainViterbi(const ArrayList<Matrix>& seqs, Matrix* g
TR.Init(M, M);
Matrix emis_prob; // to hold hmm_decodeG results
ArrayList<Matrix> emis_prob_cluster;
std::vector<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.PushBackCopy(m);
emis_prob_cluster.push_back(m);
}
double loglik = 0, oldlog;
fastlib/branches/fastlib-stl/mlpack/hmm/mixgaussHMM.h
+17 -18
View File
@@ -26,7 +26,7 @@ class MixtureofGaussianHMM {
Matrix transmission_;
/** List of Mixture of Gaussian objects corresponding to each state */
ArrayList<MixtureGauss> list_mixture_gauss_;
std::vector<MixtureGauss> list_mixture_gauss_;
OT_DEF(MixtureofGaussianHMM) {
OT_MY_OBJECT(transmission_);
@@ -35,14 +35,14 @@ class MixtureofGaussianHMM {
public:
/** Getters */
const Matrix& transmission() const { return transmission_; }
const ArrayList<MixtureGauss>& list_mixture_gauss() const { return list_mixture_gauss_; }
const std::vector<MixtureGauss>& list_mixture_gauss() const { return list_mixture_gauss_; }
/** Setter used when already initialized */
void setModel(const Matrix& transmission,
const ArrayList<MixtureGauss>& list_mixture_gauss);
const std::vector<MixtureGauss>& list_mixture_gauss);
/** Initializes from computed transmission and Mixture of Gaussian parameters */
void Init(const Matrix& transmission, const ArrayList<MixtureGauss>& list_mixture_gauss);
void Init(const Matrix& transmission, const std::vector<MixtureGauss>& list_mixture_gauss);
/** Initializes by loading from a file */
void InitFromFile(const char* profile);
@@ -50,7 +50,6 @@ class MixtureofGaussianHMM {
/** Initializes empty object */
void Init() {
transmission_.Init(0, 0);
list_mixture_gauss_.Init();
}
/** Load from file, used when already initialized */
@@ -87,7 +86,7 @@ class MixtureofGaussianHMM {
double ComputeLogLikelihood(const Matrix& data_seq) const;
/** Compute the log-likelihood of a list of sequences */
void ComputeLogLikelihood(const ArrayList<Matrix>& list_data_seq, ArrayList<double>* list_likelihood) const;
void ComputeLogLikelihood(const std::vector<Matrix>& list_data_seq, std::vector<double>* list_likelihood) const;
/** Compute the most probable sequence (Viterbi) */
void ComputeViterbiStateSequence(const Matrix& data_seq, Vector* state_seq) const;
@@ -96,18 +95,18 @@ class MixtureofGaussianHMM {
* Train the model with a list of sequences, must be already initialized
* using Baum-Welch EM algorithm
*/
void TrainBaumWelch(const ArrayList<Matrix>& list_data_seq, int max_iteration, double tolerance);
void TrainBaumWelch(const std::vector<Matrix>& list_data_seq, int max_iteration, double tolerance);
/**
* 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<Matrix>& list_data_seq, int max_iteration, double tolerance);
void TrainViterbi(const std::vector<Matrix>& list_data_seq, int max_iteration, double tolerance);
////////// Static helper functions ///////////////////////////////////////
static success_t LoadProfile(const char* profile, Matrix* trans, ArrayList<MixtureGauss>* mixs);
static success_t SaveProfile(const char* profile, const Matrix& trans, const ArrayList<MixtureGauss>& mixs);
static success_t LoadProfile(const char* profile, Matrix* trans, std::vector<MixtureGauss>* mixs);
static success_t SaveProfile(const char* profile, const Matrix& trans, const std::vector<MixtureGauss>& mixs);
/**
* Generating a sequence and states using transition and emission probabilities.
@@ -118,13 +117,13 @@ class MixtureofGaussianHMM {
* seq: generated sequence, uninitialized matrix, will have size N x L
* states: generated states, uninitialized vector, will have length L
*/
static void GenerateInit(int L, const Matrix& trans, const ArrayList<MixtureGauss>& mixs, Matrix* seq, Vector* states);
static void GenerateInit(int L, const Matrix& trans, const std::vector<MixtureGauss>& mixs, Matrix* seq, Vector* states);
/** Estimate transition and emission distribution from sequence and states */
static void EstimateInit(int NumClusters, const Matrix& seq, const Vector& states,
Matrix* trans, ArrayList<MixtureGauss>* mixs);
Matrix* trans, std::vector<MixtureGauss>* mixs);
static void EstimateInit(int numStates, int NumClusters, const Matrix& seq,
const Vector& states, Matrix* trans, ArrayList<MixtureGauss>* mixs);
const Vector& states, Matrix* trans, std::vector<MixtureGauss>* mixs);
/**
* Calculate posteriori probabilities of states at each steps
@@ -147,7 +146,7 @@ class MixtureofGaussianHMM {
static double Decode(int L, const Matrix& trans, const Matrix& emis_prob,
Matrix* pstates, Matrix* fs, Matrix* bs, Vector* scales);
static void CalculateEmissionProb(const Matrix& seq, const ArrayList<MixtureGauss>& mixs, Matrix* emis_prob);
static void CalculateEmissionProb(const Matrix& seq, const std::vector<MixtureGauss>& mixs, Matrix* emis_prob);
/**
* Calculate the most probable states for a sequence
@@ -164,9 +163,9 @@ class MixtureofGaussianHMM {
* Baum-Welch and Viterbi estimation of transition and emission
* distribution (Gaussian)
*/
static void Train(const ArrayList<Matrix>& seqs, Matrix* guessTR,
ArrayList<MixtureGauss>* guessMG, int max_iter, double tol);
static void TrainViterbi(const ArrayList<Matrix>& seqs, Matrix* guessTR,
ArrayList<MixtureGauss>* guessMG, int max_iter, double tol);
static void Train(const std::vector<Matrix>& seqs, Matrix* guessTR,
std::vector<MixtureGauss>* guessMG, int max_iter, double tol);
static void TrainViterbi(const std::vector<Matrix>& seqs, Matrix* guessTR,
std::vector<MixtureGauss>* guessMG, int max_iter, double tol);
};
#endif
fastlib/branches/fastlib-stl/mlpack/hmm/mixtureDST.cc
+21 -28
View File
@@ -11,28 +11,26 @@
using namespace hmm_support;
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.PushBackCopy(v);
means.push_back(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.PushBackCopy(m);
covs.push_back(m);
}
prior.Init(means.size());
for (int i = 0; i < prior.length(); i++) prior[i] = 1.0/K;
ACC_means.InitCopy(means);
ACC_covs.InitCopy(covs);
ACC_means.assign(means.begin(), means.end());
ACC_covs.assign(covs.begin(), covs.end());
ACC_prior.Init(K);
inv_covs.InitCopy(covs);
inv_covs.assign(covs.begin(), covs.end());
det_covs.Init(covs.size());
for (int i = 0; i < K; i++) {
double det = la::Determinant(covs[i]);
@@ -41,28 +39,26 @@ void MixtureGauss::Init(int K, int N) {
}
}
void MixtureGauss::Init(int K, const Matrix& data, const ArrayList<int>& labels) {
means.Init();
void MixtureGauss::Init(int K, const Matrix& data, const std::vector<int>& labels) {
int N = data.n_rows();
for (int i = 0; i < K; i++) {
Vector v;
v.Init(N);
means.PushBackCopy(v);
means.push_back(v);
}
covs.Init();
for (int i = 0; i < means.size(); i++) {
Matrix m;
m.Init(N, N);
covs.PushBackCopy(m);
covs.push_back(m);
}
prior.Init(means.size());
ACC_means.InitCopy(means);
ACC_covs.InitCopy(covs);
ACC_means.assign(means.begin(), means.end());
ACC_covs.assign(covs.begin(), covs.end());
ACC_prior.Init(K);
inv_covs.InitCopy(covs);
inv_covs.assign(covs.begin(), covs.end());
det_covs.Init(covs.size());
start_accumulate();
//printf("cols = %d rows = %d\n", data.n_cols(), data.n_rows());
@@ -88,12 +84,11 @@ void MixtureGauss::InitFromFile(const char* mean_fn, const char* covs_fn, const
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.PushBackCopy(m);
covs.push_back(m);
}
}
if (prior_fn != NULL) {
@@ -108,10 +103,10 @@ void MixtureGauss::InitFromFile(const char* mean_fn, const char* covs_fn, const
for (int i = 0; i < prior.length(); i++) prior[i] = 1.0/K;
}
ACC_means.InitCopy(means);
ACC_covs.InitCopy(covs);
ACC_means.assign(means.begin(), means.end());
ACC_covs.assign(covs.begin(), covs.end());
ACC_prior.Init(K);
inv_covs.InitCopy(covs);
inv_covs.assign(covs.begin(), covs.end());
det_covs.Init(covs.size());
for (int i = 0; i < K; i++) {
double det = la::Determinant(covs[i]);
@@ -120,7 +115,7 @@ void MixtureGauss::InitFromFile(const char* mean_fn, const char* covs_fn, const
}
}
void MixtureGauss::InitFromProfile(const ArrayList<Matrix>& matlst, int start, int N) {
void MixtureGauss::InitFromProfile(const std::vector<Matrix>& matlst, int start, int N) {
DEBUG_ASSERT(matlst[start].n_cols()==1);
Vector tmp;
matlst[start].MakeColumnVector(0, &tmp);
@@ -128,21 +123,19 @@ void MixtureGauss::InitFromProfile(const ArrayList<Matrix>& matlst, int start, i
// DEBUG: print_vector(prior, " prior = ");
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.PushBackCopy(m);
covs.PushBackCopy(matlst[i+1]);
means.push_back(m);
covs.push_back(matlst[i+1]);
}
ACC_means.InitCopy(means);
ACC_covs.InitCopy(covs);
ACC_means.assign(means.begin(),means.end());
ACC_covs.assign(covs.begin(),covs.end());
ACC_prior.Init(K);
inv_covs.InitCopy(covs);
inv_covs.assign(covs.begin(),covs.end());
det_covs.Init(covs.size());
for (int i = 0; i < K; i++) {
double det = la::Determinant(covs[i]);
fastlib/branches/fastlib-stl/mlpack/hmm/mixtureDST.h
+7 -7
View File
@@ -18,25 +18,25 @@ class MixtureGauss {
////////////// Member variables //////////////////////////////////////
private:
/** List of means of clusters */
ArrayList<Vector> means;
std::vector<Vector> means;
/** List of covariance matrices of clusters */
ArrayList<Matrix> covs;
std::vector<Matrix> covs;
/** Prior probabilities of the clusters */
Vector prior;
/** Inverse of covariance matrices */
ArrayList<Matrix> inv_covs;
std::vector<Matrix> inv_covs;
/** Vector of constant in normal density formula */
Vector det_covs;
/** Accumulating means */
ArrayList<Vector> ACC_means;
std::vector<Vector> ACC_means;
/** Accumulating covariance */
ArrayList<Matrix> ACC_covs;
std::vector<Matrix> ACC_covs;
/** Accumulating prior probability */
Vector ACC_prior;
@@ -64,13 +64,13 @@ class MixtureGauss {
* start with the prior vector, follows by the mean and covariance
* of each cluster.
*/
void InitFromProfile(const ArrayList<Matrix>& matlst, int start, int N);
void InitFromProfile(const std::vector<Matrix>& matlst, int start, int N);
/** Init with K clusters and dimension N */
void Init(int K, int N);
/** Init with K clusters and the data with label */
void Init(int K, const Matrix& data, const ArrayList<int>& labels);
void Init(int K, const Matrix& data, const std::vector<int>& labels);
/** Print the mixture to stdout for debugging */
void print_mixture(const char* s) const;
fastlib/branches/fastlib-stl/mlpack/hmm/support.cc
+66 -57
View File
@@ -96,19 +96,19 @@ namespace hmm_support {
return det_cov * exp(-0.5*MyMulExpert(d, inv_cov, d));
}
bool kmeans(const ArrayList<Matrix>& data, int num_clusters,
ArrayList<int> *labels_, ArrayList<Vector> *centroids_,
bool kmeans(const std::vector<Matrix>& data, int num_clusters,
std::vector<int> *labels_, std::vector<Vector> *centroids_,
int max_iter, double error_thresh)
{
ArrayList<int> counts; //number of points in each cluster
ArrayList<Vector> tmp_centroids;
std::vector<int> counts; //number of points in each cluster
std::vector<Vector> tmp_centroids;
int num_points, num_dims;
int i, j, num_iter=0;
double error, old_error;
//Assign pointers to references to avoid repeated dereferencing.
ArrayList<int> &labels = *labels_;
ArrayList<Vector> &centroids = *centroids_;
std::vector<int> &labels = *labels_;
std::vector<Vector> &centroids = *centroids_;
num_points = 0;
for (int i = 0; i < data.size(); i++)
@@ -118,11 +118,11 @@ namespace hmm_support {
num_dims = data[0].n_rows();
centroids.Init(num_clusters);
tmp_centroids.Init(num_clusters);
centroids.reserve(num_clusters);
tmp_centroids.reserve(num_clusters);
counts.Init(num_clusters);
labels.Init(num_points);
counts.reserve(num_clusters);
labels.reserve(num_points);
//Initialize the clusters to k points
for (j=0; j < num_clusters; j++) {
@@ -183,18 +183,18 @@ namespace hmm_support {
}
bool kmeans(Matrix const &data, int num_clusters,
ArrayList<int> *labels_, ArrayList<Vector> *centroids_,
std::vector<int> *labels_, std::vector<Vector> *centroids_,
int max_iter, double error_thresh)
{
ArrayList<int> counts; //number of points in each cluster
ArrayList<Vector> tmp_centroids;
std::vector<int> counts; //number of points in each cluster
std::vector<Vector> tmp_centroids;
int num_points, num_dims;
int i, j, num_iter=0;
double error, old_error;
//Assign pointers to references to avoid repeated dereferencing.
ArrayList<int> &labels = *labels_;
ArrayList<Vector> &centroids = *centroids_;
std::vector<int> &labels = *labels_;
std::vector<Vector> &centroids = *centroids_;
if (data.n_cols() < num_clusters)
return false;
@@ -202,11 +202,11 @@ namespace hmm_support {
num_points = data.n_cols();
num_dims = data.n_rows();
centroids.Init(num_clusters);
tmp_centroids.Init(num_clusters);
centroids.reserve(num_clusters);
tmp_centroids.reserve(num_clusters);
counts.Init(num_clusters);
labels.Init(num_points);
counts.reserve(num_clusters);
labels.reserve(num_points);
//Initialize the clusters to k points
for (i=0, j=0; j < num_clusters; i+=num_points/num_clusters, j++) {
@@ -263,32 +263,31 @@ namespace hmm_support {
}
void mat2arrlst(Matrix& a, ArrayList<Vector> * seqs) {
void mat2arrlst(Matrix& a, std::vector<Vector> * seqs) {
int n = a.n_cols();
ArrayList<Vector> & s_ = *seqs;
s_.Init();
std::vector<Vector> & s_ = *seqs;
for (int i = 0; i < n; i++) {
Vector seq;
a.MakeColumnVector(i, &seq);
s_.PushBackCopy(seq);
s_.push_back(seq);
}
}
void mat2arrlstmat(int N, Matrix& a, ArrayList<Matrix> * seqs) {
void mat2arrlstmat(int N, Matrix& a, std::vector<Matrix> * seqs) {
int n = a.n_cols();
ArrayList<Matrix>& s_ = *seqs;
s_.Init();
std::vector<Matrix>& s_ = *seqs;
for (int i = 0; i < n; i+=N) {
Matrix b;
a.MakeColumnSlice(i, N, &b);
s_.PushBackCopy(b);
s_.push_back(b);
}
}
bool skip_blank(TextLineReader& reader) {
for (;;){
if (!reader.MoreLines()) return false;
char* pos = reader.Peek().begin();
// char* pos = reader.Peek().begin();
std::string::iterator pos = reader.Peek().begin();
while (*pos == ' ' || *pos == ',' || *pos == '\t')
pos++;
if (*pos == '\0' || *pos == '%') reader.Gobble();
@@ -307,25 +306,29 @@ namespace hmm_support {
int n_cols = 0;
bool is_done;
{// How many columns ?
ArrayList<String> num_str;
num_str.Init();
reader.Peek().Split(", \t", &num_str);
std::vector<std::string> num_str;
tokenizeString(reader.Peek(), ", \t", num_str);
n_cols = num_str.size();
}
ArrayList<double> num_double;
num_double.Init();
std::vector<double> num_double;
for(;;) { // read each rows
n_rows++;
// deprecated: double* point = num_double.AddBack(n_cols);
ArrayList<String> num_str;
num_str.Init();
reader.Peek().Split(", \t", &num_str);
std::vector<std::string> num_str;
tokenizeString(reader.Peek(), ", \t", num_str);
// reader.Peek().Split(", \t", &num_str);
DEBUG_ASSERT(num_str.size() == n_cols);
for (int i = 0; i < n_cols; i++)
num_double.PushBackCopy(strtod(num_str[i], NULL));
std::istringstream is;
for (int i = 0; i < n_cols; i++) {
double d;
is.str(num_str[i]);
if( !(is >> d ) )
abort();
num_double.push_back(d);
}
is_done = false;
@@ -336,7 +339,7 @@ namespace hmm_support {
is_done = true;
break;
}
char* pos = reader.Peek().begin();
std::string::iterator pos = reader.Peek().begin();
while (*pos == ' ' || *pos == '\t')
pos++;
if (*pos == '\0') reader.Gobble();
@@ -348,8 +351,12 @@ namespace hmm_support {
}
if (is_done) {
num_double.Trim();
matrix->Own(num_double.ReleasePtr(), n_cols, n_rows);
{
std::vector<double> empty;
num_double.swap(empty);
}
//FIXME
// matrix->Own(num_double.ReleasePtr(), n_cols, n_rows);
return SUCCESS_PASS;
}
}
@@ -362,20 +369,23 @@ namespace hmm_support {
return SUCCESS_FAIL;
}
else {
ArrayList<double> num_double;
num_double.Init();
std::vector<double> num_double;
for(;;) { // read each rows
bool is_done = false;
ArrayList<String> num_str;
num_str.Init();
reader.Peek().Split(", \t", &num_str);
std::vector<std::string> num_str;
tokenizeString(reader.Peek(), ", \t", num_str);
// deprecated: double* point = num_double.AddBack(num_str.size());
for (int i = 0; i < num_str.size(); i++)
num_double.PushBackCopy(strtod(num_str[i], NULL));
std::istringstream is;
for (int i = 0; i < num_str.size(); i++) {
double d;
is.str(num_str[i]);
if( !(is >> d) )
abort();
}
reader.Gobble();
@@ -384,7 +394,7 @@ namespace hmm_support {
is_done = true;
break;
}
char* pos = reader.Peek().begin();
std::string::iterator pos = reader.Peek().begin();
while (*pos == ' ' || *pos == '\t')
pos++;
if (*pos == '\0') reader.Gobble();
@@ -396,37 +406,36 @@ namespace hmm_support {
}
if (is_done) {
num_double.Trim();
//num_double.Trim();
int length = num_double.size();
vec->Own(num_double.ReleasePtr(), length);
// FIXME
// vec->Own(num_double.ReleasePtr(), length);
return SUCCESS_PASS;
}
}
}
}
success_t load_matrix_list(const char* filename, ArrayList<Matrix> *matlst) {
success_t load_matrix_list(const char* filename, std::vector<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->PushBackCopy(tmp);
matlst->push_back(tmp);
}
else break;
} while (1);
return SUCCESS_PASS;
}
success_t load_vector_list(const char* filename, ArrayList<Vector> *veclst) {
success_t load_vector_list(const char* filename, std::vector<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->PushBackCopy(vec);
veclst->push_back(vec);
}
else break;
} while (1);
fastlib/branches/fastlib-stl/mlpack/hmm/support.h
+7 -7
View File
@@ -38,31 +38,31 @@ namespace hmm_support {
void print_vector(TextWriter& writer, const Vector& a, const char* msg, const char* format = "%f,");
/** Compute the centroids and label the samples by K-means algorithm */
bool kmeans(const ArrayList<Matrix>& data, int num_clusters,
ArrayList<int> *labels_, ArrayList<Vector> *centroids_,
bool kmeans(const std::vector<Matrix>& data, int num_clusters,
std::vector<int> *labels_, std::vector<Vector> *centroids_,
int max_iter = 1000, double error_thresh = 1e-3);
bool kmeans(Matrix const &data, int num_clusters,
ArrayList<int> *labels_, ArrayList<Vector> *centroids_,
std::vector<int> *labels_, std::vector<Vector> *centroids_,
int max_iter=1000, double error_thresh=1e-04);
/** Convert a matrix in to an array list of vectors of its column */
void mat2arrlst(Matrix& a, ArrayList<Vector> * seqs);
void mat2arrlst(Matrix& a, std::vector<Vector> * seqs);
/** Convert a matrix in to an array list of matrices of slice of its columns */
void mat2arrlstmat(int N, Matrix& a, ArrayList<Matrix> * seqs);
void mat2arrlstmat(int N, Matrix& a, std::vector<Matrix> * seqs);
/**
* Load an array list of matrices from file where the matrices
* are seperated by a line start with %
*/
success_t load_matrix_list(const char* filename, ArrayList<Matrix> *matlst);
success_t load_matrix_list(const char* filename, std::vector<Matrix> *matlst);
/**
* Load an array list of vectors from file where the vectors
* are seperated by a line start with %
*/
success_t load_vector_list(const char* filename, ArrayList<Vector> *veclst);
success_t load_vector_list(const char* filename, std::vector<Vector> *veclst);
};
#endif
fastlib/branches/fastlib-stl/mlpack/quicsvd/cosine_tree.cc
+21 -21
View File
@@ -16,8 +16,8 @@ bool isZero(double d) {
// Init a root cosine node from a matrix
CosineNode::CosineNode(const Matrix& A) {
this->A_.Alias(A);
origIndices_.Init(A_.n_cols());
norms_.Init(n_cols());
origIndices_.reserve(A_.n_cols());
norms_.reserve(n_cols());
for (int i_col = 0; i_col < origIndices_.size(); i_col++) {
origIndices_[i_col] = i_col;
norms_[i_col] = columnNormL2(A_, i_col);
@@ -32,10 +32,10 @@ CosineNode::CosineNode(const Matrix& A) {
// Init a child cosine node from its parent and a set of the parent's columns
CosineNode::CosineNode(CosineNode& parent,
const ArrayList<int>& indices, bool isLeft) {
const std::vector<int>& indices, bool isLeft) {
A_.Alias(parent.A_);
origIndices_.Init(indices.size());
norms_.Init(n_cols());
origIndices_.reserve(indices.size());
norms_.reserve(n_cols());
for (int i_col = 0; i_col < origIndices_.size(); i_col++) {
origIndices_[i_col] = parent.origIndices_[indices[i_col]];
norms_[i_col] = parent.norms_[indices[i_col]];
@@ -53,7 +53,7 @@ CosineNode::CosineNode(CosineNode& parent,
void CosineNode::CalStats() {
// Calculate cummlulative sum square of L2 norms
cum_norms_.Init(origIndices_.size());
cum_norms_.reserve(origIndices_.size());
for (int i_col = 0; i_col < origIndices_.size(); i_col++)
cum_norms_[i_col] = ((i_col > 0) ? cum_norms_[i_col-1]:0)
+ math::Sqr(norms_[i_col]);
@@ -77,8 +77,8 @@ void CosineNode::ChooseCenter(Vector* center) {
}
void CosineNode::CalCosines(const Vector& center,
ArrayList<double>* cosines) {
cosines->Init(n_cols());
std::vector<double>* cosines) {
cosines->reserve(n_cols());
double centerL2 = la::LengthEuclidean(center);
for (index_t i_col = 0; i_col < n_cols(); i_col++)
// if col is a zero vector then push it to the left node
@@ -93,14 +93,14 @@ void CosineNode::CalCosines(const Vector& center,
}
}
void CosineNode::CreateIndices(ArrayList<int>* indices) {
indices->Init(n_cols());
void CosineNode::CreateIndices(std::vector<int>* indices) {
indices->reserve(n_cols());
for (index_t i_col = 0; i_col < n_cols(); i_col++)
(*indices)[i_col] = i_col;
}
// Quicksort partitioning procedure
index_t qpartition(ArrayList<double>& key, ArrayList<int>& data,
index_t qpartition(std::vector<double>& key, std::vector<int>& data,
index_t left, index_t right) {
index_t j = left;
double x = key[left];
@@ -117,7 +117,7 @@ index_t qpartition(ArrayList<double>& key, ArrayList<int>& data,
// Quicksort on the cosine values
void qsort(ArrayList<double>& key, ArrayList<int>& data,
void qsort(std::vector<double>& key, std::vector<int>& data,
index_t left, index_t right) {
if (left >= right) return;
index_t middle = qpartition(key, data, left, right);
@@ -125,13 +125,13 @@ void qsort(ArrayList<double>& key, ArrayList<int>& data,
qsort(key, data, middle+1, right);
}
void Sort(ArrayList<double>& key, ArrayList<int>& data) {
void Sort(std::vector<double>& key, std::vector<int>& data) {
qsort(key, data, 0, key.size()-1);
}
// Calculate the split point where the cosines values are closer
// to the minimum cosine value than the maximum cosine value
index_t calSplitPoint(const ArrayList<double>& key) {
index_t calSplitPoint(const std::vector<double>& key) {
double leftKey = key[0];
double rightKey = key[key.size()-1];
index_t i = 0;
@@ -141,16 +141,16 @@ index_t calSplitPoint(const ArrayList<double>& key) {
}
// Init a subcopy of an array list
void InitSubCopy(const ArrayList<int>& src, index_t pos, index_t size,
ArrayList<int>* dst) {
dst->Init(size);
void InitSubCopy(const std::vector<int>& src, index_t pos, index_t size,
std::vector<int>* dst) {
dst->reserve(size);
for (index_t i = 0; i < size; i++)
(*dst)[i] = src[pos+i];
}
// Split the indices at the split point
void splitIndices(ArrayList<int>& indices, int leftSize,
ArrayList<int>* leftIdx, ArrayList<int>* rightIdx) {
void splitIndices(std::vector<int>& indices, int leftSize,
std::vector<int>* leftIdx, std::vector<int>* rightIdx) {
InitSubCopy(indices, 0, leftSize, leftIdx);
InitSubCopy(indices, leftSize, indices.size()-leftSize, rightIdx);
}
@@ -163,8 +163,8 @@ void splitIndices(ArrayList<int>& indices, int leftSize,
void CosineNode::Split() {
if (n_cols() < 2) return;
Vector center;
ArrayList<double> cosines;
ArrayList<int> indices, leftIdx, rightIdx;
std::vector<double> cosines;
std::vector<int> indices, leftIdx, rightIdx;
ChooseCenter(&center);
//ot::Print(center, "center", stdout);
fastlib/branches/fastlib-stl/mlpack/quicsvd/cosine_tree.h
+6 -6
View File
@@ -20,13 +20,13 @@ class CosineNode {
Matrix A_;
/** Indices of columns of matrix A in this node */
ArrayList<int> origIndices_;
std::vector<int> origIndices_;
/** L2 norms of the columns in this node */
ArrayList<double> norms_;
std::vector<double> norms_;
/** Cummulative sum of L2 norm squares to use in column sampling */
ArrayList<double> cum_norms_;
std::vector<double> cum_norms_;
/** Mean vector to be added to the basis when this node is chosen */
Vector mean_;
@@ -48,7 +48,7 @@ class CosineNode {
CosineNode(const Matrix& A);
/** Constructor of the childen */
CosineNode(CosineNode& parent, const ArrayList<int>& indices,
CosineNode(CosineNode& parent, const std::vector<int>& indices,
bool isLeft);
/** Get a column in this node, we have to use origIndices to
@@ -122,10 +122,10 @@ class CosineNode {
void ChooseCenter(Vector* center);
/** Calculate cosine values of all column with respect to the center */
void CalCosines(const Vector& center, ArrayList<double>* cosines);
void CalCosines(const Vector& center, std::vector<double>* cosines);
/** Create an array list of indices 0..n_cols()-1 */
void CreateIndices(ArrayList<int>* indices);
void CreateIndices(std::vector<int>* indices);
/** Friend unit test class */
friend class CosineNodeTest;
fastlib/branches/fastlib-stl/mlpack/quicsvd/quicsvd.cc
+7 -9
View File
@@ -11,9 +11,7 @@
QuicSVD::QuicSVD(const Matrix& A, double targetRelErr) : root_(A) {
A_.Alias(A);
basis_.Init(); // init empty basis
UTA_.Init(); // and empty projection, too
projMagSq_.InitRepeat(0, n_cols()); // projected magnitude squares are zero's
projMagSq_.resize(n_cols(), 0); // projected magnitude squares are zero's
sumProjMagSq_ = 0;
targetRelErr_ = targetRelErr;
dataNorm2_ = root_.getSumL2(); // keep the Frobenius norm of A
@@ -22,7 +20,7 @@ QuicSVD::QuicSVD(const Matrix& A, double targetRelErr) : root_(A) {
// Modified Gram-Schmidt method, calculate the orthogonalized
// new basis vector
bool MGS(const ArrayList<Vector>& basis, const Vector& newVec,
bool MGS(const std::vector<Vector>& basis, const Vector& newVec,
Vector* newBasisVec) {
//ot::Print(basis, "basis", stdout);
//ot::Print(newVec, "newVec", stdout);
@@ -49,9 +47,9 @@ void QuicSVD::addBasisFrom(const CosineNode& node) {
// check if new vector are independent and orthogonalize it
if (MGS(basis_, node.getMean(), &nodeBasis)) {
Vector av;
basis_.PushBackCopy(nodeBasis); // add to the basis
basis_.push_back(nodeBasis); // add to the basis
la::MulInit(nodeBasis, A_, &av); // calculate projection of A
UTA_.PushBackCopy(av); // on the new basis vector and save
UTA_.push_back(av); // on the new basis vector and save
for (index_t i_col = 0; i_col < n_cols(); i_col++) {
double magSq = math::Sqr(av[i_col]); // magnitude square of the i-th
projMagSq_[i_col] += magSq; // column of A in the new subspace
@@ -120,7 +118,7 @@ void QuicSVD::ComputeSVD(Vector* s, Matrix* U, Matrix* VT) {
// squares of singular values of UTA.
// SVD on UTA2 is more efficient as it is a square matrix
// with smaller dimension
void createUTA2(const ArrayList<Vector>& UTA, Matrix* UTA2) {
void createUTA2(const std::vector<Vector>& UTA, Matrix* UTA2) {
UTA2->Init(UTA.size(), UTA.size());
for (int i = 0; i < UTA.size(); i++)
for (int j = 0; j < UTA.size(); j++)
@@ -129,7 +127,7 @@ void createUTA2(const ArrayList<Vector>& UTA, Matrix* UTA2) {
// Matrix multiplication of a list of vector and a matrix
// C = AB
void MulInit(const ArrayList<Vector>& A, const Matrix& B,
void MulInit(const std::vector<Vector>& A, const Matrix& B,
Matrix* C) {
index_t m = A[0].length();
index_t n = A.size();
@@ -144,7 +142,7 @@ void MulInit(const ArrayList<Vector>& A, const Matrix& B,
}
// Matrix multiplication C = B' A'
void MulTransCInit(const ArrayList<Vector>& A, const Matrix& B,
void MulTransCInit(const std::vector<Vector>& A, const Matrix& B,
Matrix* C) {
index_t m = A[0].length();
index_t n = A.size();
fastlib/branches/fastlib-stl/mlpack/quicsvd/quicsvd.h
+3 -3
View File
@@ -41,13 +41,13 @@ class QuicSVD {
/** Orthonomal basis of the subspace */
ArrayList<Vector> basis_;
std::vector<Vector> basis_;
/** Projection coordinates of columns of A_ onto the subspace */
ArrayList<Vector> UTA_;
std::vector<Vector> UTA_;
/** Projection magnitude of columns of A_ in the subspace */
ArrayList<double> projMagSq_;
std::vector<double> projMagSq_;
/** Frobenius norm of A_ being projected to the subspace */
double sumProjMagSq_;
fastlib/branches/fastlib-stl/mlpack/series_expansion/series_expansion_aux.cc
+16 -16
View File
@@ -8,23 +8,23 @@ int SeriesExpansionAux::get_max_total_num_coeffs() const {
return list_total_num_coeffs_[max_order_];
}
const ArrayList < short int > *SeriesExpansionAux::get_lower_mapping_index()
const std::vector < short int > &SeriesExpansionAux::get_lower_mapping_index()
const {
return lower_mapping_index_.begin();
return lower_mapping_index_.front();
}
int SeriesExpansionAux::get_max_order() const {
return max_order_;
}
const ArrayList < short int > &SeriesExpansionAux::get_multiindex(int pos)
const std::vector < short int > &SeriesExpansionAux::get_multiindex(int pos)
const {
return multiindex_mapping_[pos];
}
const ArrayList < short int > *SeriesExpansionAux::get_multiindex_mapping()
const std::vector < short int > &SeriesExpansionAux::get_multiindex_mapping()
const {
return multiindex_mapping_.begin();
return multiindex_mapping_.front();
}
const Vector& SeriesExpansionAux::get_neg_inv_multiindex_factorials() const {
@@ -44,14 +44,14 @@ int SeriesExpansionAux::get_total_num_coeffs(int order) const {
return list_total_num_coeffs_[order];
}
const ArrayList < short int > *SeriesExpansionAux::get_upper_mapping_index()
const std::vector < short int > &SeriesExpansionAux::get_upper_mapping_index()
const {
return upper_mapping_index_.begin();
return upper_mapping_index_.front();
}
int SeriesExpansionAux::ComputeMultiindexPosition
(const ArrayList<short int> &multiindex) const {
(const std::vector<short int> &multiindex) const {
int dim = multiindex.size();
int mapping_sum = 0;
@@ -92,8 +92,8 @@ double SeriesExpansionAux::DirectLocalAccumulationCost(int order) const {
void SeriesExpansionAux::Init(int max_order, int dim) {
int p, k, t, tail, i, j;
ArrayList<int> heads;
ArrayList<int> cinds;
std::vector<int> heads;
std::vector<int> cinds;
// initialize max order and dimension
dim_ = dim;
@@ -101,7 +101,7 @@ void SeriesExpansionAux::Init(int max_order, int dim) {
// compute the list of total number of coefficients for p-th order expansion
int limit = 2 * max_order + 1;
list_total_num_coeffs_.Init(limit);
list_total_num_coeffs_.reserve(limit);
for(p = 0; p < limit; p++) {
list_total_num_coeffs_[p] = (int) math::BinomialCoefficient(p + dim, dim);
}
@@ -114,8 +114,8 @@ void SeriesExpansionAux::Init(int max_order, int dim) {
// and multiindex_combination precomputed factors
inv_multiindex_factorials_.Init(list_total_num_coeffs_[limit - 1]);
neg_inv_multiindex_factorials_.Init(list_total_num_coeffs_[limit - 1]);
multiindex_mapping_.Init(list_total_num_coeffs_[limit - 1]);
(multiindex_mapping_[0]).Init(dim_);
multiindex_mapping_.reserve(list_total_num_coeffs_[limit - 1]);
(multiindex_mapping_[0]).reserve(dim_);
for(j = 0; j < dim; j++) {
(multiindex_mapping_[0])[j] = 0;
}
@@ -123,8 +123,8 @@ void SeriesExpansionAux::Init(int max_order, int dim) {
n_choose_k_.SetZero();
// initialization of temporary variables for computation...
heads.Init(dim + 1);
cinds.Init(list_total_num_coeffs_[limit - 1]);
heads.reserve(dim + 1);
cinds.reserve(list_total_num_coeffs_[limit - 1]);
for(i = 0; i < dim; i++) {
heads[i] = 0;
@@ -147,7 +147,7 @@ void SeriesExpansionAux::Init(int max_order, int dim) {
neg_inv_multiindex_factorials_[t] =
-neg_inv_multiindex_factorials_[j] / cinds[t];
(multiindex_mapping_[t]).InitCopy(multiindex_mapping_[j]);
(multiindex_mapping_[t]).assign(multiindex_mapping_[j].begin(), multiindex_mapping_[j].end());
(multiindex_mapping_[t])[i] = (multiindex_mapping_[t])[i] + 1;
}
}
fastlib/branches/fastlib-stl/mlpack/series_expansion/series_expansion_aux.h
+23 -25
View File
@@ -19,7 +19,7 @@ class SeriesExpansionAux {
Vector factorials_;
ArrayList<int> list_total_num_coeffs_;
std::vector<int> list_total_num_coeffs_;
Vector inv_multiindex_factorials_;
@@ -27,21 +27,21 @@ class SeriesExpansionAux {
Matrix multiindex_combination_;
ArrayList< ArrayList<short int> > multiindex_mapping_;
std::vector< std::vector<short int> > multiindex_mapping_;
/**
* for each i-th multiindex m_i, store the positions of the j-th
* multiindex mapping such that m_i - m_j >= 0 (the difference in
* all coordinates is nonnegative).
*/
ArrayList< ArrayList<short int> > lower_mapping_index_;
std::vector< std::vector<short int> > lower_mapping_index_;
/**
* for each i-th multiindex m_i, store the positions of the j-th
* multiindex mapping such that m_i - m_j <= 0 (the difference in
* all coordinates is nonpositive).
*/
ArrayList< ArrayList<short int> > upper_mapping_index_;
std::vector< std::vector<short int> > upper_mapping_index_;
/** row index is for n, column index is for k */
Matrix n_choose_k_;
@@ -73,20 +73,19 @@ class SeriesExpansionAux {
void ComputeLowerMappingIndex() {
ArrayList<int> diff;
diff.Init(dim_);
std::vector<int> diff;
diff.reserve(dim_);
int limit = 2 * max_order_;
// initialize the index
lower_mapping_index_.Init(list_total_num_coeffs_[limit]);
lower_mapping_index_.reserve(list_total_num_coeffs_[limit]);
for(index_t i = 0; i < list_total_num_coeffs_[limit]; i++) {
const ArrayList<short int> &outer_mapping = multiindex_mapping_[i];
lower_mapping_index_[i].Init();
const std::vector<short int> &outer_mapping = multiindex_mapping_[i];
for(index_t j = 0; j < list_total_num_coeffs_[limit]; j++) {
const ArrayList<short int> &inner_mapping = multiindex_mapping_[j];
const std::vector<short int> &inner_mapping = multiindex_mapping_[j];
int flag = 0;
for(index_t d = 0; d < dim_; d++) {
@@ -99,7 +98,7 @@ class SeriesExpansionAux {
}
if(flag == 0) {
(lower_mapping_index_[i]).PushBackCopy(j);
(lower_mapping_index_[i]).push_back(j);
}
} // end of j-loop
} // end of i-loop
@@ -114,12 +113,12 @@ class SeriesExpansionAux {
for(index_t j = 0; j < list_total_num_coeffs_[limit]; j++) {
// beta mapping
const ArrayList<short int> &beta_mapping = multiindex_mapping_[j];
const std::vector<short int> &beta_mapping = multiindex_mapping_[j];
for(index_t k = 0; k < list_total_num_coeffs_[limit]; k++) {
// alpha mapping
const ArrayList<short int> &alpha_mapping = multiindex_mapping_[k];
const std::vector<short int> &alpha_mapping = multiindex_mapping_[k];
// initialize the factor to 1
multiindex_combination_.set(j, k, 1);
@@ -139,18 +138,17 @@ class SeriesExpansionAux {
void ComputeUpperMappingIndex() {
int limit = 2 * max_order_;
ArrayList<int> diff;
diff.Init(dim_);
std::vector<int> diff;
diff.reserve(dim_);
// initialize the index
upper_mapping_index_.Init(list_total_num_coeffs_[limit]);
upper_mapping_index_.reserve(list_total_num_coeffs_[limit]);
for(index_t i = 0; i < list_total_num_coeffs_[limit]; i++) {
const ArrayList<short int> &outer_mapping = multiindex_mapping_[i];
upper_mapping_index_[i].Init();
const std::vector<short int> &outer_mapping = multiindex_mapping_[i];
for(index_t j = 0; j < list_total_num_coeffs_[limit]; j++) {
const ArrayList<short int> &inner_mapping = multiindex_mapping_[j];
const std::vector<short int> &inner_mapping = multiindex_mapping_[j];
int flag = 0;
for(index_t d = 0; d < dim_; d++) {
@@ -163,7 +161,7 @@ class SeriesExpansionAux {
}
if(flag == 0) {
(upper_mapping_index_[i]).PushBackCopy(j);
(upper_mapping_index_[i]).push_back(j);
}
} // end of j-loop
} // end of i-loop
@@ -180,13 +178,13 @@ class SeriesExpansionAux {
const Vector& get_inv_multiindex_factorials() const;
const ArrayList< short int > * get_lower_mapping_index() const;
const std::vector< short int > & get_lower_mapping_index() const;
int get_max_order() const;
const ArrayList< short int > & get_multiindex(int pos) const;
const std::vector< short int > & get_multiindex(int pos) const;
const ArrayList< short int > * get_multiindex_mapping() const;
const std::vector< short int > & get_multiindex_mapping() const;
const Vector& get_neg_inv_multiindex_factorials() const;
@@ -194,14 +192,14 @@ class SeriesExpansionAux {
double get_n_multichoose_k_by_pos(int n, int k) const;
const ArrayList< short int > * get_upper_mapping_index() const;
const std::vector< short int > & get_upper_mapping_index() const;
// interesting functions
/**
* Computes the position of the given multiindex
*/
int ComputeMultiindexPosition(const ArrayList<short int> &multiindex) const;
int ComputeMultiindexPosition(const std::vector<short int> &multiindex) const;
/** @brief Computes the computational cost of evaluating a far-field
* expansion of order p at a single query point.