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