#include <DataFrame/DataFrameMLVisitors.h>

template<arithmetic T, typename I = unsigned long>
struct  HWESForecastVisitor;

// -------------------------------------

template<typename T, typename I = unsigned long>
using hwes_v = HWESForecastVisitor<T, I>;

explicit
HWESForecastVisitor(long periods = 3,
                    size_type season_length = 0,
                    value_type alpha = 0.25,
                    value_type beta = 0.1,
                    value_type gamma = 0.1,
                    decompose_type season_type = decompose_type::additive);

periods: Number of periods ahead to forecast.
season_length: season_length tells the algorithm how many data points make up one full seasonal cycle — that is,
               how long it takes for a pattern to repeat.
alpha: Level smoothing, Alpha controls how much weight the algorithm gives to the most recent observation when
       updating the level (the estimated base value around which data fluctuates).

       High (close to 1), the model reacts quickly to recent changes.
       Low (close to 0), the level is smoothed heavily, slow to react (more stable).
beta: Trend smoothing, Beta controls how quickly the trend (slope) adapts over time. In the additive model

      High (close to 1), trend changes rapidly (useful if the slope changes often).
      Low (close to 0), trend is stable (changes only slowly).
gamma: Seasonal smoothing, Gamma controls how fast the seasonal pattern e.g., monthly or weekly fluctuations)
       is updated.

       High (close to 1), seasonal effects can shift quickly (good if the pattern changes over time).
       Low (close to 0), seasonality assumed stable.
season_type: season_type tells the algorithm how seasonal effects combine with the base signal (level + trend).

             Additive: seasonality is added on top of the baseline.
             Multiplicative: seasonality scales the baseline.
        

