GEIM

import numpy as np
import skfem  # for Finite Element Method
from skfem import Basis, ElementTriP0
from skfem import MeshTri
import matplotlib.pyplot as plt
import time
# ============================================================
# Param
# ============================================================
Nx = 50
Ny = 50

# Maillage triangulaire du carré (0,1)^2
m = MeshTri.init_tensor(
    np.linspace(0.0, 1.0, Nx + 1),
    np.linspace(0.0, 1.0, Ny + 1)
)


dx=1./Nx
dy=1./Ny

cell_areas=0.5*dx*dy

NumberOfModes = 15
NumberOfSnapshots = 20

MuMin = -2.0
MuMax = -1.0

rw = 0.15   # width of gaussian sensors

basis = Basis(m, ElementTriP0())   # P0 : one value for each element
dof_locs = basis.doflocs           # barycenters of elements for P0
ndof = dof_locs.shape[1]

# training set
mu1_train = np.linspace(MuMin, MuMax, NumberOfSnapshots)
mu2_train = np.linspace(MuMin, MuMax, NumberOfSnapshots)

# ============================================================
# Function g(x; mu1, mu2)
# ============================================================

def g_fun(x, mu1, mu2):
    return 1.0 / ((x[0] - mu1)**2 + (x[1] - mu2))

def project_g(mu1, mu2):
    return basis.project(lambda x: g_fun(x, mu1, mu2))

# ============================================================
# Gaussian kernel ell_m(x)
# ============================================================

def gaussian_weight_noC(x, xm, rw):
    """
    Gaussian centered at xm, without normalization constant.
    x: shape (2, npts)
    xm: shape (2,)
    """
    return np.exp(-0.5 * np.sum((x - xm[:, None])**2, axis=0) / rw**2)


# ============================================================
# Construction of continuous-inspired candidate sensors l_m
# ============================================================

def sigma_m_from_center(xm, rw):
    """
    Build the discrete representation of the continuous functional

        sigma_m(g) = \int_\Omega l_m(x) g dx

    for P0 functions, using barycenter quadrature:
        \int_K l_m(x) dx ~ |K| l_m(x_K)=w[K]

    Returns a weight vector w such that
        sigma_m(g_h) ~ sum_K w[K] * g_h[K]
    """
    l_vals = gaussian_weight_noC(dof_locs, xm, rw)   # l(x_K)
    w_K = cell_areas * l_vals                          # |K| * l(x_K)

    s = np.sum(w_K)
    if s > 1e-14:
        w_K = w_K / s   # normalize so that l_m(1) = 1 approximately
    return w_K

def apply_sigma_m(w, u):
    """
    Apply the functional represented by weights w to the discrete field u.
    """
    return np.dot(w, u)

# candidate sensors: one candidate centered at each barycenter
L_candidates = []
L_centers = []

for k in range(ndof):
    xm = dof_locs[:, k]
    wm = sigma_m_from_center(xm, rw)
    L_candidates.append(wm)
    L_centers.append((xm[0], xm[1]))

# ============================================================
# OFFLINE : initialization
# ============================================================

ns = NumberOfSnapshots
ntrain = ns * ns

GGrid = np.zeros((ndof, ntrain))
GinftyGrid = np.zeros(ntrain)
mu_pairs = []

col = 0
for a in mu1_train:
    for b in mu2_train:
        u = project_g(a, b)
        GGrid[:, col] = u
        GinftyGrid[col] = np.linalg.norm(u, np.inf)
        mu_pairs.append((a, b))
        col += 1

# first parameter
muIndex = np.argmax(GinftyGrid)
muIndexlist = [muIndex]

mu1_star, mu2_star = mu_pairs[muIndex]
g1 = GGrid[:, muIndex].copy()

# first magic sensor l_1: maximize |l(g1)|
scores = np.array([np.abs(apply_sigma_m(wm, g1)) for wm in L_candidates])
LIndex = np.argmax(scores)

