U-shaped Scaling Law

import numpy as np
from scipy.optimize import minimize

# EVOLVE-BLOCK-START

def scaling_law_func(data_points, params):
    """
    6-parameter U-shaped scaling law:
      y(x) = offset + α·(x−δ)^2 + β·exp(−γ·(x−δ)^2) + θ·(x−δ)
    with α>0, β<0, γ>0 enforced via exponentials.
    params = [offset, log α, δ, log|β|, log γ, θ]
    """
    X = np.asarray(data_points).ravel()
    p = np.asarray(params).ravel()
    offset = p[0]
    alpha  = np.exp(p[1])
    delta  = p[2]
    beta   = -np.exp(p[3])
    gamma  = np.exp(p[4]) + 1e-12
    tilt   = p[5]
    Z = X - delta
    return offset + alpha * Z**2 + beta * np.exp(-gamma * Z**2) + tilt * Z

def fit_scaling_law(data_points, loss_values):
    """
    Fit the 6-parameter law by minimizing MSE with data-informed
    initialization plus a few restarts.
    Returns the best parameter vector of length 6.
    """
    X = np.asarray(data_points).ravel()
    y = np.asarray(loss_values).ravel()

    # 1) Data-informed initial guess
    offset0 = np.median(y)
    i_peak  = np.argmax(y)                 # worst loss → center of dip
    delta0  = X[i_peak]
    beta0   = max(y[i_peak] - offset0, 1e-3)
    half_h  = offset0 + beta0 * 0.5
    left    = X[(X < delta0) & (y > half_h)]
    right   = X[(X > delta0) & (y > half_h)]
    x1      = left[-1]  if left.size  else X.min()
    x2      = right[0]  if right.size else X.max()
    span    = abs(x2 - x1)
    sigma0  = span/2 if span>1e-8 else (X.max() - X.min())/4
    gamma0  = 1.0 / (2 * sigma0**2 + 1e-6)
    span_sq = max((X.max()-delta0)**2, (delta0-X.min())**2)
    alpha0  = beta0 / (span_sq + 1e-6)
    tilt0   = 0.5 * (y.max() - y.min()) / (X.max() - X.min() + 1e-6)

    init0 = np.array([
        offset0,
        np.log(alpha0),
        delta0,
        np.log(beta0),
        np.log(gamma0),
        tilt0
    ])
    # 2) Zero fallback
    init1 = np.zeros(6)
    # 3) A few small random restarts
    rng = np.random.RandomState(0)
    inits = [init0, init1] + [init0 + rng.randn(6)*0.1 for _ in range(3)]

    # Bound δ within the observed X range
    bounds = [
        (None, None),         # offset
        (None, None),         # log α
        (X.min(), X.max()),   # δ
        (None, None),         # log|β|
        (None, None),         # log γ
        (None, None)          # tilt
    ]

    def mse_obj(p):
        pred = scaling_law_func(X, p)
        return np.mean((pred - y)**2)

    best_p, best_loss = None, np.inf
    for init in inits:
        res = minimize(
            mse_obj, init, method='L-BFGS-B',
            bounds=bounds,
            options={'maxiter': 500, 'ftol': 1e-12}
        )
        if res.success and res.fun < best_loss:
            best_loss, best_p = res.fun, res.x

    return best_p if best_p is not None else init1

# EVOLVE-BLOCK-END

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

def scaling_law_func(data_points, params):
    """
    6-parameter U-shaped law:
      y(x) = a*(x - h)^2 + B / (1 + ((x - mu)/s)^2) + c

    params:
      p[0] = log(a)      (a > 0)
      p[1] = h           (parabola center)
      p[2] = B           (bump amplitude, can be + or -)
      p[3] = mu          (bump center)
      p[4] = log(s)      (s > 0, bump width)
      p[5] = c           (vertical offset)
    """
    X = np.asarray(data_points).ravel()
    p = np.asarray(params).ravel()
    a  = np.exp(p[0])
    h  = p[1]
    B  = p[2]
    mu = p[3]
    s  = np.exp(p[4]) + 1e-12
    c  = p[5]
    z  = (X - mu) / s
    return a * (X - h)**2 + B / (1.0 + z*z) + c

