SLD - MoE Scaling Law - SLDAgent + GPT-5

All Runs (sorted by R²)

Best Run 1 R² = 0.961586

▼

Python

# EVOLVE-BLOCK-START
import numpy as np
from scipy.optimize import minimize

_D0 = 1e8
_EPS = 1e-12

def _softplus(x):
    return np.log1p(np.exp(-np.abs(x))) + np.maximum(x, 0.0)

def _isp(y):
    y = np.maximum(y, 1e-12)
    return np.log(np.expm1(y))

def scaling_law_func(data_points, params):
    # data_points: (N,2) -> [num_experts, dense_parameter_count]
    X = np.atleast_2d(np.asarray(data_points, dtype=np.float64))
    if X.shape[1] != 2:
        raise ValueError("Expected data_points with 2 features: [num_experts, dense_parameter_count].")
    E = X[:, 0].astype(np.float64)
    D = (X[:, 1].astype(np.float64) / _D0)

    p = np.asarray(params, dtype=np.float64)
    if p.ndim == 1: p = p[None, :]
    T, K = p.shape
    if K < 6:
        pad = np.zeros((T, 6), dtype=np.float64)
        pad[:, :K] = p
        if K < 1: pad[:, 0] = 2.0
        if K < 2: pad[:, 1] = _isp(1.0)
        if K < 3: pad[:, 2] = _isp(0.5)
        if K < 4: pad[:, 3] = _isp(0.6)
        if K < 5: pad[:, 4] = _isp(8.0)
        if K < 6: pad[:, 5] = _isp(1.0)
        p = pad

    L_inf = p[:, 0]
    A     = _softplus(p[:, 1])
    a     = _softplus(p[:, 2])
    b     = _softplus(p[:, 3])
    E0    = _softplus(p[:, 4]) + 1e-8
    gamma = _softplus(p[:, 5])

    lnD = np.log(D + _EPS)[:, None]
    lnE = np.log(E + _EPS)[:, None]

    # P_eff = D^a + gamma * (E^b / (1 + E/E0)), computed stably in log-space
    t1 = a[None, :] * lnD
    t2 = np.log(gamma[None, :]) + b[None, :] * lnE - np.log1p(E[:, None] / E0[None, :])
    log_Peff = np.logaddexp(t1, t2)

    pred = L_inf[None, :] + A[None, :] * np.exp(-log_Peff)  # L = L_inf + A / P_eff
    return pred[:, 0] if pred.shape[1] == 1 else pred

