{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "26959ea4-522b-48ca-85ed-87cf5cd1ca76",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# An introduction to self-adaptive cooperative enhanced scatter search (saCeSS) in PyScat\n",
    "\n",
    "Goals:\n",
    "* Introduce the concepts of self-adaptive cooperative enhanced scatter search (saCeSS)\n",
    "* Show how to use the `pyscat.SacessOptimizer` and introduce its hyperparameters\n",
    "\n",
    "It is recommended to read the [eSS introduction](ess_intro.ipynb) first."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7112da22-10b7-4bc4-92f6-95dc4e14645c",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "The PyScat scatter search implementations are based on the following publications:\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.\n",
    "* David R. Penas, Patricia González, Jose A. Egea, Ramón Doallo, and Julio R. Banga. **Parameter estimation in large-scale systems biology models: a parallel and self-adaptive cooperative strategy**. BMC Bioinformatics, January 2017. doi:10.1186/s12859-016-1452-4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca2a845c1b45edd7",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import logging\n",
    "from pprint import pprint\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from pyscat import SacessOptimizer, get_default_ess_options\n",
    "from pyscat.plot import plot_sacess_history\n",
    "\n",
    "np.random.seed(1337)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd98a3c65ef0aa18",
   "metadata": {},
   "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": "4f41cfa1224f7555",
   "metadata": {},
   "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": "36ff97320e9c3c7",
   "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": "2d7df6fd26a7f02c",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Self-Adaptive Cooperative Enhanced Scatter Search (saCeSS)—`SacessOptimizer`\n",
    "\n",
    "**Motivation**\n",
    "* eSS makes it difficult to balance exploration and intensification\n",
    "* eSS hyperparameters are difficult to tune\n",
    "* eSS itself offers limited room for parallelization to reduce walltime\n",
    "\n",
    "**Approach**\n",
    "\n",
    "* **cooperative**: eSS (`ESSOptimizer`) instances with different degrees of exploration versus intensification are running concurrently and exchange promising solutions\n",
    "* **asynchronous communication**: non-blocking exchange of solutions\n",
    "* **self-adaptive**: concurrent eSS instances exchange hyperparameters\n",
    "\n",
    "**Implementation**\n",
    "\n",
    "Several `SacessWorker` are running in parallel, controlled by a `SacessManager` that handles global state, each running an instance of `ESSOptimizer`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4e36cc81d08fac6",
   "metadata": {},
   "source": [
    "### Optimization with default options\n",
    "\n",
    "For `SacessOptimizer` there are plenty of hyperparameters that can be tuned, but for a start we just use the default settings.\n",
    "The only options that have to be specified are:\n",
    "\n",
    "* the number of workers `num_workers` (i.e., the number `ESSOptimizer` instances that will run in parallel)\n",
    "* the walltime limit `max_walltime_s`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b62ede89931c03ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "optimizer = SacessOptimizer(\n",
    "    problem=problem,\n",
    "    num_workers=6,\n",
    "    max_walltime_s=2,\n",
    "    sacess_loglevel=logging.WARNING,\n",
    ")\n",
    "result = optimizer.minimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b02fee9c8f65508",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate default options for the individual eSS instances\n",
    "ess_options = get_default_ess_options(\n",
    "    num_workers=6, dim=problem.dim, local_optimizer=False\n",
    ")\n",
    "print(\"Options for the individual eSS instances:\")\n",
    "pprint(ess_options)\n",
    "\n",
    "# Initialize and run the optimizer\n",
    "sacess = SacessOptimizer(\n",
    "    problem=problem,\n",
    "    ess_init_args=ess_options,\n",
    "    max_walltime_s=2,\n",
    "    mp_start_method=\"fork\",\n",
    "    sacess_loglevel=logging.WARNING,\n",
    ")\n",
    "result = sacess.minimize()\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f40ed3a8c6252cc4",
   "metadata": {},
   "source": [
    "Visualize the optimization trajectory across iterations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdefe354632ed996",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_sacess_history(sacess.histories)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "plot_f()\n",
    "for i, history in enumerate(sacess.histories):\n",
    "    h = np.vstack(history.get_x_trace())\n",
    "    plt.plot(h[:, 0], h[:, 1], marker=\".\", label=f\"Worker {i}\")\n",
    "plt.legend(loc=\"center left\", bbox_to_anchor=(1.3, 0.5))\n",
    "\n",
    "r = np.vstack(result.optimize_result.x)\n",
    "plt.scatter(\n",
    "    r[:, 0], r[:, 1], c=\"white\", marker=\"*\", label=\"Reported optimum\", zorder=5\n",
    ")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d87834a6ad09c33",
   "metadata": {},
   "source": "Visualize exploration of the parameter space (the API is experimental, and recording all visited points is usually quite memory intensive):\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbc43d51b4287d31",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscat.eval_logger import EvalLogger\n",
    "\n",
    "el = EvalLogger()\n",
    "with el.attach(problem):\n",
    "    sacess = SacessOptimizer(\n",
    "        problem=problem,\n",
    "        ess_init_args=ess_options,\n",
    "        max_walltime_s=1,\n",
    "        mp_start_method=\"fork\",\n",
    "        sacess_loglevel=logging.WARNING,\n",
    "    )\n",
    "    result = sacess.minimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54cf46f958c356d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract traces\n",
    "x_trace = np.array([x for x, _ in el.evals])\n",
    "fx_trace = np.array([fx for _, fx in el.evals])\n",
    "\n",
    "plot_f()\n",
    "# plot all visited points\n",
    "plt.scatter(x_trace[:, 0], x_trace[:, 1], marker=\".\", alpha=0.5, c=\"w\")\n",
    "# plot optimum\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",
    "plt.legend(loc=\"center left\", bbox_to_anchor=(1.2, 0.5))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2b5b3ffff2608a6",
   "metadata": {},
   "source": [
    "### saCeSS hyperparameters\n",
    "\n",
    "**eSS settings**\n",
    "* `num_workers` / `ess_init_args`\n",
    "\n",
    "  saCeSS runs several ESSOptimizers where each instance can be configured independently.\n",
    "  To use (different) default configurations for each worker, `num_workers` can be passed; for full control over the eSS instances, worker-specific hyperparameters can be passed via `ess_init_args`.\n",
    "\n",
    "**Thresholds for propagating promising solutions**\n",
    "\n",
    "The current best parameters are not perfectly synchronized across the different workers. Promising new solutions will only be exchanged if they exceed a certain relative improvement threshold.\n",
    "\n",
    "* `manager_initial_rejection_threshold`: Initial threshold for rejecting solutions. This threshold will be halved every time `num_workers` solutions have been rejected in a row.\n",
    "* `manager_minimum_rejection_threshold`: Minimum threshold for rejecting solutions\n",
    "* `worker_acceptance_threshold`: Threshold for accepting solutions\n",
    "\n",
    "**Adaptation settings** (`adaptation_min_evals`, `adaptation_sent_coeff`, `adaptation_sent_offset`)\n",
    "\n",
    "Worker hyperparameters are updated if one of the following conditions is met:\n",
    "\n",
    "* The number of function evaluations since the last solution was sent\n",
    "  to the manager times the number of optimization parameters is greater\n",
    "  than ``adaptation_min_evals``.\n",
    "\n",
    "* The number of solutions received by the worker since the last\n",
    "  solution it sent to the manager is greater than\n",
    "  ``adaptation_sent_coeff * n_sent_solutions + adaptation_sent_offset``,\n",
    "  where ``n_sent_solutions`` is the number of solutions sent to the\n",
    "  manager by the given worker.\n",
    "\n",
    "**Exit criteria**\n",
    "\n",
    "  * `max_walltime_s`: So far, only a time limit is considered on top of the per-scatter-search function evaluation limit\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2505e3c6e738300a",
   "metadata": {},
   "source": [
    "\n",
    "### Parallelization within `SacessOptimizer`\n",
    "\n",
    "For parallelization of `pypesto.SacessOptimizer` optimizations, the following options are available:\n",
    "\n",
    "* **Parallelization of the individual eSS instances**: This is required and is based on `multiprocessing`. The number of processes is the number of workers.\n",
    "* **Parallelization of different objective function evaluations**: This is optional and can be controlled by `n_procs` and `n_threads` in the `ess_init_args` dictionary.\n",
    "* **Parallelization inside a single objective evaluation**: This is independent of `SacessOptimizer`"
   ]
  }
 ],
 "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
}