C++ HMM类代码示例

int main_run(char * args[]){
	vector<vector<T>> ufile = read(args[1]);
		vector<vector<T>> tfile = read(args[2]);
		vector<vector<T>> logs = read(args[3]);
		HMM<T> hmm;
		/*for (auto i = transformed.begin(); i != transformed.end(); ++i) {
		  for (auto j = i->begin(); j != i->end(); ++j) {
		  cout<<j->first<<"\t"<<j->second<<"\n";
		  }
		  }*/
		hmm.Train(ufile, tfile);
		map<int, string> tagcodemap;
		vector<vector<T>> taglist = hmm.TagLogs(logs, tagcodemap);
		vector<vector<T>> transformed = transform2(taglist);
		//FP GROWTH PART FOR CLUSTERING after classification
		int s = atoi(args[4]);
		FPTree<T> tree(s);
		//consider the position of token in the line
		tree.build(transformed);
		std::vector<std::vector<T>> r;
		clock_t t = clock();
		tree.mine(r);
		std::cout<<"Initial pattern no: "<<r.size()<<"\n";
		Cluster<T> data;
		data.AssociatePatterns(transformed, r);
		cout << "Time: " << double(clock() - t) / CLOCKS_PER_SEC << "\n";
		data.DisplayCluster(logs, transformed, r, tagcodemap);
		Analyze<T> analyze;
		vector<vector<T>> manual = read("manual");
		cout<<"\n"<<analyze.Efficiency(taglist, manual);
}

void HMMModel::loadHMMs(QString strPath)
{
	if (!strPath.endsWith('/') && !strPath.endsWith('\\'))
	{
		strPath += "/";
	}
	QStringList filters;
	filters << "*.hmm";
	QDir dir(strPath);
	QStringList list = dir.entryList(filters);

	cleanHMMs();
	QString strFileName;
	int nPos;
	std::string strSignName, strPathFile;
	HMM* phmm;
	const char* pCharCom;
	char* tmp;
	for (int l=0; l<list.count(); l++)
	{
		strFileName = list.at(l);
		nPos = strFileName.lastIndexOf(".hmm");
		strSignName = strFileName.left(nPos).toStdString();
		strPathFile = strPath.toLocal8Bit()+strFileName.toLocal8Bit();
		pCharCom = strPathFile.c_str();
		tmp = new char[strlen(pCharCom)+1];
		strcpy(tmp, pCharCom);

		phmm = new HMM();
		phmm->Load(tmp);
		m_mHMM.insert(std::pair<std::string, HMM*>(strSignName, phmm));
	}
}

void LoadHMM(HMM<distribution::DiscreteDistribution>& hmm,
             util::SaveRestoreUtility& sr)
{
  std::string type;
  size_t states;

  sr.LoadParameter(type, "hmm_type");
  if (type != "discrete")
  {
    Rcpp::Rcout << "Cannot load non-discrete HMM (of type " << type << ") as "
        << "discrete HMM!" << std::endl;
  }

  sr.LoadParameter(states, "hmm_states");

  // Load transition matrix.
  sr.LoadParameter(hmm.Transition(), "hmm_transition");

  // Now each emission distribution.
  hmm.Emission().resize(states);
  for (size_t i = 0; i < states; ++i)
  {
    std::stringstream s;
    s << "hmm_emission_distribution_" << i;
    sr.LoadParameter(hmm.Emission()[i].Probabilities(), s.str());
  }

  hmm.Dimensionality() = 1;
}

int 
Clusters::
reclusterNeyEssen()
{
  cerr << "Calculating MLE for prior probabilities" << endl;
  vector<double> prior(numberClasses,0.0l);
  double xxx = 1.0l/(double)numberTypes;
  for (int i = 0; i < numberTypes; i++){
    int c = classVector[i];
    prior[c] += xxx;
  }
  for (int i = 0; i < numberClasses;i++)
    cerr << i << " " << prior[i] << endl;
  if (numberStates > 0){
    cerr << "Training all the HMMs" << endl;
    for (int c = 0; c < numberClasses; c++){
      //      cerr << "Training HMM " << c << endl;
      HMM* hmmPtr = hmms[c];
      HMM* newHmmPtr = new HMM(numberStates, alphabetSize);
      for (int i = 0; i < numberTypes; i++){
	if (classVector[i] == c){
	  // then this word is in the right class
	  // so train it on word i
	  const string & word = *(corpus.wordArray[i]);
	  vector<int> v;
	  hmmPtr->convertString(word,v);
	  // FIXME 
	  double weight = 1.0l;
	  if (USE_TRUE_WEIGHT){
	    weight = corpus.countArray[i];
	  }
	  hmmPtr->emSingle(*newHmmPtr, weight, v);
	}
      }
      newHmmPtr->normalise();
      hmms[c] = newHmmPtr;
      delete hmmPtr;
    }
  }
  int something = 0;
  for (int i = 0; i < numberTypes; i++){
    //    cerr << "Word " << i;
    int w = sortedWords[i];
    //cerr << *(corpus.wordArray[w]) << endl;
    if (counts[w] > FREQ_CUTOFF){
      //cerr << "Doing " << w << endl;
      if (bestCluster(w, prior)){
	something++;
      }
    }
  }
 
  return something;
}