def fit_scaling_law(data_points, loss_values):
    """
    Fit the 6-parameter U-shape law by:
      1) Robust global search (Tukey's biweight via Differential Evolution)
      2) Local robust refinement (L-BFGS-B on Tukey loss)
      3) Final MSE refinement (L-BFGS-B on MSE)
    """
    X = np.asarray(data_points).ravel()
    y = np.asarray(loss_values).ravel()
    y_std = np.std(y) + 1e-12

    # Tukey's biweight loss for robustness
    def tukey_obj(p):
        r    = scaling_law_func(X, p) - y
        c0   = 4.685 * y_std
        absr = np.abs(r)
        t    = r / c0
        rho  = np.empty_like(r)
        mask = absr <= c0
        rho[mask]  = (c0**2 / 6.0) * (1.0 - (1.0 - t[mask]**2)**3)
        rho[~mask] = (c0**2 / 6.0)
        return np.mean(rho)

    # Mean-squared-error for final polishing
    def mse_obj(p):
        r = scaling_law_func(X, p) - y
        return np.mean(r*r)

    # parameter bounds based on data ranges
    xmin, xmax = X.min(), X.max()
    ymin, ymax = y.min(), y.max()
    b = 2 * y_std
    bounds = [
        (-10.0, 10.0),               # log(a)
        (xmin - 1.0, xmax + 1.0),    # h
        (ymin - b, ymax + b),        # B
        (xmin - 1.0, xmax + 1.0),    # mu
        (-10.0,  5.0),               # log(s)
        (ymin - b, ymax + b)         # c
    ]

    # 1) Global search with Tukey loss
    de_res = differential_evolution(
        tukey_obj, bounds,
        maxiter=80, popsize=20,
        tol=1e-6, polish=True
    )

    # 2) Local robust refinement
    local_robust = minimize(
        tukey_obj, de_res.x,
        method='L-BFGS-B', bounds=bounds,
        options={'ftol':1e-12, 'maxiter':2000}
    )
    best_p = local_robust.x if local_robust.success else de_res.x

    # 3) Final MSE polishing
    final = minimize(
        mse_obj, best_p,
        method='L-BFGS-B', bounds=bounds,
        options={'ftol':1e-12, 'maxiter':1000}
    )
    return final.x if final.success else best_p
# EVOLVE-BLOCK-END

import numpy as np
from scipy.optimize import minimize

# EVOLVE-BLOCK-START
def scaling_law_func(data_points, params):
    """
    6-parameter U-shaped model: convex parabola + Cauchy bump
      f(x) = a*(x - mu)^2 + d*x + e + b / (1 + (|x - mu|/c)^2)
    params = [log_a, mu, b, log_c, d, e]
    """
    X = np.asarray(data_points).ravel()
    p = np.asarray(params).ravel()
    a  = np.exp(p[0])
    mu = p[1]
    b  = p[2]
    c  = np.exp(p[3]) + 1e-12
    d  = p[4]
    e  = p[5]
    dx = X - mu
    # convex parabola + Cauchy-style bump
    return a * dx * dx + d * X + e + b / (1.0 + (np.abs(dx) / c)**2)

