// g++ *.cpp -o main
// g++ -g -ffast-math -O3 -march=native -Wall -pedantic -Wextra -Weffc++ -Woverloaded-virtual -Wfloat-equal -Wshadow -Wredundant-decls -fmax-errors=1 *.cpp -o main


#include "mylib.h"
#include "timing.h"

#include <cassert>          // assert
#include <vector>
#include <iostream>
#include <cmath>
#include <tuple>
#include <string>
#include <algorithm>
#include <fstream>
#include <list>
#include <stdexcept>
using namespace std;
using namespace std::chrono;  // timing

static void task_a() {
    printf("\n\n-------------- Task A --------------\n\n");

    auto [a,b,c] = means0(1,4,16);
    auto [d,e,f] = means0(2,3,5);
    auto [g,h,i] = means0(1000,4000,16000);
    printf("means(1,4,16)          = (%f, %f, %f)\n", a, b, c);
    printf("means(2,3,5)           = (%f, %f, %f)\n", d, e, f);
    printf("means(1000,4000,16000) = (%f, %f, %f)\n", g, h, i);

    vector<double> v = {4,8,15,16,23,42};
    auto [j,k,l] = means(v);
    printf("means(4,8,15,16,23,42) = (%f, %f, %f)\n", j, k, l);
    

}

static void task_b() {
    printf("\n\n-------------- Task B --------------\n\n");
    
    // Read vector
    vector<double> a;
    read_vector_from_file("data_1.txt", a);

    // Print numbers
    // for (unsigned int k=0; k<a.size(); ++k)
    // {
    //     cout << "  " << a.at(k);
    // }
    // cout << endl;

    // min and max
    auto min = min_element(a.begin(), a.end());
    auto max = max_element(a.begin(), a.end());
    printf("Minimum: %f\n", *min);
    printf("Maximum: %f\n", *max);

    // means
    auto [x,y,z] = means(a);
    printf("Arithmetic: %f\n", x);
    printf("Geometric: %f\n", y);
    printf("Harmonic: %f\n", z);

    // deviation
    double deviation(0.0);
    for (long unsigned int i=0; i<a.size(); i++){
        deviation += pow(x - a.at(i),2);
    }
    deviation = sqrt(deviation/static_cast<double>(a.size()));
    printf("Deviation: %f\n", deviation);

    // write results to file
    vector<double> b = {*min,*max,x,y,z,deviation};
    write_vector_to_file("out_1.txt", b);

}

static void task_c() {
    printf("\n\n-------------- Task C --------------\n\n");

    vector<int> n_values = {15, 1001, 1432987};
    
    for (int n : n_values) {
        printf("n = %d\n", n);
        long int sum = 0;
        double loops = 1000;

        // Timing first function
        tic();
        for (int i=0; i<loops; i++){
            sum = sum_of_spec(n);
        }
        double sec1 = toc();
        printf("For-loop funtion: result = %ld | time = %f milliseconds\n", sum, sec1*1000);
        
        // Timing second function
        tic();
        for (int i=0; i<loops; i++){
            sum = static_cast<long int>(formula(n));
        }
        double sec2 = toc();
        printf("Formula funtion:  result = %ld | time = %f milliseconds\n", sum, sec2*1000);
        }
}

#ifdef __GNUC__
#pragma GCC push_options
#pragma GCC optimize("O1")
#endif

