//
// Created by Elie Baier on 23.02.24.
//

#include <cmath>
#include "iostream"
#include "string"
#include "fstream"
#include "sstream"
#include "matplot/matplot.h"
#include <filesystem>

using namespace std;
namespace plt = matplot;
namespace fs = std::filesystem;

int MIN_ZERO_SEPARATION = 100;
int POWERED_MIN = 3;
int POWERED_MAX = 3;

long double average(vector<long double> a)
{
    return reduce(a.begin(), a.end(), 0.0) / a.size();
}

vector<long double> findZeroCrossings(const vector<long double>& data, long double absZero) {
    vector<long double> zeroCrossings;
    for (int i = 1; i < data.size(); ++i) {
        if ((data[i - 1] - absZero) * (data[i] - absZero) < 0 || (data[i - 1] - absZero) * (data[i] - absZero) == 0) {

            // Cleaning the duplicate zero which are +-MIN_ZERO_SEPARATION
            if(zeroCrossings.size() >= 1) {
                if(zeroCrossings[zeroCrossings.size() - 1] - MIN_ZERO_SEPARATION < i && zeroCrossings[zeroCrossings.size() - 1] + MIN_ZERO_SEPARATION > i) {
                    // cout << "Found duplicate zero: " << i << endl;
                    continue;
                }
            }

            zeroCrossings.push_back(i);
        }
    }
    return zeroCrossings;
}

pair<long double, bool> poweredMax(const vector<long double>& data, int n) {
    if (data.empty()) {
        return {0.0, false}; // If vector is empty, return false
    }

    unordered_map<long double, int> freqMap;
    for (long double value : data) {
        freqMap[value]++;
    }

    long double maxVal = data[0]; // Initialize maxVal with the first element
    int maxFreq = freqMap[maxVal];

    for (const auto& entry : freqMap) {
        if (entry.second >= n && entry.first > maxVal) {
            maxVal = entry.first;
            maxFreq = entry.second;
        }
    }

    bool found = maxFreq >= n;

    return {maxVal, found};
}

std::pair<long double, bool> poweredMin(const std::vector<long double>& data, int n) {
    if (data.empty()) {
        return {0.0, false}; // If vector is empty, return false
    }

    std::unordered_map<long double, int> freqMap;
    for (long double value : data) {
        freqMap[value]++;
    }

    long double minVal = std::numeric_limits<long double>::max(); // Initialize minVal with maximum double value
    int minFreq = freqMap[minVal];

    for (const auto& entry : freqMap) {
        if (entry.second >= n && entry.first < minVal) {
            minVal = entry.first;
            minFreq = entry.second;
        }
    }

    bool found = minFreq >= n;

    return {minVal, found};
}

std::vector<long double> slice(const std::vector<long double>& data, int start, int end) {
    if (start < 0 || end >= static_cast<int>(data.size()) || start > end) {
        std::cerr << "Invalid slice indices." << std::endl;
        return {};
    }

    return {data.begin() + start, data.begin() + end + 1};
}

ifstream openFile(string path) {
    ifstream file;
    file.open(path);
    if(!file.is_open()) {
        std::cerr << "Error opening file." << std::endl;
    }
    return file;
}