def fit_scaling_law(data_points, loss_values):
    X = np.atleast_2d(np.asarray(data_points, dtype=np.float64))
    y = np.asarray(loss_values, dtype=np.float64)
    if X.shape[1] != 2:
        raise ValueError("Expected data_points with 2 features: [num_experts, dense_parameter_count].")
    y2d = y[:, None] if y.ndim == 1 else y
    T = y2d.shape[1]

    E = X[:, 0].astype(np.float64)
    D = (X[:, 1].astype(np.float64) / _D0)
    lnD = np.log(D + _EPS)
    lnE = np.log(E + _EPS)

    def solve_LA(Peff, yt):
        # Ridge-stabilized least squares for [L_inf, A] in yt ≈ L_inf + A / Peff
        Xd = np.stack([np.ones_like(Peff), 1.0 / (Peff + 1e-12)], axis=1)
        lam = 1e-6
        XtX = Xd.T @ Xd
        XtX[0,0] += lam; XtX[1,1] += lam
        sol = np.linalg.solve(XtX, Xd.T @ yt)
        return float(sol[0]), float(sol[1])

    def obj_core(q, yt, delta):
        a = _softplus(q[0]); b = _softplus(q[1]); E0 = _softplus(q[2]) + 1e-8; g = _softplus(q[3])
        t1 = a * lnD
        t2 = np.log(g) + b * lnE - np.log1p(E / E0)
        log_Peff = np.logaddexp(t1, t2)
        Peff = np.exp(log_Peff)
        L_opt, A_opt = solve_LA(Peff, yt)
        r = (L_opt + A_opt / (Peff + 1e-12)) - yt
        # pseudo-Huber
        loss = np.mean(delta**2 * (np.sqrt(1.0 + (r / delta)**2) - 1.0))
        # priors to keep parameters physically plausible
        reg = 1e-4 * ((a - 0.5)**2 + (b - 0.6)**2 + (np.log(E0) - np.log(8.0))**2 + (g - 1.0)**2)
        reg += 5e-5 * (L_opt - (float(np.min(yt)) - 0.05))**2
        return loss + reg, L_opt, A_opt

    rng = np.random.default_rng(42)
    out = np.zeros((T, 6), dtype=np.float64)

    for t in range(T):
        yt = y2d[:, t]
        delta = max(0.08, 0.15 * float(np.std(yt)))
        base_q = np.array([_isp(0.5), _isp(0.6), _isp(8.0), _isp(1.0)], dtype=np.float64)

        def objective(q):
            v, _, _ = obj_core(q, yt, delta)
            return v

        best_q = base_q.copy()
        best_v = np.inf
        inits = [base_q.copy()]
        for _ in range(8):
            init = base_q + rng.normal(0.0, 0.3, size=4)
            inits.append(init)
        for E0_try in (4.0, 12.0, 16.0, 24.0):
            for g_try in (0.7, 1.0, 1.4):
                for a_try in (0.4, 0.5, 0.6):
                    for b_try in (0.5, 0.6, 0.7):
                        init = np.array([_isp(a_try), _isp(b_try), _isp(E0_try), _isp(g_try)], dtype=np.float64)
                        inits.append(init)

        for init in inits:
            res = minimize(objective, init, method='L-BFGS-B', options={'maxiter': 800})
            q_cand = res.x if res.success else init
            v, L_opt, A_opt = obj_core(q_cand, yt, delta)
            if v < best_v:
                best_v = v
                best_q = q_cand
                best_L = L_opt
                best_A = A_opt

        # Compose full raw parameter vector: [L_inf_raw, A_raw, qa, qb, qE0, qgamma]
        out[t] = np.array([best_L, _isp(best_A), best_q[0], best_q[1], best_q[2], best_q[3]], dtype=np.float64)

    return out[0] if T == 1 else out
# EVOLVE-BLOCK-END

#2 Run 3 R² = 0.960711

▼

Python

# EVOLVE-BLOCK-START
import numpy as np
from scipy.optimize import minimize

def scaling_law_func(data_points, params):
    X = np.atleast_2d(np.asarray(data_points, dtype=float))
    if X.shape[1] != 2:
        raise ValueError("data_points must be (N,2): [num_experts, dense_parameter_count]")
    # Normalize and smooth features
    E = np.clip(X[:, 0].astype(float), 1.0, None)
    Pn = np.clip(X[:, 1].astype(float), 1.0, None) / 1e8  # roughly in [1,8]
    Es = E + 0.5  # smooth low-expert regime

    p = np.atleast_2d(np.asarray(params, dtype=float))
    if p.shape[1] < 6:
        p = np.concatenate([p, np.zeros((p.shape[0], 6 - p.shape[1]))], 1)
    p = p[:, :6]
    # Parameters: [L0, A, uP, a, uE, b]
    L0 = p[:, 0]
    A  = np.clip(p[:, 1], 0.0, None)
    uP = np.clip(p[:, 2], 1e-12, None)
    a  = np.clip(p[:, 3], 1e-8, None)
    uE = np.clip(p[:, 4], 1e-12, None)
    b  = np.clip(p[:, 5], 1e-8, None)

    lnP = np.log(Pn + 1e-12)[None, :]
    lnE = np.log(Es + 1e-12)[None, :]
    Ppow = np.exp(a[:, None] * lnP)
    Epow = np.exp(b[:, None] * lnE)

    # Multiplicative saturation: synergistic yet stable
    sP = 1.0 + uP[:, None] * Ppow
    sE = 1.0 + uE[:, None] * Epow
    S  = sP * sE

    pred = (L0[:, None] + A[:, None] / S).T
    return pred[:, 0] if pred.shape[1] == 1 else pred