def fit_scaling_law(data_points, loss_values):
    """
    Fit the 6-parameter Cauchy-bump U-shaped model:
      1) Quadratic baseline (a, d, e) via least squares on [X^2, X, 1]
      2) Bump center mu & amplitude b from largest positive residual
      3) Bump width c from IQR of X
      4) Multi-start L-BFGS-B to refine all parameters
    """
    X = np.asarray(data_points).ravel()
    y = np.asarray(loss_values).ravel()

    # 1) Quadratic baseline init: y ≈ a*X^2 + d*X + e
    A = np.vstack([X**2, X, np.ones_like(X)]).T
    coef, *_ = np.linalg.lstsq(A, y, rcond=None)
    a0 = max(coef[0], 1e-8)
    d0 = coef[1]
    e0 = coef[2]

    # 2) Residuals for bump init
    baseline = a0 * X**2 + d0 * X + e0
    resid    = y - baseline
    idx      = np.argmax(resid)
    mu0      = X[idx]
    b0       = resid[idx]

    # 3) Width init via interquartile range
    iqr = np.percentile(X, 75) - np.percentile(X, 25)
    c0  = max(iqr, 1e-6)

    # Pack into params: [log_a, mu, b, log_c, d, e]
    p0 = np.array([np.log(a0), mu0, b0, np.log(c0), d0, e0])

    # Parameter bounds
    bounds = [
        (-20, 20),             # log_a
        (X.min(), X.max()),    # mu
        (None, None),          # b
        (-20, 20),             # log_c
        (None, None),          # d
        (None, None)           # e
    ]

    # Objective: mean squared error
    def mse(p):
        return np.mean((scaling_law_func(X, p) - y)**2)

    # 4) Multi-start L-BFGS-B refinement
    best_p   = p0.copy()
    best_val = mse(best_p)
    # one baseline start + several small perturbations
    for init in [p0] + [p0 + np.random.randn(6) * 0.1 for _ in range(6)]:
        res = minimize(mse, init, method='L-BFGS-B',
                       bounds=bounds,
                       options={'maxiter': 1000, 'ftol': 1e-9})
        if res.success and res.fun < best_val:
            best_val, best_p = res.fun, res.x

    return best_p
# EVOLVE-BLOCK-END

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

def scaling_law_func(data_points, params):
    """
    6-parameter U-shaped law: quadratic baseline + skewed Gaussian dip/hump.
    y(x) = C + A*(x - M)^2 + B*(x - μ)*exp(-0.5*((x - μ)/σ)^2)
    params = [C, logA, M, B, μ, logσ]
    """
    X = np.asarray(data_points).ravel()
    C, logA, M, B, mu, logS = params
    A = np.exp(logA)
    S = np.exp(logS) + 1e-8
    # Quadratic baseline plus skewed Gaussian term for U-shape flexibility
    return C + A * (X - M)**2 + B * (X - mu) * np.exp(-0.5 * ((X - mu) / S)**2)

def fit_scaling_law(data_points, loss_values):
    """
    Fit the 6-parameter U-shaped law via robust least-squares with multiple starts.
    Uses a soft-L1 loss (pseudo-Huber) for outlier resilience.
    """
    X = np.asarray(data_points).ravel()
    y = np.asarray(loss_values).ravel()
    # Basic statistics for initialization and bounds
    XR = max(X.max() - X.min(), 1e-6)
    yr = max(y.max() - y.min(), 1e-6)
    C0     = np.median(y)
    logA0  = np.log(yr / (XR**2 + 1e-6) + 1e-6)
    M0     = X[np.argmax(y)]
    B0     = 0.0
    mu0    = X.mean()
    logS0  = np.log(XR / 4 + 1e-6)
    p0     = np.array([C0, logA0, M0, B0, mu0, logS0], dtype=float)

    # Parameter bounds
    lower = [y.min(),   -20, X.min() - XR, -2*yr, X.min() - XR, -10]
    upper = [y.max(),    20, X.max() + XR,  2*yr, X.max() + XR,  10]
    bounds = (lower, upper)

    # Residual function for least_squares
    def residuals(p):
        return scaling_law_func(X, p) - y

    # Multi-start to avoid local minima
    best_p, best_cost = p0, np.inf
    rng = np.random.default_rng(42)
    # try the base init plus a few noisy restarts
    for noise in np.vstack([np.zeros(6), rng.normal(scale=0.1, size=(4, 6))]):
        p_init = p0 + noise
        try:
            res = least_squares(
                residuals,
                p_init,
                bounds=bounds,
                loss='soft_l1',
                f_scale=yr * 0.5,
                max_nfev=1000
            )
            if res.success and res.cost < best_cost:
                best_cost, best_p = res.cost, res.x
        except Exception:
            pass

    return best_p
# EVOLVE-BLOCK-END