* MSC-P-035 saturation and adstock. MIT License.
BEGIN PROGRAM Python3.
import csv, math
required = ["week","media_spend_keur","revenue_keur"]
path = "public/datasets/msc-p035-media-series.csv"
with open(path, newline="", encoding="utf-8") as handle:
    reader = csv.DictReader(handle)
    if reader.fieldnames != required: raise ValueError("exact schema required")
    rows = list(reader)
if len(rows) < 104: raise ValueError("too few weeks for the declared grid search")
if any(set(r) != set(required) or any(r[k] is None for k in required) for r in rows): raise ValueError("exactly three cells required per row")
if [r["week"].strip() for r in rows] != [str(i) for i in range(1, len(rows) + 1)]: raise ValueError("weeks must be numbered 1, 2, ... without gaps")
x, y = [], []
for r in rows:
    a = float(r["media_spend_keur"]); b = float(r["revenue_keur"])
    if not math.isfinite(a) or a <= 0: raise ValueError("media spend must be finite and strictly positive")
    if not math.isfinite(b) or b <= 0: raise ValueError("revenue must be finite and strictly positive")
    x.append(a); y.append(b)
PERIOD, HARM = 52, 2
LAMBDAS = [round(0.1 * s, 1) for s in range(10)]
SHAPES = [0.6, 1.0, 1.4, 1.8, 2.2, 2.6, 3.0]
GAMMAS = [25.0, 35.0, 45.0, 55.0, 65.0, 75.0, 85.0, 95.0]
TRUTH = (0.6, 1.8, 55.0)
TOL, LAM_RANGE_MAX, RATIO_MAX = 0.01, 0.20, 1.50
n = len(y)
base = [[1.0] * n, [float(w) for w in range(1, n + 1)]]
for k in range(1, HARM + 1):
    angles = [2.0 * math.pi * k * w / PERIOD for w in range(1, n + 1)]
    base.append([math.sin(a) for a in angles]); base.append([math.cos(a) for a in angles])
def adstock(lam):
    out, carried = [], 0.0
    for v in x:
        carried = v + lam * carried
        out.append(carried)
    return out
def sat(series, b, g):
    d = g ** b
    return [v ** b / (d + v ** b) for v in series]
def slope(level, b, g):
    d = g ** b
    return b * d * level ** (b - 1.0) / (d + level ** b) ** 2
def solve(matrix, vector):
    m = len(vector)
    aug = [list(r) + [v] for r, v in zip(matrix, vector)]
    for c in range(m):
        p = max(range(c, m), key=lambda i: abs(aug[i][c]))
        if abs(aug[p][c]) < 1e-12: raise ValueError("the declared design matrix is singular for this grid point")
        aug[c], aug[p] = aug[p], aug[c]
        for i in range(c + 1, m):
            f = aug[i][c] / aug[c][c]
            for j in range(c, m + 1): aug[i][j] -= f * aug[c][j]
    out = [0.0] * m
    for c in range(m - 1, -1, -1):
        out[c] = (aug[c][m] - sum(aug[c][j] * out[j] for j in range(c + 1, m))) / aug[c][c]
    return out
def fit(columns):
    rowset = list(zip(*columns)); size = len(columns)
    normal = [[sum(r[i] * r[j] for r in rowset) for j in range(size)] for i in range(size)]
    right = [sum(r[i] * v for r, v in zip(rowset, y)) for i in range(size)]
    beta = solve(normal, right)
    rss = sum((v - sum(a * b for a, b in zip(r, beta))) ** 2 for r, v in zip(rowset, y))
    return beta, rss
carried = {lam: adstock(lam) for lam in LAMBDAS}
grid = []
for lam in LAMBDAS:
    series = carried[lam]; level = sum(series) / n
    for b in SHAPES:
        for g in GAMMAS:
            beta, rss = fit(base + [sat(series, b, g)])
            grid.append({"lam": lam, "b": b, "g": g, "rss": rss, "amp": beta[-1],
                         "marg": beta[-1] * slope(level, b, g)})
best = min(grid, key=lambda e: (e["rss"], e["lam"], e["b"], e["g"]))
region = [e for e in grid if e["rss"] <= best["rss"] * (1.0 + TOL)]
lams = [e["lam"] for e in region]; shapes = [e["b"] for e in region]
gammas = [e["g"] for e in region]; margs = [e["marg"] for e in region]
if min(margs) <= 0: raise ValueError("a grid point in the region implies a non-positive marginal return")
lam_range = max(lams) - min(lams); ratio = max(margs) / min(margs)
lam_flag = "PASS" if lam_range <= LAM_RANGE_MAX else "FAIL"
marg_flag = "PASS" if ratio <= RATIO_MAX else "FAIL"
truth = next(e for e in grid if (e["lam"], e["b"], e["g"]) == TRUTH)
mean_y = sum(y) / n
total = sum((v - mean_y) ** 2 for v in y)
print("weeks=%d" % n)
print("grid_points=%d" % len(grid))
print("best_carryover=%.6f" % best["lam"])
print("best_shape=%.6f" % best["b"])
print("best_half_saturation=%.6f" % best["g"])
print("best_amplitude=%.6f" % best["amp"])
print("r_squared=%.6f" % (1.0 - best["rss"] / total))
print("residual_sd=%.6f" % math.sqrt(best["rss"] / (n - len(base) - 1)))
print("region_tolerance=%.2f" % TOL)
print("region_points=%d" % len(region))
print("carryover_min=%.6f" % min(lams))
print("carryover_max=%.6f" % max(lams))
print("carryover_range=%.6f" % lam_range)
print("carryover_flag=%s" % lam_flag)
print("shape_min=%.6f" % min(shapes))
print("shape_max=%.6f" % max(shapes))
print("half_saturation_min=%.6f" % min(gammas))
print("half_saturation_max=%.6f" % max(gammas))
print("marginal_min=%.6f" % min(margs))
print("marginal_max=%.6f" % max(margs))
print("marginal_ratio=%.6f" % ratio)
print("marginal_flag=%s" % marg_flag)
print("best_marginal_return=%.6f" % best["marg"])
print("true_parameters_inside_region=%s" % ("YES" if truth in region else "NO"))
print("true_parameters_fit_excess_percent=%.6f" % (100.0 * (truth["rss"] / best["rss"] - 1.0)))
print("true_marginal_return=%.6f" % truth["marg"])
print("verdict=%s" % ("RESPONSE_CURVE_IDENTIFIED_FOR_REALLOCATION" if lam_flag == "PASS" and marg_flag == "PASS" else "DIAGNOSTIC_BLOCKS_RESPONSE_CURVE"))
END PROGRAM.