def fit_scaling_law(data_points, loss_values):
    X = np.atleast_2d(np.asarray(data_points, dtype=float))
    y = np.asarray(loss_values, dtype=float)
    if X.shape[1] != 2:
        raise ValueError("data_points must be (N,2): [num_experts, dense_parameter_count]")
    y2d = y[:, None] if y.ndim == 1 else y
    N, T = y2d.shape

    E = np.clip(X[:, 0].astype(float), 1.0, None)
    Pn = np.clip(X[:, 1].astype(float), 1.0, None) / 1e8
    Es = E + 0.5
    lnP = np.log(Pn + 1e-12)
    lnE = np.log(Es + 1e-12)

    # Parameter bounds: [L0, A, uP, a, uE, b]
    lower = np.array([0.8, 0.05, 1e-12, 0.05, 1e-12, 0.05], dtype=float)
    upper = np.array([5.0, 10.0, 200.0, 2.5, 200.0, 2.5], dtype=float)
    bnds = list(zip(lower, upper))
    rng = np.random.default_rng(123)

    def pseudo_huber_grad(r, d):
        s = r / d
        return r / np.sqrt(1.0 + s * s)

    def obj_grad(p, yt, d):
        L0, A, uP, a, uE, b = p
        A = max(A, 0.0)
        uP = max(uP, 1e-12)
        uE = max(uE, 1e-12)
        a  = max(a, 1e-8)
        b  = max(b, 1e-8)

        Ppow = np.exp(a * lnP)
        Epow = np.exp(b * lnE)
        sP = 1.0 + uP * Ppow
        sE = 1.0 + uE * Epow
        S  = sP * sE

        pred = L0 + A / S
        r = pred - yt
        g = pseudo_huber_grad(r, d)

        invS = 1.0 / S
        invS2 = invS * invS
        dL0 = np.mean(g)
        dA  = np.mean(g * invS)

        # dS terms
        dS_uP = Ppow * sE
        dS_uE = Epow * sP
        dS_a  = uP * Ppow * lnP * sE
        dS_b  = uE * Epow * lnE * sP

        G = -A * invS2
        duP = np.mean(g * G * dS_uP)
        duE = np.mean(g * G * dS_uE)
        da  = np.mean(g * G * dS_a)
        db  = np.mean(g * G * dS_b)

        # Regularization (stability and parsimony)
        lam_exp = 1e-4
        lam_u   = 1e-5
        loss = np.mean(d**2 * (np.sqrt(1.0 + (r / d) ** 2) - 1.0)) \
             + lam_exp * (a * a + b * b) + lam_u * (uP * uP + uE * uE)
        grad = np.array([dL0, dA, duP + 2 * lam_u * uP, da + 2 * lam_exp * a,
                         duE + 2 * lam_u * uE, db + 2 * lam_exp * b], dtype=float)
        return loss, grad

    params_opt = np.zeros((T, 6), dtype=float)

    # Small grid for exponent initialization
    a_grid = np.array([0.4, 0.7, 1.0, 1.3])
    b_grid = np.array([0.4, 0.7, 1.0, 1.3])

    for t in range(T):
        yt = y2d[:, t]
        yr = max(np.ptp(yt), 1e-3)
        d  = 0.25 * np.std(yt) + 0.03

        # Base initialization
        L0 = float(np.clip(np.min(yt) - 0.05 * yr, lower[0], upper[0]))
        A0 = float(np.clip(0.8 * yr, lower[1], upper[1]))

        # Choose a,b via small grid using linearized S ≈ 1 + uP P^a + uE E^b
        best_seed = None
        best_val  = np.inf
        for a0 in a_grid:
            for b0 in b_grid:
                Ppow0 = np.exp(a0 * lnP)
                Epow0 = np.exp(b0 * lnE)
                rhs = np.maximum(A0 / np.maximum(yt - L0, 0.1) - 1.0, 0.0)
                M = np.stack([Ppow0, Epow0], 1)
                try:
                    u_est, *_ = np.linalg.lstsq(M, rhs, rcond=None)
                except np.linalg.LinAlgError:
                    u_est = np.array([0.5, 0.1])
                uP0 = float(np.clip(u_est[0], lower[2], upper[2]))
                uE0 = float(np.clip(u_est[1], lower[4], upper[4]))
                seed = np.array([L0, A0, uP0, a0, uE0, b0], dtype=float)
                val = obj_grad(seed, yt, d)[0]
                if val < best_val:
                    best_val, best_seed = val, seed

        seed = best_seed

        # Multi-start local refinement
        best_p, best_v = seed.copy(), obj_grad(seed, yt, d)[0]
        inits = [seed]
        for _ in range(6):
            init = seed.copy()
            init[0] = rng.uniform(lower[0], upper[0])
            init[1] *= rng.uniform(0.7, 1.4)
            init[2] *= rng.uniform(0.5, 1.8)
            init[4] *= rng.uniform(0.5, 1.8)
            init[3] = rng.uniform(lower[3], upper[3])
            init[5] = rng.uniform(lower[5], upper[5])
            inits.append(np.minimum(np.maximum(init, lower), upper))

        for init in inits:
            res = minimize(lambda p: obj_grad(p, yt, d)[0], init,
                           jac=lambda p: obj_grad(p, yt, d)[1],
                           method='L-BFGS-B', bounds=bnds,
                           options={'maxiter': 600, 'ftol': 1e-9})
            val = res.fun if res.success else obj_grad(init, yt, d)[0]
            cand = res.x if res.success else init
            if val < best_v:
                best_v, best_p = val, cand

        params_opt[t] = best_p

    return params_opt[0] if T == 1 else params_opt
