{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  },
  "colab": {
   "name": "week-07-h1-sequences-lstm.ipynb",
   "provenance": []
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# W07 H1: Sequential Data and LSTM\n\nMachine Learning in Engineering Science | NCKU + LAiMM | Fall 2026\n\n**Teaching date:** Monday 19 October 2026. **Guided notebook time:** 25 minutes within a 60-minute module.\n\n### What you will be able to explain\n- Define a sequence, a forecast origin, and a prediction horizon.\n- Trace how the four LSTM gates update cell and hidden state.\n- Compare a trained LSTM with persistence and ridge baselines without future leakage.\n\n### Start here, even if this is your first notebook\n1. In Colab choose **Runtime > Change runtime type > CPU**. In Jupyter choose a Python 3 kernel.\n2. Run cells from top to bottom with **Shift+Enter**. Read the Markdown before each code cell.\n3. The setup installs only missing free packages. A network connection is needed for first installation; the data itself is generated locally.\n4. If a variable is undefined, restart the kernel and run all earlier cells. Do not skip a setup cell.\n5. Save a personal copy before changing parameters. Every core cell is complete: there are no hidden files, tokens, paid services, or unfinished code blocks.\n\n**Data:** synthetic teaching data, not engineering measurements. Record your actual output; a model is not required to outperform a baseline.\n**Practice policy:** in-class exercise and optional Homework extension, not a new weekly graded assignment.\n",
   "id": "cell-000"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## Mechanism reference: the complete LSTM cell\n\nLet $u_t=[h_{t-1};x_t]$ concatenate the previous hidden vector and current input. A standard LSTM computes\n\n$$f_t=\\sigma(W_fu_t+b_f),\\quad i_t=\\sigma(W_iu_t+b_i),$$\n$$g_t=\\tanh(W_gu_t+b_g),\\quad o_t=\\sigma(W_ou_t+b_o),$$\n$$c_t=f_t\\odot c_{t-1}+i_t\\odot g_t,\\quad h_t=o_t\\odot\\tanh(c_t).$$\n\nSigmoid gives coefficients between 0 and 1; $\\odot$ multiplies matching coordinates. Forget $f$ retains old memory, input $i$ scales a signed write $g$, and output $o$ controls exposed state. Weights are shared across time. Each gate has one coefficient per hidden dimension.\n\nThe additive path can improve gradient flow but does not guarantee infinite memory. With gates held fixed, the direct retained contribution includes a product of forget coefficients; full gradients also include gate dependencies. Gradient clipping limits explosions but does not restore vanished gradients.\n\n",
   "id": "cell-001"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 1. Record the environment\n\nRun setup once. The random seed fixes the synthetic data; recording package versions makes differences between runs traceable.\n",
   "id": "cell-002"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "metadata": {},
   "source": "import sys, subprocess, importlib.util, platform, json, random, hashlib\nfrom pathlib import Path\nrequired = {\"numpy\": \"numpy\", \"matplotlib\": \"matplotlib\", \"sklearn\": \"scikit-learn\"}\nmissing = [package for module, package in required.items() if importlib.util.find_spec(module) is None]\nif missing:\n    subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", *missing])\nimport numpy as np\nimport matplotlib.pyplot as plt\nSEED = 2026\nrandom.seed(SEED)\nnp.random.seed(SEED)\nrng = np.random.default_rng(SEED)\nARTIFACTS = Path(\"artifacts\")\nARTIFACTS.mkdir(exist_ok=True)\nprint(\"Python:\", platform.python_version(), \"NumPy:\", np.__version__, \"Seed:\", SEED)\n",
   "id": "cell-003"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**Checkpoint.** You should see the Python version, NumPy version, and seed 2026. If installation fails, check internet access and rerun this cell.\n",
   "id": "cell-004"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 2. Enable a small CPU model\n\nNo GPU is required. We use one CPU thread and deterministic operations so this small example is inexpensive to rerun.\n",
   "id": "cell-005"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "metadata": {},
   "source": "if importlib.util.find_spec(\"torch\") is None:\n    subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", \"torch\"])\nimport torch\nfrom torch import nn\nimport copy, time\ntorch.set_num_threads(1)\ntorch.manual_seed(SEED)\ntorch.use_deterministic_algorithms(True)\ndevice = torch.device(\"cpu\")\nprint(\"PyTorch:\", torch.__version__, \"Device:\", device)\n",
   "id": "cell-006"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**Checkpoint.** The printed device must be cpu. Small numerical differences can still occur across package versions or hardware.\n",
   "id": "cell-007"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 3. Define the physical question before the model\n\nImagine a sensor sampled at equal time intervals. We predict the next sample from the preceding 24 samples. The signal below is synthetic: a slow trend, two oscillations, and noise. Index t represents a sample, not a calibrated second.\n",
   "id": "cell-008"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "metadata": {},
   "source": "N, WINDOW = 1600, 24\nt = np.arange(N)\nsignal = (0.0004*t + np.sin(2*np.pi*t/40) + 0.35*np.sin(2*np.pi*t/13)\n          + rng.normal(0, 0.12, N)).astype(np.float32)\nTRAIN_END, VAL_END = 960, 1280\nplt.figure(figsize=(11,3))\nplt.plot(t, signal, lw=1, label=\"synthetic sensor\")\nplt.axvline(TRAIN_END,color=\"tab:orange\",label=\"train/validation boundary\")\nplt.axvline(VAL_END,color=\"tab:red\",label=\"validation/test boundary\")\nplt.xlabel(\"Sample index\"); plt.ylabel(\"Synthetic sensor units\"); plt.legend(); plt.show()\nprint(\"Dataset SHA256:\", hashlib.sha256(signal.tobytes()).hexdigest())\n",
   "id": "cell-009"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**Checkpoint.** Training occurs on earlier samples. Validation and test occur later. Do not shuffle the raw timeline before splitting.\n",
   "id": "cell-010"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 4. Build windows using only available history\n\nA row is a complete history window, not a single sensor value. With length L=24 and target index t, input is x[t-24:t] and the target is x[t]. Thus every input index is strictly smaller than its target index. A test window may use previously observed validation/test sensor history: this is **rolling one-step forecasting**, not a forecast of the entire future from one fixed origin.\n",
   "id": "cell-011"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "metadata": {},
   "source": "target_index = np.arange(WINDOW, N)\nX_raw = np.stack([signal[i-WINDOW:i] for i in target_index])\ny_raw = signal[target_index]\ntrain = target_index < TRAIN_END\nvalid = (target_index >= TRAIN_END) & (target_index < VAL_END)\ntest = target_index >= VAL_END\nmu, sigma = signal[:TRAIN_END].mean(), signal[:TRAIN_END].std()\nX = ((X_raw-mu)/sigma).astype(np.float32)[..., None]\ny = ((y_raw-mu)/sigma).astype(np.float32)\nassert X.shape == (N-WINDOW, WINDOW, 1)\nassert target_index[train].max() < target_index[valid].min()\nassert target_index[valid].max() < target_index[test].min()\nprint(\"X shape (batch,time,features):\", X.shape)\nprint(\"Counts:\", train.sum(), valid.sum(), test.sum())\nprint(\"First window uses indices\", target_index[0]-WINDOW, \"to\", target_index[0]-1,\n      \"and predicts\", target_index[0])\n",
   "id": "cell-012"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**Checkpoint.** The shape is (1576, 24, 1), with 936 training, 320 validation, and 320 test targets. Mean and standard deviation use only raw training samples.\n",
   "id": "cell-013"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 5. Establish two baselines\n\nPersistence predicts that the next value equals the last observation. Ridge regression can exploit the full window without a recurrent state. Choose ridge regularization on validation only, then evaluate the selected baseline on test once.\n",
   "id": "cell-014"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "metadata": {},
   "source": "from sklearn.linear_model import Ridge\nfrom sklearn.metrics import mean_absolute_error, mean_squared_error\nridge_candidates = []\nfor alpha in [0.01, 0.1, 1.0, 10.0]:\n    candidate = Ridge(alpha=alpha).fit(X[train,:,0], y[train])\n    val_mae = mean_absolute_error(y[valid], candidate.predict(X[valid,:,0]))\n    ridge_candidates.append((val_mae, alpha, candidate))\n_, best_alpha, ridge = min(ridge_candidates, key=lambda item:item[0])\npredictions = {\"Persistence\": X[test,-1,0]*sigma+mu,\n               \"Ridge\": ridge.predict(X[test,:,0])*sigma+mu}\nprint(\"Validation-selected ridge alpha:\", best_alpha)\n",
   "id": "cell-015"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**Checkpoint.** Persistence has no fitted parameters. Both baselines predict exactly the same target indices as the LSTM.\n",
   "id": "cell-016"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 6. Compute one LSTM update by hand\n\nThe forget gate f retains old cell memory; input gate i controls writing; candidate g supplies new content; output gate o exposes the updated cell. All products are elementwise. Here the gate values are **chosen**, not trained. For c_previous=0.8, f=0.9, i=0.2, g=-0.5, and o=0.7, compute c_new=f*c_previous+i*g before running.\n",
   "id": "cell-017"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "metadata": {},
   "source": "c_previous, f, i, g, o = 0.8, 0.9, 0.2, -0.5, 0.7\nc_new = f*c_previous + i*g\nh_new = o*np.tanh(c_new)\nprint(\"New cell:\", c_new, \"New hidden:\", h_new)\nfig, ax = plt.subplots(figsize=(8,3))\nax.bar([\"retained old memory\", \"new written content\", \"updated cell\", \"exposed hidden\"],\n       [f*c_previous, i*g, c_new, h_new], color=[\"#1e8798\",\"#e4002b\",\"#12233f\",\"#f3b23a\"])\nax.axhline(0,color=\"gray\"); ax.set_ylabel(\"Illustrative scalar state\"); plt.show()\nassert np.isclose(c_new, 0.62)\n",
   "id": "cell-018"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**Checkpoint.** The new cell is 0.62 and hidden state is about 0.386. The cell and hidden state are different quantities; an LSTM does not simply copy the last input.\n",
   "id": "cell-019"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 7. Implement the prediction model\n\nWith batch_first=True, input and output are (B,T,F) and (B,T,H). Hidden state has shape (layers,B,H), regardless of batch_first. A linear head turns the final hidden output into one forecast. Model weights are shared across all 24 time steps.\n",
   "id": "cell-020"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "metadata": {},
   "source": "class SensorLSTM(nn.Module):\n    def __init__(self, hidden=16):\n        super().__init__()\n        self.rnn = nn.LSTM(input_size=1, hidden_size=hidden, batch_first=True)\n        self.head = nn.Linear(hidden, 1)\n    def forward(self, x):\n        states, (h_n, c_n) = self.rnn(x)\n        return self.head(states[:, -1, :]).squeeze(-1)\n\ntorch.manual_seed(SEED)\nmodel = SensorLSTM().to(device)\nX_t = torch.from_numpy(X)\ny_t = torch.from_numpy(y)\nprint(\"Parameters:\", sum(p.numel() for p in model.parameters()))\nprint(\"Output shape:\", tuple(model(X_t[:4]).shape))\nassert model(X_t[:4]).shape == (4,)\n",
   "id": "cell-021"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**Checkpoint.** There must be one prediction per window. A bidirectional model is unnecessary here; every input window already contains only known past values.\n",
   "id": "cell-022"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 8. Train on training data; select on validation\n\nEach epoch performs gradient updates on training batches. Validation is evaluated with model.eval() and no gradients. We retain the state with the lowest validation loss. The test set is not consulted during model selection.\n",
   "id": "cell-023"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "metadata": {},
   "source": "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)\nloss_fn = nn.MSELoss()\ntrain_ids = np.flatnonzero(train)\nhistory, best_val, best_state = [], float(\"inf\"), None\nstart = time.perf_counter()\nfor epoch in range(24):\n    model.train()\n    order = rng.permutation(train_ids)\n    epoch_losses = []\n    for batch in np.array_split(order, 12):\n        optimizer.zero_grad()\n        loss = loss_fn(model(X_t[batch]), y_t[batch])\n        loss.backward()\n        nn.utils.clip_grad_norm_(model.parameters(), 1.0)\n        optimizer.step()\n        epoch_losses.append(loss.item())\n    model.eval()\n    with torch.no_grad():\n        val_loss = loss_fn(model(X_t[valid]), y_t[valid]).item()\n    history.append((float(np.mean(epoch_losses)), val_loss))\n    if val_loss < best_val:\n        best_val, best_state = val_loss, copy.deepcopy(model.state_dict())\nmodel.load_state_dict(best_state)\nprint(\"Training seconds:\", round(time.perf_counter()-start,2), \"Best validation MSE:\", best_val)\nplt.plot(np.array(history), label=None)\nplt.legend([\"training MSE\", \"validation MSE\"]); plt.xlabel(\"Epoch\"); plt.ylabel(\"Standardized MSE\"); plt.show()\n",
   "id": "cell-024"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**Checkpoint.** Your actual loss curve is the evidence. If validation worsens while training improves, discuss overfitting rather than increasing epochs blindly.\n",
   "id": "cell-025"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 9. Compare held-out errors and inspect residuals\n\nReport errors in original synthetic sensor units. MAE is the average absolute error; RMSE penalizes large residuals more strongly. A lower score on this toy signal is not evidence of reliability on a real sensor.\n",
   "id": "cell-026"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "metadata": {},
   "source": "model.eval()\nwith torch.no_grad():\n    predictions[\"LSTM\"] = model(X_t[test]).numpy()*sigma+mu\nmetrics = {}\nfor name, prediction in predictions.items():\n    mae = mean_absolute_error(y_raw[test], prediction)\n    rmse = np.sqrt(mean_squared_error(y_raw[test], prediction))\n    metrics[name] = {\"MAE\": float(mae), \"RMSE\": float(rmse)}\n    print(f\"{name:12s} MAE={mae:.4f} RMSE={rmse:.4f}\")\nfig, axes = plt.subplots(2,1,figsize=(10,5),sharex=True)\naxes[0].plot(target_index[test], y_raw[test], label=\"observed\",color=\"black\")\nfor name, prediction in predictions.items(): axes[0].plot(target_index[test],prediction,label=name,alpha=.7)\naxes[0].legend(); axes[0].set_ylabel(\"Sensor units\")\naxes[1].plot(target_index[test], predictions[\"LSTM\"]-y_raw[test])\naxes[1].axhline(0,color=\"black\"); axes[1].set_ylabel(\"LSTM residual\"); axes[1].set_xlabel(\"Target sample index\")\nplt.tight_layout(); plt.show()\n",
   "id": "cell-027"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**Checkpoint.** Explain the comparison even if persistence or ridge wins. Inspect whether errors cluster at peaks; one aggregate score hides this structure.\n",
   "id": "cell-028"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 10. Inspect real gates in the trained model\n\nPyTorch stores gate chunks in i,f,g,o order. Reconstruct one cell from learned weights and compare it to nn.LSTM. This connects the equation to an actual implementation, not an illustrative diagram.\n",
   "id": "cell-029"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "metadata": {},
   "source": "x_first = X_t[test][:1,:1,:]\nH = model.rnn.hidden_size\nh0 = torch.zeros(1,H); c0 = torch.zeros(1,H)\nwith torch.no_grad():\n    pre = (x_first[:,0,:] @ model.rnn.weight_ih_l0.T + model.rnn.bias_ih_l0\n           + h0 @ model.rnn.weight_hh_l0.T + model.rnn.bias_hh_l0)\n    pi, pf, pg, po = pre.chunk(4,dim=-1)\n    gate_i, gate_f, candidate, gate_o = pi.sigmoid(), pf.sigmoid(), pg.tanh(), po.sigmoid()\n    cell = gate_f*c0 + gate_i*candidate\n    hidden = gate_o*cell.tanh()\n    reference, (_, reference_cell) = model.rnn(x_first)\n    assert torch.allclose(hidden, reference[:,0,:], atol=1e-6)\n    assert torch.allclose(cell, reference_cell[0], atol=1e-6)\nprint(\"Manual learned-cell update matches PyTorch.\")\nplt.imshow(torch.stack([gate_i[0],gate_f[0],candidate[0],gate_o[0]]).numpy(),aspect=\"auto\",cmap=\"coolwarm\",vmin=-1,vmax=1)\nplt.yticks(range(4),[\"input i\",\"forget f\",\"candidate g\",\"output o\"])\nplt.xlabel(\"Hidden dimension\"); plt.colorbar(); plt.show()\n",
   "id": "cell-030"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**Checkpoint.** The gates are vectors with one value per hidden dimension, not a single universal decision. Their values do not prove causal importance of input samples.\n",
   "id": "cell-031"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 11. Save a rerunnable evidence record\n\nSave data, boundaries, normalization, versions, model state, and measured results. These files stay in your notebook runtime until downloaded. Colab files disappear after the session ends.\n",
   "id": "cell-032"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "metadata": {},
   "source": "record = {\"seed\":SEED,\"window\":WINDOW,\"train_end\":TRAIN_END,\"validation_end\":VAL_END,\n          \"forecast_protocol\":\"rolling one-step, observed history available\",\n          \"data_sha256\":hashlib.sha256(signal.tobytes()).hexdigest(),\n          \"numpy\":np.__version__,\"torch\":torch.__version__,\"python\":platform.python_version(),\n          \"normalization\":{\"mean\":float(mu),\"std\":float(sigma)},\"metrics\":metrics}\n(ARTIFACTS/\"w07-h1-evidence.json\").write_text(json.dumps(record,indent=2))\nnp.savez(ARTIFACTS/\"w07-h1-data.npz\",signal=signal,target_index=target_index)\ntorch.save(model.state_dict(),ARTIFACTS/\"w07-h1-lstm.pt\")\nprint(json.dumps(record,indent=2))\n",
   "id": "cell-033"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**Checkpoint.** Download the artifacts folder or zip it. A classmate should be able to rebuild the dataset and obtain a comparable baseline from the recorded configuration.\n",
   "id": "cell-034"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 12. Independent practice and explain-back\n1. Predict how WINDOW=12 versus 48 changes the information available. Rerun from Step 3 with one change only; keep the original result record.\n2. Explain why fitting normalization on all 1600 samples leaks future information even though no labels are directly passed to a model.\n3. What does the forget gate do when it is 0, 0.5, or 1? Does f near 1 guarantee unlimited memory? Explain the additive path and remaining optimization limits.\n4. Why might a linear model beat an LSTM on this signal? State model capacity, data size, and signal structure.\n\n**Stop when:** you have a held-out comparison table, a residual plot, an evidence JSON, and a two-sentence limitation.\n\n**Bridge to H2:** an encoder that passes only its final state to a decoder compresses the entire history. Can we retrieve relevant earlier states directly for each output?\n\n### Sources\n- [Hochreiter & Schmidhuber (1997), Long Short-Term Memory](https://doi.org/10.1162/neco.1997.9.8.1735)\n- [Dive into Deep Learning: LSTM](https://d2l.ai/chapter_recurrent-modern/lstm.html)\n- [PyTorch sequence models tutorial](https://docs.pytorch.org/tutorials/beginner/nlp/sequence_models_tutorial.html)\n",
   "id": "cell-035"
  }
 ]
}