test_gradient_root.cpp

#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this
                          // in one cpp file
#include "catch.hpp"
#include "differentiator.h"
#include "optimizer.h"
#include "sample_functions.h"
#define _USE_MATH_DEFINES
#include <cmath>
#include <cstdint>
#include <fstream>
#include <iomanip>
#include <string>
#include <vector>

TEST_CASE("Minimization via gradient root finding on functions with "
          "multi-dimensional argument",
          "[gradient_root_finding]") {
    using namespace numerics;
    using std::cerr;
    using std::cout;
    using std::endl;
    FivePointDifferentiator fp;
    double diff_precision = 1e-6;
    GradientRootfinder optimizer(fp, diff_precision);

    Coordinate<double> position1(1), position2(2), ref2(2), res1(1), res2(2);

    double precision = 1e-9;
    double error = 1e-5;

    using funct = Function<Coordinate<double>, double>;

    SECTION("Testing 1D parabola") {
        funct& parabola = functions::parabola_1d;
        cerr << "Tested function: " << functions::parabola_1d << endl;
        position1[0] = 10;
        res1 = optimizer.optimize(parabola, position1, precision);
        REQUIRE(res1.l2Norm() < error);
        REQUIRE(std::abs(parabola(res1)) < error);
    }
    SECTION("Testing 1D higher-order parabola") {
        double ref_pos = sqrt(4.6 / 4);

        funct& higher_parabola = functions::higher_parabola_1d;
        cerr << "Tested function: " << functions::higher_parabola_1d << endl;
        position1[0] = 10;
        res1 = optimizer.optimize(higher_parabola, position1, precision);
        REQUIRE(std::abs(res1[0] - ref_pos) < error);

        position1[0] = -ref_pos + 0.5;
        res1 = optimizer.optimize(higher_parabola, position1, precision);
        REQUIRE(std::abs(res1[0] - ref_pos) < error);
    }
    SECTION("Testing 1D reference Polynomial") {
        funct& reference_polynomial = functions::ref_polynomial_1d;
        cerr << "Tested function: " << functions::ref_polynomial_1d << endl;
        position1[0] = 20;
        res1 = optimizer.optimize(reference_polynomial, position1, precision);
        REQUIRE(std::abs(reference_polynomial(res1)) < error);
    }

    SECTION("Testing narrow 1d Lennard Jones potential") {
        double rm = functions::narrow_LJ.get_rm();

        funct& narrow_LJ = functions::narrow_LJ;
        position1[0] = rm * 3;
        res1 = optimizer.optimize(narrow_LJ, position1, precision);
        REQUIRE(std::abs(res1[0] - rm) < error);

        position1[0] = rm / 2.0;
        res1 = optimizer.optimize(narrow_LJ, position1, precision);
        REQUIRE(std::abs(res1[0] - rm) < error);
    }

    SECTION("Testing wide 1d Lennard Jones potential") {
        double rm = functions::wide_LJ.get_rm();
        funct& wide_LJ = functions::wide_LJ;
        position1[0] = rm * 3;
        res1 = optimizer.optimize(wide_LJ, position1, precision);
        REQUIRE(std::abs(res1[0] - rm) < error);

        position1[0] = rm / 2.0;
        res1 = optimizer.optimize(wide_LJ, position1, precision);
        REQUIRE(std::abs(res1[0] - rm) < error);
    }

    SECTION("Testing 2D parabola") {
        funct& parabola2d = functions::parabola_2d;
        position2[0] = 10;
        position2[1] = -123.0;
        res2 = optimizer.optimize(parabola2d, position2, precision);
        REQUIRE(res2.l2Norm() < error);
        REQUIRE(std::abs(parabola2d(res2)) < error);
        position2[0] = -133;
        position2[1] = -0.1;
        res2 = optimizer.optimize(parabola2d, position2, precision);
        REQUIRE(res2.l2Norm() < error);
        REQUIRE(std::abs(parabola2d(res2)) < error);
    }
    SECTION("Testing 2D parabola") {
        funct& periodic_wave = functions::periodic_wave_2d;
        uint32_t k = 13;
        uint32_t l = 31;

        ref2[0] = (2 * k - 0.5) * M_PI;
        ref2[1] = (2 * l + 1.0) * M_PI;
        position2[0] = ref2[0] + 0.5;
        position2[1] = ref2[1] - 0.2;
        res2 = optimizer.optimize(periodic_wave, position2, precision);
        REQUIRE((res2 - ref2).l2Norm() < error);
        k = -4;
        l = 5;

        ref2[0] = (2 * k - 0.5) * M_PI;
        ref2[1] = (2 * l + 1.0) * M_PI;
        position2[0] = ref2[0] - 0.3;
        position2[1] = ref2[1] + 0.4;
        res2 = optimizer.optimize(periodic_wave, position2, precision);
        REQUIRE((res2 - ref2).l2Norm() < error);
    }
}