# EVOLVE-BLOCK-END

#3 Run 5 R² = 0.957837

▼

Python

# EVOLVE-BLOCK-START
import numpy as np
from scipy.optimize import least_squares

def scaling_law_func(data_points, params):
    X = np.atleast_2d(np.asarray(data_points))
    e = np.maximum(X[:, 0].astype(float) / 8.0, 1e-12)   # experts ~[0.125, 8]
    p = np.maximum(X[:, 1].astype(float) / 1e8, 1e-12)   # dense params ~[1, 8]

    params = np.asarray(params)
    if params.ndim == 1:
        params = params[None, :]
    T = params.shape[0]

    c = params[:, 0]
    A = params[:, 1]
    a = params[:, 2]
    s = params[:, 3]
    u = params[:, 4]
    v = params[:, 5] if params.shape[1] > 5 else np.ones(T)

    # Additive capacity: dense term plus saturating expert contribution
    e_pow = e[:, None] ** v[None, :]
    sat_e = e_pow / (1.0 + e_pow)                   # bounded in [0,1)
    cap = (p[:, None] ** u[None, :]) + s[None, :] * sat_e
    cap = np.maximum(cap, 1e-12)

    pred = c[None, :] + A[None, :] * (cap ** (-a[None, :]))
    return pred[:, 0] if pred.shape[1] == 1 else pred