static void task_d() {
    printf("\n\n-------------- Task D --------------\n\n");
    int const NLOOPS = 25;        // chose a value such that the benchmark runs at least 10 sec.
    unsigned int N = 50000001;
    //##########################################################################
    //  Memory allocation
        cout << "Memory allocation\n";

        vector<double> x(N), y(N);

        cout.precision(2);
        cout << 2.0 * N *sizeof(x[0]) / 1024 / 1024 / 1024 << " GByte Memory allocated\n";
        cout.precision(6);

    //##########################################################################
    //  Data initialization
    //  Special:  x_i = i+1;  y_i = 1/x_i  ==> <x,y> == N
        for (unsigned int i = 0; i < N; ++i)
        {
            x[i] = i + 1;
            y[i] = 1.0 / pow(x[i], 2);
        }

    //##########################################################################
        cout << "\nStart Benchmarking Normal sum\n";

    // Do calculation
        auto t1 = system_clock::now(); // start timer
        double sk1(0.0),ss(0.0);
        for (int i = 0; i < NLOOPS; ++i)
        {
            sk1 = normal_sum(y);
            ss += sk1;                   // prevents the optimizer from removing unused calculation results.
        }

        auto t2 = system_clock::now();  // stop timer
        auto duration = duration_cast<microseconds>(t2 - t1);        // duration in microseconds
        double t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds
        t_diff = t_diff/NLOOPS;     


    // Print result
        printf("\nSum = %.16f\n", sk1);
        
    //##########################################################################

    // Timings  and Performance
        cout << endl;
        cout.precision(2);
        cout << "Timing in sec. : " << t_diff << endl;
        cout << "GFLOPS         : " << 2.0 * N / t_diff / 1024 / 1024 / 1024 << endl;
        cout << "GiByte/s       : " << 2.0 * N / t_diff / 1024 / 1024 / 1024 * sizeof(x[0]) << endl;

    //##########################################################################

        cout << "\nStart Benchmarking Kahan summation\n";

    // Do calculation
        t1 = system_clock::now(); // start timer
        double sk2(0.0),sss(0.0);
        for (int i = 0; i < NLOOPS; ++i)
        {
            sk2 = Kahan_skalar(y);
            sss += sk2;                   // prevents the optimizer from removing unused calculation results.
        }

        t2 = system_clock::now();  // stop timer
        duration = duration_cast<microseconds>(t2 - t1);        // duration in microseconds
        t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds
        t_diff = t_diff/NLOOPS;                                      // duration per loop seconds
                                    // duration per loop seconds

    // Print result
        printf("\nSum = %.16f\n", sk2);
        

    //##########################################################################

    // Timings  and Performance
        cout << endl;
        cout.precision(2);
        cout << "Timing in sec. : " << t_diff << endl;
        cout << "GFLOPS         : " << 2.0 * N / t_diff / 1024 / 1024 / 1024 << endl;
        cout << "GiByte/s       : " << 2.0 * N / t_diff / 1024 / 1024 / 1024 * sizeof(x[0]) << endl;

    //##########################################################################

    // Print limit
    printf("\nLimit = %.16f\n\n", pow(M_PI,2) / 6.0f);
}

#ifdef __GNUC__
#pragma GCC pop_options
#endif

static void task_e() {
    printf("\n\n-------------- Task E --------------\n\n");

    for (int n : {100, 1000, 10000}) {
        vector<int> vec(n);
        list<int> lst(n);

        // Initialize
        for (int i = 1; i < n+1; ++i) {
            vec.push_back(i);
            lst.push_back(i);
        }

        // Insert into vector
        tic();
        insert_into_vector(vec, n);
        double sec1 = toc();
        printf("Vector insertion time for n = %d: %.f microseconds.\n", n, sec1*1000*1000);

        // Insert into list
        tic();
        insert_into_list(lst, n);
        double sec2 = toc();
        printf("List insertion time for   n = %d: %.f microseconds.\n", n, sec2*1000*1000);
    }
}

