changeset 242:320efced24c1 subsequence

Fix indexing bug, + add optional debug output

--- a/src/FullDTW.cpp	Thu Jul 16 15:17:02 2020 +0100
+++ b/src/FullDTW.cpp	Fri Jul 17 14:25:36 2020 +0100
@@ -19,6 +19,10 @@
 
 #include <stdexcept>
 #include <algorithm>
+#include <iostream>
+#include <cmath>
+
+//#define DEBUG_DTW 1
 
 FullDTW::FullDTW(Parameters params, DistanceMetric::Parameters dparams) :
     m_params(params),
@@ -47,25 +51,43 @@
 
     pathcostmat_t costs(s1.size(), pathcostvec_t(s2.size(), 0));
 
+#ifdef DEBUG_DTW
+    std::cerr << "Distance matrix:" << std::endl;
+#endif
+
     for (size_t j = 0; j < s1.size(); ++j) {
         for (size_t i = 0; i < s2.size(); ++i) {
             distance_t dist = m_metric.calcDistance(s1[j], s2[i]);
+#ifdef DEBUG_DTW
+            std::cerr << pathcost_t(dist) << " "; // dist may be char, want a number!
+#endif
+            
             if (i == 0 && m_params.subsequence) {
                 costs[j][i] = dist;
             } else {
                 costs[j][i] = choose
                     (
                         { j > 0,
-                          j > 0 ? dist + costs[j-1][i] : 0
+                          j > 0 ?
+                          (dist + costs[j-1][i])
+                          : 0
                         },
                         { i > 0,
-                          i > 0 ? dist + costs[j][i-1] : 0
+                          i > 0 ?
+                          (dist + costs[j][i-1])
+                          : 0
                         },
-                        { j > 0 && i > 0,
-                          j > 0 && i > 0 ? dist + costs[j-1][i-1] : 0
+                        { (j > 0 && i > 0),
+                          (j > 0 && i > 0) ?
+                          pathcost_t(round(m_params.diagonalWeight * dist)
+                                     + costs[j-1][i-1])
+                          : 0
                         });
             }
         }
+#ifdef DEBUG_DTW
+        std::cerr << "\n";
+#endif
     }
 
     return costs;
@@ -85,12 +107,12 @@
     auto costs = costSequences(s1, s2);
 
 #ifdef DEBUG_DTW
-    SVCERR << "Cost matrix:" << endl;
-    for (auto v: cost) {
+    std::cerr << "Path cost matrix:" << std::endl;
+    for (auto v: costs) {
         for (auto x: v) {
-            SVCERR << x << " ";
+            std::cerr << x << " ";
         }
-        SVCERR << "\n";
+        std::cerr << "\n";
     }
 #endif
 
@@ -108,15 +130,30 @@
         }
         j = minidx;
 #ifdef DEBUG_DTW
-        SVCERR << "Lowest cost at end of subsequence = " << min
-               << " at index " << j << ", tracking back from there" << endl;
+        std::cerr << "Lowest cost at end of subsequence = " << min
+                  << " at index " << j << ", tracking back from there"
+                  << std::endl;
 #endif
     }
-        
-    while (i > 0 && (j > 0 || m_params.subsequence)) {
+
+    while (i > 0 || j > 0) {
 
         alignment[i] = j;
-            
+        
+        if (i == 0) {
+            if (m_params.subsequence) {
+                break;
+            } else {
+                --j;
+                continue;
+            }
+        }
+
+        if (j == 0) {
+            --i;
+            continue;
+        }
+
         pathcost_t a = costs[j-1][i];
         pathcost_t b = costs[j][i-1];
         pathcost_t both = costs[j-1][i-1];
@@ -137,6 +174,34 @@
     if (m_params.subsequence) {
         alignment[0] = j;
     }
-        
+
+#ifdef DEBUG_DTW
+    std::cerr << "Alignment:" << std::endl;
+    pathcost_t prevcost = 0;
+    int indent = 0;
+    size_t prevj = 0;
+    for (size_t i = 0; i < alignment.size(); ++i) {
+        size_t j = alignment[i];
+        pathcost_t cost = costs[j][i];
+        if (prevcost > 1000) indent += 5;
+        else if (prevcost > 100) indent += 4;
+        else if (prevcost > 10) indent += 3;
+        else if (prevcost > 0) indent += 2;
+        if (j > prevj) {
+            while (j > prevj) {
+                std::cerr << "\n";
+                ++prevj;
+            }
+            for (int k = 0; k < indent; ++k) std::cerr << " ";
+        } else {
+            std::cerr << " ";
+        }
+        std::cerr << cost;
+        prevcost = cost;
+        if (prevcost == 0) prevcost = 1;
+    }
+    std::cerr << "\n";
+#endif
+    
     return alignment;
 }