def fit_scaling_law(data_points, loss_values):
    X = np.atleast_2d(np.asarray(data_points))
    y = np.asarray(loss_values)
    if y.ndim == 1:
        y = y[:, None]
    T = y.shape[1]

    def softplus(x):
        return np.log1p(np.exp(-np.abs(x))) + np.maximum(x, 0.0)
    def inv_softplus(z):
        z = np.maximum(z, 1e-12)
        return np.log(np.expm1(z))

    def run_fit(yvec):
        yvec = yvec.ravel()
        y_min = float(np.min(yvec))
        y_max = float(np.max(yvec))

        # Seeds: [c, A, alpha, s, u, v]
        base = [
            np.array([y_min - 0.05, max(y_max - y_min, 0.12), 0.60, 1.00, 0.95, 0.85]),
            np.array([y_min,        max(y_max - y_min, 0.20), 0.50, 0.80, 0.85, 1.00]),
            np.array([y_min + 0.10, max(y_max - y_min, 0.18), 0.70, 1.30, 1.10, 0.75]),
            np.array([y_min - 0.08, max(y_max - y_min, 0.25), 0.40, 0.65, 0.70, 1.20]),
        ]
        rng = np.random.RandomState(0)
        for _ in range(4):
            jitter = np.array([
                rng.normal(0, 0.05),                    # c
                rng.lognormal(mean=-1.9, sigma=0.35),   # A
                rng.lognormal(mean=-0.7, sigma=0.25),   # alpha
                rng.lognormal(mean=-0.3, sigma=0.3),    # s
                rng.lognormal(mean=-0.1, sigma=0.25),   # u
                rng.lognormal(mean=-0.2, sigma=0.25),   # v
            ])
            base.append(np.array([y_min, 0.22, 0.55, 0.95, 0.95, 0.90]) + jitter)

        def pack(theta):
            c, A, a, s, u, v = theta
            return np.array([c, inv_softplus(A), inv_softplus(a),
                             inv_softplus(s), inv_softplus(u), inv_softplus(v)])
        def unpack(w):
            c = w[0]
            return np.array([c, softplus(w[1]), softplus(w[2]),
                             softplus(w[3]), softplus(w[4]), softplus(w[5])])

        def residuals(w):
            p = unpack(w)
            pred = scaling_law_func(X, p)
            r = pred - yvec
            # Regularization as residuals (sqrt-lambda scaling)
            c, A, a, s, u, v = p
            reg = []
            lam = 1e-4
            reg.extend([np.sqrt(lam)*A, np.sqrt(lam)*a, np.sqrt(lam)*s, np.sqrt(lam)*u, np.sqrt(lam)*v])
            lam_bias = 1e-4
            reg.extend([np.sqrt(lam_bias)*(u - 1.0), np.sqrt(lam_bias)*(v - 0.9)])
            lam_c = 5e-4
            reg.extend([np.sqrt(lam_c)*np.maximum(1.25 - c, 0.0),
                        np.sqrt(lam_c)*np.maximum(c - 5.0, 0.0)])
            return np.hstack([r, np.array(reg, dtype=float)])

        best_w, best_cost = None, np.inf
        lb = np.array([1.0, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf])
        ub = np.array([5.5,  np.inf,  np.inf,  np.inf,  np.inf,  np.inf])
        for th in base:
            w0 = pack(th)
            res = least_squares(residuals, w0, method='trf', loss='soft_l1', f_scale=0.22, bounds=(lb, ub))
            if res.cost < best_cost:
                best_cost, best_w = res.cost, res.x
        return unpack(best_w)

    params = [run_fit(y[:, t]) for t in range(T)]
    return params[0] if T == 1 else np.vstack(params)
# EVOLVE-BLOCK-END

#4 Run 2 R² = 0.955857

▼

Python

# EVOLVE-BLOCK-START
import numpy as np
from scipy.optimize import minimize

def scaling_law_func(data_points, params):
    # L(P,E) = L_inf + A * (P/1e8)^(-alpha) * (1 + k * s(E))^(-gamma)
    # s(E) = E^eta / (1 + E^eta), smooth, monotone, saturating in [0,1)
    X = np.atleast_2d(np.asarray(data_points, dtype=float))
    E = np.clip(X[:, 0], 1.0, None)
    Pn = np.clip(X[:, 1] / 1e8, 1e-12, None)

    p = np.asarray(params, dtype=float).reshape(-1)
    if p.size < 6:
        pad = np.array([0.0, 1.0, 0.5, 0.0, 0.8, 1.0], dtype=float)
        pad[:p.size] = p
        p = pad
    L_inf, A, alpha, k, eta, gamma = p

    A = max(A, 1e-12)
    alpha = np.clip(alpha, 1e-3, 5.0)
    eta = np.clip(eta, 0.05, 5.0)
    gamma = np.clip(gamma, 0.0, 6.0)
    k = np.clip(k, -0.95, 50.0)  # ensure 1 + k*s(E) stays positive

    t = np.clip(eta * np.log(E), -60.0, 60.0)
    sE = np.exp(t) / (1.0 + np.exp(t))
    fac = np.clip(1.0 + k * sE, 1e-8, None)

    pred = L_inf + A * np.exp(-alpha * np.log(Pn)) * np.exp(-gamma * np.log(fac))
    return pred


