/** @file monte_carlo.cpp
 *
 *  @author Cory Alexander Balaton (coryab)
 *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
 *
 *  @version 1.0
 *
 *  @brief Implementation of the monte carlo functions
 *
 *  @bug No known bugs
 * */
#include "monte_carlo.hpp"

#include <cmath>
#include <cstdint>

void monte_carlo_progression(double T, int L, int cycles,
                             const std::string filename)
{
    // Set some variables
    data_t data, tmp;
    int n_spins = L * L;

    // File stuff
    std::string directory = utils::dirname(filename);
    std::ofstream ofile;

    // Create random engine using the mersenne twister
    std::random_device rd;
    uint32_t rd_sample = rd();
    std::mt19937 engine(rd_sample);

    std::cout << "Seed: " << rd_sample << std::endl;

    IsingModel ising(L, T);

    // Create path and open file
    utils::mkpath(directory);
    ofile.open(filename);

    // Loop through cycles
    for (size_t i = 1; i <= cycles; i++) {
        data += ising.Metropolis();
        tmp = data / (i * n_spins);
        ofile << i << ',' << tmp.E << ',' << tmp.E2 << ',' << tmp.M << ','
              << tmp.M2 << ',' << tmp.M_abs << '\n';
    }
    ofile.close();
}

void monte_carlo_progression(double T, int L, int cycles, int value,
                             const std::string filename)
{
    // Set some variables
    data_t data, tmp;
    int n_spins = L * L;

    // File stuff
    std::string directory = utils::dirname(filename);
    std::ofstream ofile;

    // Create random engine using the mersenne twister
    std::random_device rd;
    uint32_t rd_sample = rd();
    std::mt19937 engine(rd());

    std::cout << "Seed: " << rd_sample << std::endl;

    IsingModel ising(L, T, value);

    // Create path and open file
    utils::mkpath(directory);
    ofile.open(filename);

    // Loop through cycles
    for (size_t i = 1; i <= cycles; i++) {
        data += ising.Metropolis();
        tmp = data / (i * n_spins);
        ofile << i << ',' << tmp.E << ',' << tmp.E2 << ',' << tmp.M << ','
              << tmp.M2 << ',' << tmp.M_abs << '\n';
    }
    ofile.close();
}

void pd_estimate(double T, int L, int cycles, const std::string filename)
{
    // Set some variables
    data_t data, tmp;
    int n_spins = L * L;

    // File stuff
    std::string directory = utils::dirname(filename);
    std::ofstream ofile;

    IsingModel ising(L, T);

    // Create path and open file
    utils::mkpath(directory);
    ofile.open(filename);

    // Loop through cycles
    for (size_t i = 1; i <= cycles; i++) {
        data = ising.Metropolis() / n_spins;
        ofile << data.E << ',' << data.E2 << ',' << data.M << ',' << data.M2
              << ',' << data.M_abs << '\n';
    }
    ofile.close();
}

// Code for seeing phase transitions.
data_t monte_carlo_serial(int L, double T, int cycles)
{
    data_t data;
    IsingModel model(L, T);

    for (size_t i = 0; i < BURN_IN_TIME; i++) {
        model.Metropolis();
    }

    double E, M;
    for (size_t i = 0; i < (uint)cycles; i++) {
        data += model.Metropolis();
    }

    return data;
}

data_t monte_carlo_parallel(int L, double T, int cycles)
{
    data_t data;
#pragma omp parallel
    {
        // Each thread creates an instance of IsingModel.
        IsingModel model(L, T);

        // Each thread runs the Metropolis algorithm before starting to collect
        // samples
        for (size_t i = 0; i < BURN_IN_TIME; i++) {
            model.Metropolis();
        }

        // Now each thread work on one loop together, but using their own
        // instances of things, but the total of cycles add up.
        // static ensure that each thread gets the same amount of iterations
#pragma omp for schedule(static) reduction(+ : data)
        for (size_t i = 0; i < (uint)cycles; i++) {
            data += model.Metropolis();
        }
    }

    double norm = 1. / (double)cycles;

    return data * norm;
}

void phase_transition(int L, double start, double end, int points,
                      std::function<data_t(int, double, int)> monte_carlo,
                      std::string outfile)
{
    double dt = (end - start) / (double)points;
    int cycles = 10000;
    int N = L * L;
    std::ofstream ofile;

    data_t data;

    utils::mkpath(utils::dirname(outfile));
    ofile.open(outfile);

    double temp, CV, X, E_var, M_var;

    using utils::scientific_format;
    for (size_t i = 0; i < points; i++) {
        temp = start + dt * i;
        data = monte_carlo(L, temp, cycles);
        E_var = (data.E2 - data.E * data.E) / (double)N;
        M_var = (data.M2 - data.M_abs * data.M_abs) / (double)N;

        ofile << scientific_format(temp) << ','
              << scientific_format(data.E / (double)N) << ','
              << scientific_format(data.M_abs / N) << ','
              << scientific_format(E_var / (temp * temp)) << ','
              << scientific_format(M_var / temp) << '\n';
    }
    ofile.close();
}