static void task_f() {
    printf("\n\n-------------- Task F --------------\n\n");

    // single_goldbach(k)
    int k = 694;
    printf("single_goldbach(k = %d) = %d\n", k, single_goldbach(k));

    // Prints decompositions
    print_decomps(k);
    
    // count_goldbach(n)
    printf("\nNOTE: For n=2'000'000 it will take ~30 seconds.\n");
    for (int n : {10'000, 100'000, 400'000, 1'000'000}) {       //, 2'000'000}) {
        tic();
        vector<int> counts = count_goldbach(n);
        double sec = toc();

        auto max = max_element(counts.begin(), counts.end());
        printf("count_goldbach(n = %d): k = %ld, decompositions = %d, time elapsed: %f milliseconds\n", n, max-counts.begin(), *max, sec*1000);
    }

    // Results
    // count_goldbach(n =     10'000): k =    9240, decompositions =    329, time elapsed:      1.235096 milliseconds
    // count_goldbach(n =    100'000): k =   99330, decompositions =   2168, time elapsed:     39.003922 milliseconds
    // count_goldbach(n =    400'000): k =  390390, decompositions =   7094, time elapsed:    497.282572 milliseconds
    // count_goldbach(n =  1'000'000): k =  990990, decompositions =  15594, time elapsed:   3236.044944 milliseconds
    // count_goldbach(n =  2'000'000): k = 1981980, decompositions =  27988, time elapsed:  29864.384370 milliseconds
    // count_goldbach(n = 10'000'000): k = 9699690, decompositions = 124180, time elapsed: 825392.110981 milliseconds

}

static void task_g() {
    printf("\n\n-------------- Task G --------------\n\n");

    DenseMatrix const M(5,3);

    vector<double> const u{{1,2,3}};
    vector<double> f1 = M.Mult(u);

    vector<double> const v{{-1,2,-3,4,-5}};
    vector<double> f2 = M.MultT(v);


    cout << "M = " << endl;
    M.print();
    cout << endl << "u       = ";
    for (size_t i=0; i<u.size(); i++){ cout << u[i] << " ";}
    cout << endl << "M * u   = ";
    for (size_t i=0; i<f1.size(); i++){cout << f1[i] << " ";}
    cout << endl << "v       = ";
    for (size_t i=0; i<v.size(); i++){cout << v[i] << " ";}
    cout << endl << "M^T * v = ";
    for (size_t i=0; i<f2.size(); i++){cout << f2[i] << " ";}
    cout << endl << endl;

    // #######################################################

    int const NLOOPS = 100;
    int const n = 3000;

    // Time initialization
    tic();
    DenseMatrix const M2(n,n);
    vector<double> w(n,0);
    for (int i=0; i<n; i++){w[i]=i+1;}
    double t1 = toc();

    // Time Mult
    tic();
    vector<double> f3 = M2.Mult(w);
    for (size_t k=1; k<NLOOPS; ++k){
        f3 = M2.Mult(w);
    }
    double t2 = toc();

    // Time MultT
    tic();
    vector<double> f4 = M2.MultT(w);
    for (size_t k=1; k<NLOOPS; ++k){
        f4 = M2.MultT(w);
    }
    double t3 = toc();

    // Print results
    printf("Results for %dx%d matrix vector multiplication doing %d loops\n", n, n, NLOOPS);
    printf("Time for initialization: %f seconds.\n", t1);
    printf("Time for Mult          : %f seconds, %f per loop.\n", t2, t2/NLOOPS);
    printf("Time for MultT         : %f seconds, %f per loop.\n", t3, t3/NLOOPS);

    // Check if resulting vectors are equal
    for (size_t i = 0; i < n; i++)
    {
        double err = f3[i] - f4[i];
        if(abs(err) > 1e-4)
        {
            cout << "Resulting vectors are not equal" << endl;
        }
    }
// ################################################

    // Time initialization
    tic();
    vector<double> x(n,0);
    for (int i=0; i<n; i++){x[i] = sigmoid( (10.0*static_cast<double>(i))/(n-1) - 5 );}
    DenseMatrix2 M3(x,x);
    double t4 = toc();

    // Time Mult
    tic();
    vector<double> f5 = M3.Mult(w);
    for (size_t k=1; k<NLOOPS; ++k){
        f5 = M3.Mult(w);
    }
    double t5 = toc();

    // Time MultT
    tic();
    vector<double> f6 = M3.MultT(w);
    for (size_t k=1; k<NLOOPS; ++k){
        f6 = M3.MultT(w);
    }
    double t6 = toc();
    
     // Print results
    printf("\nResults for %dx%d matrix vector multiplication doing %d loops taking advantage of tensor product structure of the matrix\n", n, n, NLOOPS);
    printf("Time for initialization: %f seconds.\n", t4);
    printf("Time for Mult          : %f seconds, %f per loop.\n", t5, t5/NLOOPS);
    printf("Time for MultT         : %f seconds, %f per loop.\n", t6, t6/NLOOPS);

    // Check if resulting vectors are equal
    for (size_t i = 0; i < n; i++)
    {
        double err = f5[i] - f6[i];
        if(abs(err) > 1e-4)
        {
            cout << "Resulting vectors are not equal" << endl;
        }
    }

}


int main(){

    task_a();
    task_b();
    task_c();
    task_d();
    task_e();
    task_f();
    task_g();

    return 0;
}