Project-4/src/test_suite.cpp

/** @file test_suite.cpp
 *
 *  @author Cory Alexander Balaton (coryab)
 *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
 *
 *  @version 1.0
 *
 *  @brief Sweep over different temperatures and generate data.
 *
 *  @bug No known bugs
 * */
#include "IsingModel.hpp"
#include "testlib.hpp"

#include <fstream>

#define EPS_2 (-2 * std::sinh(8.)) / (std::cosh(8.) + 3)

#define MAG_2 (std::exp(8.) + 1) / (2 * (cosh(8.) + 3))

#define CV_2                                                                   \
    16 * (3 * std::cosh(8.) + 1) / ((std::cosh(8.) + 3) * (std::cosh(8.) + 3))

#define X_2                                                                    \
    (3 * std::exp(8.) + std::exp(-8.) + 3)                                     \
        / ((std::cosh(8.) + 3) * (std::cosh(8.) + 3))

/** @brief Test class for the Ising model
 * */
class IsingModelTest {
public:
    /** @brief Test That initializing works as intended.
     * */
    void test_init_functions()
    {
        IsingModel test;
        test.L = 3;
        test.T = 1.;

        // Test that initializing the lattice only yields 1s and -1s.
        test.initialize_lattice();
        std::function<bool(int)> f = [](int x) { return x == 1 || x == -1; };
        ASSERT(testlib::assert_each(f, test.lattice),
               "Test lattice initialization.");

        test.initialize_neighbors();
        arma::Mat<int> neighbor_matrix("2, 1 ; 0, 2 ; 1, 0");
        ASSERT(testlib::is_equal(neighbor_matrix, test.neighbors),
               "Test neighbor matrix.");

        // Fill the lattice with 1s to be able to test the next functions.
        test.lattice.fill(1);

        // Test the initial magnetization.
        test.initialize_magnetization();
        ASSERT(std::fabs(test.M - 9.) < 1e-8, "Test intial magnetization");

        // Test that the initial energy is correct
        test.initialize_energy();
        ASSERT(std::fabs(test.E - (-18)) < 1e-8, "Test initial energy.");
    }

    /** @brief Test numerical data with analytical data.
     *
     *  @param tol The tolerance between the analytical and numerical solution.
     *  @param max_cycles The max number of Monte Carlo cycles.
     *
     *  return int
     * */
    int test_2x2_lattice(double tol, int max_cycles)
    {
        data_t data, tmp;
        size_t L = 2;
        size_t n_spins = L * L;
        double T = 1.;
        size_t cycles = 0;

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

        IsingModel test(L, T);

        int arr[]{0, 0, 0, 0};

        // Loop through cycles
        //std::ofstream ofile;
        //ofile.open("output/test_2x2.txt");
        while (cycles++ < max_cycles) {
            data += test.Metropolis();
            tmp = data / cycles;
            //ofile << cycles << ',' << tmp.E / n_spins << ','
                  //<< tmp.M_abs / n_spins << ','
                  //<< (tmp.E2 - tmp.E * tmp.E) / (T * T) / n_spins << ','
                  //<< (tmp.M2 - tmp.M_abs * tmp.M_abs) / T / n_spins << '\n';
            if (testlib::close_to(EPS_2, tmp.E / n_spins, tol)
                && testlib::close_to(MAG_2, tmp.M_abs / n_spins, tol)
                && testlib::close_to(CV_2, (tmp.E2 - tmp.E * tmp.E) / (T * T)
                                               / n_spins, tol)
                && testlib::close_to(X_2, (tmp.M2 - tmp.M_abs * tmp.M_abs) / T
                                              / n_spins, tol)) {
                return cycles;
            }
        }
        //std::cout << EPS_2 << ',' << MAG_2 << ',' << CV_2 << ',' << X_2
                  //<< std::endl;
        //ofile.close();
        // cycles = 0;
        // data = 0;
        // IsingModel test_mag(L, T);
        // while (cycles++ < max_cycles) {
        // data += test.Metropolis();
        // tmp = data / (cycles * n_spins);
        // if (testlib::close_to(MAG_2, tmp.M, tol)) {
        // arr[1] = cycles;
        // break;
        //}
        //}
        // cycles = 0;
        // data = 0;
        // IsingModel test_CV(L, T);
        // while (cycles++ < max_cycles) {
        // data += test.Metropolis();
        // tmp = data / (cycles * n_spins);
        // if (testlib::close_to(CV_2, (tmp.E2 - tmp.E * tmp.E) / (T * T),
        // tol)) {
        // arr[2] = cycles;
        // break;
        //}
        //}
        // cycles = 0;
        // data = 0;
        // IsingModel test_X(L, T);
        // while (cycles++ < max_cycles) {
        // data += test.Metropolis();
        // tmp = data / (cycles * n_spins);
        // if (testlib::close_to(X_2, (tmp.M2 - tmp.M_abs * tmp.M_abs) / T,
        // tol)) {
        // arr[3] = cycles;
        // break;
        //}
        //}
        return 0;
    }
};

/** @brief The main function.*/
int main()
{
    IsingModelTest test;

    test.test_init_functions();
    int res = 0;
    int tmp;
    for (size_t i=0; i < 1000; i++) {
        tmp = test.test_2x2_lattice(1e-2, 1e5);
        if (tmp == 0) {
            std::cout << "not enough cycles\n";
            break;
        }
        res += tmp;
    }
    
    std::cout << "Res: " << res / 1000 << std::endl;

    return 0;
}