# HG changeset patch
# User Chris Cannam
# Date 1593614086 -3600
# Node ID 8fa98f89eda8f7616246cae2c51148dc23e26870
# Parent  5de2b710cfae265e2253547a641c77be66cf1977
Add subsequence DTW (not yet in use)

diff -r 5de2b710cfae -r 8fa98f89eda8 align/DTW.h
--- a/align/DTW.h	Wed Jul 01 11:41:07 2020 +0100
+++ b/align/DTW.h	Wed Jul 01 15:34:46 2020 +0100
@@ -27,54 +27,22 @@
     DTW(std::function<double(const Value &, const Value &)> distanceMetric) :
         m_metric(distanceMetric) { }
 
-    std::vector<size_t> alignSeries(std::vector<Value> s1,
-                                    std::vector<Value> s2) {
+    /**
+     * Align the sequence s2 against the whole of the sequence s1,
+     * returning the index into s1 for each element in s2.
+     */
+    std::vector<size_t> alignSequences(std::vector<Value> s1,
+                                       std::vector<Value> s2) {
+        return align(s1, s2, false);
+    }
 
-        // Return the index into s2 for each element in s1
-        
-        std::vector<size_t> alignment(s1.size(), 0);
-
-        if (s1.empty() || s2.empty()) {
-            return alignment;
-        }
-
-        auto costs = costSeries(s1, s2);
-
-#ifdef DEBUG_DTW
-        SVCERR << "Cost matrix:" << endl;
-        for (auto v: costs) {
-            for (auto x: v) {
-                SVCERR << x << " ";
-            }
-            SVCERR << "\n";
-        }
-#endif
-        
-        size_t j = s1.size() - 1;
-        size_t i = s2.size() - 1;
-
-        while (j > 0 && i > 0) {
-
-            alignment[j] = i;
-            
-            cost_t a = costs[j-1][i];
-            cost_t b = costs[j][i-1];
-            cost_t both = costs[j-1][i-1];
-
-            if (a < b) {
-                --j;
-                if (both <= a) {
-                    --i;
-                }
-            } else {
-                --i;
-                if (both <= b) {
-                    --j;
-                }
-            }
-        }
-
-        return alignment;
+    /**
+     * Align the sequence sub against the best-matching subsequence of
+     * s, returning the index into s for each element in sub.
+     */
+    std::vector<size_t> alignSubsequence(std::vector<Value> s,
+                                         std::vector<Value> sub) {
+        return align(s, sub, true);
     }
 
 private:
@@ -102,8 +70,9 @@
         }
     }
 
-    std::vector<std::vector<cost_t>> costSeries(std::vector<Value> s1,
-                                                std::vector<Value> s2) {
+    std::vector<std::vector<cost_t>> costSequences(std::vector<Value> s1,
+                                                   std::vector<Value> s2,
+                                                   bool subsequence) {
 
         std::vector<std::vector<cost_t>> costs
             (s1.size(), std::vector<cost_t>(s2.size(), 0.0));
@@ -111,22 +80,97 @@
         for (size_t j = 0; j < s1.size(); ++j) {
             for (size_t i = 0; i < s2.size(); ++i) {
                 cost_t c = m_metric(s1[j], s2[i]);
-                costs[j][i] = choose
-                    (
-                        { j > 0,
-                          j > 0 ? c + costs[j-1][i] : 0.0
-                        },
-                        { i > 0,
-                          i > 0 ? c + costs[j][i-1] : 0.0
-                        },
-                        { j > 0 && i > 0,
-                          j > 0 && i > 0 ? c + costs[j-1][i-1] : 0.0
-                        });
+                if (i == 0 && subsequence) {
+                    costs[j][i] = c;
+                } else {
+                    costs[j][i] = choose
+                        (
+                            { j > 0,
+                              j > 0 ? c + costs[j-1][i] : 0.0
+                            },
+                            { i > 0,
+                              i > 0 ? c + costs[j][i-1] : 0.0
+                            },
+                            { j > 0 && i > 0,
+                              j > 0 && i > 0 ? c + costs[j-1][i-1] : 0.0
+                            });
+                }
             }
         }
 
         return costs;
     }
+
+    std::vector<size_t> align(const std::vector<Value> &s1,
+                              const std::vector<Value> &s2,
+                              bool subsequence) {
+
+        // Return the index into s1 for each element in s2
+        
+        std::vector<size_t> alignment(s2.size(), 0);
+
+        if (s1.empty() || s2.empty()) {
+            return alignment;
+        }
+
+        auto costs = costSequences(s1, s2, subsequence);
+
+#ifdef DEBUG_DTW
+        SVCERR << "Cost matrix:" << endl;
+        for (auto v: cost) {
+            for (auto x: v) {
+                SVCERR << x << " ";
+            }
+            SVCERR << "\n";
+        }
+#endif
+
+        size_t j = s1.size() - 1;
+        size_t i = s2.size() - 1;
+
+        if (subsequence) {
+            cost_t min = 0.0;
+            size_t minidx = 0;
+            for (size_t j = 0; j < s1.size(); ++j) {
+                if (j == 0 || costs[j][i] < min) {
+                    min = costs[j][i];
+                    minidx = j;
+                }
+            }
+            j = minidx;
+#ifdef DEBUG_DTW
+            SVCERR << "Lowest cost at end of subsequence = " << min
+                   << " at index " << j << ", tracking back from there" << endl;
+#endif
+        }
+        
+        while (i > 0 && (j > 0 || subsequence)) {
+
+            alignment[i] = j;
+            
+            cost_t a = costs[j-1][i];
+            cost_t b = costs[j][i-1];
+            cost_t both = costs[j-1][i-1];
+
+            if (a < b) {
+                --j;
+                if (both <= a) {
+                    --i;
+                }
+            } else {
+                --i;
+                if (both <= b) {
+                    --j;
+                }
+            }
+        }
+
+        if (subsequence) {
+            alignment[0] = j;
+        }
+        
+        return alignment;
+    }
 };
 
 class MagnitudeDTW
