Chris@9: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@9: 
Chris@9: /*
Chris@9:     pYIN - A fundamental frequency estimator for monophonic audio
Chris@9:     Centre for Digital Music, Queen Mary, University of London.
Chris@9:     
Chris@9:     This program is free software; you can redistribute it and/or
Chris@9:     modify it under the terms of the GNU General Public License as
Chris@9:     published by the Free Software Foundation; either version 2 of the
Chris@9:     License, or (at your option) any later version.  See the file
Chris@9:     COPYING included with this distribution for more information.
Chris@9: */
Chris@9: 
matthiasm@0: #include "SparseHMM.h"
matthiasm@0: #include <vector>
matthiasm@0: #include <cstdio>
matthiasm@0: #include <iostream>
matthiasm@0: 
matthiasm@0: using std::vector;
matthiasm@0: using std::pair;
matthiasm@0: 
matthiasm@0: const vector<double>
matthiasm@0: SparseHMM::calculateObsProb(const vector<pair<double, double> > data)
matthiasm@0: {
matthiasm@0:     // dummy (virtual?) implementation to be overloaded
matthiasm@0:     return(vector<double>());
matthiasm@0: }
matthiasm@0: 
matthiasm@0: const std::vector<int> 
matthiasm@0: SparseHMM::decodeViterbi(std::vector<vector<double> > obsProb,
matthiasm@0:                          vector<double> *scale) 
matthiasm@0: {
Chris@38:     if (obsProb.size() < 1) {
Chris@38:         return vector<int>();
Chris@38:     }
Chris@38: 
matthiasm@0:     size_t nState = init.size();
matthiasm@0:     size_t nFrame = obsProb.size();
matthiasm@0:     
matthiasm@0:     // check for consistency    
matthiasm@0:     size_t nTrans = transProb.size();
matthiasm@0:     
matthiasm@0:     // declaring variables
matthiasm@0:     std::vector<double> delta = std::vector<double>(nState);
matthiasm@0:     std::vector<double> oldDelta = std::vector<double>(nState);
matthiasm@0:     vector<vector<int> > psi; //  "matrix" of remembered indices of the best transitions
matthiasm@0:     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@0: 
matthiasm@0:     double deltasum = 0;
matthiasm@0: 
matthiasm@0:     // initialise first frame
matthiasm@0:     for (size_t iState = 0; iState < nState; ++iState)
matthiasm@0:     {
matthiasm@0:         oldDelta[iState] = init[iState] * obsProb[0][iState];
matthiasm@0:         // std::cerr << iState << " ----- " << init[iState] << std::endl;
matthiasm@0:         deltasum += oldDelta[iState];
matthiasm@0:     }
matthiasm@0: 
matthiasm@0:     for (size_t iState = 0; iState < nState; ++iState)
matthiasm@0:     {
matthiasm@0:         oldDelta[iState] /= deltasum; // normalise (scale)
matthiasm@0:         // std::cerr << oldDelta[iState] << std::endl;
matthiasm@0:     }
matthiasm@0: 
matthiasm@0:     scale->push_back(1.0/deltasum);
matthiasm@0:     psi.push_back(vector<int>(nState,0));
matthiasm@0: 
matthiasm@0:     // rest of forward step
matthiasm@0:     for (size_t iFrame = 1; iFrame < nFrame; ++iFrame)
matthiasm@0:     {
matthiasm@0:         deltasum = 0;
matthiasm@0:         psi.push_back(vector<int>(nState,0));
matthiasm@0: 
matthiasm@0:         // calculate best previous state for every current state
matthiasm@0:         size_t fromState;
matthiasm@0:         size_t toState;
matthiasm@0:         double currentTransProb;
matthiasm@0:         double currentValue;
matthiasm@0:         
matthiasm@0:         // this is the "sparse" loop
matthiasm@0:         for (size_t iTrans = 0; iTrans < nTrans; ++iTrans)
matthiasm@0:         {
matthiasm@0:             fromState = from[iTrans];
matthiasm@0:             toState = to[iTrans];
matthiasm@0:             currentTransProb = transProb[iTrans];
matthiasm@0:             
matthiasm@0:             currentValue = oldDelta[fromState] * currentTransProb;
matthiasm@0:             if (currentValue > delta[toState])
matthiasm@0:             {
matthiasm@0:                 delta[toState] = currentValue; // will be multiplied by the right obs later!
matthiasm@0:                 psi[iFrame][toState] = fromState;
matthiasm@0:             }            
matthiasm@0:         }
matthiasm@0:         
matthiasm@0:         for (size_t jState = 0; jState < nState; ++jState)
matthiasm@0:         {
matthiasm@0:             delta[jState] *= obsProb[iFrame][jState];
matthiasm@0:             deltasum += delta[jState];
matthiasm@0:         }
matthiasm@0: 
matthiasm@0:         if (deltasum > 0)
matthiasm@0:         {
matthiasm@0:             for (size_t iState = 0; iState < nState; ++iState)
matthiasm@0:             {
matthiasm@0:                 oldDelta[iState] = delta[iState] / deltasum; // normalise (scale)
matthiasm@0:                 delta[iState] = 0;
matthiasm@0:             }
matthiasm@0:             scale->push_back(1.0/deltasum);
matthiasm@0:         } else
matthiasm@0:         {
matthiasm@27:             std::cerr << "WARNING: Viterbi has been fed some zero probabilities, at least they become zero at frame " <<  iFrame << " in combination with the model." << std::endl;
matthiasm@0:             for (size_t iState = 0; iState < nState; ++iState)
matthiasm@0:             {
matthiasm@0:                 oldDelta[iState] = 1.0/nState;
matthiasm@0:                 delta[iState] = 0;
matthiasm@0:             }
matthiasm@0:             scale->push_back(1.0);
matthiasm@0:         }
matthiasm@0:     }
matthiasm@0: 
matthiasm@0:     // initialise backward step
matthiasm@0:     double bestValue = 0;
matthiasm@0:     for (size_t iState = 0; iState < nState; ++iState)
matthiasm@0:     {
matthiasm@0:         double currentValue = oldDelta[iState];
matthiasm@0:         if (currentValue > bestValue)
matthiasm@0:         {
matthiasm@0:             bestValue = currentValue;            
matthiasm@0:             path[nFrame-1] = iState;
matthiasm@0:         }
matthiasm@0:     }
matthiasm@0: 
matthiasm@0:     // rest of backward step
matthiasm@0:     for (int iFrame = nFrame-2; iFrame != -1; --iFrame)
matthiasm@0:     {
matthiasm@0:         path[iFrame] = psi[iFrame+1][path[iFrame+1]];
matthiasm@0:     }
matthiasm@0:     
matthiasm@27:     // for (size_t iState = 0; iState < nState; ++iState)
matthiasm@27:     // {
matthiasm@27:     //     // std::cerr << psi[2][iState] << std::endl;
matthiasm@27:     // }
matthiasm@0:     
matthiasm@0:     return path;
matthiasm@0: }