#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;

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);
    

}

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;
    for (long unsigned int i=0; i<a.size(); i++){
        deviation += pow(x - a.at(i),2);
    }
    deviation = sqrt(deviation/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);

}

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);
        double 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 = %.f | time = %f milliseconds\n", sum, sec1*1000);
        
        // Timing second function
        tic();
        for (int i=0; i<loops; i++){
            sum = formula(n);
        }
        double sec2 = toc();
        printf("Formula funtion:  result = %.f | time = %f milliseconds\n", sum, sec2*1000);
        }
}

void task_d() {
    printf("\n\n-------------- Task D --------------\n\n");
    printf("See folder D.\n");
}

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);
    }
}

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

}

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 = 1000;

    // 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 (int k=1; k<NLOOPS; ++k){
        f3 = M2.Mult(w);
    }
    double t2 = toc();

    // Time MultT
    tic();
    vector<double> f4 = M2.MultT(w);
    for (int k=1; k<NLOOPS; ++k){
        f4 = M2.Mult(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 (int 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*i)/(n-1) - 5 );}
    DenseMatrix2 M3(x,x);
    double t4 = toc();

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

    // Time MultT
    tic();
    vector<double> f6 = M3.MultT(w);
    for (int k=1; k<NLOOPS; ++k){
        f6 = M3.Mult(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 (int i = 0; i < n; i++)
    {
        double err = f5[i] - f6[i];
        if(abs(err) > 1e-4)
        {
            cout << "Resulting vectors are not equal" << endl;
        }
    }

    printf("\nNOTE: difference in runtime noticable with n=10.000 (~30 seconds)\n");


}


int main(){

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

    return 0;
};