Chris@235: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@235: 
Chris@235: /*
Chris@235:     Vamp feature extraction plugin using the MATCH audio alignment
Chris@235:     algorithm.
Chris@235: 
Chris@235:     Centre for Digital Music, Queen Mary, University of London.
Chris@236:     Copyright (c) 2007-2020 Simon Dixon, Chris Cannam, and Queen Mary
Chris@235:     University of London, Copyright (c) 2014-2015 Tido GmbH.
Chris@235:     
Chris@235:     This program is free software; you can redistribute it and/or
Chris@235:     modify it under the terms of the GNU General Public License as
Chris@235:     published by the Free Software Foundation; either version 2 of the
Chris@235:     License, or (at your option) any later version.  See the file
Chris@235:     COPYING included with this distribution for more information.
Chris@235: */
Chris@235: 
Chris@235: #include "FullDTW.h"
Chris@235: 
Chris@235: #include <stdexcept>
Chris@238: #include <algorithm>
Chris@242: #include <iostream>
Chris@242: #include <cmath>
Chris@242: 
Chris@242: //#define DEBUG_DTW 1
Chris@235: 
Chris@235: FullDTW::FullDTW(Parameters params, DistanceMetric::Parameters dparams) :
Chris@235:     m_params(params),
Chris@235:     m_metric(dparams)
Chris@235: {
Chris@235: }
Chris@235: 
Chris@235: pathcost_t
Chris@235: FullDTW::choose(CostOption x, CostOption y, CostOption d) {
Chris@235:     if (x.present && y.present) {
Chris@235:         if (!d.present) {
Chris@235:             throw std::logic_error("if x & y both exist, so must diagonal");
Chris@235:         }
Chris@235:         return std::min(std::min(x.cost, y.cost), d.cost);
Chris@235:     } else if (x.present) {
Chris@235:         return x.cost;
Chris@235:     } else if (y.present) {
Chris@235:         return y.cost;
Chris@235:     } else {
Chris@235:         return 0.0;
Chris@235:     }
Chris@235: }
Chris@235: 
Chris@235: pathcostmat_t
Chris@235: FullDTW::costSequences(const featureseq_t &s1, const featureseq_t &s2) {
Chris@235: 
Chris@235:     pathcostmat_t costs(s1.size(), pathcostvec_t(s2.size(), 0));
Chris@235: 
Chris@242: #ifdef DEBUG_DTW
Chris@242:     std::cerr << "Distance matrix:" << std::endl;
Chris@242: #endif
Chris@242: 
Chris@235:     for (size_t j = 0; j < s1.size(); ++j) {
Chris@235:         for (size_t i = 0; i < s2.size(); ++i) {
Chris@235:             distance_t dist = m_metric.calcDistance(s1[j], s2[i]);
Chris@242: #ifdef DEBUG_DTW
Chris@244:             std::cerr << distance_print_t(dist) << " ";
Chris@242: #endif
Chris@242:             
Chris@235:             if (i == 0 && m_params.subsequence) {
Chris@235:                 costs[j][i] = dist;
Chris@235:             } else {
Chris@235:                 costs[j][i] = choose
Chris@235:                     (
Chris@235:                         { j > 0,
Chris@242:                           j > 0 ?
Chris@242:                           (dist + costs[j-1][i])
Chris@242:                           : 0
Chris@235:                         },
Chris@235:                         { i > 0,
Chris@242:                           i > 0 ?
Chris@242:                           (dist + costs[j][i-1])
Chris@242:                           : 0
Chris@235:                         },
Chris@242:                         { (j > 0 && i > 0),
Chris@242:                           (j > 0 && i > 0) ?
Chris@242:                           pathcost_t(round(m_params.diagonalWeight * dist)
Chris@242:                                      + costs[j-1][i-1])
Chris@242:                           : 0
Chris@235:                         });
Chris@235:             }
Chris@235:         }
Chris@242: #ifdef DEBUG_DTW
Chris@242:         std::cerr << "\n";
Chris@242: #endif
Chris@235:     }
Chris@235: 
Chris@235:     return costs;
Chris@235: }
Chris@235: 
Chris@235: std::vector<size_t>
Chris@235: FullDTW::align(const featureseq_t &s1, const featureseq_t &s2) {
Chris@235: 
Chris@235:     // Return the index into s1 for each element in s2
Chris@235:         
Chris@235:     std::vector<size_t> alignment(s2.size(), 0);
Chris@235: 
Chris@235:     if (s1.empty() || s2.empty()) {
Chris@235:         return alignment;
Chris@235:     }
Chris@235: 
Chris@235:     auto costs = costSequences(s1, s2);
Chris@235: 
Chris@235: #ifdef DEBUG_DTW
Chris@242:     std::cerr << "Path cost matrix:" << std::endl;
Chris@242:     for (auto v: costs) {
Chris@235:         for (auto x: v) {
Chris@242:             std::cerr << x << " ";
Chris@235:         }
Chris@242:         std::cerr << "\n";
Chris@235:     }
Chris@235: #endif
Chris@235: 
Chris@235:     size_t j = s1.size() - 1;
Chris@235:     size_t i = s2.size() - 1;
Chris@235: 
Chris@235:     if (m_params.subsequence) {
Chris@235:         pathcost_t min = 0.0;
Chris@235:         size_t minidx = 0;
Chris@235:         for (size_t j = 0; j < s1.size(); ++j) {
Chris@235:             if (j == 0 || costs[j][i] < min) {
Chris@235:                 min = costs[j][i];
Chris@235:                 minidx = j;
Chris@235:             }
Chris@235:         }
Chris@235:         j = minidx;
Chris@235: #ifdef DEBUG_DTW
Chris@242:         std::cerr << "Lowest cost at end of subsequence = " << min
Chris@242:                   << " at index " << j << ", tracking back from there"
Chris@242:                   << std::endl;
Chris@235: #endif
Chris@235:     }
Chris@242: 
Chris@242:     while (i > 0 || j > 0) {
Chris@235: 
Chris@235:         alignment[i] = j;
Chris@242:         
Chris@242:         if (i == 0) {
Chris@242:             if (m_params.subsequence) {
Chris@242:                 break;
Chris@242:             } else {
Chris@242:                 --j;
Chris@242:                 continue;
Chris@242:             }
Chris@242:         }
Chris@242: 
Chris@242:         if (j == 0) {
Chris@242:             --i;
Chris@242:             continue;
Chris@242:         }
Chris@242: 
Chris@235:         pathcost_t a = costs[j-1][i];
Chris@235:         pathcost_t b = costs[j][i-1];
Chris@235:         pathcost_t both = costs[j-1][i-1];
Chris@235: 
Chris@235:         if (a < b) {
Chris@235:             --j;
Chris@235:             if (both <= a) {
Chris@235:                 --i;
Chris@235:             }
Chris@235:         } else {
Chris@235:             --i;
Chris@235:             if (both <= b) {
Chris@235:                 --j;
Chris@235:             }
Chris@235:         }
Chris@235:     }
Chris@235: 
Chris@235:     if (m_params.subsequence) {
Chris@235:         alignment[0] = j;
Chris@235:     }
Chris@242: 
Chris@242: #ifdef DEBUG_DTW
Chris@246:     std::cerr << "Costed path:" << std::endl;
Chris@242:     pathcost_t prevcost = 0;
Chris@242:     int indent = 0;
Chris@242:     size_t prevj = 0;
Chris@242:     for (size_t i = 0; i < alignment.size(); ++i) {
Chris@242:         size_t j = alignment[i];
Chris@242:         pathcost_t cost = costs[j][i];
Chris@242:         if (prevcost > 1000) indent += 5;
Chris@242:         else if (prevcost > 100) indent += 4;
Chris@242:         else if (prevcost > 10) indent += 3;
Chris@242:         else if (prevcost > 0) indent += 2;
Chris@242:         if (j > prevj) {
Chris@242:             while (j > prevj) {
Chris@242:                 std::cerr << "\n";
Chris@242:                 ++prevj;
Chris@242:             }
Chris@242:             for (int k = 0; k < indent; ++k) std::cerr << " ";
Chris@242:         } else {
Chris@242:             std::cerr << " ";
Chris@242:         }
Chris@242:         std::cerr << cost;
Chris@242:         prevcost = cost;
Chris@242:         if (prevcost == 0) prevcost = 1;
Chris@242:     }
Chris@242:     std::cerr << "\n";
Chris@246: 
Chris@246:     std::cerr << "Alignment:" << std::endl;
Chris@246:     for (size_t i = 0; i < alignment.size(); ++i) {
Chris@246:         std::cerr << i << " -> " << alignment[i] << "\n";
Chris@246:     }
Chris@242: #endif
Chris@242:     
Chris@235:     return alignment;
Chris@235: }