Dynamic Time Warping

I've read about Dynamic Time Warping yesterday and here is how I would do it in C++.

TimeSeriesUtils.hpp

#ifndef TSUTIL_H
#define TSUTIL_H

#include <vector>
#include <cstdint>

namespace TimeSeriesUtils {

    double CalculateEuclideanDistance(double x, double y);

    double CalculateDynamicTimeWarpedDistance(std::vector<double> t0, std::vector<double> t1);

    int CalculateLongestCommonSubsequence(std::vector<double> t0, std::vector<double> t1);

    int CalculateTimeWarpedLongestCommonSubsequence(std::vector<double> t0, std::vector<double> t1);

}

#endif

TimeSeriesUtils.cpp

#include <vector>
#include <algorithm>
#include "../include/TimeSeriesUtils.hpp"

namespace TimeSeriesUtils {

    double CalculateEuclideanDistance(double x, double y) {
        return std::sqrt(std::pow((x - y), 2));
    }   

    double CalculateDynamicTimeWarpedDistance(std::vector<double> t0, std::vector<double> t1) {

        size_t m = t0.size();
        size_t n = t1.size();

        // Allocate the Matrix to work on:
        std::vector<std::vector<double>> cost(m, std::vector<double>(n));

        cost[0][0] = CalculateEuclideanDistance(t0[0], t1[0]);

        // Calculate the first row:
        for (int i = 1; i < m; i++) {
            cost[i][0] = cost[i - 1][0] + CalculateEuclideanDistance(t0[i], t1[0]);
        }

        // Calculate the first column:
        for (int j = 1; j < n; j++) {
            cost[0][j] = cost[0][j - 1] + CalculateEuclideanDistance(t0[0], t1[j]);
        }

        // Fill the matrix:
        for (int i = 1; i < m; i++) {
            for (int j = 1; j < n; j++) {
                cost[i][j] = std::min(cost[i - 1][j], std::min(cost[i][j - 1], cost[i - 1][j - 1])) 
                    + CalculateEuclideanDistance(t0[i], t1[j]);
            }
        }

        return cost[m-1][n-1];
    }

    int CalculateLongestCommonSubsequence(std::vector<double> t0, std::vector<double> t1) {

        size_t m = t0.size() + 1;
        size_t n = t1.size() + 1;

        // Allocate the Matrix to work on:
        std::vector<std::vector<int>> cost(m, std::vector<int>(n));

        cost[0][0] = 0;

        // Calculate the first row:
        for (int i = 1; i < m; i++) {
            cost[i][0] = 0;
        }

        // Calculate the first column:
        for (int j = 1; j < n; j++) {
            cost[0][j] = 0;
        }


        for(int i = 1; i < m; i++) {
            for(int j = 1; j < n; j++) {
                 if(std::abs(t0[i-1] - t1[j-1]) <= 1e-10) {
                    cost[i][j] = 1 + cost[i-1][j-1];
                } else {
                    cost[i][j] = std::max(cost[i][j-1], cost[i-1][j]);
                }
            }
        }

        return cost[m-1][n-1];
    }

    int CalculateTimeWarpedLongestCommonSubsequence(std::vector<double> t0, std::vector<double> t1) {

        size_t m = t0.size() + 1;
        size_t n = t1.size() + 1;

        // Allocate the Matrix to work on:
        std::vector<std::vector<int>> cost(m, std::vector<int>(n));

        cost[0][0] = 0;

        // Calculate the first row:
        for (int i = 1; i < m; i++) {
            cost[i][0] = 0;
        }

        // Calculate the first column:
        for (int j = 1; j < n; j++) {
            cost[0][j] = 0;
        }


        for(int i = 1; i < m; i++) {
            for(int j = 1; j < n; j++) {
                 if(std::abs(t0[i-1] - t1[j-1]) <= 1e-10) {
                    cost[i][j] = 1 + std::max(cost[i][j-1], std::max(cost[i-1][j],  cost[i-1][j-1]));
                } else {
                    cost[i][j] = std::max(cost[i][j-1], cost[i-1][j]);
                }
            }
        }

        return cost[m-1][n-1];
    }
}

usage

#include "../include/TimeSeriesUtils.hpp"
#include <vector>
#include <cstdint>
#include <catch.hpp>

namespace {

    TEST_CASE("DynamicTimeWarpedDistanceTest")
    {
        std::vector<double> timeSeries0;
        std::vector<double> timeSeries1;

        timeSeries0.push_back(4);
        timeSeries0.push_back(4);
        timeSeries0.push_back(5);
        timeSeries0.push_back(5);
        timeSeries0.push_back(6);
        timeSeries0.push_back(6);
        timeSeries0.push_back(7);
        timeSeries0.push_back(7);

        timeSeries1.push_back(23);
        timeSeries1.push_back(4);
        timeSeries1.push_back(5);
        timeSeries1.push_back(6);
        timeSeries1.push_back(7);

        double result = TimeSeriesUtils::CalculateDynamicTimeWarpedDistance(timeSeries0, timeSeries1);

        REQUIRE(result == Approx(19).epsilon(0.01));
    }

    TEST_CASE("LongestCommonSubsequenceTest")
    {
        std::vector<double> timeSeries0;
        std::vector<double> timeSeries1;

        timeSeries0.push_back(4);
        timeSeries0.push_back(4);
        timeSeries0.push_back(5);
        timeSeries0.push_back(5);
        timeSeries0.push_back(6);
        timeSeries0.push_back(6);
        timeSeries0.push_back(7);
        timeSeries0.push_back(7);

        timeSeries1.push_back(4);
        timeSeries1.push_back(5);
        timeSeries1.push_back(6);
        timeSeries1.push_back(7);

        int result = TimeSeriesUtils::CalculateLongestCommonSubsequence(timeSeries0, timeSeries1);

        REQUIRE(result == 4);
    }
}

notes

You can obtain the sourcecode from my github page at:

To find the minimum warp path between two sequences, you can use a greedy strategy from cost[m][n] to cost[0][0] by picking the minimum of the three adjacenct values.

How to contribute

One of the easiest ways to contribute is to participate in discussions. You can also contribute by submitting pull requests.

General feedback and discussions?

Do you have questions or feedback on this article? Please create an issue on the GitHub issue tracker.

Something is wrong or missing?

There may be something wrong or missing in this article. If you want to help fixing it, then please make a Pull Request to this file on GitHub.