// [[Rcpp::export]]
List launch_pmmh_cpp(List inputs, List modellist,
                 List algoparameters){
  string strmodel = as<string>(modellist["model"]);
  string strprior = as<string>(modellist["prior"]);
  NumericMatrix current = inputs["current"];
  NumericMatrix observations = inputs["observations"];
  List theta = inputs["theta"];
  int nparticles = as<int>(algoparameters["nparticles"]);
  int niterations = as<int>(algoparameters["niterations"]);
  NumericMatrix cholesky_proposal = algoparameters["cholesky_proposal"];
  NumericVector initial_theta = algoparameters["initial_theta"];
  HMM* model;
  Prior* prior;
  if (strmodel.compare(string("model1")) == 0){
    model = new BatteryModel();
    if (strprior.compare(string("uniform")) == 0){
      prior = new BatteryModelUniformPrior();
    }
    if (strprior.compare(string("normal")) == 0){
      prior = new BatteryModelNormalPrior();
    }
  }
  if (strmodel.compare(string("model2")) == 0){
    model = new BatteryModel2();
    if (strprior.compare(string("uniform")) == 0){
      prior = new BatteryModel2UniformPrior();
    }
    if (strprior.compare(string("normal")) == 0){
      prior = new BatteryModel2NormalPrior();
    }
  }
  model->set_input(current);
  model->set_parameters(theta);
  model->set_observations(observations);
  PMMH pmmh(nparticles, niterations, model->dim_states);
  pmmh.set_prior(prior);
  pmmh.init(model, initial_theta);
  pmmh.set_proposal_cholesky(cholesky_proposal);
  pmmh.run();

  delete model;
  delete prior;
  return Rcpp::List::create(
                            Rcpp::Named("chain")= pmmh.chain_parameters,
                            Rcpp::Named("naccepts") = pmmh.naccepts,
                            Rcpp::Named("loglikelihood") = pmmh.loglikelihood,
                            Rcpp::Named("loglikelihood_proposal") = pmmh.loglikelihood_proposal,
                            Rcpp::Named("proposals") = pmmh.proposals,
                            Rcpp::Named("nparticles") = nparticles,
                            Rcpp::Named("niterations") = niterations,
                            Rcpp::Named("cholesky_proposal") = cholesky_proposal);
}

PosteriorViterbi::PosteriorViterbi(HMM &hmm,
				   bool shouldAdd)
  : numStates(hmm.countStates()), hmmGraph(hmm),
    addRatherThanMultiply(shouldAdd)
{
  // ctor
}

void SaveHMM(const HMM<gmm::GMM<> >& hmm,
             util::SaveRestoreUtility& sr)
{
  std::string type = "gmm";
  size_t states = hmm.Transition().n_rows;

  sr.SaveParameter(type, "hmm_type");
  sr.SaveParameter(states, "hmm_states");
  sr.SaveParameter(hmm.Transition(), "hmm_transition");

  // Now the emissions.
  for (size_t i = 0; i < states; ++i)
  {
    // Generate name.
    std::stringstream s;
    s << "hmm_emission_" << i << "_gaussians";
    sr.SaveParameter(hmm.Emission()[i].Gaussians(), s.str());

    s.str("");
    s << "hmm_emission_" << i << "_weights";
    sr.SaveParameter(hmm.Emission()[i].Weights(), s.str());

    for (size_t g = 0; g < hmm.Emission()[i].Gaussians(); ++g)
    {
      s.str("");
      s << "hmm_emission_" << i << "_gaussian_" << g << "_mean";
      sr.SaveParameter(hmm.Emission()[i].Means()[g], s.str());

      s.str("");
      s << "hmm_emission_" << i << "_gaussian_" << g << "_covariance";
      sr.SaveParameter(hmm.Emission()[i].Covariances()[g], s.str());
    }
  }
}

