experiments/tests/test_minimisation.py

import numpy as np
from dolfinx import fem, mesh
from mpi4py import MPI

from dolfinx_access import eval_Function_f64, min_Function_f64_p2


# A simple test that the finite-dimensional subproblem is solved correctly
def test_minimisation(domain_dims=[[-5, 5], [-5, 5]], degree=2, nx=256, ny=256, n=512):
    domain = mesh.create_rectangle(
        MPI.COMM_SELF,
        [
            np.array([domain_dims[0][0], domain_dims[1][0]]),
            np.array([domain_dims[0][1], domain_dims[1][1]]),
        ],
        [nx, ny],
        mesh.CellType.triangle,
    )
    V = fem.functionspace(domain, ("Lagrange", degree))
    test = fem.Function(V)

    # Beale function
    def beale(x):
        return (
            (1.5 - x[0] + x[0] * x[1]) ** 2
            + (2.25 - x[0] + x[0] * (x[1] ** 2)) ** 2
            + (2.625 - x[0] + x[0] * (x[1] ** 3)) ** 2
        )

    def quadratic(x):
        return x[0] ** 2 - 5 * x[1] ** 2 + x[0] * x[1] - x[0] + x[1]

    def sphere(x):
        return x[0] ** 2 + x[1] ** 2

    def booth(x):
        return (x[0] + 2 * x[1] - 7) ** 2 + (2 * x[0] + x[1] - 5) ** 2

    def himmelblau(x):
        return (x[0] ** 2 + x[1] - 11) ** 2 + (x[0] + x[1] ** 2 - 7) ** 2

    for testf, name in [
        (beale, "Beale"),
        (quadratic, "Quadratic"),
        (sphere, "Sphere"),
        (booth, "Booth"),
        (himmelblau, "Himmelblau"),
    ]:
        test.interpolate(testf)
        m = min_Function_f64_p2(test._cpp_object)
        print(
            name,
            "minimizer:",
            m.v,
            "(verification ",
            eval_Function_f64(test._cpp_object, [m.x[0], m.x[1], 0]),
            " original ",
            testf(m.x),
            ") @ ",
            m.x,
        )

        dist = 0.0
        m = float("inf")
        m0 = float("inf")
        mp = None
        for x in np.linspace(domain_dims[0][0], domain_dims[0][1], n):
            for y in np.linspace(domain_dims[1][0], domain_dims[1][1], n):
                p = [x, y, 0.0]
                v = eval_Function_f64(test._cpp_object, p)
                v0 = testf(p)
                if v0 < m0:
                    m = v
                    m0 = v0
                    mp = p
                d = np.linalg.norm(v - v0)
                dist = max(d, dist)
        print(name, "evaluation max discrepancy:", dist)
        print(
            name,
            "rough scan minimiser ",
            m,
            "(model",
            m0,
            ") @ [",
            float(mp[0]),
            ", ",
            float(mp[1]),
            "]",
        )
        print()


test_minimisation()
exit(0)