def fit_scaling_law(data_points, loss_values):
    X = np.atleast_2d(np.asarray(data_points, dtype=float))
    y = np.asarray(loss_values, dtype=float).reshape(-1)
    E = np.clip(X[:, 0], 1.0, None)
    Pn = np.clip(X[:, 1] / 1e8, 1e-12, None)

    def pseudo_huber(r, d):
        z = r / d
        return d * d * (np.sqrt(1.0 + z * z) - 1.0)

    ymin, ymax = float(np.min(y)), float(np.max(y))
    L0 = ymin - 0.05
    f = np.clip(y - L0, 1e-6, None)
    msk = np.isfinite(Pn) & np.isfinite(f)
    if msk.sum() >= 3:
        c1, c0 = np.polyfit(np.log(Pn[msk]), np.log(f[msk]), 1)
        alpha0 = float(np.clip(-c1, 0.05, 3.0))
        A0 = float(np.exp(c0))
    else:
        alpha0, A0 = 0.6, max(ymax - L0, 0.2)

    bounds = [
        (ymin - 0.8, ymax + 0.2),  # L_inf
        (1e-6, 100.0),             # A
        (0.05, 3.5),               # alpha
        (-0.95, 50.0),             # k
        (0.05, 3.0),               # eta
        (0.0, 6.0)                 # gamma
    ]

    delta = max(0.06, 0.2 * float(np.std(y)))
    lam = 1e-6

    def obj(theta):
        pred = scaling_law_func(X, theta)
        r = pred - y
        reg = lam * np.sum(theta[1:] ** 2)
        return float(np.mean(pseudo_huber(r, delta)) + reg)

    seeds = []
    L_seeds = [L0, 0.5 * (ymin + ymax), ymin + 0.1]
    k_seeds = [-0.5, 0.0, 0.5, 1.5]
    eta_seeds = [0.6, 0.9, 1.2]
    gam_seeds = [0.5, 1.0, 2.0]
    for Ls in L_seeds:
        Aseed = np.clip(max(ymax - Ls, 0.15), 0.1, 20.0)
        for kk in k_seeds:
            for et in eta_seeds:
                for gm in gam_seeds:
                    seeds.append([Ls, Aseed, alpha0, kk, et, gm])

    rng = np.random.default_rng(123)
    base = np.array([L0, A0, alpha0, 0.6, 0.9, 1.2])
    for _ in range(8):
        j = rng.normal(0, [0.06, 0.3, 0.08, 0.5, 0.25, 0.8])
        s = np.clip(base + j, [b[0] for b in bounds], [b[1] for b in bounds])
        seeds.append(s)

    best_val = np.inf
    best_theta = np.array([L0, A0, alpha0, 0.0, 0.8, 1.0])
    for s in seeds:
        try:
            res = minimize(obj, s, method='L-BFGS-B', bounds=bounds,
                           options={'maxiter': 600, 'ftol': 1e-9})
            th = res.x if res.success else np.asarray(s, dtype=float)
            val = obj(th)
        except Exception:
            th = np.asarray(s, dtype=float)
            val = obj(th)
        if val < best_val:
            best_val, best_theta = val, th

    best_theta[0] = np.clip(best_theta[0], ymin - 1.0, ymax + 1.0)
    return best_theta
# EVOLVE-BLOCK-END

#5 Run 4 R² = 0.952085

▼

Python

# EVOLVE-BLOCK-START
import numpy as np
from scipy.optimize import minimize

_P0 = 1e8

def scaling_law_func(data_points, params):
    X = np.atleast_2d(np.asarray(data_points, dtype=float))
    e, p = X[:, 0], X[:, 1]
    p_norm = np.clip(p / _P0, 1e-12, None)
    e_eff  = np.clip(np.log1p(e), 1e-12, None)

    pr = np.asarray(params, dtype=float)
    if pr.ndim == 1: pr = pr[None, :]
    T, K = pr.shape
    if K < 6: pr = np.pad(pr, ((0, 0), (0, 6 - K)), mode='constant')
    elif K > 6: pr = pr[:, :6]

    a, b, d, alpha, beta, m = (pr[:, i] for i in range(6))
    alpha = np.maximum(alpha, 1e-12)
    beta  = np.maximum(beta,  1e-12)
    m     = np.maximum(m,     0.0)

    Xp = p_norm[:, None] ** (-alpha[None, :])
    G  = (e_eff[:,  None] + m[None, :]) ** (-beta[None, :])
    pred = a[None, :] + b[None, :]*Xp + d[None, :]*Xp*G
    return pred[:, 0] if pred.shape[1] == 1 else pred