HMM
buildHMM(const ESBTL::Default_system & system, GMMCreator creator, ClusteringFunctor func) {
  arma::mat data;
  arma::urowvec labels;
  HMM hmm;
  unsigned int offset = 0;
  for (ESBTL::Default_system::Models_const_iterator it_model = system.models_begin();
      it_model != system.models_end();
      ++it_model) {
    const ESBTL::Default_system::Model & model = *it_model;
    arma::mat tempdata;
    arma::urowvec templabels = func(model, tempdata); 
    templabels += offset;
    data = arma::join_rows(data, tempdata);
    labels = arma::join_rows(labels, templabels);
    offset = arma::max(labels);
  }
  hmm.baumWelchCached(data, creator(data, labels));
  return hmm;
}

int main() {
    kmeans("./dataset", "./result");
    Forecast forecast("./result");
    HMM hmm = HMM("./dataset", "./dictionnary.txt", forecast);

    hmm.print();
    {
        cout << "Save" << endl;
        std::ofstream ofs("hmm_save");
        boost::archive::text_oarchive oa(ofs);
        oa << hmm;
    }
    
    HMM hmm2 = HMM();
    {
        cout << "Load" << endl;
        std::ifstream ifs("hmm_save");
        boost::archive::text_iarchive ia(ifs);
        ia >> hmm2;
    }
    hmm2.print();    
    return (EXIT_SUCCESS);
}

int main(int argc, char* argv[])
{
  HMM h;

  h.init();
  std::cout << h;

  forward_viterbi(h.get_observations(), 
                  h.get_states(), 
                  h.get_start_probability(), 
                  h.get_transition_probability(), 
                  h.get_emission_probability());
}

// Calculate secondary structure for given HMM and return prediction
void CalculateSS(HMM& q, char *ss_pred, char *ss_conf)
{
  char tmpfile[]="/tmp/hhCalcSSXXXXXX";
  if (mkstemp(tmpfile) == -1) {
    cerr << "ERROR! Could not create tmp-file!\n"; 
    exit(4);
  }
  
  // Write log-odds matrix from q to tmpfile.mtx
  char filename[NAMELEN];
  FILE* mtxf = NULL;

  strcpy(filename,tmpfile);
  strcat(filename,".mtx");
  mtxf = fopen(filename,"w");
  if (!mtxf) OpenFileError(filename);

  fprintf(mtxf,"%i\n",q.L);
  fprintf(mtxf,"%s\n",q.seq[q.nfirst]+1);
  fprintf(mtxf,"2.670000e-03\n4.100000e-02\n-3.194183e+00\n1.400000e-01\n2.670000e-03\n4.420198e-02\n-3.118986e+00\n1.400000e-01\n3.176060e-03\n1.339561e-01\n-2.010243e+00\n4.012145e-01\n");
  
  for (int i = 1; i <= q.L; ++i) 
    {
      fprintf(mtxf,"-32768 ");
      for (int a = 0; a < 20; ++a)
	{
	  int tmp = iround(50*flog2(q.p[i][s2a[a]]/pb[s2a[a]]));
	  fprintf(mtxf,"%5i ",tmp);
	  if (a == 0) {   // insert logodds value for B
	    fprintf(mtxf,"%5i ",-32768);
	  } else if (a == 18) {   // insert logodds value for X
	    fprintf(mtxf,"%5i ",-100);
	  } else if (a == 19) {   // insert logodds value for Z
	    fprintf(mtxf,"%5i ",-32768);
	  }
	}
      fprintf(mtxf,"-32768 -400\n");
    }
  fclose(mtxf);

  // Calculate secondary structure
  CalculateSS(ss_pred, ss_conf, tmpfile);
  
  q.AddSSPrediction(ss_pred, ss_conf);

  // Remove temp-files
  std::string command = "rm " + (std::string)tmpfile + "*";
  runSystem(command,v);
}

/////////////////////////////////////////////////////////////////////////////////////
// Do precalculations for q and t to prepare comparison
/////////////////////////////////////////////////////////////////////////////////////
void PrepareTemplate(HMM& q, HMM& t, int format)
{
    if (format==0) // HHM format
    {
        // Add transition pseudocounts to template
        t.AddTransitionPseudocounts();

	// Don't use CS-pseudocounts because of runtime!!!
	// Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
	t.PreparePseudocounts();
	
	// Add amino acid pseudocounts to query:  p[i][a] = (1-tau)*f[i][a] + tau*g[i][a]
	t.AddAminoAcidPseudocounts(par.pcm, par.pca, par.pcb, par.pcc);

        t.CalculateAminoAcidBackground();
    }
    else // HHMER format
    {
        // Don't add transition pseudocounts to template
        // t.AddTransitionPseudocounts(par.gapd, par.gape, par.gapf, par.gapg, par.gaph, par.gapi, 0.0);

        // Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
        // t.PreparePseudocounts();

        // DON'T ADD amino acid pseudocounts to temlate: pcm=0!  t.p[i][a] = t.f[i][a]
        t.AddAminoAcidPseudocounts(0, par.pca, par.pcb, par.pcc);
        t.CalculateAminoAcidBackground();
    }

    if (par.forward>=1) t.Log2LinTransitionProbs(1.0);

    // Factor Null model into HMM t
    // ATTENTION! t.p[i][a] is divided by pnul[a] (for reasons of efficiency) => do not reuse t.p
    t.IncludeNullModelInHMM(q,t);  // Can go BEFORE the loop if not dependent on template

    return;
}

