FunctionDynamicTimeWarping
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © RicardoSantos

//@version=5

// @description 
//  "In time series analysis, dynamic time warping (DTW) is an algorithm for 
//   measuring similarity between two temporal sequences, which may vary in 
//   speed. For instance, similarities in walking could be detected using DTW, 
//   even if one person was walking faster than the other, or if there were 
//   accelerations and decelerations during the course of an observation. 
//   DTW has been applied to temporal sequences of video, audio, and graphics 
//   data — indeed, any data that can be turned into a linear sequence can be 
//   analyzed with DTW. A well-known application has been automatic speech 
//   recognition, to cope with different speaking speeds. Other applications 
//   include speaker recognition and online signature recognition. 
//   It can also be used in partial shape matching applications." 
//
//   "Dynamic time warping is used in finance and econometrics to assess the 
//   quality of the prediction versus real-world data."
//
//    ~~ wikipedia
//
// reference:
//      https://en.wikipedia.org/wiki/Dynamic_time_warping
//      https://towardsdatascience.com/dynamic-time-warping-3933f25fcdd
//      https://github.com/shunsukeaihara/pydtw/blob/master/pydtw/dtw.pyx
//
library(title='FunctionDynamicTimeWarping')


d_argmin (float a, float b, float c) => a <= b or a <= c ? 0 : (b <= c ? 1 : 2)
distance (float a, float b)          => math.abs(a - b)


// @function Dynamic Time Warping procedure.
// @param a array<float>, data series.
// @param b array<float>, data series.
// @param w  int         , minimum window size.
// @returns
//      matrix<float> optimum match matrix.
export cost_matrix (
  array<float> a, 
  array<float> b, 
  int           w
  ) => //{
    int _size_a = array.size(a) + 1
    int _size_b = array.size(b) + 1
    //
    float         _inf = 1.0e308 // near enough to infinity :o)
    int           _w   = math.max(w, math.abs(_size_a - _size_b)) // adapt window size
    matrix<float> _dtw = matrix.new<float>(_size_a, _size_b, _inf)
    matrix.set(_dtw, 0, 0, 0.0)
    for _i = 1 to _size_a - 1
        for _j = math.max(1, _i - _w) to math.min(_size_b - 1, _i + _w)
            matrix.set(_dtw, _i, _j, 0)
    for _i = 1 to _size_a - 1
        for _j = math.max(1, _i - _w) to math.min(_size_b - 1, _i + _w)
            float _cost = distance(array.get(a, _i-1), array.get(b, _j-1))
            matrix.set(_dtw, _i, _j, _cost + math.min(matrix.get(_dtw, _i - 1, _j), matrix.get(_dtw, _i, _j - 1), matrix.get(_dtw, _i - 1, _j - 1)))
    matrix.remove_col(_dtw, 0)
    matrix.remove_row(_dtw, 0)
    _dtw
// example:
a = array.from(1.0, 2, 3, 3, 5)
b = array.from(1.0, 2, 2, 2, 2, 2, 2, 4)
c = cost_matrix(a, b, 3)
var la = label.new(bar_index, 0.0, '')
// if barstate.islast
//     label.set_x(la, bar_index)
//     label.set_text(la, str.tostring(c))
//}

// @function perform a backtrace on the cost matrix and retrieve optimal paths and cost between arrays.
// @param M matrix<float>, cost matrix.
// @returns
//      tuple:
//      array<int> aligned 1st array of indices.
//      array<int> aligned 2nd array of indices.
//      float      final cost.
// reference:
//      https://github.com/shunsukeaihara/pydtw/blob/master/pydtw/dtw.pyx
export traceback (
  matrix<float> M
  ) => //{
    int _i =    matrix.rows(M) - 1
    int _j = matrix.columns(M) - 1
    float _cost = matrix.get(M, _i, _j)
    array<int> _a = array.new<int>(1, _i)
    array<int> _b = array.new<int>(1, _j)
    int _match = int(na)
    while _i > 0 and _j > 0
        _match := d_argmin(matrix.get(M, _i - 1, _j - 1), matrix.get(M, _i - 1, _j), matrix.get(M, _i, _j - 1))
        switch _match
            0 => _i -= 1 , _j -= 1
            1 => _i -= 1
            =>   _j -= 1
        array.push(_a, _i), array.push(_b, _j)
    array.reverse(_a), array.reverse(_b)
    [_a, _b, _cost]
//}

// @function report ordered arrays, cost and cost matrix.
// @param a array<float>, data series.
// @param b array<float>, data series.
// @param w  int        , minimum window size.
// @returns
//      string report.
export report (
  array<float> a, 
  array<float> b, 
  int w
  ) => //{
    matrix<float> _cost_mat = cost_matrix(a, b, w)
    [_align_a, _align_b, _cost] = traceback(_cost_mat)
    str.format('cost: {0}\nalign a: {1}\nalign b: {2}\ncost matrix:\n{3}', str.tostring(_cost), str.tostring(_align_a), str.tostring(_align_b), str.tostring(_cost_mat))
//
s = report(a, b, 3)
if barstate.islast
    label.set_x(la, bar_index)
    label.set_text(la, s)
//}