int analyseOne(int i, ofstream &out) {
    cout << setprecision(13);
    out << setprecision(13);
    vector<ifstream> channels;
    if(i < 10) {
        channels.push_back(openFile("/Users/eliebaier/Workspace/Physique/TP8/DATA/ALL000" + to_string(i) + "/F000" + to_string(i) + "CH1.CSV"));
        channels.push_back(openFile("/Users/eliebaier/Workspace/Physique/TP8/DATA/ALL000" + to_string(i) + "/F000" + to_string(i) + "CH2.CSV"));
    } else {
        channels.push_back(openFile("/Users/eliebaier/Workspace/Physique/TP8/DATA/ALL00" + to_string(i) + "/F00" + to_string(i) + "CH1.CSV"));
        channels.push_back(openFile("/Users/eliebaier/Workspace/Physique/TP8/DATA/ALL00" + to_string(i) + "/F00" + to_string(i) + "CH2.CSV"));
    }

    vector<vector<long double>> channelsXData, channelsYData, channelsZeros;
    vector<long double> channelsP2P, channelsAbsoluteZero, channelsTimings, channelsFrequency;
    long double channelOneT = 0;

    for(int channel = 0; channel < channels.size(); channel++) {
        string line;

        cout << "========== CHANNEL " << channel << " ==========" << endl;
        cout << "[INFO] Starting data fetching for channel " << channel << endl;

        // Skipping the TEKTRONICS DEFINITION lines
        cout << "[INFO] Skipping lines used by TEKTRONICS to store device data" << endl;
        for(auto i = 0; i < 18; i++) {
            getline(channels[channel], line);
        }

        cout << "[INFO] Reading channel data" << endl;
        vector<long double> xs, ys;
        while(getline(channels[channel], line)) {
            // Cleaning the lines
            line.erase(0, 3);
            line.erase(line.size() - 1, line.size());

            std::string cleanedLine;
            for (char c : line) {
                if (std::isdigit(c) || c == '.' || c == ',') {
                    cleanedLine += c;
                }
            }

            long double x, y;
            istringstream iss(line);
            char comma;
            if (!(iss >> x >> comma >> y)) {
                std::cerr << "Error parsing line: " << line << std::endl;
                continue;
            }
            xs.push_back(x);
            ys.push_back(y);
        }
        channelsXData.push_back(xs);
        channelsYData.push_back(ys);
        cout << "[INFO] Closing channel " << channel << " file" << endl;
        channels[channel].close();

        // Finding the zeros and finding the average minimum and the average maximum (also finding the reference zero)
        // Finding the absolute max (powered)
        cout << "[INFO] Finding the absolute max and min" << endl;
        double absoluteMax = poweredMax(channelsYData[channel], POWERED_MAX).first;
        cout << "[INFO] Absolute max is " << absoluteMax << endl;
        double absoluteMin = poweredMin(channelsYData[channel], POWERED_MAX).first;
        cout << "[INFO] Absolute min is " << absoluteMin << endl;

        // Find the reference zero and the offset to the measured one
        channelsAbsoluteZero.push_back((absoluteMax + absoluteMin)/2);
        cout << "[INFO] Reference zero is " << channelsAbsoluteZero[channel] << endl;

        cout << "[INFO] Finding zeros" << endl;
        channelsZeros.push_back(findZeroCrossings(channelsYData[channel], channelsAbsoluteZero[channel]));
        vector<long double> mins, maxs;
        for(int i = 1; i < channelsZeros[channel].size(); i++) {
            // Checking if it's a min or a max by taking the middle and looking at the sign
            // Add subtraction with channelsAbsoluteZero[channel]
            if(channelsYData[channel][int (channelsZeros[channel][i - 1] + channelsZeros[channel][i])/2] >= 0) {
                // This is a max
                // cout << "Found max between " << xData[zeroCrossings[i - 1]] << " and " << xData[zeroCrossings[i]] << " ";
                // cout << "max is: " << poweredMax(slice(yData, zeroCrossings[i - 1], zeroCrossings[i]), 5).first << endl;
                long double newlyMax = poweredMax(slice(channelsYData[channel], channelsZeros[channel][i - 1], channelsZeros[channel][i]), POWERED_MAX).first;
                cout << "Max found " << newlyMax << endl;
                maxs.push_back(newlyMax);
            } else {
                mins.push_back(poweredMin(slice(channelsYData[channel], channelsZeros[channel][i - 1], channelsZeros[channel][i]), POWERED_MIN).first);
            }
        }

        vector<long double> frequency;
        for(int j = 0; j + 2 < channelsZeros[channel].size(); j=j+2) {
            // cout << "[INFO] Processing zeros " << j << " and " << j+1 << endl;
            frequency.push_back(((channelsXData[channel][channelsZeros[channel][j + 2]]) - channelsXData[channel][channelsZeros[channel][j]]));
            cout << "[INFO] Distance between zero " << j << " (" << channelsZeros[channel][j] << ") and " << j+2 << " (" << channelsZeros[channel][j + 2] << ") is " << ((channelsXData[channel][channelsZeros[channel][j + 2]]) - channelsXData[channel][channelsZeros[channel][j]]) << " (" << channelsZeros[channel][j+2] - channelsZeros[channel][j] << ")" << endl;
        }

        long double avgFrequency = average(frequency);
        // avgFrequency *= 2;
        cout << "[INFO] Average peak time for channel " << channel << " is " << avgFrequency << "s" << endl;
        channelsTimings.push_back(avgFrequency);
        cout << "[INFO] Average frequency for channel " << channel << " is " << 1/avgFrequency << "Hz" << endl;
        channelsFrequency.push_back(1/avgFrequency);
        if(channel == 0) {
            channelOneT = avgFrequency;
            // cout << "Set chan 1 freq: " << channelOneT << endl;
        }

        long double avgMax = 0, avgMin = 0;
        for(auto max : maxs) {
            avgMax += max;
        }
        for(auto min : mins) {
            avgMin += min;
        }
        avgMax /= maxs.size();
        avgMin /= mins.size();
        cout << "[INFO] Average Peak for channel " << channel << " max is " << avgMax << endl;
        cout << "[INFO] Average Peak for channel " << channel << " min is " << avgMin << endl;
        cout << "[INFO] Average Peak-to-Peak for channel " << channel << " is " << abs(avgMin) + abs(avgMax) << endl;
        // Change this to subtraction for cases which are all negative
        channelsP2P.push_back(abs(avgMin) + abs(avgMax));
    }

    cout << "[INFO] Starting the calculate the phase shift" << endl;
    long double deltas = 0;
    int ndeltas = 0;
    for(int i = 0; i < channelsZeros[0].size(); i++) {
        if(channelsZeros[1].size() <= i) continue;
        long double d = channelsXData[1][channelsZeros[1][i]] - channelsXData[0][channelsZeros[0][i]];
        cout << "[INFO] Phase shift for zero " << i << " is " << ((d) / channelOneT) * (-360) << endl;
        deltas += (d / channelOneT) * (-360);
        ndeltas++;
    }
    deltas /= ndeltas;
    // cout << ndeltas << deltas << channelOneT << endl;
    cout << "[INFO] Average phase shift (CHA0 to CHA1) " << deltas << endl;

    cout << "[INFO] Writing to CSV" << endl;
    out << i << "," << deltas << "," << channelsP2P[0] << "," << channelsP2P[1] << "," << channelsTimings[0] << ","  << channelsTimings[1] << ","  << channelsFrequency[0] << ","  << channelsFrequency[1] << "\n";

    plt::tiledlayout(1, 1);

    auto ax1 = plt::nexttile();

    plt::plot(channelsXData[0], channelsYData[0]);
    plt::xlabel("Temps [s]");
    plt::ylabel("Tension [V]");
    plt::title("Tension en fonction du temps (Mesure 00" + to_string(i) + ")");

    plt::hold(matplot::on);

    plt::plot(channelsXData[1], channelsYData[1])->use_y2(true);

    plt::grid(matplot::on);
    plt::gca()->minor_grid(true);

    auto lgd =  ::matplot::legend(ax1, {"Générateur", "Bornes de la résistance"});
    lgd->location(matplot::legend::general_alignment::bottomleft);
    lgd->box(false);

    /* auto ax2 = plt::nexttile();
    plt::plot(ax2, channelsXData[1], channelsYData[1]);
    plt::xlabel(ax2, "Temps [s]");
    plt::ylabel(ax2, "Tension [V]");
    plt::title(ax2, "Tension aux bornes de la résistance en fonction du temps (Mesure 00" + to_string(i) + ")");

    plt::grid(matplot::on);
    plt::gca()->minor_grid(true); */

    /* vector<double> zeroCrossings = findZeroCrossings(yData);

    // Drawing lines on the found zeros
    for(auto zero : zeroCrossings) {
        cout << "Zero: " << zero << endl;
        plt::line(xData[zero], -2, xData[zero], 2);
    } */

    // Show the plot
    // plt::show();

    /* cout << "[INFO] Saving svg to " << "/Users/eliebaier/Workspace/Physique/TP8/build/" + to_string(i) + "/ch1.svg" << endl;
    fs::create_directories("./build/" + to_string(i));
    plt::save("/Users/eliebaier/Workspace/Physique/TP8/build/" + to_string(i) + "/ch1.svg"); */

    /* plt::tiledlayout(2, 1);

    auto ax2 = plt::nexttile();

    plt::plot(ax2, channelsXData[1], channelsYData[1]);
    plt::xlabel(ax2, "Temps [s]");
    plt::ylabel(ax2, "Tension [V]");
    plt::title(ax2, "Tension aux bornes de la résistance en fonction du temps (Mesure 00" + to_string(i) + ")"); */

    /* for(auto zero : channelsZeros[1]) {
        cout << "Zero: " << zero << endl;
        plt::line(channelsXData[1][zero], -2, channelsXData[1][zero], 2);
    } */

    // plt::show();

    cout << "[INFO] Saving svg to " << "/Users/eliebaier/Workspace/Physique/TP8/build/" + to_string(i) + "/both-with-label.svg" << endl;
    fs::create_directories("./build/" + to_string(i));
    plt::save("/Users/eliebaier/Workspace/Physique/TP8/build/" + to_string(i) + "/both-with-label.svg");

    return 0;

}

int mains() {
    std::ofstream out;
    out.open("/Users/eliebaier/Workspace/Physique/TP8/build/out.csv");
    out << "Mesure, Average Phase Shift, Ch1 Peak-to-Peak, Ch2 Peak-to-Peak, Ch1 Time, Ch2 Time, Ch1 Frequency, Ch2 Frequency" << endl;
    for(int i = 0; i < 53; i++) {
        analyseOne(i, out);
    }
    out.close();
}