{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Response surfaces and steepest ascent in two factors\n",
    "\n",
    "**Source worksheets:** [yint.org/w10](https://yint.org/w10) and [yint.org/w11](https://yint.org/w11) - weeks 10 and 11 of the applied DoE short course.\n",
    "\n",
    "Module 6 finished the 1-D version of optimization: pick a factor, take\n",
    "small steps, fit a linear model, swap to a quadratic when the surface\n",
    "curves, predict the peak, confirm.  Module 7 generalizes that loop to\n",
    "**two factors**, using the classical response-surface trio of\n",
    "*factorial points*, a *center point* to detect curvature, and *axial\n",
    "points* to fit a quadratic.\n",
    "\n",
    "This module walks the loop end-to-end on a 2-factor process whose\n",
    "response surface is quadratic in both factors plus an interaction."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. tip::\n",
    "\n",
    "   Two factors is the largest case you can *see* on a flat page.\n",
    "   Use the contour plot as the main output of every step: it\n",
    "   shows the *direction* of steepest ascent at a glance, and where\n",
    "   the next experiment should go."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q1 from w11 - picking sensible coded ranges\n",
    "\n",
    "The ecommerce team has narrowed eight factors down to two:\n",
    "\n",
    "- **A = items in the menu/navigation bar** (1 to 10)\n",
    "- **B = products shown per web page** (1 to 50)\n",
    "\n",
    "The worksheet asks: what are reasonable starting values for the ``-1``\n",
    "and ``+1`` levels?\n",
    "\n",
    "The principle: pick a *narrow* initial range that you can change\n",
    "safely on the live site, centered near your current operating point.\n",
    "You can always widen later.  A *wide* range invites the linear model\n",
    "to badly approximate a curving surface."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. admonition:: Solution\n",
    "\n",
    "   Sensible starting choices:\n",
    "\n",
    "   - For A (menu items): center at the current site's value, say 5,\n",
    "     and use ``-1 = 4`` and ``+1 = 6``.  A 20% perturbation is enough\n",
    "     to learn the slope without surprising users.\n",
    "   - For B (products per page): center at 12 (a common pagination\n",
    "     length), with ``-1 = 8`` and ``+1 = 16``.\n",
    "\n",
    "   Run the 2x2 plus a center point.  If the center is consistent with\n",
    "   the linear fit, take a step in the direction of steepest ascent.\n",
    "   If the center disagrees with the fit, the surface is already\n",
    "   curving and you need axial points (CCD) right away."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q2 - a 2-D optimization, end to end\n",
    "\n",
    "Below is a concrete 2-factor system.  Think of two coded process\n",
    "inputs ``x1`` and ``x2`` (temperature and reaction time, for example)\n",
    "driving some yield-like response ``y``.  The hidden surface is\n",
    "\n",
    "```\n",
    "y(x1, x2) = 70 + 8*x1 + 4*x2\n",
    "            - 3*x1**2 - 2*x2**2\n",
    "            - 2*x1*x2\n",
    "            + noise(0, 0.5)\n",
    "```\n",
    "\n",
    "The true peak sits at ``(x1, x2) = (1.22, 0.33)`` with response\n",
    "``y = 75.6`` - but the optimizer does not know that.  Only the function\n",
    "call ``observe(x1, x2)`` is allowed, and each call returns a *noisy*\n",
    "measurement."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import plotly.graph_objects as go\n",
    "import plotly.io as pio\n",
    "\n",
    "from process_improve.experiments import c, gather, lm, predict\n",
    "\n",
    "pio.renderers.default = \"notebook_connected\"\n",
    "\n",
    "rng = np.random.default_rng(seed=7)\n",
    "\n",
    "\n",
    "def observe(x1: float, x2: float) -> float:\n",
    "    \"\"\"Return a noisy measurement of the hidden 2-D surface.\"\"\"\n",
    "    truth = (70 + 8 * x1 + 4 * x2\n",
    "             - 3 * x1 ** 2 - 2 * x2 ** 2\n",
    "             - 2 * x1 * x2)\n",
    "    return float(truth + rng.normal(scale=0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: 2x2 factorial plus one center run\n",
    "\n",
    "Spend five runs first: the four corners of a 2x2 plus the center.\n",
    "Fit a linear-plus-interaction model and check whether the center's\n",
    "measured response matches the model's prediction at ``(0, 0)``.  A\n",
    "big mismatch is the surface's way of telling you it is curving."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "x1 = c(-1, +1, -1, +1, 0, name=\"x1\")\n",
    "x2 = c(-1, -1, +1, +1, 0, name=\"x2\")\n",
    "y_obs = [observe(a, b) for a, b in zip([-1, 1, -1, 1, 0], [-1, -1, 1, 1, 0], strict=True)]\n",
    "y = c(*y_obs, name=\"y\")\n",
    "step1 = gather(x1=x1, x2=x2, y=y)\n",
    "step1"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "m_lin = lm(\"y ~ x1 + x2 + x1:x2\", step1)\n",
    "print(\"Linear-with-interaction fit on 2x2+center:\")\n",
    "print(m_lin.get_parameters(drop_intercept=False).to_string())\n",
    "centre_pred = float(predict(m_lin, x1=0, x2=0).iloc[0])\n",
    "print()\n",
    "print(f\"Predicted center y at (0, 0): {centre_pred:.2f}\")\n",
    "print(f\"Observed center y at (0, 0):  {y_obs[4]:.2f}\")\n",
    "print(f\"Difference: {y_obs[4] - centre_pred:+.2f}\")"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. admonition:: Solution\n",
    "\n",
    "   The intercept fits at ``~65.9`` from the corners alone, but the\n",
    "   center run measured ``~69.8`` - about four units above the linear\n",
    "   prediction.  Four units is much larger than the measurement noise\n",
    "   (sigma about 0.5), so the surface is **curving** inside the design\n",
    "   region.  The linear model is no longer enough; we need at least\n",
    "   a quadratic, which means adding **axial points**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: axial points -> central composite design (CCD)\n",
    "\n",
    "A central composite design adds four axial runs at\n",
    "``(+/-alpha, 0)`` and ``(0, +/-alpha)``.  Rotatable choice is\n",
    "``alpha = sqrt(2) ~= 1.414`` for two factors.\n",
    "\n",
    "Combining the original five runs with the four new ones gives a 9-run\n",
    "CCD that supports a full quadratic model\n",
    "``y = b_0 + b_1 x_1 + b_2 x_2 + b_11 x_1^2 + b_22 x_2^2 + b_12 x_1 x_2``."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "alpha = np.sqrt(2)\n",
    "x1_axial = [-alpha, alpha, 0, 0]\n",
    "x2_axial = [0, 0, -alpha, alpha]\n",
    "y_axial = [observe(a, b) for a, b in zip(x1_axial, x2_axial, strict=True)]\n",
    "\n",
    "x1_ccd = c(-1, +1, -1, +1, 0, *x1_axial, name=\"x1\")\n",
    "x2_ccd = c(-1, -1, +1, +1, 0, *x2_axial, name=\"x2\")\n",
    "y_ccd  = c(*y_obs, *y_axial, name=\"y\")\n",
    "ccd = gather(x1=x1_ccd, x2=x2_ccd, y=y_ccd)\n",
    "ccd"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "m_quad = lm(\"y ~ x1 + x2 + I(x1**2) + I(x2**2) + x1:x2\", ccd)\n",
    "print(\"Quadratic fit on the 9-run CCD:\")\n",
    "print(m_quad.get_parameters(drop_intercept=False).to_string())"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. admonition:: Guidance\n",
    "\n",
    "   The fitted intercept (``~69.8``) now sits at the *measured* center\n",
    "   value, not the linear-model extrapolation.  The quadratic\n",
    "   coefficients ``b_11 ~= -3.0`` and ``b_22 ~= -1.8`` confirm the\n",
    "   surface curves *downward* away from the center in both directions\n",
    "   - the classic peak shape.  ``b_12 ~= -2.1`` says the two factors\n",
    "   *interact*: pushing both high together gives diminishing returns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: find the predicted optimum\n",
    "\n",
    "For the quadratic ``y = b_0 + b_1 x_1 + b_2 x_2 + b_11 x_1^2 + b_22 x_2^2 + b_12 x_1 x_2``,\n",
    "the stationary point satisfies the gradient equation\n",
    "``[2 b_11, b_12; b_12, 2 b_22] * x = -[b_1; b_2]``.\n",
    "A two-by-two solve gives the predicted peak in closed form."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "p = m_quad.get_parameters(drop_intercept=False)\n",
    "A = np.array([[2 * p[\"I(x1 ** 2)\"], p[\"x1:x2\"]],\n",
    "              [p[\"x1:x2\"], 2 * p[\"I(x2 ** 2)\"]]])\n",
    "b = np.array([-p[\"x1\"], -p[\"x2\"]])\n",
    "x_star = np.linalg.solve(A, b)\n",
    "y_star = float(predict(m_quad, x1=x_star[0], x2=x_star[1]).iloc[0])\n",
    "print(f\"Predicted optimum: x1 = {x_star[0]:.3f}, x2 = {x_star[1]:.3f}\")\n",
    "print(f\"Predicted response at the optimum: y = {y_star:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Confirm the optimum with one more run.\n",
    "\n",
    "y_confirm = observe(x_star[0], x_star[1])\n",
    "print(f\"Confirmation run at the predicted optimum: y = {y_confirm:.2f}\")\n",
    "print(f\"Predicted: {y_star:.2f} | Observed: {y_confirm:.2f} | Diff: {y_confirm - y_star:+.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Contour plot of the fitted surface, overlaid with the runs.\n",
    "\n",
    "grid = np.linspace(-1.8, 1.8, 50)\n",
    "G1, G2 = np.meshgrid(grid, grid)\n",
    "Z = (p[\"Intercept\"]\n",
    "     + p[\"x1\"] * G1 + p[\"x2\"] * G2\n",
    "     + p[\"I(x1 ** 2)\"] * G1 ** 2\n",
    "     + p[\"I(x2 ** 2)\"] * G2 ** 2\n",
    "     + p[\"x1:x2\"] * G1 * G2)\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Contour(x=grid, y=grid, z=Z, contours_coloring=\"lines\",\n",
    "                         line_width=1, showscale=False))\n",
    "fig.add_trace(go.Scatter(x=ccd[\"x1\"], y=ccd[\"x2\"], mode=\"markers\",\n",
    "                         marker={\"size\": 9, \"color\": \"black\"},\n",
    "                         name=\"CCD runs\"))\n",
    "fig.add_trace(go.Scatter(x=[x_star[0]], y=[x_star[1]], mode=\"markers\",\n",
    "                         marker={\"size\": 14, \"symbol\": \"star\", \"color\": \"red\"},\n",
    "                         name=\"Predicted optimum\"))\n",
    "fig.update_layout(width=560, height=480,\n",
    "                  xaxis_title=\"x1 (coded)\", yaxis_title=\"x2 (coded)\",\n",
    "                  title=\"Fitted response surface with predicted optimum\")\n",
    "fig"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. admonition:: Solution\n",
    "\n",
    "   The fitted quadratic places the optimum near\n",
    "   ``(x1, x2) ~= (1.22, 0.33)`` with predicted response ``~75.6``.\n",
    "   The confirmation run at that point should land within a noise\n",
    "   width of the prediction.\n",
    "\n",
    "   **When to stop.**  The quadratic surface has a single stationary\n",
    "   point; the contour plot is concentric ellipses around it; and the\n",
    "   confirmation run matches the prediction.  All three signals agree -\n",
    "   you are at the optimum *of this model*.  Whether the optimum of the\n",
    "   model is also the optimum of the true system is a separate\n",
    "   question; a sensitivity analysis around the predicted point (run a\n",
    "   couple of nearby points and confirm none beat it) is the final\n",
    "   check.\n",
    "\n",
    ".. admonition:: Check yourself\n",
    "\n",
    "   In coded units the predicted optimum sits at about\n",
    "   ``(1.22, 0.33)``.  How would you convert that back to real-world\n",
    "   ecommerce parameters if the design used ``A: 4 to 6`` for menu\n",
    "   items and ``B: 8 to 16`` for products per page?"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. admonition:: Solution\n",
    "\n",
    "   The real-world value for any factor is\n",
    "   ``real = center + coded * half_range``.\n",
    "\n",
    "   - For A (menu items), center = 5, half-range = 1, so\n",
    "     ``A* = 5 + 1.22 * 1 = 6.22`` items.  In practice you would round\n",
    "     to 6 (since you cannot have fractional menu items) and run a\n",
    "     final confirmation at A = 6.\n",
    "   - For B (products per page), center = 12, half-range = 4, so\n",
    "     ``B* = 12 + 0.33 * 4 = 13.3`` products - round to 13."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wrap-up\n",
    "\n",
    "The 2-D response-surface workflow in one paragraph:\n",
    "\n",
    "1. **Factorial + center point** (5 runs): if the center disagrees with\n",
    "   the linear fit, you have curvature.\n",
    "2. **Add axial points** to make a CCD (9 runs total).  Fit a full\n",
    "   quadratic.\n",
    "3. **Solve the gradient = 0 equation** for the predicted peak.\n",
    "4. **Confirm** with one or two runs at and near the predicted peak.\n",
    "5. **Sensitivity check**: explore a small region around the optimum\n",
    "   to confirm no nearby point beats it.\n",
    "\n",
    "**Next:** Module 8 wraps up the whole course - we collect the\n",
    "vocabulary in one place, map it back to the modules, and point to\n",
    "the parts of `process_improve` that handle each step."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}