@@ -134,9 +178,14 @@
 public:
     MagnitudeDTW() : m_dtw(metric) { }
 
-    std::vector<size_t> alignSeries(std::vector<double> s1,
-                                    std::vector<double> s2) {
-        return m_dtw.alignSeries(s1, s2);
+    std::vector<size_t> alignSequences(std::vector<double> s1,
+                                       std::vector<double> s2) {
+        return m_dtw.alignSequences(s1, s2);
+    }
+
+    std::vector<size_t> alignSubsequence(std::vector<double> s,
+                                         std::vector<double> sub) {
+        return m_dtw.alignSubsequence(s, sub);
     }
 
 private:
@@ -163,9 +212,14 @@
 
     RiseFallDTW() : m_dtw(metric) { }
 
-    std::vector<size_t> alignSeries(std::vector<Value> s1,
-                                    std::vector<Value> s2) {
-        return m_dtw.alignSeries(s1, s2);
+    std::vector<size_t> alignSequences(std::vector<Value> s1,
+                                       std::vector<Value> s2) {
+        return m_dtw.alignSequences(s1, s2);
+    }
+
+    std::vector<size_t> alignSubsequence(std::vector<Value> s,
+                                         std::vector<Value> sub) {
+        return m_dtw.alignSubsequence(s, sub);
     }
 
 private:
diff -r 5de2b710cfae -r 8fa98f89eda8 align/TransformDTWAligner.cpp
--- a/align/TransformDTWAligner.cpp	Wed Jul 01 11:41:07 2020 +0100
+++ b/align/TransformDTWAligner.cpp	Wed Jul 01 15:34:46 2020 +0100
@@ -297,19 +297,19 @@
     
     for (int i = 0; in_range_for(alignment, i); ++i) {
 
-        // DTW returns "the index into s2 for each element in s1"
-        sv_frame_t refFrame = refFrames[i];
-
-        if (!in_range_for(otherFrames, alignment[i])) {
+        // DTW returns "the index into s1 for each element in s2"
+        sv_frame_t alignedFrame = otherFrames[i];
+        
+        if (!in_range_for(refFrames, alignment[i])) {
             SVCERR << "TransformDTWAligner::makePath: Internal error: "
-                   << "DTW maps index " << i << " in reference frame vector "
-                   << "(size " << refFrames.size() << ") onto index "
-                   << alignment[i] << " in other frame vector "
-                   << "(only size " << otherFrames.size() << ")" << endl;
+                   << "DTW maps index " << i << " in other frame vector "
+                   << "(size " << otherFrames.size() << ") onto index "
+                   << alignment[i] << " in ref frame vector "
+                   << "(only size " << refFrames.size() << ")" << endl;
             continue;
         }
             
-        sv_frame_t alignedFrame = otherFrames[alignment[i]];
+        sv_frame_t refFrame = refFrames[alignment[i]];
         path.add(PathPoint(alignedFrame, refFrame));
     }
 
@@ -361,7 +361,7 @@
                << "]: serialising DTW to avoid over-allocation" << endl;
 #endif
         QMutexLocker locker(&m_dtwMutex);
-        alignment = dtw.alignSeries(s1, s2);
+        alignment = dtw.alignSequences(s1, s2);
     }
 
 #ifdef DEBUG_TRANSFORM_DTW_ALIGNER
@@ -450,7 +450,7 @@
                << "]: serialising DTW to avoid over-allocation" << endl;
 #endif
         QMutexLocker locker(&m_dtwMutex);
-        alignment = dtw.alignSeries(s1, s2);
+        alignment = dtw.alignSequences(s1, s2);
     }
 
 #ifdef DEBUG_TRANSFORM_DTW_ALIGNER