int main (int argc, const char * argv[])
{
    TimeSeriesClassificationData trainingData;      //This will store our training data

    GestureRecognitionPipeline pipeline;               //This is a wrapper for our classifier and any pre/post processing modules

	string dirPath = "/home/vlad/AndroidStudioProjects/DataCapture/dataSetGenerator/build";


	if (!trainingData.loadDatasetFromFile(dirPath + "/acc-training-set-segmented.data")) {
		printf("Cannot open training set\n");
		return 0;
	}

	printf("Successfully opened training data set ...\n");

    HMM hmm;
    hmm.setHMMType( HMM_CONTINUOUS );
    hmm.setDownsampleFactor( 5 );
    hmm.setAutoEstimateSigma( true );
    hmm.setSigma( 20.0 );

    hmm.setModelType( HMM_LEFTRIGHT );
    hmm.setDelta( 1 );
//    LowPassFilter lpf(0.1, 1, 3);
//    pipeline.setPreProcessingModule(lpf);

    pipeline.setClassifier( hmm );

    pipeline.train(trainingData, 20);

    //You can then get then get the accuracy of how well the pipeline performed during the k-fold cross validation testing
    double accuracy = pipeline.getCrossValidationAccuracy();

    printf("Accuracy: %f\n", accuracy);

}

tuple<vector<double>, Matrix<double>, Matrix<double>> forward_backward(string obs, const HMM& model)
{
    if (!model.isFinalized())
        throw runtime_error("Model should be finalized!");
    
    // Forward algorithm
    Matrix<double> forward(obs.length(), model.numStates(), 0);
    vector<double> cs(obs.length(), 0);
    
    // Calculate c1
    for (size_t state = 0; state < model.numStates(); state++)
        cs[0] += model.startProb(state) * model.emissionProb(state, obs.substr(0,1));
    // Base case
    for (size_t state = 0; state < model.numStates(); state++)
        forward(0, state) = model.startProb(state) * model.emissionProb(state, obs.substr(0,1)) / cs[0];

    // Recursion
    for (size_t i = 1; i < obs.length(); i++) {
        vector<double> delta(model.numStates(), 0);
        for (size_t state = 0; state < model.numStates(); state++) {
            if (i < model.stateArity(state))
                continue;
            
            for (auto prevState : model.incommingStates(state)) {
                double val = forward(i - model.stateArity(state), prevState) * model.transitionProb(prevState, state);
                for (size_t k = 1; k < model.stateArity(state); k++)
                    val /= cs[i - k];
                
                delta[state] += val;
            }
            delta[state] *= model.emissionProb(state, obs.substr(i - model.stateArity(state) + 1, model.stateArity(state)));
            
            cs[i] += delta[state];
        }
        
        for (size_t state = 0; state < model.numStates(); state++) {
            forward(i, state) = delta[state] / cs[i];
        }
    }
    
    // Backward algorithm
    Matrix<double> backward(obs.length(), model.numStates(), 0);
    const size_t N = obs.length() - 1;
    for (size_t state = 0; state < model.numStates(); state++)
        backward(N, state) = 1;
    
    for (long i = N - 1; i >= 0; i--) {
        for (size_t state = 0; state < model.numStates(); state++) {
            double prob = 0;
            
            for (auto nextState : model.outgoingStates(state)) {
                if (i + model.stateArity(nextState) > N)
                    continue;
                
                double val = backward(i + model.stateArity(nextState), nextState) * model.transitionProb(state, nextState)
                               * model.emissionProb(nextState, obs.substr(i + 1, model.stateArity(nextState)));
                
                for (size_t k = 0; k < model.stateArity(nextState); k++)
                    val /= cs[i + 1 + k];
                
                prob += val;
            }
            backward(i, state) = prob;
        }
    }
    
    return make_tuple(cs, forward, backward);
}

FastViterbi::FastViterbi(HMM &hmm,bool posterior)
  : numStates(hmm.countStates()), hmmGraph(hmm), posterior(posterior)
{
  // ctor
}