#include "utk/metrics/PCF.hpp"
#include "utk/metrics/RDF_Heck.hpp"
#include <filesystem>
#include <format>
#include <utk/metrics/RadialSpectrum.hpp>
#include <utk/samplers/SamplerStep.hpp>
#include <utk/utils/Pointset.hpp>
#include <utk/utils/PointsetIO.hpp>

#define DIMENSION 1025
#define NSAMPLES 4096
#define NBINS 1000
#define SEED 8970

template<typename T>
void writeSpectrumToFile(const std::string &filename, std::vector<T> spectrum) {
    std::ofstream file;
    file.open(filename);

    file << DIMENSION << std::endl;

    for (auto freq: spectrum) {
        file << std::setprecision(std::numeric_limits<long double>::digits10 + 2)
             << std::fixed;
        file << freq << std::endl;
    }
}

template<typename T>
void writePCFToFile(const std::string &filename, std::vector<T> spectrum) {
    std::ofstream file;
    file.open(filename);

    for (auto freq: spectrum) {
        file << std::setprecision(std::numeric_limits<long double>::digits10 + 2)
             << std::fixed;
        file << freq << std::endl;
    }
}

template<typename T>
void writeRadspecToFile(const std::string &filename,
                        std::pair<std::vector<T>, std::vector<T>> radspec) {
    std::ofstream file;
    file.open(filename);

    auto xs = radspec.first;
    auto ys = radspec.second;

    if (xs.size() != ys.size()) {
        std::cerr << "Dimensions of radial spactrum are unequal: xDim: "
                  << xs.size() << " yDim: " << ys.size() << std::endl;
        std::terminate();
    }

    for (int i = 0; i < xs.size(); ++i) {
        file << std::setprecision(std::numeric_limits<long double>::digits10 + 2)
             << std::fixed;
        file << xs[i] << ", " << ys[i] << std::endl;
    }
}

template<typename T>
void writeRDFtoFile(const std::string &filename,
                    std::pair<std::vector<T>, std::vector<T>> radspec) {
    std::ofstream file;
    file.open(filename);

    auto xs = radspec.first;
    auto ys = radspec.second;

    if (xs.size() != ys.size()) {
        std::cerr << "Dimensions of radial spectrum are unequal: xDim: "
                  << xs.size() << " yDim: " << ys.size() << std::endl;
        std::terminate();
    }

    for (int i = 0; i < xs.size(); ++i) {
        file << std::setprecision(std::numeric_limits<long double>::digits10 + 2)
             << std::fixed;
        file << xs[i] << ", " << ys[i] << std::endl;
    }
}

int main() {

    auto resultDataPath =
            std::filesystem::path{"../result_data/pointset_series/pointsets/"};
    auto outputPath =
            std::filesystem::path{"../result_data/pointset_series/radspecs/"};

    std::vector<std::pair<std::string, utk::Pointset<long double>>> pointsets;

    for (const auto &file: std::filesystem::directory_iterator(resultDataPath)) {

        auto loadedPoints =
                utk::read_pointsets<long double>(file.path().c_str())[0];

        auto copiedPointset = utk::Pointset<long double>(std::move(loadedPoints));

        pointsets.emplace_back(file.path().stem().string(), copiedPointset);
    }

    auto radSpec =
            utk::RadialSpectrum(
                    0,  // Number of bins
                    0.5,// Scale for distances (compress / expands bins size)
                    0,  // Resolution for underlying fourier spectrum, only odd number.
                    // If even will use res - 1
                    true// Cancel DC
            );


    for (auto &pointset: pointsets) {
        std::cout << "computing radial spectrum for " << pointset.first << std::endl;

        auto result = radSpec.compute(pointset.second);

        writeRadspecToFile((outputPath / pointset.first).string() + ".rdf", result);
    }

    return 0;
}