{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "26959ea4-522b-48ca-85ed-87cf5cd1ca76",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# An introduction to enhanced scatter search (ESS) in PyScat\n",
    "\n",
    "Goals:\n",
    "* Introduce the concepts of enhanced scatter search (eSS)\n",
    "* Show how to use the `pyscat.ESSOptimizer` and introduce its hyperparameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7112da22-10b7-4bc4-92f6-95dc4e14645c",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "The PyScat scatter search implementations is based on:\n",
    "\n",
    "Jose A. Egea, Eva Balsa-Canto, María-Sonia G. García, and Julio R. Banga. **Dynamic optimization of nonlinear processes with an enhanced scatter search method**. Industrial & Engineering Chemistry Research, 48(9):4388–4401, April 2009. doi:[10.1021/ie801717t](https://doi.org/10.1021/ie801717t)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec84235f52ca52a9",
   "metadata": {},
   "source": [
    "## What is scatter search?\n",
    "\n",
    "Scatter search is a **meta-heuristic for global optimization**. It is based on the idea of exploring the parameter space by evolving a population of diverse candidate solutions, i.e., an evolutionary algorithm. The fundamental challenge is to balance exploration and exploitation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca2a845c1b45edd7",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from itertools import product\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from IPython.display import display_markdown\n",
    "from pypesto.history import MemoryHistory, NoHistory\n",
    "\n",
    "# Note that this demo uses some private API,\n",
    "#  not relevant to regular users, that may change without notice.\n",
    "from pyscat import ESSOptimizer\n",
    "from pyscat.function_evaluator import FunctionEvaluator\n",
    "from pyscat.refset import RefSet\n",
    "\n",
    "np.random.seed(1337)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c18bb58bee4a36dc",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Set up problem\n",
    "\n",
    "To run any optimization, we first need to specify the optimization problem. PyScat currently heavily relies on the pyPESTO framework and requires a `pypesto.Problem`.\n",
    "For this demo, we use the Schwefel function which is one of the examples included in PyScat:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "507488411355f0b2",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from pyscat.examples import plot_problem, problem_info, xyz\n",
    "\n",
    "cur_problem_info = problem_info[\"Schwefel\"]\n",
    "\n",
    "problem = cur_problem_info[\"problem\"]\n",
    "\n",
    "plot_problem(problem, title=\"Schwefel function\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad74797e442bf656",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate data for plotting\n",
    "X, Y, Z = xyz(problem)\n",
    "\n",
    "\n",
    "# plotting function for our objective landscape\n",
    "def plot_f(ax=None):\n",
    "    \"\"\"contour plot\"\"\"\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "\n",
    "    c = ax.contourf(X, Y, Z, cmap=\"viridis\")\n",
    "    plt.colorbar(c, ax=ax, label=\"fval\")\n",
    "    ax.set_xlabel(\"$x_1$\")\n",
    "    ax.set_ylabel(\"$x_2$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e84a290114d55e3c",
   "metadata": {},
   "source": [
    "## Enhanced Scatter Search (eSS) --- `ESSOptimizer`\n",
    "\n",
    "The idea of ESS is to maintain some reference set (RefSet) comprising a constant number of points (`dim_refset`) that (a) explores the parameter space and (b) approaches minima.\n",
    "\n",
    "The basic steps of ESS are:\n",
    "\n",
    "* **Initialization**: Generate a diverse set of points in the parameter space.\n",
    "* **Recombination**: Generate new points by recombining the RefSet.\n",
    "* **Improvement**: Improve the RefSet by replacing points with better ones.\n",
    "\n",
    "The steps are repeated until a stopping criterion is met.\n",
    "\n",
    "ESS itself is gradient-free, but if gradient information is available, a gradient-based local optimizer can be used during the *Improvement* step (see below).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b0f29cb7e52c5e1",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "### Initialize ESS\n",
    "\n",
    "Create the initial RefSet:\n",
    "* Sample `n_diverse` points\n",
    "* Fill half of the RefSet with the best points\n",
    "* Fill the other half with random points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "786d335c448b7b6a",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# number of points in the RefSet\n",
    "# (chosen for visualization, not a general recommendation)\n",
    "REFSET_SIZE = 8\n",
    "N_DIVERSE = 2 * REFSET_SIZE\n",
    "\n",
    "# just some object that will evaluate\n",
    "# the objective function and help us sample random points\n",
    "evaluator = FunctionEvaluator(problem)\n",
    "# create initial population\n",
    "x, fx = evaluator.multiple_random(N_DIVERSE)\n",
    "order = np.argsort(fx)\n",
    "# the first half of the refset is the best points\n",
    "# the second half is randomly selected\n",
    "order[int(REFSET_SIZE / 2) :] = np.random.permutation(\n",
    "    order[int(REFSET_SIZE / 2) :]\n",
    ")\n",
    "x = x[order]\n",
    "fx = fx[order]\n",
    "\n",
    "# initialize RefSet\n",
    "refset = RefSet(\n",
    "    x=x[:REFSET_SIZE, :],\n",
    "    fx=fx[:REFSET_SIZE],\n",
    "    evaluator=evaluator,\n",
    ")\n",
    "refset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee09c19b4067c5ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualize initialization\n",
    "fig, axs = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(16, 6))\n",
    "\n",
    "ax = axs[0]\n",
    "plot_f(ax=ax)\n",
    "ax.scatter(\n",
    "    x[:REFSET_SIZE, 0],\n",
    "    x[:REFSET_SIZE, 1],\n",
    "    c=\"yellow\",\n",
    "    marker=\"*\",\n",
    "    label=\"Initial RefSet\",\n",
    ")\n",
    "ax.scatter(\n",
    "    x[REFSET_SIZE:, 0],\n",
    "    x[REFSET_SIZE:, 1],\n",
    "    c=\"red\",\n",
    "    marker=\".\",\n",
    "    label=\"Dismissed\",\n",
    ")\n",
    "ax.legend(loc=\"center left\", bbox_to_anchor=(1.2, 1.1))\n",
    "ax.set_title(\"Initialization – n_diverse\")\n",
    "\n",
    "ax = axs[1]\n",
    "plot_f(ax=ax)\n",
    "ax.scatter(\n",
    "    x[:REFSET_SIZE, 0],\n",
    "    x[:REFSET_SIZE, 1],\n",
    "    c=\"yellow\",\n",
    "    marker=\"*\",\n",
    "    label=\"refset\",\n",
    ")\n",
    "ax.set_title(\"Initialization – RefSet\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33fd4676593c0e1a",
   "metadata": {},
   "source": [
    "Hyperparameter controlling initialization:\n",
    "\n",
    "* `dim_refset`: number of points in the RefSet\n",
    "* `n_diverse`: number of initial random points to generate (default: `10 * dim_refset`)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7360c980a7df19d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "### Recombination --- generate new points based on the current RefSet\n",
    "\n",
    "Every iteration generates $dim\\_refset^2 - dim\\_refset$ new points from pairwise recombination of the RefSet members. Only the best offspring from each RefSet member will be retained. A new point can only replace its parent, but not any other RefSet member. This ensures that the RefSet remains diverse.\n",
    "\n",
    "Currently, there are no hyperparameters to control recombination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fc7d6add68ce157",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    f\"RefSet size if {refset.dim}. \"\n",
    "    f\"Thus, recombination will generate {refset.dim**2 - refset.dim} \"\n",
    "    f\"new points in each iteration.\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f47dfe9d2a545003",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recombination scheme\n",
    "from pyscat.ess import DefaultRecombination\n",
    "\n",
    "recombinator = DefaultRecombination()\n",
    "# for i, j in [(0, 8), (8, 0)]:\n",
    "all_pairs = list(product(range(refset.dim), range(refset.dim)))\n",
    "for i, j in all_pairs[:3] + all_pairs[-4:-1]:\n",
    "    if i == j:\n",
    "        continue\n",
    "    c1, c2 = recombinator.get_hyper_rect(refset, evaluator, i, j)\n",
    "    new_x = np.random.uniform(low=c1, high=c2, size=problem.dim)\n",
    "\n",
    "    marker_i, marker_j = (\"o\", \".\") if i < j else (\".\", \"o\")\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(6, 4))\n",
    "    plot_f(ax=ax)\n",
    "    ax.scatter(x[i, 0], x[i, 1], c=\"red\", marker=marker_i, label=\"refset\")\n",
    "    ax.text(x[i, 0], x[i, 1], str(i), color=\"red\")\n",
    "    ax.scatter(x[j, 0], x[j, 1], c=\"red\", marker=marker_j, label=\"refset\")\n",
    "    ax.text(x[j, 0], x[j, 1], str(j), color=\"red\")\n",
    "    # draw the rectangle\n",
    "    ax.plot(\n",
    "        [c1[0], c2[0], c2[0], c1[0], c1[0]],\n",
    "        [c1[1], c1[1], c2[1], c2[1], c1[1]],\n",
    "        c=\"white\",\n",
    "    )\n",
    "    # draw the new point\n",
    "    ax.scatter(new_x[0], new_x[1], c=\"white\", marker=\"*\", label=\"new\")\n",
    "    # label the new point\n",
    "    ax.text(new_x[0], new_x[1], \"new\", color=\"white\")\n",
    "    ax.title.set_text(f\"Recombination of {i} and {j}\")\n",
    "    plt.show()\n",
    "    print(f\"#{i}: {problem.objective(x[i])}\")\n",
    "    print(f\"#{j}: {problem.objective(x[j])}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30bab9a0cfb9a1fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "plot_f(ax=ax)\n",
    "\n",
    "all_children = []\n",
    "x_best_children = []\n",
    "for i in range(refset.dim):\n",
    "    children_i: list[np.ndarray] = []\n",
    "    for j in range(refset.dim):\n",
    "        if i == j:\n",
    "            continue\n",
    "        c1, c2 = recombinator.get_hyper_rect(refset, evaluator, i, j)\n",
    "        new_x = np.random.uniform(low=c1, high=c2, size=problem.dim)\n",
    "        children_i.append(new_x)\n",
    "        # plot children with markers and plot line from parent to child\n",
    "        ax.plot(\n",
    "            [refset.x[i, 0], new_x[0]],\n",
    "            [refset.x[i, 1], new_x[1]],\n",
    "            c=\"red\",\n",
    "            linestyle=\"--\",\n",
    "            linewidth=1,\n",
    "            zorder=1,\n",
    "        )\n",
    "        ax.scatter(new_x[0], new_x[1], c=\"white\", marker=\"*\", label=\"new\")\n",
    "\n",
    "    best_child_idx = np.array(\n",
    "        [problem.objective(x) for x in children_i]\n",
    "    ).argmin()\n",
    "    ax.scatter(\n",
    "        children_i[best_child_idx][0],\n",
    "        children_i[best_child_idx][1],\n",
    "        c=\"yellow\",\n",
    "        marker=\"*\",\n",
    "        label=\"new_best\",\n",
    "        zorder=10,\n",
    "    )\n",
    "    all_children.append(children_i)\n",
    "    x_best_children.append(children_i[best_child_idx])\n",
    "\n",
    "    ax.scatter(\n",
    "        refset.x[i, 0],\n",
    "        refset.x[i, 1],\n",
    "        c=\"red\",\n",
    "        marker=\".\",\n",
    "        label=\"refset\",\n",
    "        zorder=5,\n",
    "    )\n",
    "\n",
    "plt.legend(\n",
    "    [\n",
    "        plt.Line2D([0], [0], marker=\".\", c=\"red\", linestyle=\"None\"),\n",
    "        plt.Line2D([0], [0], marker=\"*\", c=\"white\", linestyle=\"None\"),\n",
    "        plt.Line2D([0], [0], marker=\"*\", c=\"yellow\", linestyle=\"None\"),\n",
    "    ],\n",
    "    [\"parent\", \"child\", \"best child\"],\n",
    ")\n",
    "plt.title(\"Recombination – parent & children\")\n",
    "\n",
    "\n",
    "# new plot with parents and best children connected by lines\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "plot_f(ax=ax)\n",
    "for i in range(refset.dim):\n",
    "    ax.plot(\n",
    "        [refset.x[i, 0], x_best_children[i][0]],\n",
    "        [refset.x[i, 1], x_best_children[i][1]],\n",
    "        c=\"red\",\n",
    "        linestyle=\"--\",\n",
    "        linewidth=1,\n",
    "        zorder=1,\n",
    "    )\n",
    "    ax.scatter(\n",
    "        x_best_children[i][0],\n",
    "        x_best_children[i][1],\n",
    "        c=\"yellow\",\n",
    "        marker=\"*\",\n",
    "        label=\"new_best\",\n",
    "        zorder=10,\n",
    "    )\n",
    "    ax.scatter(\n",
    "        refset.x[i, 0],\n",
    "        refset.x[i, 1],\n",
    "        c=\"red\",\n",
    "        marker=\".\",\n",
    "        label=\"refset\",\n",
    "        zorder=5,\n",
    "    )\n",
    "plt.title(\"Recombination – parent & best child\")\n",
    "\n",
    "x_best_children = np.array(x_best_children)\n",
    "fx_best_children = np.array([problem.objective(x) for x in x_best_children])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f4b055f54878938",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "### Go-beyond\n",
    "\n",
    "The *go-beyond strategy* aims at improving the best children from recombination. If the offspring is better than the parent, the offspring will be used as the new parent for the next iteration. The offspring will be generated by sampling from a hyper-rectangle around the parent. The size of the hyper-rectangle is determined by the relative improvement of the offspring over the parent. The process is repeated until no further improvement is possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec7bef9932937b25",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Re-implementation of the go-beyond strategy\n",
    "# because we need some internal state\n",
    "def go_beyond(\n",
    "    x_best_children: np.ndarray,\n",
    "    fx_best_children: np.ndarray,\n",
    "    refset: RefSet,\n",
    "    evaluator: FunctionEvaluator,\n",
    ") -> tuple[list[np.ndarray], list[np.ndarray], list[list[tuple[np.ndarray]]]]:\n",
    "    trials_x = []\n",
    "    trials_fx = []\n",
    "    rects = []\n",
    "    for i in range(refset.dim):\n",
    "        cur_trials_x = [x_best_children[i][np.newaxis, :]]\n",
    "        cur_trials_fx = [fx_best_children[i]]\n",
    "        cur_rects = []\n",
    "\n",
    "        if fx_best_children[i] >= refset.fx[i]:\n",
    "            # include child before go-beyond,\n",
    "            # since x_best_children will be updated here\n",
    "            trials_x.append(cur_trials_x)\n",
    "            trials_fx.append(cur_trials_fx)\n",
    "            continue\n",
    "\n",
    "        # offspring is better than parent\n",
    "        x_parent = refset.x[[i]].copy()\n",
    "        fx_parent = refset.fx[i]\n",
    "        x_child = x_best_children[[i]].copy()\n",
    "        fx_child = fx_best_children[i]\n",
    "        improvement = 1\n",
    "        # Multiplier used in determining the hyper-rectangle from which to\n",
    "        # sample children. Will be increased in case of 2 consecutive\n",
    "        # improvements.\n",
    "        # (corresponds to 1/\\Lambda in [Egea2009]_ algorithm 1)\n",
    "        go_beyond_factor = 1\n",
    "        while fx_child < fx_parent:\n",
    "            # update best child\n",
    "            x_best_children[i] = x_child\n",
    "            fx_best_children[i] = fx_child\n",
    "\n",
    "            # create new solution, child becomes parent\n",
    "            # hyper-rectangle for sampling child\n",
    "            box_lb = x_child - (x_parent - x_child) * go_beyond_factor  # if\n",
    "            box_ub = x_child\n",
    "            # clip to bounds\n",
    "            ub, lb = evaluator.problem.ub, evaluator.problem.lb\n",
    "            box_lb = np.fmax(np.fmin(box_lb, ub), lb)\n",
    "            box_ub = np.fmax(np.fmin(box_ub, ub), lb)\n",
    "            # sample parameters\n",
    "            x_new = np.random.uniform(low=box_lb, high=box_ub)\n",
    "            cur_rects.append(\n",
    "                (box_lb, box_ub),\n",
    "            )\n",
    "            x_parent = x_child\n",
    "            fx_parent = fx_child\n",
    "            x_child = x_new\n",
    "            fx_child = evaluator.single(x_child)\n",
    "            cur_trials_x.append(x_child)\n",
    "            cur_trials_fx.append(fx_child)\n",
    "            improvement += 1\n",
    "            if improvement == 2:\n",
    "                go_beyond_factor *= 2\n",
    "                improvement = 0\n",
    "        trials_x.append(cur_trials_x)\n",
    "        trials_fx.append(cur_trials_fx)\n",
    "        rects.append(cur_rects)\n",
    "    trials_x = list(map(np.vstack, trials_x))\n",
    "    trials_fx = list(map(np.array, trials_fx))\n",
    "    return trials_x, trials_fx, rects\n",
    "\n",
    "\n",
    "trials_x, trials_fx, go_beyond_rects = go_beyond(\n",
    "    x_best_children, fx_best_children, refset, evaluator\n",
    ")\n",
    "assert len(trials_x) == len(trials_fx) == refset.dim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ff44d48039de01d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# new plot with parents and best children connected by lines\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "plot_f(ax=ax)\n",
    "for i in range(refset.dim):\n",
    "    # line from old refset member to the best child from recombination\n",
    "    ax.plot(\n",
    "        [refset.x[i, 0], trials_x[i][0][0]],\n",
    "        [refset.x[i, 1], trials_x[i][0][1]],\n",
    "        c=\"grey\",\n",
    "        linestyle=\"-\",\n",
    "        linewidth=1,\n",
    "        zorder=1,\n",
    "    )\n",
    "\n",
    "    # final points after go-beyond\n",
    "    ax.scatter(\n",
    "        x_best_children[i][0],\n",
    "        x_best_children[i][1],\n",
    "        c=\"yellow\",\n",
    "        marker=\"*\",\n",
    "        label=\"new_best\",\n",
    "        zorder=10,\n",
    "    )\n",
    "\n",
    "    # before go-beyond\n",
    "    ax.scatter(\n",
    "        refset.x[i, 0],\n",
    "        refset.x[i, 1],\n",
    "        c=\"red\",\n",
    "        marker=\".\",\n",
    "        label=\"refset\",\n",
    "        zorder=5,\n",
    "    )\n",
    "    ax.text(refset.x[i, 0], refset.x[i, 1], str(i))\n",
    "\n",
    "    # the go-beyond line\n",
    "    ax.plot(\n",
    "        trials_x[i][:-1, 0],\n",
    "        trials_x[i][:-1, 1],\n",
    "        linestyle=\"-\",\n",
    "        c=\"white\",\n",
    "        marker=\".\",\n",
    "        zorder=1,\n",
    "        linewidth=1,\n",
    "    )\n",
    "    for x, y in trials_x[i]:\n",
    "        if problem.objective([x, y]) < refset.fx[i]:\n",
    "            # ax.scatter(x, y)\n",
    "            ax.text(x, y, str(i))\n",
    "\n",
    "plt.title(\"Go-beyond – parent & best child\")\n",
    "plt.legend(\n",
    "    [\n",
    "        plt.Line2D([0], [0], marker=\".\", c=\"red\", linestyle=\"None\"),\n",
    "        plt.Line2D([0], [0], marker=\"*\", c=\"yellow\", linestyle=\"None\"),\n",
    "        plt.Line2D([0], [0], marker=\".\", c=\"white\", linestyle=\"-\"),\n",
    "    ],\n",
    "    [\"parent\", \"best child\", \"go-beyond steps\"],\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed8ab98e800f6992",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "### Local search (optional)\n",
    "\n",
    "If a local optimizer is provided (`local_optmizer`), *one* of the current solutions *may* be improved further during each eSS iteration.\n",
    "This is heavily controlled by hyperparameters:\n",
    "\n",
    "* whether a local search will take place at all:\n",
    "  * `local_n1`: The number of eSS iterations after which the *first* local search will take place\n",
    "  * `local_n2`: The number of eSS iterations that have to pass after the previous local optimization\n",
    "* which solution will be optimized:\n",
    "  * `balance` $\\in [0, 1]$: A priority score will be calculated for each RefSet member based on 1) its function value rank 2) its rank-based distance to previous *local optimization results*.\n",
    "    `balance` determines the relative weight of (1) and (2) --- a lower value favors quality, a higher value favors diversity.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "924cbb37f9ce7c56",
   "metadata": {},
   "source": [
    "### Replacement of stuck RefSet members\n",
    "\n",
    "If there was no improvement of a specific RefSet member during recombination, go-beyond, or local optimization over `n_change` (default: `20`) iterations, the respective RefSet member will be replaced by a random point."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a867dc0cc4b1e4a4",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### Full example\n",
    "\n",
    "Run a `ESSOptimizer` optimization and show the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee138156cefc2233",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize and run the optimizer\n",
    "ess = ESSOptimizer(max_eval=5000, dim_refset=10)\n",
    "result = ess.minimize(problem)\n",
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95f81a5d7103eaaf",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Show the optimization result\n",
    "display_markdown(result.summary().replace(\"#\", \"##\"), raw=True)\n",
    "\n",
    "# Visualize the optimum\n",
    "plot_f()\n",
    "plt.scatter(\n",
    "    result.optimize_result.x[0][0],\n",
    "    result.optimize_result.x[0][1],\n",
    "    c=\"magenta\",\n",
    "    marker=\"*\",\n",
    "    label=\"Reported optimum\",\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aeb82a0564120437",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert (\n",
    "    abs(cur_problem_info[\"global_best\"] - result.optimize_result.fval[0])\n",
    "    < 1e-3\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d79eb3d4bf2a0aa",
   "metadata": {},
   "source": [
    "Visualize optimizer trajectory across iterations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8294ab552defb297",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscat.plot import plot_ess_history\n",
    "\n",
    "plot_ess_history(ess.history)\n",
    "plt.axhline(\n",
    "    cur_problem_info[\"global_best\"],\n",
    "    linestyle=\"dotted\",\n",
    "    label=\"known global optimum\",\n",
    ")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "plot_f()\n",
    "h = np.vstack(ess.history.get_x_trace())\n",
    "plt.plot(h[:, 0], h[:, 1], marker=\".\", c=\"white\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91f442dc21530999",
   "metadata": {},
   "source": [
    "Let's look at the exploration of the parameter space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2231ebdb1d5906a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show exploration of the parameter space\n",
    "problem.objective.history = MemoryHistory()\n",
    "ess = ESSOptimizer(max_eval=5000, dim_refset=10)\n",
    "result = ess.minimize(problem)\n",
    "x_trace = np.vstack(problem.objective.history.get_x_trace())\n",
    "fval_trace = np.vstack(problem.objective.history.get_fval_trace())\n",
    "problem.objective.history = NoHistory()\n",
    "\n",
    "plot_f()\n",
    "plt.scatter(\n",
    "    x_trace[:, 0],\n",
    "    x_trace[:, 1],\n",
    "    c=np.arange(len(x_trace)),\n",
    "    cmap=\"Greys\",\n",
    "    marker=\".\",\n",
    "    label=\"function evaluation\",\n",
    "    alpha=1,\n",
    ")\n",
    "plt.scatter(\n",
    "    result.optimize_result.x[0][0],\n",
    "    result.optimize_result.x[0][1],\n",
    "    c=\"magenta\",\n",
    "    marker=\"*\",\n",
    "    label=\"Reported optimum\",\n",
    ")\n",
    "plt.gcf().set_size_inches(12, 8)\n",
    "\n",
    "# add color bar for function evaluations, normalize to length\n",
    "cbar = plt.colorbar(\n",
    "    plt.cm.ScalarMappable(\n",
    "        cmap=\"Greys\", norm=plt.Normalize(vmin=0, vmax=len(x_trace))\n",
    "    ),\n",
    "    ax=plt.gca(),\n",
    ")\n",
    "cbar.set_label(\"function evaluation index\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39af3f690114caf4",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "### General hyperparameters for `ESSOptimizer`\n",
    "\n",
    "* parallelization of objective evaluation (mutually exclusive):\n",
    "  * `n_procs`: Number of processes for `multiprocessing`-based parallelization\n",
    "  * `n_threads`: Number of threads for threading-based parallelization\n",
    "* exit criteria:\n",
    "  * `max_eval`: Maximum number of objective evaluations\n",
    "  * `max_walltime_s`: Maximum walltime (seconds)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}