def fit_scaling_law(data_points, loss_values):
    X = np.atleast_2d(np.asarray(data_points, dtype=float))
    y = np.asarray(loss_values, dtype=float)
    if y.ndim == 1: y = y[:, None]
    N, T = y.shape

    p_norm = np.clip(X[:, 1] / _P0, 1e-12, None)
    e_eff  = np.clip(np.log1p(X[:, 0]), 1e-12, None)

    def pseudo_huber(r, d): return d*d * (np.sqrt(1.0 + (r/d)**2) - 1.0)
    def mad(x): 
        m = np.median(x)
        return np.median(np.abs(x - m))

    def solve_abd(alpha, beta, m, yt, lam=1e-6, delta=None, iters=2):
        Xp  = p_norm ** (-alpha)
        G   = (e_eff + m) ** (-beta)
        Xpg = Xp * G
        F   = np.column_stack([np.ones(N), Xp, Xpg])
        w   = np.ones(N, dtype=float)
        for _ in range(max(1, iters)):
            W   = np.diag(w)
            C   = F.T @ W @ F + lam * np.eye(3)
            rhs = F.T @ W @ yt
            try: a_, b_, d_ = map(float, np.linalg.solve(C, rhs))
            except np.linalg.LinAlgError: a_, b_, d_ = map(float, np.linalg.lstsq(C, rhs, rcond=None)[0])
            b_ = max(b_, 0.0); d_ = max(d_, 0.0)
            r  = yt - (a_ + b_*Xp + d_*Xpg)
            di = max(0.12, 1.4826 * mad(r)) if delta is None else delta
            w  = 1.0 / np.maximum(1.0, np.abs(r) / di)
        a_ = float(np.clip(np.average(yt - (b_*Xp + d_*Xpg), weights=w), 1.2, 5.0))
        return a_, b_, d_, Xp, G, Xpg

    def make_obj(yt):
        base = max(0.15, 0.5 * 1.4826 * mad(yt))
        def obj(x):
            alpha = np.clip(x[0], 1e-6, 2.5)
            beta  = np.clip(x[1], 1e-6, 2.5)
            m     = np.clip(x[2], 0.0, 2.0)
            a, b, d, Xp, G, Xpg = solve_abd(alpha, beta, m, yt, lam=1e-6, delta=base, iters=2)
            r = (a + b*Xp + d*Xpg) - yt
            reg = 1e-4*(alpha*alpha + beta*beta) + 1e-5*m + 1e-7*(b*b + d*d)
            return float(np.mean(pseudo_huber(r, base))) + reg
        return obj

    grid_a = np.array([0.1, 0.2, 0.3, 0.5, 0.8, 1.2])
    grid_b = np.array([0.1, 0.2, 0.3, 0.5, 0.8, 1.2])
    grid_m = np.array([0.0, 0.2, 0.5, 1.0])

    params_out = np.zeros((T, 6), dtype=float)
    rng = np.random.default_rng(123)

    for t in range(T):
        yt = y[:, t]
        obj = make_obj(yt)
        best = np.inf; x_best = np.array([0.5, 0.5, 0.25], dtype=float)
        for aa in grid_a:
            for bb in grid_b:
                for mm in grid_m:
                    x0 = np.array([aa, bb, mm], dtype=float)
                    v = obj(x0)
                    if v < best: best, x_best = v, x0

        bounds = [(1e-6, 2.5), (1e-6, 2.5), (0.0, 2.0)]
        starts = [x_best, [0.3,0.4,0.2], [0.8,0.6,0.5], [1.2,0.8,0.3]] + \
                 [rng.uniform([0.1,0.1,0.0],[1.5,1.5,1.2]) for _ in range(6)]

        best = np.inf; xb = x_best.copy()
        for s in starts:
            x0 = np.array(s, dtype=float)
            res = minimize(obj, x0, method='L-BFGS-B', bounds=bounds,
                           options={'maxiter': 400, 'ftol': 1e-9})
            val = res.fun if res.success else obj(x0)
            if val < best: best, xb = val, (res.x if res.success else x0)

        alpha = float(np.clip(xb[0], 1e-6, 2.5))
        beta  = float(np.clip(xb[1], 1e-6, 2.5))
        m     = float(np.clip(xb[2], 0.0, 2.0))
        a, b, d, _, _, _ = solve_abd(alpha, beta, m, yt, lam=1e-6, iters=3)
        params_out[t] = np.array([a, b, d, alpha, beta, m], dtype=float)

    return params_out[0] if T == 1 else params_out
# EVOLVE-BLOCK-END