#include "utk/metrics/PCF.hpp"
#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;
  }
}

int main() {

  // load Points from file and generate metrics
  auto loadedPoints =
      utk::read_pointsets<long double>("../result_data/sampled.txt")[0];
  auto pointView =
      utk::Pointset<long double>::View(loadedPoints.Data(), NSAMPLES, 2);

  // Radial spectrum
  std::cout << "Calculating radial spectrum with " << NBINS
            << " bins and resolution " << DIMENSION << std::endl;
  auto radSpec =
      utk::RadialSpectrum(NBINS, 0.5, DIMENSION, true).compute(pointView);

  writeRadspecToFile("radSpec.txt", radSpec);

  for (float i = 0.015625f; i <= 0.5f; i += 0.015625f) {
    // PCF
    std::cout << "Calculating PCF with " << NBINS << " bins for max radius " << i << "…" << std::endl;
    auto pcf = utk::PCF{true, 0.01, i, NBINS, 0.001}.compute(pointView);

    writePCFToFile("pcffiles/pcf_" + std::to_string(i) + ".txt", pcf);
  }

  // FFT based spectrum
  std::cout << "Calculating spectrum with resolution " << DIMENSION << "…"
            << std::endl;

  auto spec = utk::Spectrum{DIMENSION, true}.compute(pointView);
  writeSpectrumToFile("spectrum.txt", spec);
}