dtw.cpp 5.34 KB
Newer Older
mzed's avatar
mzed committed
1
2
3
4
5
6
7
8
/**
 * @file dtw.cpp
 * RapidLib
 *
 * @author Michael Zbyszynski
 * @date 07 Jun 2017
 * @copyright Copyright © 2017 Goldsmiths. All rights reserved.
 */
mzed's avatar
mzed committed
9
10
11

#include <vector>
#include <cmath>
mzed's avatar
mzed committed
12
#include <cassert>
13
#include <limits>
mzed's avatar
mzed committed
14
#include "dtw.h"
mzed's avatar
mzed committed
15

mzed's avatar
mzed committed
16
17
template<typename T>
dtw<T>::dtw() {};
mzed's avatar
mzed committed
18

mzed's avatar
mzed committed
19
20
template<typename T>
dtw<T>::~dtw() {};
mzed's avatar
mzed committed
21

mzed's avatar
mzed committed
22
23
template<typename T>
inline T dtw<T>::distanceFunction(const std::vector<T> &x, const std::vector<T> &y) {
mzed's avatar
mzed committed
24
    assert(x.size() == y.size());
mzed's avatar
mzed committed
25
    T euclidianDistance = 0;
mzed's avatar
mzed committed
26
    for(int j = 0; j < x.size() ; ++j){
mzed's avatar
mzed committed
27
28
29
30
31
32
        euclidianDistance = euclidianDistance + pow((x[j] - y[j]), 2);
    }
    euclidianDistance = sqrt(euclidianDistance);
    return euclidianDistance;
};

mzed's avatar
mzed committed
33
/* Just returns the cost, doesn't calculate the path */
mzed's avatar
mzed committed
34
35
template<typename T>
T dtw<T>::getCost(const std::vector<std::vector<T> > &seriesX, const std::vector<std::vector<T> > &seriesY) {
mzed's avatar
mzed committed
36
37
38
39
    if (seriesX.size() < seriesY.size()) {
        return getCost(seriesY, seriesX);
    }
    
mzed's avatar
mzed committed
40
41
    costMatrix.clear();
    for (int i = 0; i < seriesX.size(); ++i) {
mzed's avatar
mzed committed
42
        std::vector<T> tempVector;
mzed's avatar
mzed committed
43
44
45
46
47
48
49
        for (int j = 0; j < seriesY.size(); ++j) {
            tempVector.push_back(0);
        }
        costMatrix.push_back(tempVector);
    }
    int maxX = int(seriesX.size()) - 1;
    int maxY = int(seriesY.size()) - 1;
mzed's avatar
mzed committed
50
51
    
    //Calculate values for the first column
mzed's avatar
mzed committed
52
    costMatrix[0][0] = distanceFunction(seriesX[0], seriesY[0]);
mzed's avatar
mzed committed
53
    for (int j = 1; j <= maxY; ++j) {
mzed's avatar
mzed committed
54
        costMatrix[0][j] = costMatrix[0][j - 1] + distanceFunction(seriesX[0], seriesY[j]);
mzed's avatar
mzed committed
55
56
    }
    
mzed's avatar
mzed committed
57
    for (int i = 1; i <= maxX; ++i) {
mzed's avatar
mzed committed
58
        //Bottom row of current column
mzed's avatar
mzed committed
59
        costMatrix[i][0] = costMatrix[i - 1][0] + distanceFunction(seriesX[i], seriesY[0]);
mzed's avatar
mzed committed
60
        for (int j = 1; j <= maxY; ++j) {
mzed's avatar
mzed committed
61
            T minGlobalCost = fmin(costMatrix[i-1][j-1], costMatrix[i][j-1]);
mzed's avatar
mzed committed
62
            costMatrix[i][j] = minGlobalCost + distanceFunction(seriesX[i], seriesY[j]);
mzed's avatar
mzed committed
63
64
        }
    }
mzed's avatar
mzed committed
65
    return costMatrix[maxX][maxY];
mzed's avatar
mzed committed
66
67
};

