DCT Laboratory — Volume II, Chapter 10

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi']=110

import numpy as np
SEED = 26210
# --- memberships: Low tri(0,0,5), Med tri(0,5,10), High tri(5,10,10) on [0,10] ---
def mu_L(x): return max(0.0, (5-x)/5) if x <= 5 else 0.0
def mu_M(x): return max(0.0, x/5) if x <= 5 else max(0.0, (10-x)/5)
def mu_H(x): return max(0.0, (x-5)/5) if x >= 5 else 0.0
# --- hand-written (expert) Sugeno consequents ---
EXPERT = {"L": (1.0, 0.2), "M": (3.0, 0.6), "H": (2.0, 1.0)}
def sugeno(x, params):
    ws = np.array([mu_L(x), mu_M(x), mu_H(x)])
    ys = np.array([a + b*x for a, b in (params["L"], params["M"], params["H"])])
    return float(ws @ ys / ws.sum())
# --- crisp switch policy (threshold at x = 5 between M and H rules) ---
def crisp(x, params):
    a, b = params["M"] if x < 5 else params["H"]
    return a + b*x
# --- the ANFIS least-squares step against an enterprise cost curve ---
def target(x): return 20 - 3*x + 0.35*x*x
XS = np.arange(0, 11, 1.0)
def fit_lse():
    rows = []
    for x in XS:
        ws = np.array([mu_L(x), mu_M(x), mu_H(x)]); ws = ws/ws.sum()
        rows.append([ws[0], ws[0]*x, ws[1], ws[1]*x, ws[2], ws[2]*x])
    A = np.array(rows); y = target(XS)
    theta, *_ = np.linalg.lstsq(A, y, rcond=None)
    theta = np.round(theta, 4)          # the rounded params ARE the deliverable model
    return {"L": (theta[0], theta[1]), "M": (theta[2], theta[3]), "H": (theta[4], theta[5])}
def rmse(params):
    return float(np.sqrt(np.mean([(sugeno(x, params)-target(x))**2 for x in XS])))

def reference_values():
    fit = fit_lse()
    r0, r1 = rmse(EXPERT), rmse(fit)
    return {
        "mu_L_3": round(mu_L(3),4), "mu_M_3": round(mu_M(3),4), "mu_H_3": round(mu_H(3),4),
        "partition_sum_3": round(mu_L(3)+mu_M(3)+mu_H(3),4),
        "yhat_3": round(sugeno(3, EXPERT),4), "yhat_7": round(sugeno(7, EXPERT),4),
        "crisp_jump": round(abs(crisp(5.1, EXPERT)-crisp(4.9, EXPERT)),4),
        "fuzzy_delta": round(abs(sugeno(5.1, EXPERT)-sugeno(4.9, EXPERT)),4),
        "robustness_ratio": round(abs(crisp(5.1, EXPERT)-crisp(4.9, EXPERT))/abs(sugeno(5.1, EXPERT)-sugeno(4.9, EXPERT)),4),
        "rmse_expert": round(r0,4), "rmse_fitted": round(r1,4),
        "yfit_3": round(sugeno(3, fit_lse()),4), "yfit_7": round(sugeno(7, fit_lse()),4),
    }
if __name__ == "__main__":
    fit = fit_lse()
    print("fitted consequents:", {k:(round(a,4),round(b,4)) for k,(a,b) in fit.items()})
    [print(f"{k:20s} {v}") for k,v in reference_values().items()]

xs = np.linspace(0, 10, 400)
fig, ax = plt.subplots(figsize=(7.8,3.8))
for f, lab, c in ((mu_L,"Low","#8A8F8B"),(mu_M,"Medium","#C8A24B"),(mu_H,"High","#0B3D2E")):
    ax.plot(xs, [f(x) for x in xs], lw=2.2, label=lab, c=c)
ax.axvline(3, c="#B0532F", ls=":", lw=1.4)
ax.set(xlabel="demand x", ylabel="membership μ", title="Low / Medium / High — a partition of unity (seed 26210)")
ax.legend(frameon=False); ax.grid(alpha=.25); plt.tight_layout(); plt.show()
print(f"at x=3: mu_L={mu_L(3):.4f}  mu_M={mu_M(3):.4f}  mu_H={mu_H(3):.4f}  sum={mu_L(3)+mu_M(3)+mu_H(3):.4f}")

xs = np.linspace(0, 10, 400)
fig, ax = plt.subplots(figsize=(7.8,4.0))
ax.plot(xs, [crisp(x, EXPERT) for x in xs], c="#8A8F8B", lw=2, label="crisp switch at x = 5")
ax.plot(xs, [sugeno(x, EXPERT) for x in xs], c="#C8A24B", lw=2.4, label="fuzzy Sugeno blend")
ax.axvspan(4.9, 5.1, color="#B0532F", alpha=.15, label="boundary ±0.1")
ax.set(xlabel="demand x", ylabel="decision y", title="The blend has no cliff — seed 26210")
ax.legend(frameon=False, fontsize=9); ax.grid(alpha=.25); plt.tight_layout(); plt.show()
print(f"yhat(3) = {sugeno(3,EXPERT):.4f}   yhat(7) = {sugeno(7,EXPERT):.4f}")
print(f"crisp jump across x=5: {abs(crisp(5.1,EXPERT)-crisp(4.9,EXPERT)):.4f}   fuzzy change: {abs(sugeno(5.1,EXPERT)-sugeno(4.9,EXPERT)):.4f}")

fit = fit_lse()
print("fitted consequents (a, b) per rule:", {k:(round(a,4),round(b,4)) for k,(a,b) in fit.items()})
xs = np.linspace(0, 10, 400)
fig, ax = plt.subplots(figsize=(7.8,4.0))
ax.plot(xs, [target(x) for x in xs], c="#17231F", lw=2, ls="--", label="true cost curve")
ax.plot(xs, [sugeno(x, EXPERT) for x in xs], c="#8A8F8B", lw=2, label=f"expert rules (RMSE {rmse(EXPERT):.2f})")
ax.plot(xs, [sugeno(x, fit) for x in xs], c="#C8A24B", lw=2.4, label=f"ANFIS-calibrated (RMSE {rmse(fit):.4f})")
ax.scatter(XS, target(XS), c="#0B3D2E", s=28, zorder=5, label="data")
ax.set(xlabel="x", ylabel="cost", title="Same rule grammar, learned coefficients (seed 26210)")
ax.legend(frameon=False, fontsize=9); ax.grid(alpha=.25); plt.tight_layout(); plt.show()
print(f"yfit(3) = {sugeno(3,fit):.4f}   yfit(7) = {sugeno(7,fit):.4f}")

ref = reference_values()
expected = {"mu_L_3":0.4,"mu_M_3":0.6,"mu_H_3":0.0,"partition_sum_3":1.0,
 "yhat_3":3.52,"yhat_7":7.92,"crisp_jump":1.16,"fuzzy_delta":0.22,
 "robustness_ratio":5.2727,"rmse_expert":11.777,"rmse_fitted":0.0003,
 "yfit_3":14.1499,"yfit_7":16.1497}
for k,v in expected.items():
    assert abs(ref[k]-v)<5e-4, f"MISMATCH {k}"
    print(f"PASS  {k:20s} {ref[k]}")
print("\nAll checkpoints agree — seed 26210.")