//***************************************************************************
//* Copyright (c) 2015 Saint Petersburg State University
//* Copyright (c) 2011-2014 Saint Petersburg Academic University
//* All Rights Reserved
//* See file LICENSE for details.
//***************************************************************************

#pragma once

#include <cmath>

namespace math {

template<typename T>
static T MedianOf3(T u, T v, T w) {
    /* Median(u,v,w): */
    if ((u <= v && v <= w) ||
        (u >= v && v >= w))
        return v;
    if ((u <= w && w <= v) ||
        (u >= w && w >= v))
        return w;

    /* else */ return u;
}

/* Return (Index-1) of  median(u,v,w) , i.e.,
-1 : u
0 : v
1 : w
*/
template<typename T>
static int IndexOfMedianOf3(T u, T v, T w) {
    if ((u <= v && v <= w) ||
        (u >= v && v >= w))
        return 0;
    if ((u <= w && w <= v) ||
        (u >= w && w >= v))
        return 1;

    /* else */ return -1;
}

enum class SmoothEndRule {
    No,
    Copy,
    Tukey
};

template<typename T>
static bool SmoothEndStep(const T *x, T *y, size_t n, SmoothEndRule end_rule) {
    switch (end_rule) {
        default:
        case SmoothEndRule::No:
            return false;
        case SmoothEndRule::Copy:
            y[0] = x[0];
            y[n - 1] = x[n - 1];
            return false;
        case SmoothEndRule::Tukey: {
            bool chg = false;
            y[0] = MedianOf3(3 * y[1] - 2 * y[2], x[0], y[1]);
            chg = chg || (y[0] != x[0]);
            y[n - 1] = MedianOf3(y[n - 2], x[n - 1], 3 * y[n - 2] - 2 * y[n - 3]);
            chg = chg || (y[n - 1] != x[n - 1]);
            return chg;
        }
    }
}

template<typename T>
static bool Smooth3(const T *x, T *y, size_t n, SmoothEndRule end_rule) {
    // y[] := Running Median of three (x) = "3 (x[])" with "copy ends"
    // ---  return chg := ( y != x )
    bool chg = false;

    for (size_t i = 1; i < n - 1; i++) {
        int j = IndexOfMedianOf3(x[i - 1], x[i], x[i + 1]);
        y[i] = x[(int) i + j];
        chg = chg || j;
    }

    chg |= SmoothEndStep(x, y, n, end_rule);

    return chg;
}

template<typename T>
static size_t Smooth3R(const T *x, T *y, T *z, size_t n, SmoothEndRule end_rule) {
    // y[] := "3R"(x) ; 3R = Median of three, repeated until convergence
    size_t iter;
    bool chg;

    iter = chg = Smooth3(x, y, n, SmoothEndRule::Copy);

    while (chg) {
        if ((chg = Smooth3(y, z, n, SmoothEndRule::No))) {
            iter += 1;
            for (size_t i = 1; i < n - 1; i++)
                y[i] = z[i];
        }
    }

    chg |= SmoothEndStep(x, y, n, end_rule);

    return (iter ? iter : chg);
    /* = 0   <==>  only one "3" w/o any change
       = 1   <==>  either ["3" w/o change + endchange]
       or   [two "3"s, 2nd w/o change  ] */
}

template<typename T>
static bool SplitTest(const T *x, size_t i) {
    // Split test:
    //  Are we at a /-\ or \_/ location => split should be made ?

    if (x[i] != x[i + 1])
        return false;

    if ((x[i - 1] <= x[i] && x[i + 1] <= x[i + 2]) ||
        (x[i - 1] >= x[i] && x[i + 1] >= x[i + 2]))
        return false;

    /* else */ return true;
}

template<typename T>
static bool SmoothSplit3(const T *x, T *y, size_t n, bool do_ends) {
    // y[] := S(x[])  where S() = "sm_split3"
    bool chg = false;

    for (size_t i = 0; i < n; i++)
        y[i] = x[i];

    if (do_ends && SplitTest(x, 1)) {
        chg = true;
        y[1] = x[0];
        y[2] = MedianOf3(x[2], x[3], 3 * x[3] - 2 * x[4]);
    }

    for (size_t i = 2; i < n - 3; i++) {
        if (SplitTest(x, i)) {
            int j;
            // plateau at x[i] == x[i+1]

            // at left:
            if (-1 < (j = IndexOfMedianOf3(x[i], x[i - 1], 3 * x[i - 1] - 2 * x[i - 2]))) {
                y[i] = (j == 0 ? x[i - 1] : 3 * x[i - 1] - 2 * x[i - 2]);
                chg = (y[i] != x[i]);
            }

            // at right:
            if (-1 < (j = IndexOfMedianOf3(x[i + 1], x[i + 2], 3 * x[i + 2] - 2 * x[i + 3]))) {
                y[i + 1] = (j == 0 ? x[i + 2] : 3 * x[i + 2] - 2 * x[i + 3]);
                chg = (y[i + 1] != x[i + 1]);
            }
        }
    }

    if (do_ends && SplitTest(x, n - 3)) {
        chg = true;
        y[n - 2] = x[n - 1];
        y[n - 3] = MedianOf3(x[n - 3], x[n - 4], 3 * x[n - 4] - 2 * x[n - 5]);
    }

    return chg;
}

template<typename T>
size_t Smooth3RS3R(std::vector <T> &y, const std::vector <T> &x,
                   SmoothEndRule end_rule = SmoothEndRule::Tukey, bool split_ends = false) {
    // y[1:n] := "3R S 3R"(x[1:n]);  z = "work";
    size_t iter;
    bool chg;
    size_t n = x.size();

    y.resize(n);
    std::vector <T> z(n), w(n);

    iter = Smooth3R(&x[0], &y[0], &z[0], n, end_rule);
    chg = SmoothSplit3(&y[0], &z[0], n, split_ends);
    if (chg)
        iter += Smooth3R(&z[0], &y[0], &w[0], n, end_rule);

    /* else y == z already */
    return (iter + chg);
}

}