monte_carlo.cpp

Go to the documentation of this file.

/** @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"
 
namespace montecarlo {
void progression(double T, int L, int cycles, const std::string filename,
                 int burn_in_time)
{
    // 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);
 
    for (size_t i = 0; i < burn_in_time; i++) {
        ising.Metropolis();
    }
 
    // 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 progression(double T, int L, int cycles, int value,
                 const std::string filename, int burn_in_time)
{
    // 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);
 
    for (size_t i = 0; i < burn_in_time; i++) {
        ising.Metropolis();
    }
 
    // 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,
                 int burn_in_time)
{
    // 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);
 
    for (size_t i = 0; i < burn_in_time; i++) {
        ising.Metropolis();
    }
 
    // 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 mcmc_serial(int L, double T, int cycles, int burn_in_time)
{
    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();
    }
 
    double norm = 1. / (double)cycles;
 
    return data * norm;
}
 
data_t mcmc_parallel(int L, double T, int cycles, int burn_in_time)
{
    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, int cycles,
                      std::function<data_t(int, double, int, int)> monte_carlo,
                      std::string outfile, int burn_in_time)
{
    double dt = (end - start) / (double)points;
    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, burn_in_time);
        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();
}
} // namespace montecarlo