LIndexlist = [LIndex]
Llist = [L_candidates[LIndex]]

print("First magic sensor center:", L_centers[LIndex])

# first basis function q1 normalized by l_1(q1)=1

q1 = g1 /  apply_sigma_m(Llist[0], g1)
qlist = [q1]

# ============================================================
# OFFLINE : GEIM iterations
# ============================================================

for M in range(1, NumberOfModes):
    print(f"M = {M}")

    # Interpolation matrix B_M with entries l_i(q_j)
    Qmat = np.zeros((M, M))
    for i in range(M):
        for j in range(M):
            Qmat[i, j] = apply_sigma_m(Llist[i], qlist[j])

    ResidualGrid = np.zeros((ndof, ntrain))
    ResidualInf = np.zeros(ntrain)

    # residual for each training parameter
    for col, (a, b) in enumerate(mu_pairs):
        u = GGrid[:, col].copy()

        rhs = np.array([apply_sigma_m(Llist[i], u) for i in range(M)])

        alpha = np.linalg.solve(Qmat, rhs)

        residual = u.copy()
        for i in range(M):
            residual -= alpha[i] * qlist[i]

        ResidualGrid[:, col] = residual
        ResidualInf[col] = np.linalg.norm(residual, np.inf)

    # avoid already selected parameters
    for old_idx in muIndexlist:
        ResidualInf[old_idx] = -1e30

    # new parameter
    muIndex = np.argmax(ResidualInf)
    muIndexlist.append(muIndex)

    mu1_star, mu2_star = mu_pairs[muIndex]
    rM = ResidualGrid[:, muIndex].copy()

    print("Selected parameter:", (mu1_star, mu2_star))
    print("Residual max norm :", np.max(ResidualInf))

    # choose new magic sensor maximizing |l(rM)|
    scores = np.array([abs(apply_sigma_m(w, rM)) for w in L_candidates])

    # avoid already selected sensors
    for old_L in LIndexlist:
        scores[old_L] = -np.inf

    LIndex = np.argmax(scores)
    LIndexlist.append(LIndex)
    Llist.append(L_candidates[LIndex])

    print("New magic sensor center:", L_centers[LIndex])

    # normalize new basis function by l_M(q_M)=1
    qM = rM / apply_sigma_m(Llist[-1], rM)
    qlist.append(qM)

# ============================================================
# ONLINE
# ============================================================

mu1_target = -0.5
mu2_target = -0.5

u_true = project_g(mu1_target, mu2_target)

# complete interpolation matrix
Qmat = np.zeros((NumberOfModes, NumberOfModes))
for i in range(NumberOfModes):
    for j in range(NumberOfModes):
        Qmat[i, j] = apply_sigma_m(Llist[i], qlist[j])

# rhs = [l_i(u_true)]
rhs = np.array([apply_sigma_m(Llist[i], u_true) for i in range(NumberOfModes)])

# coefficients GEIM
alpha = np.linalg.solve(Qmat, rhs)

# reconstruction GEIM
u_geim = np.zeros(ndof)
for j in range(NumberOfModes):
    u_geim += alpha[j] * qlist[j]

# ============================================================
# Test in one point
# ============================================================

xseek, yseek = 0.5, 0.5

element_finder = m.element_finder()
element_index = element_finder(np.array([xseek]), np.array([yseek]))[0]

approx_value = u_geim[element_index]
true_value_fem = u_true[element_index]
true_value_exact = g_fun(np.array([[xseek], [yseek]]), mu1_target, mu2_target)[0]

print()
print(f"Tested point : ({xseek}, {yseek})")
print(f"Approximation GEIM  : {approx_value}")
print(f"Projected FEM       : {true_value_fem}")
print(f"Exact value g(x,mu) : {true_value_exact}")
print(f"Error |GEIM - FEM|  : {abs(approx_value - true_value_fem)}")
print(f"Error |GEIM - exact|: {abs(approx_value - true_value_exact)}")