#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <algorithm>
 #include <execution>   // only include if the compiler supports it
#include <omp.h>

using namespace std;

// OpenMP parallel computations (already implemented)
void computeMeansVectorParallel(const vector<int> &v, double &arith, double &geom, double &harm) {
    int n = v.size();
    double sum = 0.0, logSum = 0.0, recSum = 0.0;

    #pragma omp parallel for reduction(+:sum,logSum,recSum)
    for (int i = 0; i < n; ++i) {
        int x = v[i];
        sum += x;
        logSum += log(x);
        recSum += 1.0 / x;
    }

    arith = sum / n;
    geom  = exp(logSum / n);
    harm  = n / recSum;
}

double computeStdDevParallel(const vector<int> &v, double mean) {
    int n = v.size();
    double sumSq = 0.0;

    #pragma omp parallel for reduction(+:sumSq)
    for (int i = 0; i < n; ++i) {
        double diff = v[i] - mean;
        sumSq += diff * diff;
    }

    return sqrt(sumSq / n);
}

int main() {
    cout << "Exercise Sheet 5: Task 2\n";
    // Read data
    vector<int> v;
    ifstream fin("data_1.txt");
    int x;
    while (fin >> x) v.push_back(x);
    fin.close();

    cout << "Data size: " << v.size() << "\n";

    double t_start = omp_get_wtime();

    // OpenMP parallel computations
    #pragma omp parallel
    {
        #pragma omp single
        cout << "Using " << omp_get_num_threads() << " OpenMP threads\n";
    }

    int minVal = v[0], maxVal = v[0];
    #pragma omp parallel
    {
        int local_min = v[0], local_max = v[0];
        #pragma omp for nowait
        for (int i = 0; i < v.size(); ++i) {
            if (v[i] < local_min) local_min = v[i];
            if (v[i] > local_max) local_max = v[i];
        }
        #pragma omp critical
        {
            if (local_min < minVal) minVal = local_min;
            if (local_max > maxVal) maxVal = local_max;
        }
    }

    double arith, geom, harm;
    computeMeansVectorParallel(v, arith, geom, harm);

    double stddev = computeStdDevParallel(v, arith);

    double t_end = omp_get_wtime();

    cout << "OpenMP results:\n";
    cout << "Min: " << minVal << ", Max: " << maxVal << "\n";
    cout << "Arithmetic mean: " << arith
         << ", Geometric mean: " << geom
         << ", Harmonic mean: " << harm << "\n";
    cout << "Standard deviation: " << stddev << "\n";
    cout << "Elapsed time (s): " << t_end - t_start << "\n\n";

    cout << "\nNote: The C++17 execution::par version is implemented in the main file, "
        "but my Mac's Clang compiler does not support parallel execution policies.\n";

    // C++17 execution policy version
   
    double t_ep_start = omp_get_wtime();

    auto p = minmax_element(execution::par, v.begin(), v.end());
    double min_ep = *p.first;
    double max_ep = *p.second;

    double sum_ep = reduce(execution::par, v.begin(), v.end(), 0.0);
    double ar_ep = sum_ep / v.size();

    double geom_ep = transform_reduce(
        execution::par,
        v.begin(), v.end(),
        0.0,
        plus<>(),
        [n=v.size()](int val){ return log(val); }
    );
    geom_ep = exp(geom_ep / v.size());

    double sum_inv_ep = transform_reduce(
        execution::par,
        v.begin(), v.end(),
        0.0,
        plus<>(),
        [](int val){ return 1.0 / val; }
    );
    double harm_ep = v.size() / sum_inv_ep;

    double sum_sq_ep = transform_reduce(
        execution::par,
        v.begin(), v.end(),
        0.0,
        plus<>(),
        [ar_ep](int val){ double d = val - ar_ep; return d*d; }
    );
    double std_ep = sqrt(sum_sq_ep / v.size());

    double t_ep_end = omp_get_wtime();

    cout << "Execution policy results:\n";
    cout << "Min: " << min_ep << ", Max: " << max_ep << "\n";
    cout << "Arithmetic mean: " << ar_ep
         << ", Geometric mean: " << geom_ep
         << ", Harmonic mean: " << harm_ep << "\n";
    cout << "Standard deviation: " << std_ep << "\n";
    cout << "Elapsed time (s): " << t_ep_end - t_ep_start << "\n";
   
    return 0;
}