mzed's avatar
mzed committed
68
template<typename T>
mzed's avatar
mzed committed
69
warpPath dtw<T>::calculatePath(int seriesXsize, int seriesYsize) const {
mzed's avatar
mzed committed
70
    warpPath warpPath;
mzed's avatar
mzed committed
71
72
    int i = seriesXsize - 1;
    int j = seriesYsize - 1;
mzed's avatar
mzed committed
73
74
75
    warpPath.add(i, j);
    
    while ((i > 0) || (j > 0)) {
mzed's avatar
mzed committed
76
77
78
        T diagonalCost = std::numeric_limits<T>::infinity();
        T leftCost = std::numeric_limits<T>::infinity();
        T downCost = std::numeric_limits<T>::infinity();
mzed's avatar
mzed committed
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
        if ((i > 0) && (j > 0)) {
            diagonalCost = costMatrix[i - 1][j - 1];
        }
        if (i > 0) {
            leftCost = costMatrix[i - 1][j];
        }
        if (j > 0) {
            downCost = costMatrix[i][j - 1];
        }
        if ((diagonalCost <= leftCost) && (diagonalCost <= downCost)) {
            i--;
            j--;
        } else if ((leftCost < diagonalCost) && (leftCost < downCost)){
            i--;
        } else if ((downCost < diagonalCost) && (downCost < leftCost)) {
            j--;
        } else if (i <= j) {
            j--;
        } else {
            i--;
        }
        warpPath.add(i, j);
    }
mzed's avatar
mzed committed
102
    return warpPath;
mzed's avatar
mzed committed
103
104
105
};

/* calculates both the cost and the warp path*/
mzed's avatar
mzed committed
106
107
template<typename T>
warpInfo dtw<T>::dynamicTimeWarp(const std::vector<std::vector<T> > &seriesX, const std::vector<std::vector<T> > &seriesY) {
108
    warpInfo info;
mzed's avatar
mzed committed
109
110
    //calculate cost matrix
    info.cost = getCost(seriesX, seriesY);
mzed's avatar
mzed committed
111
    info.path = calculatePath(int(seriesX.size()), int(seriesY.size()));
112
    return info;
mzed's avatar
mzed committed
113
114
}

mzed's avatar
mzed committed
115
/* calculates warp info based on window */
mzed's avatar
mzed committed
116
template<typename T>
mzed's avatar
mzed committed
117
warpInfo dtw<T>::constrainedDTW(const std::vector<std::vector<T> > &seriesX, const std::vector<std::vector<T> > &seriesY, searchWindow<T> window) {
mzed's avatar
mzed committed
118
119
120
121
    
    //initialize cost matrix
    costMatrix.clear();
    for (int i = 0; i < seriesX.size(); ++i) { //TODO: this could be smaller, since most cells are unused
mzed's avatar
mzed committed
122
        std::vector<T> tempVector;
mzed's avatar
mzed committed
123
        for (int j = 0; j < seriesY.size(); ++j) {
mzed's avatar
mzed committed
124
            tempVector.push_back(std::numeric_limits<T>::max());
mzed's avatar
mzed committed
125
126
127
128
129
130
131
        }
        costMatrix.push_back(tempVector);
    }
    int maxX = int(seriesX.size()) - 1;
    int maxY = int(seriesY.size()) - 1;
    
    //fill cost matrix cells based on window
132
133
    for (int currentX = 0; currentX < window.minMaxValues.size(); ++currentX) {
        for (int currentY = window.minMaxValues[currentX].first; currentY <= window.minMaxValues[currentX].second; ++currentY) { //FIXME: should be <= ?
mzed's avatar
mzed committed
134
135
136
            if (currentX == 0 && currentY == 0) { //bottom left cell
                costMatrix[0][0] = distanceFunction(seriesX[0], seriesY[0]);
            } else if (currentX == 0) { //first column
137
                costMatrix[0][currentY] = distanceFunction(seriesX[0], seriesY[currentY]) + costMatrix[0][currentY - 1];
mzed's avatar
mzed committed
138
139
140
            } else if (currentY == 0) { //first row
                costMatrix[currentX][0] = distanceFunction(seriesX[currentX], seriesY[0]) + costMatrix[currentX - 1][0];
            } else {
mzed's avatar
mzed committed
141
                T minGlobalCost = fmin(costMatrix[currentX - 1][currentY], fmin(costMatrix[currentX-1][currentY-1], costMatrix[currentX][currentY-1]));
mzed's avatar
mzed committed
142
143
144
145
146
147
                costMatrix[currentX][currentY] = distanceFunction(seriesX[currentX], seriesY[currentY]) + minGlobalCost;
            }
        }
    }
    warpInfo info;
    info.cost = costMatrix[maxX][maxY];
mzed's avatar
mzed committed
148
    info.path = calculatePath(int(seriesX.size()), int(seriesY.size()));
mzed's avatar
mzed committed
149
150
    return info;
}
mzed's avatar
mzed committed
151
152
153
154

//explicit instantiation
template class dtw<double>;
template class dtw<float>;