#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <limits>
#include <chrono>
#include <thread>
#include <future>
#include <omp.h>

using namespace std;

double round_to(double value, double precision = 1.0)
{
    return round(value / precision) * precision;
}

// Forward fill function to fill missing data points with the last valid value
void forward_fill(vector<double>& data) {
    for (size_t i = 1; i < data.size(); ++i) {
        // If the current value is zero (or NaN), replace it with the last valid value
        if (data[i] == 0 || std::isnan(data[i])) {
            data[i] = data[i - 1];
        }
    }
}

// Replace infinity values with a specified fallback value (e.g., 0 or last valid value)
void replace_infinity(vector<double>& data, double fallback_value = 90) {
    for (size_t i = 1; i < data.size(); ++i) {
        if (std::isinf(data[i])) {
            // If the value is +∞ or -∞, replace it with fallback_value
         if(data[i]<std::numeric_limits<int>::lowest())
            data[i] = -fallback_value;
         else
            data[i] = fallback_value;
        }
    }
}
// Function to calculate Exponential Moving Average (EMA)
vector<double> ema(const vector<double>& data, size_t n) {
    vector<double> result;
    double multiplier = 2.0 / (n + 1);
    
    // Initialize the first EMA value as the first data point
    result.push_back(data[0]);
    for (size_t i = 1; i < data.size(); ++i) {
        result.push_back((data[i] - result[i - 1]) * multiplier + result[i - 1]);
    }

    return result;
}

// Function to calculate the rolling standard deviation (std)
vector<double> rolling_std(const vector<double>& data, size_t n) {
    vector<double> result;
    
    for (size_t i = 0; i < data.size(); ++i) {
        if (i >= n - 1) {
            double sum = 0, mean = 0, sq_sum = 0;
            for (size_t j = i; j > i - n; --j) {

                sum += data[j];
                sq_sum += data[j] * data[j];
                
            }
            mean = sum / n;
            result.push_back(sqrt(sq_sum / n - mean * mean));
            
        }
    }

    return result;
}

// Function to calculate the Bollinger Bands and return the result as a vector
vector<double> bollinger(const vector<double>& data, const vector<double>& ema_data, const vector<double>& std_data, size_t n) {
    
    vector<double> result(data.size(), 0);
    vector<double> result_diff(data.size(), 0);
    for (size_t i = 0; i < data.size(); ++i) {
        if (i >= n - 1) {
            result[i] = (data[i] - ema_data[i]) / std_data[i]; // Standardized difference
        }
    }
    for (size_t i = 1; i < data.size(); ++i) {
        
            result_diff[i] = round_to(result[i]-result[i-1], 0.000001); // Standardized difference
    }
    return result_diff;
}


void bollinger_exec(const vector<double>& data, const vector<double>& ema_data, const vector<double>& std_data, size_t n, vector<double>& result_diff) {
    vector<double> result(data.size(), 0);
    
    #pragma omp parallel num_threads(10)
    {
    #pragma omp parallel for  num_threads(10)
    for (size_t i = 0; i < data.size(); ++i) {
        if (i >= n - 1) {
            result[i] = (data[i] - ema_data[i]) / std_data[i]; // Standardized difference
        }
    }
    }
    
    result_diff.emplace_back(0);
    for (size_t i = 1; i < data.size(); ++i) 
    {        
        //result_diff[i] = round_to(result[i]-result[i-1], 0.000001); // Standardized difference
        result_diff.emplace_back(round_to(result[i]-result[i-1], 0.000001));
    }
    
}

// Function to transform data (calculate EMA, std, Bollinger) and apply forward fill and infinity replacement
vector<vector<double>> transform(const vector<double>& data, size_t n) {
    // Calculate EMA
    vector<double> bollinger_data;
    vector<double> ema_data;
    vector<double> std_data;
    ema_data = ema(data, n);

    // Calculate rolling standard deviation
    std_data = rolling_std(data, n);
    
    // Calculate Bollinger Bands
    //future<vector<double>> ema_async = async(launch::deferred, ema, ref(data), n);
    //future<vector<double>> std_async = async(launch::deferred, rolling_std, ref(data), n);

    //vector<double> ema_data = ema_async.get();
    //vector<double> std_data = std_async.get();
    //vector<double> bollinger_data = bollinger(data, ema_data, std_data, n);
    bollinger_exec(data, ema_data, std_data, n, bollinger_data);

    //future<void> boll_async = async(launch::deferred,bollinger_exec, ref(data), ref(ema_data), ref(std_data), n,ref(bollinger_data));
    //boll_async.wait();
    //thread t = thread(bollinger_exec, ref(data), ref(ema_data), ref(std_data), n, ref(bollinger_data));
    //t.join();
    // Forward fill the Bollinger Bands, EMA, and std data to handle any missing values
    //forward_fill(ema_data);
    //forward_fill(std_data);
    forward_fill(bollinger_data);
    //thread t_forward_fill = thread(forward_fill, ref(bollinger_data));
    
    //thread t_replace_infinity = thread(replace_infinity, ref(bollinger_data), n);

    // Replace infinity values with fallback_value (0 for simplicity)
    //replace_infinity(ema_data);
    //replace_infinity(std_data);
    replace_infinity(bollinger_data);
    
    //t_forward_fill.join();
    //t_replace_infinity.join();

    //future<void> t_forward_fill = async(launch::deferred,forward_fill, ref(bollinger_data));
    //t_forward_fill.wait();
    //future<void> t_replace_infinity = async(launch::deferred, replace_infinity, ref(bollinger_data), n);    
    //t_replace_infinity.wait();
    // Combine the results into a 2D vector (for simplicity)
    //vector<vector<double>> result;
    vector<vector<double>> result(data.size(), vector<double>(4, 0));
    
    for (size_t i = 0; i < data.size(); ++i) {
            result[i][0] = data[i];        // Original Data
            result[i][1] = ema_data[i];    // EMA
            result[i][2] = std_data[i];    // Rolling STD
            result[i][3] = bollinger_data[i]; // Bollinger Bands
        }
    return result;
    
}

int main() {
    // Sample data (example: closing prices)
    //vector<double> data = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};

    
       
    // Introduce some infinity values for testing
    //data[3] = std::numeric_limits<double>::infinity();  // +∞ at index 3
    //data[5] = -std::numeric_limits<double>::infinity(); // -∞ at index 5
    
    for (size_t i = 0; i < 100; i++)
    {
    size_t n = 3; // Period for EMA and Bollinger
    //replace_infinity(data);
    // Call the transform function
    auto t1 = chrono::high_resolution_clock::now();
    vector<double> data;
    
    for (size_t i = 1; i <= 50; i++)
    {
        double ii(i);
        data.emplace_back(ii);
    }
    
    vector<vector<double>> result = transform(data, n);
    auto t2 = chrono::high_resolution_clock::now();

    chrono::duration<double, std::milli> ms_double = t2 - t1;
    std::cout << ms_double.count() << "ms\n"<<endl;
    std::cout << result.size() <<"n"<<endl;
    
    // Output the results
    
    /*for (size_t i = n+1; i < result.size(); ++i) {
        cout << "Data: " << result[i][0] 
             << " EMA: " << result[i][1] 
             << " Std: " << result[i][2] 
             << " Bollinger: " << result[i][3] << endl;
    }*/
    }
    
    system("pause");
    return 0;
}