static void test_HWESForecastVisitor()  {

    std::cout << "\nTesting HWESForecastVisitor{ } ..." << std::endl;

    std::vector<unsigned long>  idxvec = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 };
    std::vector<double>         col1 = { 266.0, 145.9, 183.1, 119.3, 180.3, 168.5, 231.8, 224.5, 192.8, 122.9, 336.5, 185.9, 194.3 };
    std::vector<double>         oscil = { 1.5, 1.8, 1.62, 1.78, 1.5, 1.68, 1.6, 1.8, 1.71, 1.9, 1.78, 1.84, 1.69 };
    std::vector<double>         constant = { 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56 };
    std::vector<double>         increasing = { 10.56, 10.68, 10.78, 10.90, 11.01, 11.45, 11.99, 12.01, 12.21, 12.35, 12.67, 13.89, 13.01 };
    std::vector<double>         decreasing = { 10.56, 10.30, 10.12, 10.01, 9.80, 9.74, 9.41, 9.03, 9.0, 8.20, 8.01, 7.9, 7.55 };
    ULDataFrame                 df;

    df.load_data(std::move(idxvec),
                 std::make_pair("col1", col1),
                 std::make_pair("oscil", oscil),
                 std::make_pair("constant", constant),
                 std::make_pair("increasing", increasing),
                 std::make_pair("decreasing", decreasing));

    HWESForecastVisitor<double>   hwes;

    df.single_act_visit<double>("col1", hwes);

    const auto  result1 = hwes.get_result();

    assert(result1.size() == 3);
    assert(std::fabs(result1[0] - 208.351) < 0.001);
    assert(std::fabs(result1[1] - 208.69) < 0.001);
    assert(std::fabs(result1[2] - 209.03) < 0.001);

    HWESForecastVisitor<double>   hwes2 { 3, 2 };

    df.single_act_visit<double>("oscil", hwes2);

    const auto  result2 = hwes2.get_result();

    assert(result2.size() == 3);
    assert(std::fabs(result2[0] - 1.73499) < 0.00001);
    assert(std::fabs(result2[1] - 1.9216) < 0.00001);
    assert(std::fabs(result2[2] - 1.76383) < 0.00001);

    df.single_act_visit<double>("constant", hwes);

    const auto  result3 = hwes.get_result();

    assert(result3.size() == 3);
    assert(std::fabs(result3[0] - 10.56) < 0.00001);
    assert(std::fabs(result3[1] - 10.56) < 0.00001);
    assert(std::fabs(result3[2] - 10.56) < 0.00001);

    df.single_act_visit<double>("increasing", hwes);

    const auto  result4 = hwes.get_result();

    assert(result4.size() == 3);
    assert(std::fabs(result4[0] - 13.4675) < 0.0001);
    assert(std::fabs(result4[1] - 13.6977) < 0.0001);
    assert(std::fabs(result4[2] - 13.928) < 0.0001);

    df.single_act_visit<double>("decreasing", hwes);

    const auto  result5 = hwes.get_result();

    assert(result5.size() == 3);
    assert(std::fabs(result5[0] - 7.41545) < 0.00001);
    assert(std::fabs(result5[1] - 7.16314) < 0.00001);
    assert(std::fabs(result5[2] - 6.91084) < 0.00001);

    // Now some real data
    //
    StrDataFrame    df2;

    try  {
        df2.read("IBM.csv", io_format::csv2);
    }
    catch (const DataFrameError &ex)  {
        std::cout << ex.what() << std::endl;
        ::exit(-1);
    }

    HWESForecastVisitor<double, std::string>    hwes3 { 4 };

    df2.single_act_visit<double>("IBM_Close", hwes3);

    const auto  result6 = hwes3.get_result();

    assert(result6.size() == 4);
    assert(std::fabs(result6[0] - 109.264) < 0.001);
    assert(std::fabs(result6[1] - 108.313) < 0.001);
    assert(std::fabs(result6[2] - 107.361) < 0.001);
    assert(std::fabs(result6[3] - 106.41) < 0.001);

    // Now multidimensional data
    //
    constexpr std::size_t   dim { 3 };
    constexpr std::size_t   n { 24 };

    using ary_col_t = std::array<double, dim>;
    using vec_col_t = std::vector<double>;

    std::vector<vec_col_t>  vec_col1(n, vec_col_t(dim));

    for (std::size_t t { 0 }; t < n; ++t)  {
        vec_col1[t][0] = 10.0 + 2.0 * t;   // dim 0: upward
        vec_col1[t][1] = 100.0 - 2.0 * t;  // dim 1: downward
        vec_col1[t][2] = 50.0;             // dim 2: flat
    }

    // Add a repeating additive seasonal pattern of period 12 to all dims
    // so the seasonal path is exercised meaningfully.
    //
    //
    const std::vector<double>   pattern { 3, -3, 2, -2, 1, -1, 3, -3, 2, -2, 1, -1, };

    for (size_t t { 0 }; t < n; ++t)
        for (size_t d { 0 }; d < dim; ++d)
            vec_col1[t][d] += pattern[t % pattern.size()];

    std::vector<ary_col_t>  ary_col1(n);

    for (size_t t { 0 }; t < n; ++t)
        for (size_t d { 0 }; d < dim; ++d)
            ary_col1[t][d] = vec_col1[t][d];

    df2.load_column<vec_col_t>("VEC OBSV", std::move(vec_col1), nan_policy::dont_pad_with_nans);
    df2.load_column<ary_col_t>("ARY OBSV", std::move(ary_col1), nan_policy::dont_pad_with_nans);

    // No seasons
    //
    HWESForecastVisitor<vec_col_t, std::string> vec_hwes1 { 5, 0, 0.3, 0.1 };
    HWESForecastVisitor<ary_col_t, std::string> ary_hwes1 { 5, 0, 0.3, 0.1 };

    df2.single_act_visit<vec_col_t>("VEC OBSV", vec_hwes1);
    df2.single_act_visit<ary_col_t>("ARY OBSV", ary_hwes1);

    assert(vec_hwes1.get_result().size() == 5);
    for (const auto &vec : vec_hwes1.get_result())
        assert(vec.size() == dim);
    assert(ary_hwes1.get_result().size() == 5);
    for (const auto &ary : ary_hwes1.get_result())
        assert(ary.size() == dim);

    // dim 0 is trending up — each forecast step should be increasing
    //
    for (std::size_t i { 1 }; i < vec_hwes1.get_result().size(); ++i)
        assert(vec_hwes1.get_result()[i][0] > vec_hwes1.get_result()[i - 1][0]);
    // dim 1 is trending down — each step should be decreasing
    //
    for (std::size_t i { 1 }; i < ary_hwes1.get_result().size(); ++i)
        assert(ary_hwes1.get_result()[i][1] < ary_hwes1.get_result()[i - 1][1]);

    // With seasons 12 periods
    //
    HWESForecastVisitor<vec_col_t, std::string> vec_hwes2 { 6, 12, 0.3, 0.1, 0.2, decompose_type::additive };

    df2.single_act_visit<vec_col_t>("VEC OBSV", vec_hwes2);
    assert(vec_hwes2.get_result().size() == 6);
    assert(vec_hwes2.get_fitted().size() == n);
    assert(vec_hwes2.get_seasonal_factors().size() == pattern.size());
    for (const auto &vec : vec_hwes2.get_result())
        assert(vec.size() == dim);
    assert(std::abs(vec_hwes2.get_result()[0][0] - 45.1255) < 0.0001);
    assert(std::abs(vec_hwes2.get_result()[2][2] - 52.0) < 0.0001);
    assert(std::abs(vec_hwes2.get_result()[4][1] - 59.49) < 0.0001);
    assert(std::abs(vec_hwes2.get_result()[5][0] - 53.3878) < 0.0001);
}