matthiasm@43: 
matthiasm@43: #include "viterbi.h"
matthiasm@43: #include <iostream>
matthiasm@43: 
matthiasm@106: std::vector<int> ViterbiPath(std::vector<double> init, std::vector<vector<double> > trans, std::vector<vector<double> > obs, double *delta, vector<double> *scale) {
matthiasm@43:     
matthiasm@43:     int nState = init.size();
matthiasm@43:     int nFrame = obs.size();
matthiasm@43:     // cerr << nState << " " << nFrame << endl;
matthiasm@43:     
matthiasm@43:     // check for consistency
matthiasm@43:     if (trans[0].size() != nState || trans.size() != nState || obs[0].size() != nState) {
matthiasm@43:         cerr << "ERROR: matrix sizes inconsistent." << endl;
matthiasm@43:     }
matthiasm@43:         
matthiasm@50:     // vector<vector<double> > delta; // "matrix" of conditional probabilities    
matthiasm@43:     vector<vector<int> > psi; //  "matrix" of remembered indices of the best transitions
matthiasm@43:     vector<int> path = vector<int>(nFrame, nState-1); // the final output path (current assignment arbitrary, makes sense only for Chordino, where nChord-1 is the "no chord" label)
matthiasm@106:     // vector<double> scale = vector<double>(nFrame, 0); // remembers by how much the vectors in delta are scaled.
matthiasm@43:     
matthiasm@43:     double deltasum = 0;
matthiasm@43:     
matthiasm@43:     /* initialise first frame */
matthiasm@50:     // delta.push_back(init);    
matthiasm@43:     for (int iState = 0; iState < nState; ++iState) {
matthiasm@94:         delta[iState] = init[iState] * obs[0][iState];
matthiasm@50:         deltasum += delta[iState];
matthiasm@43:     }
matthiasm@50:     for (int iState = 0; iState < nState; ++iState) delta[iState] /= deltasum; // normalise (scale)
matthiasm@106:     scale->push_back(1.0/deltasum);
matthiasm@43:     psi.push_back(vector<int>(nState,0));
matthiasm@43:     
matthiasm@43:     /* rest of the forward step */
matthiasm@43:     for (int iFrame = 1; iFrame < nFrame; ++iFrame) {
matthiasm@50:         // delta.push_back(vector<double>(nState,0));
matthiasm@43:         deltasum = 0;
matthiasm@43:         psi.push_back(vector<int>(nState,0));
matthiasm@43:         /* every state wants to know which previous state suits him best */
matthiasm@43:         for (int jState = 0; jState < nState; ++jState) {            
matthiasm@43:             int bestState = nState - 1;
matthiasm@43:             double bestValue = 0;
matthiasm@44:             if (obs[iFrame][jState] > 0) {
matthiasm@44:                 for (int iState = 0; iState < nState; ++iState) {
matthiasm@50:                     double currentValue = delta[(iFrame-1) * nState + iState] * trans[iState][jState];
matthiasm@44:                     if (currentValue > bestValue) {
matthiasm@44:                         bestValue = currentValue;
matthiasm@44:                         bestState = iState;
matthiasm@44:                     }
matthiasm@43:                 }
matthiasm@43:             }
matthiasm@43:             // cerr << bestState <<" ::: " << bestValue << endl ;
matthiasm@50:             delta[iFrame * nState + jState] = bestValue * obs[iFrame][jState];
matthiasm@50:             deltasum += delta[iFrame * nState + jState];
matthiasm@43:             psi[iFrame][jState] = bestState;
matthiasm@43:         }
matthiasm@43:         if (deltasum > 0) {
matthiasm@43:             for (int iState = 0; iState < nState; ++iState) {            
matthiasm@50:                 delta[iFrame * nState + iState] /= deltasum; // normalise (scale)
matthiasm@43:             }
matthiasm@106:             scale->push_back(1.0/deltasum);
matthiasm@43:         } else {
matthiasm@43:             for (int iState = 0; iState < nState; ++iState) {            
matthiasm@50:                 delta[iFrame * nState + iState] = 1.0/nState;
matthiasm@43:             }
matthiasm@106:             scale->push_back(1.0);
matthiasm@43:         }
matthiasm@43:         
matthiasm@43:     }
matthiasm@43:     
matthiasm@43:     /* initialise backward step */
matthiasm@43:     int bestValue = 0;
matthiasm@43:     for (int iState = 0; iState < nState; ++iState) {
matthiasm@50:         double currentValue = delta[(nFrame-1) * nState + iState];
matthiasm@43:         if (currentValue > path[nFrame-1]) {
matthiasm@43:             bestValue = currentValue;            
matthiasm@43:             path[nFrame-1] = iState;
matthiasm@43:         }
matthiasm@43:     }
matthiasm@43: 
matthiasm@43:     /* rest of backward step */
matthiasm@43:     for (int iFrame = nFrame-2; iFrame > -1; --iFrame) {
matthiasm@43:         path[iFrame] = psi[iFrame+1][path[iFrame+1]];
matthiasm@43:         // cerr << path[iFrame] << endl;
matthiasm@43:     }    
matthiasm@43:     
matthiasm@43:     return path;
matthiasm@43: }