{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "h1-ext-00",
   "metadata": {},
   "source": [
    "# Hands-On 1: optional exercises\n",
    "Run these after the six-stage session. These numerical exercises complement the\n",
    "mathematical appendices and the main session. Each exercise asks one\n",
    "additional question. The setup recomputes the full calibration and the\n",
    "short-recalibration comparison, so this notebook runs independently."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h1-ext-01",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import session_plotting as plots\n",
    "from session_plotting import tr\n",
    "plots.configure(\"en\")\n",
    "np.set_printoptions(precision=3, suppress=True)\n",
    "\n",
    "cal = np.loadtxt(\"data/calibration.csv\", delimiter=\",\", skiprows=1)\n",
    "T_cal, Y = cal[:, 1], cal[:, 2:]\n",
    "N, C = Y.shape\n",
    "X = np.column_stack([np.ones(N), T_cal])\n",
    "Theta = np.linalg.lstsq(X, Y, rcond=None)[0]\n",
    "a, b = Theta\n",
    "R = Y - X @ Theta\n",
    "E = R / b\n",
    "s = np.sqrt((R**2).sum(axis=0) / (N - 2))\n",
    "S = E.T @ E / (N - 2)\n",
    "quick = np.loadtxt(\"data/recal_quick.csv\", delimiter=\",\", skiprows=1)\n",
    "T_q, y_q = quick.T\n",
    "Xq = np.column_stack([np.ones(len(T_q)), T_q])\n",
    "theta_ls = np.linalg.lstsq(Xq, y_q, rcond=None)[0]\n",
    "sigma = s[2]\n",
    "cov_ls = sigma**2 * np.linalg.inv(Xq.T @ Xq)\n",
    "others = np.arange(C) != 2\n",
    "mu0 = Theta[:, others].mean(axis=1)\n",
    "sd0 = Theta[:, others].std(axis=1, ddof=1)\n",
    "Lambda0 = np.diag(1 / sd0**2)\n",
    "P = Xq.T @ Xq / sigma**2 + Lambda0\n",
    "theta_map = np.linalg.solve(P, Xq.T @ y_q / sigma**2 + Lambda0 @ mu0)\n",
    "S_post = np.linalg.inv(P)\n",
    "full = np.loadtxt(\"data/recal_full.csv\", delimiter=\",\", skiprows=1)\n",
    "Xf = np.column_stack([np.ones(len(full)), full[:, 0]])\n",
    "theta_full = np.linalg.lstsq(Xf, full[:, 1], rcond=None)[0]\n",
    "rng = np.random.default_rng(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-ext-02",
   "metadata": {},
   "source": [
    "## Exercise 1. Why does one matrix solve equal ten separate fits?\n",
    "Construct the block-diagonal problem, then compare its answer with the\n",
    "matrix-right-hand-side solve. The shared-design WLS proof is in appendix A.\n",
    "Also inspect coefficient uncertainty: why does the temperature spread matter?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h1-ext-03",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_big = np.kron(np.eye(C), X)\n",
    "theta_big = np.linalg.lstsq(X_big, Y.T.reshape(-1), rcond=None)[0]\n",
    "print(\"same fit:\", np.allclose(theta_big.reshape(C, 2).T, Theta))\n",
    "XtX_inv = np.linalg.inv(X.T @ X)\n",
    "sd_theta = np.sqrt(np.outer(np.diag(XtX_inv), s**2))\n",
    "print(\"sd [a; b]:\", sd_theta)\n",
    "print(\"S_TT:\", ((T_cal - T_cal.mean())**2).sum())\n",
    "\n",
    "z = (R / s).ravel()\n",
    "grid = np.linspace(-4, 4, 300)\n",
    "fig, ax = plt.subplots(figsize=(7, 3))\n",
    "ax.hist(z, bins=50, density=True, alpha=0.5)\n",
    "ax.plot(grid, stats.norm.pdf(grid))\n",
    "ax.set(xlabel=\"r / s\", ylabel=tr(\"density\", \"πυκνότητα\"))\n",
    "plt.tight_layout(); plt.show()\n",
    "print(\"P(|r/s| > 2):\", np.mean(np.abs(z) > 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-ext-04",
   "metadata": {},
   "source": [
    "## Exercise 2. How can the local fit look convincing while the gain is unreliable?\n",
    "Draw plausible lines, simulate repeated checks, and inspect the SVD weak\n",
    "direction. Then draw the ridge filter factors. Why does a wide range of\n",
    "regularisation strengths give almost the same coefficients? See appendix B."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h1-ext-05",
   "metadata": {},
   "outputs": [],
   "source": [
    "one_bath = T_q == 25.0\n",
    "theta_one = np.linalg.lstsq(Xq[one_bath], y_q[one_bath], rcond=None)[0]\n",
    "print(\"one-bath rank / minimum-norm fit:\", np.linalg.matrix_rank(Xq[one_bath]), theta_one)\n",
    "fan = rng.multivariate_normal(theta_ls, cov_ls, size=40)\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10, 3.4))\n",
    "for ax, limits in zip(axes, [(24.9, 25.3), (0, 45)]):\n",
    "    t = np.array(limits)\n",
    "    for theta in fan:\n",
    "        ax.plot(t, theta[0] + theta[1] * t, color=\"red\", alpha=0.2, lw=0.6)\n",
    "    ax.plot(t, theta_full[0] + theta_full[1] * t, \"k--\")\n",
    "    ax.plot(T_q, y_q, \"k.\")\n",
    "    ax.set(xlabel=\"T [°C]\", ylabel=\"counts\", xlim=limits)\n",
    "plt.tight_layout(); plt.show()\n",
    "\n",
    "Y_sim = (Xq @ Theta[:, 2])[:, None] + sigma * rng.standard_normal((len(T_q), 5000))\n",
    "b_sim = np.linalg.lstsq(Xq, Y_sim, rcond=None)[0][1]\n",
    "print(\"gain: mean / sd / P(b<0):\", b_sim.mean(), b_sim.std(), np.mean(b_sim < 0))\n",
    "_, sv, Vt = np.linalg.svd(Xq, full_matrices=False)\n",
    "print(\"singular values / weak direction:\", sv, Vt[1])\n",
    "\n",
    "def ridge(lam):\n",
    "    return np.linalg.solve(Xq.T @ Xq + lam * np.eye(2), Xq.T @ y_q)\n",
    "\n",
    "lams = np.logspace(-7, 5, 120)\n",
    "path = np.array([ridge(lam) for lam in lams])\n",
    "fig, axes = plt.subplots(1, 3, figsize=(11, 3.2))\n",
    "for i in range(2):\n",
    "    axes[0].semilogx(lams, sv[i]**2 / (sv[i]**2 + lams), label=f\"f{i+1}\")\n",
    "    axes[i+1].semilogx(lams, path[:, i])\n",
    "    axes[i+1].axhline(theta_full[i], color=\"black\", ls=\"--\")\n",
    "axes[0].legend(); axes[1].set_ylabel(\"a [counts]\"); axes[2].set_ylabel(\"b [counts/°C]\")\n",
    "for ax in axes:\n",
    "    ax.set_xlabel(\"λ\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-ext-06",
   "metadata": {},
   "source": [
    "## Exercise 3. What do the prior and likelihood each contribute?\n",
    "The contour plot evaluates the log likelihood, prior and posterior directly.\n",
    "The uncertainty curves and repeated-check table distinguish posterior\n",
    "uncertainty from sampling variability. Change the prior mean and scales;\n",
    "a small posterior uncertainty need not imply a small bias. Appendix B gives\n",
    "the derivations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h1-ext-07",
   "metadata": {},
   "outputs": [],
   "source": [
    "A, B = np.meshgrid(np.linspace(-700, 500, 500), np.linspace(32, 72, 500))\n",
    "residual = y_q[:, None, None] - A[None] - B[None] * T_q[:, None, None]\n",
    "loglik = -0.5 * (residual**2).sum(axis=0) / sigma**2\n",
    "logprior = -0.5 * (((A-mu0[0])/sd0[0])**2 + ((B-mu0[1])/sd0[1])**2)\n",
    "logpost = loglik + logprior\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "for values, color, label in [(loglik, \"red\", \"likelihood\"),\n",
    "                              (logprior, \"blue\", \"prior\"), (logpost, \"purple\", \"posterior\")]:\n",
    "    ax.contour(A, B, values - values.max(), levels=[-4.5, -2, -0.5], colors=color)\n",
    "    ax.plot([], [], color=color, label=label)\n",
    "ax.plot(*theta_map, \"ko\")\n",
    "ax.set(xlabel=\"a [counts]\", ylabel=\"b [counts/°C]\")\n",
    "ax.legend(); plt.tight_layout(); plt.show()\n",
    "\n",
    "T_line = np.linspace(0, 50, 201)\n",
    "XT = np.column_stack([np.ones_like(T_line), T_line])\n",
    "curves = {\"LS\": cov_ls, \"prior\": np.diag(sd0**2), \"MAP\": S_post,\n",
    "          \"full\": sigma**2 * np.linalg.inv(Xf.T @ Xf), \"rack\": sigma**2 * XtX_inv}\n",
    "fig, ax = plt.subplots(figsize=(8, 3.4))\n",
    "for label, covariance in curves.items():\n",
    "    sd_line = np.sqrt(np.einsum(\"ti,ij,tj->t\", XT, covariance, XT)) / theta_full[1]\n",
    "    ax.semilogy(T_line, sd_line, label=label)\n",
    "ax.set(xlabel=\"T [°C]\", ylabel=\"sd [°C]\"); ax.legend()\n",
    "plt.tight_layout(); plt.show()\n",
    "\n",
    "fits = {\"LS\": np.linalg.lstsq(Xq, Y_sim, rcond=None)[0],\n",
    "        \"ridge\": np.linalg.solve(Xq.T @ Xq + np.eye(2), Xq.T @ Y_sim),\n",
    "        \"MAP\": np.linalg.solve(P, Xq.T @ Y_sim / sigma**2 + (Lambda0 @ mu0)[:, None])}\n",
    "for label, fitted in fits.items():\n",
    "    err = (np.array([1,40]) @ (fitted - Theta[:, 2, None])) / b[2]\n",
    "    print(label, \"bias / sd / RMSE:\", err.mean(), err.std(), np.sqrt((err**2).mean()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-ext-08",
   "metadata": {},
   "source": [
    "## Exercise 4. Why does a two-dimensional “2-sigma ellipse” contain only 86%?\n",
    "Simulate χ² at several dimensions, then draw ellipses in the original\n",
    "coordinates. Finally compare Cholesky and symmetric whitening: different\n",
    "coordinates, identical squared distances. See appendices D and E."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h1-ext-09",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(11, 3.1))\n",
    "for ax, df in zip(axes, [1, 2, 10]):\n",
    "    draws = (rng.standard_normal((20000, df))**2).sum(axis=1)\n",
    "    ax.hist(draws, bins=70, density=True, alpha=0.5)\n",
    "    grid = np.linspace(0.01, stats.chi2.ppf(0.999, df), 300)\n",
    "    ax.plot(grid, stats.chi2.pdf(grid, df))\n",
    "    ax.axvline(stats.chi2.ppf(0.99, df), color=\"red\", ls=\"--\")\n",
    "    ax.set_title(f\"χ²({df})\")\n",
    "plt.tight_layout(); plt.show()\n",
    "\n",
    "S2 = S[:2, :2]\n",
    "values, vectors = np.linalg.eigh(S2)\n",
    "t = np.linspace(0, 2*np.pi, 200)\n",
    "fig, ax = plt.subplots(figsize=(4.5, 4))\n",
    "ax.plot(E[:, 0], E[:, 1], \".\", alpha=0.4)\n",
    "for c in [1, 2, 3]:\n",
    "    ellipse = vectors @ (c*np.sqrt(values)[:, None] * np.array([np.cos(t), np.sin(t)]))\n",
    "    ax.plot(*ellipse, label=f\"c={c}, P={stats.chi2.cdf(c*c, 2):.0%}\")\n",
    "ax.set(xlabel=\"e1 [°C]\", ylabel=\"e2 [°C]\", aspect=\"equal\")\n",
    "ax.legend(); plt.tight_layout(); plt.show()\n",
    "L = np.linalg.cholesky(S)\n",
    "values, vectors = np.linalg.eigh(S)\n",
    "W = (vectors * values**-0.5) @ vectors.T\n",
    "q_chol = np.linalg.solve(L, E.T).T\n",
    "q_sym = E @ W\n",
    "print(\"same distances:\", np.allclose((q_chol**2).sum(1), (q_sym**2).sum(1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-ext-10",
   "metadata": {},
   "source": [
    "## Exercise 5. Why can correlation make a small single-probe fault easier to detect?\n",
    "This is an **answer-key analysis** of the validation data: probe 7 has a\n",
    "+0.30 °C offset from snapshot 200. It uses reference temperatures to compare\n",
    "single-probe z-scores with the χ²(C) distance. The main session instead diagnoses\n",
    "unknown temperatures with C−1 degrees of freedom. Appendix E derives the\n",
    "conditional variance and noncentral χ² power curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h1-ext-11",
   "metadata": {},
   "outputs": [],
   "source": [
    "ops = np.loadtxt(\"data/operations.csv\", delimiter=\",\", skiprows=1)\n",
    "errors = (ops[:, 2:] - a) / b - ops[:, 1, None]\n",
    "q = np.linalg.solve(L, errors.T).T\n",
    "d2 = (q**2).sum(axis=1)\n",
    "threshold = stats.chi2.ppf(0.99, C)\n",
    "precision = np.linalg.inv(S)\n",
    "j = 6\n",
    "sd_total = np.sqrt(S[j, j])\n",
    "sd_conditional = 1 / np.sqrt(precision[j, j])\n",
    "print(\"total / conditional sd:\", sd_total, sd_conditional)\n",
    "delta = np.linspace(0, 0.8, 161)\n",
    "zcrit = stats.norm.ppf(0.995)\n",
    "power_z = stats.norm.sf(zcrit - delta/sd_total) + stats.norm.cdf(-zcrit - delta/sd_total)\n",
    "power_m = stats.ncx2.sf(threshold, C, delta**2 * precision[j, j])\n",
    "fig, ax = plt.subplots(figsize=(8, 3.4))\n",
    "ax.plot(delta, power_z, label=\"z\")\n",
    "ax.plot(delta, power_m, label=\"Mahalanobis\")\n",
    "ax.set(xlabel=\"δ [°C]\", ylabel=tr(\"alarm probability\", \"πιθανότητα συναγερμού\"))\n",
    "ax.legend(); plt.tight_layout(); plt.show()\n",
    "print(\"observed before / after:\", (d2[:200] > threshold).mean(), (d2[200:] > threshold).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-ext-12",
   "metadata": {},
   "source": [
    "## Further exercises\n",
    "- Estimate the covariance using only 25 snapshots from each calibration bath.\n",
    "  How does the held-out false-alarm rate change? Appendix C explains inverse-covariance bias.\n",
    "- Regenerate a clean operations run before replacing the step fault with a\n",
    "  gradual drift or extra noise; avoid adding a new fault on top of the existing one.\n",
    "- Try two faulty probes. Explain why the main session's single-candidate rule may\n",
    "  retain all probes and report an unresolved cause.\n",
    "- Use a separate new simulated run to assess any modified diagnostic rule."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}