ToPS
GeneralizedPairHiddenMarkovModel.cpp
00001 #include "Alphabet.hpp"
00002 #include "GeneralizedPairHiddenMarkovModel.hpp"
00003 #include "ProbabilisticModelParameter.hpp"
00004 #include "util.hpp"
00005 #include "Symbol.hpp"
00006 #include <iostream>
00007 #include <cmath>
00008 #include <sstream>
00009 #include <vector>
00010 #include <iterator>
00011 #include <stdio.h>
00012 #include <algorithm>
00013 
00014 namespace tops {
00015 
00016   double GeneralizedPairHiddenMarkovModel::forward(const Sequence & seq1, const Sequence & seq2, vector<Matrix> &a)
00017   {
00018     int nstates = _states.size();
00019     int length1 = seq1.size();
00020     int length2 = seq2.size();
00021     vector<Matrix> alpha;
00022     alpha.resize(nstates);
00023     for(int i = 0; i < nstates; i++){
00024       (alpha[i]).resize(length1+1,length2+1);
00025     }
00026 
00027     //Initialization
00028     for(int k = 0; k < nstates; k++){
00029       alpha[k](0,0) = initialProbabilities()->log_probability_of(k);
00030     }
00031 
00032     //Recursion
00033     for (int i = 0; i <= length1; i++){
00034       for(int j = 0; j <= length2; j++){
00035         for(int k = 0; k < nstates; k++){
00036           if(i == 0 && j == 0)
00037             continue;
00038     alpha[k](i,j) = -HUGE;
00039           int maxS1 = getState(k)->maxSeq1();
00040           int minS1 = getState(k)->minSeq1();
00041           int maxS2 = getState(k)->maxSeq2();
00042           int minS2 = getState(k)->minSeq2();
00043     if(i - minS1 < 0 || j - minS2 < 0){
00044       continue;
00045     }
00046     for(int ne1 = minS1; ne1 <= maxS1; ne1++){
00047       if(i - ne1 < 0)
00048         break;
00049       for(int ne2 = minS2; ne2 <= maxS2; ne2++){
00050         if(j - ne2 < 0)
00051     break;
00052         if(ne1 == 0 && ne2 == 0)
00053     continue;
00054         double aux = -HUGE;
00055         if(i - ne1 == 0 && j - ne2 == 0){
00056     aux = alpha[k](i-ne1,j-ne2);
00057     aux += getState(k)->duration()->log_probability_of_pair(ne1,ne2);
00058     if(maxS1 == 1 && maxS2 == 1)
00059       aux += getState(k)->emission()->log_probability_of_pair(seq1[i-1],seq2[j-1]);
00060     if(maxS1 == 1 && maxS2 == 0)
00061       aux += getState(k)->emission()->log_probability_of_pair(seq1[i-1],_gap_id);
00062     if(maxS1 == 0 && maxS2 == 1)
00063       aux += getState(k)->emission()->log_probability_of_pair(_gap_id,seq2[j-1]);
00064     if(maxS1 > 1 && maxS2 > 1){
00065       vector<Matrix> f;
00066       aux += getState(k)->emission()->pairDecodable()->forward(sub_seq(seq1, i-ne1, i-1), sub_seq(seq2, j-ne2, j-1), f);
00067     }
00068     alpha[k](i,j) = log_sum(alpha[k](i,j), aux);
00069     continue;
00070         }
00071         for(int l = 0; l < (int)(getState(k)->iTransitions()).size(); l++){
00072     int id = getState(k)->getITransId(l);
00073     aux = log_sum(aux, alpha[id](i-ne1,j-ne2) + getState(id)->transitions()->log_probability_of(k));
00074         }
00075         aux += getState(k)->duration()->log_probability_of_pair(ne1,ne2);
00076         if(maxS1 == 1 && maxS2 == 1)
00077     aux += getState(k)->emission()->log_probability_of_pair(seq1[i-1],seq2[j-1]);
00078         if(maxS1 == 1 && maxS2 == 0)
00079     aux += getState(k)->emission()->log_probability_of_pair(seq1[i-1],_gap_id);
00080         if(maxS1 == 0 && maxS2 == 1)
00081     aux += getState(k)->emission()->log_probability_of_pair(_gap_id,seq2[j-1]);
00082         if(maxS1 > 1 && maxS2 > 1){
00083     vector<Matrix> f;
00084     aux += getState(k)->emission()->pairDecodable()->forward(sub_seq(seq1, i-ne1, i-1), sub_seq(seq2, j-ne2, j-1), f);
00085         }
00086         alpha[k](i,j) = log_sum(alpha[k](i,j), aux);
00087       }
00088     }
00089   }
00090       }
00091     }
00092 
00093     //Termination
00094     double sum = -HUGE;
00095     for(int i = 0; i < nstates; i++){
00096       alpha[i](length1,length2) += getState(i)->transitions()->log_probability_of(_end_id);
00097       sum = log_sum(sum, alpha[i](length1,length2));
00098     }
00099     a = alpha;
00100     return sum;
00101   }
00102 
00103   double GeneralizedPairHiddenMarkovModel::viterbi(const Sequence & seq1, const Sequence & seq2, Sequence & statePath, Sequence & alignment1, Sequence & alignment2, vector<Matrix> &a)
00104   {
00105     int nstates = _states.size();
00106     int length1 = seq1.size();
00107     int length2 = seq2.size();
00108     vector<Matrix> alpha;
00109     vector<vMatrix> pathMatrix;
00110     alpha.resize(nstates);
00111     pathMatrix.resize(nstates);
00112     for(int i = 0; i < nstates; i++){
00113       (alpha[i]).resize(length1+1,length2+1);
00114       (pathMatrix[i]).resize(length1+1,length2+1);
00115     }
00116 
00117     //Initialization
00118     for(int k = 0; k < nstates; k++){
00119       alpha[k](0,0) = initialProbabilities()->log_probability_of(k);
00120     }
00121 
00122     //Recursion
00123     for (int i = 0; i <= length1; i++){
00124       for(int j = 0; j <= length2; j++){
00125         for(int k = 0; k < nstates; k++){
00126           if(i == 0 && j == 0)
00127             continue;
00128           int maxS1 = getState(k)->maxSeq1();
00129           int minS1 = getState(k)->minSeq1();
00130           int maxS2 = getState(k)->maxSeq2();
00131           int minS2 = getState(k)->minSeq2();
00132     alpha[k](i,j) = -HUGE;
00133     if(i - minS1 < 0 || j - minS2 < 0){
00134       (pathMatrix[k](i,j)).lastState = -1;
00135       continue;
00136     }
00137     for(int ne1 = minS1; ne1 <= maxS1; ne1++){
00138       if(i - ne1 < 0)
00139         break;
00140       for(int ne2 = minS2; ne2 <= maxS2; ne2++){
00141         if(j - ne2 < 0)
00142     break;
00143         if(ne1 == 0 && ne2 == 0)
00144     continue;
00145         double aux = -HUGE;
00146         Sequence saux1, saux2, saux3;
00147         vector<Matrix> f;
00148         int state;
00149         if(i - ne1 == 0 && j - ne2 == 0){
00150     aux = alpha[k](i-ne1,j-ne2);
00151     aux += getState(k)->duration()->log_probability_of_pair((sub_seq(seq1, i-ne1, i-1)).size(),(sub_seq(seq2, j-ne2, j-1)).size());
00152     if(maxS1 == 1 && maxS2 == 1){
00153       aux += getState(k)->emission()->log_probability_of_pair(seq1[i-1],seq2[j-1]);
00154       saux1.push_back(seq1[i-1]);
00155       saux2.push_back(seq2[j-1]);
00156     }
00157     if(maxS1 == 1 && maxS2 == 0){
00158       aux += getState(k)->emission()->log_probability_of_pair(seq1[i-1],_gap_id);
00159       saux1.push_back(seq1[i-1]);
00160       saux2.push_back(_gap_id);
00161     }
00162     if(maxS1 == 0 && maxS2 == 1){
00163       aux += getState(k)->emission()->log_probability_of_pair(_gap_id,seq2[j-1]);
00164       saux1.push_back(_gap_id);
00165       saux2.push_back(seq2[j-1]);
00166     }
00167     if(maxS1 > 1 && maxS2 > 1){
00168       aux += getState(k)->emission()->pairDecodable()->viterbi(sub_seq(seq1, i-ne1, i-1), sub_seq(seq2, j-ne2, j-1), saux3, saux1, saux2, f);
00169     }
00170     if(aux > alpha[k](i,j)){
00171       alpha[k](i,j) = aux;
00172       (pathMatrix[k](i,j)).lastState = -1;
00173       (pathMatrix[k](i,j)).ne1 = ne1;
00174       (pathMatrix[k](i,j)).ne2 = ne2;
00175       (pathMatrix[k](i,j)).al1 = saux1;
00176       (pathMatrix[k](i,j)).al2 = saux2;
00177     }
00178     continue;
00179         }
00180         for(int l = 0; l < (int)(getState(k)->iTransitions()).size(); l++){
00181     int id = getState(k)->getITransId(l);
00182     if(aux < alpha[id](i-ne1,j-ne2) + getState(id)->transitions()->log_probability_of(k)){
00183       aux = alpha[id](i-ne1,j-ne2) + getState(id)->transitions()->log_probability_of(k);
00184       state = id;
00185     }
00186         }
00187         aux += getState(k)->duration()->log_probability_of_pair((sub_seq(seq1, i-ne1, i-1)).size(),(sub_seq(seq2, j-ne2, j-1)).size());
00188         if(maxS1 == 1 && maxS2 == 1){
00189     aux += getState(k)->emission()->log_probability_of_pair(seq1[i-1],seq2[j-1]);
00190     saux1.push_back(seq1[i-1]);
00191     saux2.push_back(seq2[j-1]);
00192         }
00193         if(maxS1 == 1 && maxS2 == 0){
00194     aux += getState(k)->emission()->log_probability_of_pair(seq1[i-1],_gap_id);
00195     saux1.push_back(seq1[i-1]);
00196     saux2.push_back(_gap_id);
00197         }
00198         if(maxS1 == 0 && maxS2 == 1){
00199     aux += getState(k)->emission()->log_probability_of_pair(_gap_id,seq2[j-1]);
00200     saux1.push_back(_gap_id);
00201     saux2.push_back(seq2[j-1]);
00202         }
00203         if(maxS1 > 1 && maxS2 > 1){
00204     aux += getState(k)->emission()->pairDecodable()->viterbi(sub_seq(seq1, i-ne1, i-1), sub_seq(seq2, j-ne2, j-1), saux3, saux1, saux2, f);
00205         }
00206         if(aux > alpha[k](i,j)){
00207     alpha[k](i,j) = aux;
00208     (pathMatrix[k](i,j)).lastState = state;
00209     (pathMatrix[k](i,j)).ne1 = ne1;
00210     (pathMatrix[k](i,j)).ne2 = ne2;
00211     (pathMatrix[k](i,j)).al1 = saux1;
00212     (pathMatrix[k](i,j)).al2 = saux2;
00213         }
00214       }
00215     }
00216   }
00217       }
00218     }
00219 
00220     //Termination
00221     double max = -HUGE;
00222     int lastState = -1;
00223     for(int i = 0; i < nstates; i++){
00224       alpha[i](length1,length2) += getState(i)->transitions()->log_probability_of(_end_id);
00225       if(alpha[i](length1,length2) > max){
00226   max = alpha[i](length1,length2);
00227   lastState = i;
00228       }
00229     }
00230 
00231     //Traceback
00232     int i = length1;
00233     int j = length2;
00234     Sequence al1,al2,path;
00235     while(lastState != -1){
00236       alignment1.insert(alignment1.end(),((pathMatrix[lastState](i,j)).al1).begin(),((pathMatrix[lastState](i,j)).al1).end());
00237       alignment2.insert(alignment2.end(),((pathMatrix[lastState](i,j)).al2).begin(),((pathMatrix[lastState](i,j)).al2).end());
00238       int alsize = ((pathMatrix[lastState](i,j)).al1).size();
00239       for(int k = 0; k < alsize; k++)
00240   statePath.push_back(lastState);
00241       int ne1 = (pathMatrix[lastState](i,j)).ne1;
00242       int ne2 = (pathMatrix[lastState](i,j)).ne2;
00243       lastState = (pathMatrix[lastState](i,j)).lastState;
00244       i -= ne1;
00245       j -= ne2;
00246     }
00247     a = alpha;
00248     return max;
00249   }
00250 
00251   double GeneralizedPairHiddenMarkovModel::backward(const Sequence & seq1, const Sequence & seq2, vector<Matrix> &a)
00252   {
00253     int nstates = _states.size();
00254     int length1 = seq1.size();
00255     int length2 = seq2.size();
00256     vector<Matrix> alpha;
00257     alpha.resize(nstates);
00258     for(int i = 0; i < nstates; i++){
00259       (alpha[i]).resize(length1+1,length2+1);
00260     }
00261 
00262     //Initialization
00263     for(int k = 0; k < nstates; k++){
00264       alpha[k](length1,length2) = getState(k)->transitions()->log_probability_of(_end_id);
00265     }
00266 
00267     //Recursion
00268     for (int i = length1; i >= 0; i--){
00269       for(int j = length2; j >= 0; j--){
00270         for(int k = 0; k < nstates; k++){
00271 
00272           if(i == length1 && j == length2)
00273             continue;
00274 
00275     alpha[k](i,j) = -HUGE;
00276     if(i == 0 && j == 0){
00277       int maxS1 = getState(k)->maxSeq1();
00278       int minS1 = getState(k)->minSeq1();
00279       int maxS2 = getState(k)->maxSeq2();
00280       int minS2 = getState(k)->minSeq2();
00281       for(int ne1 = minS1; ne1 <= maxS1; ne1++){
00282         if(i + ne1 > length1)
00283     break;
00284         for(int ne2 = minS2; ne2 <= maxS2; ne2++){
00285     if(j + ne2 > length2)
00286       break;
00287     if(ne1 == 0 && ne2 == 0)
00288       continue;
00289     if(maxS1 == 1 && maxS2 == 1)
00290       alpha[k](i,j) = log_sum(alpha[k](i,j),
00291             alpha[k](i+ne1,j+ne2) +
00292             initialProbabilities()->log_probability_of(k) +
00293             getState(k)->emission()->log_probability_of_pair(seq1[i-1+ne1],seq2[j-1+ne2]));
00294     else if(maxS1 == 1 && maxS2 == 0)
00295       alpha[k](i,j) = log_sum(alpha[k](i,j),
00296             alpha[k](i+ne1,j+ne2) +
00297             initialProbabilities()->log_probability_of(k)+
00298             getState(k)->emission()->log_probability_of_pair(seq1[i-1+ne1],_gap_id));
00299     else if(maxS1 == 0 && maxS2 == 1)
00300       alpha[k](i,j) = log_sum(alpha[k](i,j),
00301             alpha[k](i+ne1,j+ne2) +
00302             initialProbabilities()->log_probability_of(k)+
00303             getState(k)->emission()->log_probability_of_pair(_gap_id,seq2[j-1+ne2]));
00304     else if(maxS1 > 1 && maxS2 > 1){
00305       vector<Matrix> b;
00306       alpha[k](i,j) = log_sum(alpha[k](i,j),
00307             alpha[k](i+ne1,j+ne2) +
00308             initialProbabilities()->log_probability_of(k) +
00309             getState(k)->duration()->log_probability_of_pair(ne1,ne2)+
00310             getState(k)->emission()->pairDecodable()->forward(sub_seq(seq1, i, i+ne1-1),sub_seq(seq2, j, j+ne2-1),b));
00311     }
00312         }
00313       }
00314       continue;
00315     }
00316 
00317     for(int l = 0; l < (int)(getState(k)->oTransitions()).size(); l++){
00318       int id = getState(k)->getOTransId(l);
00319       int maxS1 = getState(id)->maxSeq1();
00320       int minS1 = getState(id)->minSeq1();
00321       int maxS2 = getState(id)->maxSeq2();
00322       int minS2 = getState(id)->minSeq2();
00323       for(int ne1 = minS1; ne1 <= maxS1; ne1++){
00324         if(i + ne1 > length1)
00325     break;
00326         for(int ne2 = minS2; ne2 <= maxS2; ne2++){
00327     if(j + ne2 > length2)
00328       break;
00329     if(ne1 == 0 && ne2 == 0)
00330       continue;
00331     if(maxS1 == 1 && maxS2 == 1)
00332       alpha[k](i,j) = log_sum(alpha[k](i,j),
00333             alpha[id](i+ne1,j+ne2) +
00334             getState(k)->transitions()->log_probability_of(id) +
00335             getState(id)->emission()->log_probability_of_pair(seq1[i-1+ne1],seq2[j-1+ne2]));
00336     else if(maxS1 == 1 && maxS2 == 0)
00337       alpha[k](i,j) = log_sum(alpha[k](i,j),
00338             alpha[id](i+ne1,j+ne2) +
00339             getState(k)->transitions()->log_probability_of(id)+
00340             getState(id)->emission()->log_probability_of_pair(seq1[i-1+ne1],_gap_id));
00341     else if(maxS1 == 0 && maxS2 == 1)
00342       alpha[k](i,j) = log_sum(alpha[k](i,j),
00343             alpha[id](i+ne1,j+ne2) +
00344             getState(k)->transitions()->log_probability_of(id)+
00345             getState(id)->emission()->log_probability_of_pair(_gap_id,seq2[j-1+ne2]));
00346     else if(maxS1 > 1 && maxS2 > 1){
00347       vector<Matrix> b;
00348       alpha[k](i,j) = log_sum(alpha[k](i,j),
00349             alpha[id](i+ne1,j+ne2) +
00350             getState(k)->transitions()->log_probability_of(id)+
00351             getState(id)->duration()->log_probability_of_pair(ne1,ne2)+
00352             getState(id)->emission()->pairDecodable()->forward(sub_seq(seq1, i, i+ne1-1),sub_seq(seq2, j, j+ne2-1),b));
00353     }
00354         }
00355       }
00356     }
00357         }
00358       }
00359     }
00360 
00361     //Termination
00362     double sum = -HUGE;
00363     for(int i = 0; i < nstates; i++){
00364       sum = log_sum(sum, alpha[i](0,0));
00365     }
00366     a = alpha;
00367     return sum;
00368   }
00369 
00370   void GeneralizedPairHiddenMarkovModel::initialize(const ProbabilisticModelParameters & parameters) {
00371     ProbabilisticModelParameterValuePtr state_names = parameters.getMandatoryParameterValue("state_names");
00372     ProbabilisticModelParameterValuePtr observation_symbols = parameters.getMandatoryParameterValue("observation_symbols");
00373     ProbabilisticModelParameterValuePtr number_of_emissions = parameters.getMandatoryParameterValue("number_of_emissions");
00374     ProbabilisticModelParameterValuePtr initial_probabilities = parameters.getMandatoryParameterValue("initial_probabilities");
00375     ProbabilisticModelParameterValuePtr end_probabilities = parameters.getMandatoryParameterValue("end_probabilities");
00376     ProbabilisticModelParameterValuePtr transitions = parameters.getMandatoryParameterValue("transitions");
00377     ProbabilisticModelParameterValuePtr emissions = parameters.getMandatoryParameterValue("emission_probabilities");
00378     ProbabilisticModelParameterValuePtr durations = parameters.getOptionalParameterValue("duration_probabilities");
00379 
00380     ProbabilisticModelCreatorClient creator;
00381 
00382     AlphabetPtr states = AlphabetPtr(new Alphabet());
00383     states->initializeFromVector(state_names->getStringVector());
00384 
00385     _end_id = states->size();
00386 
00387     AlphabetPtr observations = AlphabetPtr(new Alphabet());
00388     observations->initializeFromVector(observation_symbols->getStringVector());
00389     if(observations->has("_"))
00390       _gap_id = observations->getSymbol("_")->id();
00391     else{
00392       observations->createSymbol(std::string("_"));
00393       _gap_id = observations->getSymbol("_")->id();
00394     }
00395 
00396     DiscreteIIDModelPtr init = DiscreteIIDModelPtr(new DiscreteIIDModel());
00397     init->initializeFromMap(initial_probabilities->getDoubleMap(), states);
00398 
00399     std::map<std::string,std::string> emissModels = emissions->getStringMap();
00400 
00401     std::map<std::string,std::string> durModels = durations->getStringMap();
00402 
00403     std::map<std::string,std::string> numEmiss = number_of_emissions->getStringMap();
00404 
00405     std::map<std::string,double> transpar = transitions->getDoubleMap();
00406     std::map<std::string,double> endMap = end_probabilities->getDoubleMap();
00407     std::map<std::string,double>::const_iterator transparit;
00408     std::map<std::string,DoubleVector> trans;
00409     std::map<std::string,IntVector> inTrans;
00410     std::map<std::string,IntVector> outTrans;
00411     for(transparit = transpar.begin(); transparit != transpar.end(); transparit++){
00412       std::vector<std::string> splited;
00413       boost::regex separator("\\|");
00414       split_regex(transparit->first, splited, separator);
00415       if(splited.size() == 1)
00416         splited.push_back("");
00417 
00418       std::string to(splited[0]);
00419       std::string from(splited[1]);
00420       if(trans.find(from) == trans.end())
00421         {
00422           int id = states->getSymbol(to)->id();
00423           DoubleVector probs;
00424           IntVector outgoing;
00425           outgoing.push_back(id);
00426           outTrans[from] = outgoing;
00427           probs.resize(states->size() + 1);
00428     trans[from]=probs;
00429     if(endMap.find(from) != endMap.end())
00430       trans[from][_end_id] = endMap[from];
00431     else
00432       trans[from][_end_id] = 0.0;
00433     if(id < (int)trans[from].size())
00434             trans[from][id] = transparit->second;
00435         }
00436       else
00437         {
00438           int id = states->getSymbol(to)->id();
00439           outTrans[from].push_back(id);
00440           if(id < (int)trans[from].size())
00441             trans[from][id] = transparit->second;
00442         }
00443       if(inTrans.find(to) == inTrans.end()){
00444         int id = states->getSymbol(from)->id();
00445         IntVector incoming;
00446         incoming.push_back(id);
00447         inTrans[to] = incoming;
00448       }
00449       else{
00450         int id = states->getSymbol(from)->id();
00451         inTrans[to].push_back(id);
00452       }
00453     }
00454 
00455     std::vector<PairDecodableStatePtr> states_list;
00456 
00457     for(int i = 0; i < (int)states->size(); i++){
00458       SymbolPtr state_name = states->getSymbol(i);
00459       DiscreteIIDModelPtr t;
00460       IntVector itr, otr;
00461       ProbabilisticModelPtr dur, emiss;
00462       int maxS1, minS1, maxS2, minS2;
00463 
00464       std::map<std::string,DoubleVector>::const_iterator strdvit = trans.find(state_name->name());
00465       if(strdvit != trans.end())
00466   t = DiscreteIIDModelPtr(new DiscreteIIDModel(strdvit->second));
00467       else{
00468   cerr << "State " << state_name->name() << " has no transitions." << endl;
00469   exit(-1);
00470       }
00471 
00472       std::map<std::string,IntVector>::const_iterator strivit = inTrans.find(state_name->name());
00473       if(strivit != inTrans.end())
00474   itr = strivit->second;
00475       else{
00476   cerr << "State " << state_name->name() << " has no incoming transitions." << endl;
00477   exit(-1);
00478       }
00479 
00480       strivit = outTrans.find(state_name->name());
00481       if(strivit != outTrans.end())
00482   otr = strivit->second;
00483       else{
00484   cerr << "State " << state_name->name() << " has no transitions coming out of it." << endl;
00485   exit(-1);
00486       }
00487 
00488       std::map<std::string,std::string>::const_iterator strmapit = durModels.find(state_name->name());
00489       if(strmapit != durModels.end()){
00490   std::string model_file = strmapit->second;
00491   ProbabilisticModelPtr model = creator.create(model_file);
00492   if(model == NULL)
00493     {
00494       cerr << "Could not create model from file " << model_file << endl;
00495       exit(-1);
00496     }
00497   dur = model;
00498       }
00499       else
00500   dur = ConstantModelPtr(new ConstantModel());
00501 
00502       strmapit = emissModels.find(state_name->name());
00503       if(strmapit != emissModels.end()){
00504   std::string model_file = strmapit->second;
00505   ProbabilisticModelPtr model = creator.create(model_file);
00506   if(model == NULL)
00507     {
00508       cerr << "Could not create model from file " << model_file << endl;
00509       exit(-1);
00510     }
00511   emiss = model;
00512       }
00513       else{
00514   cerr << "Emission probabilities for state " << state_name->name() << " could not be initialized." << endl;
00515   exit(-1);
00516       }
00517 
00518       strmapit = numEmiss.find(state_name->name());
00519       if(strmapit != numEmiss.end()){
00520   std::vector<std::string> spl1;
00521   std::vector<std::string> spl2;
00522   std::vector<std::string> spl3;
00523   boost::regex sep1(",");
00524   boost::regex sep2("-");
00525   split_regex(strmapit->second, spl1, sep1);
00526   split_regex(spl1[0], spl2, sep2);
00527   split_regex(spl1[1], spl3, sep2);
00528   minS1 = atoi(spl2[0].c_str());
00529   maxS1 = atoi(spl2[1].c_str());
00530   minS2 = atoi(spl3[0].c_str());
00531   maxS2 = atoi(spl3[1].c_str());
00532       }
00533       else{
00534   cerr << "Number of emissions for state " << state_name->name() << " could not be initialized." << endl;
00535   exit(-1);
00536       }
00537 
00538       GPHMMStatePtr statePtr = GPHMMStatePtr(new GPHMMState(i, state_name, emiss, t, dur, itr, otr, maxS1, minS1, maxS2, minS2));
00539       states_list.push_back(statePtr);
00540     }
00541 
00542     setStates(states_list, states);
00543     setObservationSymbols(observations);
00544     setInitialProbabilities(init);
00545   }
00546 }