Parallel simulation refactoring

docs/how_to_guide/21_running_simulations_in_parallel.ipynb

@@ -0,0 +1,248 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# How to run simulations in parallel"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Many simulators are computationally expensive. When performing Simulation-Based Inference, we often need thousands of simulations. If your simulator takes even a fraction of a second, running them sequentially can be prohibitively slow.\n",
+    "\n",
+    "**Vectorization (or batching)** is one strategy to improve speed, where the simulator processes multiple parameters at once in a single call. However, this is not always straightforward to implement (e.g., rewriting legacy code) and can lead to high memory consumption if the batches are too large.\n",
+    "\n",
+    "**Parallelization** is an alternative that allows us to run multiple simulations at the same time in separate processes, reducing the total wall-clock time without requiring the simulator to be vectorized.\n",
+    "\n",
+    "This guide illustrates how to parallelize simulations using `sbi`'s `simulate_from_theta` utility. We will use a toy simulator that returns a file path—a common scenario where the simulator is an external executable that writes results to disk, making standard in-memory batching impossible or impractical."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import time\n",
+    "from typing import Any\n",
+    "\n",
+    "import numpy as np\n",
+    "\n",
+    "\n",
+    "def simulator(theta: Any) -> str:\n",
+    "    \"\"\"\n",
+    "    A toy simulator that sleeps and returns a dummy file path.\n",
+    "    \"\"\"\n",
+    "    time.sleep(0.1)  # Simulate expensive computation\n",
+    "\n",
+    "    # Format the parameter into a string suitable for a filename\n",
+    "    # We use 4 decimal places and replace the dot with a hyphen\n",
+    "    # e.g., 0.1234 -> 0-1234\n",
+    "    theta_id = f\"{theta[0]:.4f}\".replace(\".\", \"_\")\n",
+    "\n",
+    "    return f\"/path/to/simulation/output_{theta_id}.npy\"\n",
+    "\n",
+    "thetas = np.random.uniform(0, 1, size=(100, 1)) # 100 parameter sets"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Naive for loop\n",
+    "First, let's establish a baseline by running the simulations sequentially in a simple for loop. This represents the wall-clock time it takes without any parallelization optimization."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Elapsed time: 10.022978331000104 seconds\n",
+      "['/path/to/simulation/output_0-5192.npy', '/path/to/simulation/output_0-0983.npy', '/path/to/simulation/output_0-6823.npy', '/path/to/simulation/output_0-9948.npy', '/path/to/simulation/output_0-3902.npy', '/path/to/simulation/output_0-4715.npy', '/path/to/simulation/output_0-0382.npy', '/path/to/simulation/output_0-9635.npy', '/path/to/simulation/output_0-8807.npy', '/path/to/simulation/output_0-2155.npy']\n"
+     ]
+    }
+   ],
+   "source": [
+    "start = time.perf_counter()\n",
+    "x = [simulator(np.array(theta)) for theta in thetas]\n",
+    "end = time.perf_counter()\n",
+    "print(f\"Elapsed time: {end - start} seconds\")\n",
+    "print(x[:10])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Parallel execution\n",
+    "\n",
+    "`sbi` provides the `simulate_from_theta` utility to easily parallelize simulations. It uses `joblib` under the hood.\n",
+    "\n",
+    "We can use `joblib.parallel_config` context manager to specify the number of workers (`n_jobs`)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "3f2e3e21b42048b7b727bfc75756f4ba",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "  0%|          | 0/100 [00:00<?, ?it/s]"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Elapsed time: 4.3353425840000455 seconds\n",
+      "['/path/to/simulation/output_0-5192.npy', '/path/to/simulation/output_0-0983.npy', '/path/to/simulation/output_0-6823.npy', '/path/to/simulation/output_0-9948.npy', '/path/to/simulation/output_0-3902.npy', '/path/to/simulation/output_0-4715.npy', '/path/to/simulation/output_0-0382.npy', '/path/to/simulation/output_0-9635.npy', '/path/to/simulation/output_0-8807.npy', '/path/to/simulation/output_0-2155.npy']\n"
+     ]
+    }
+   ],
+   "source": [
+    "import joblib\n",
+    "\n",
+    "from sbi.utils.simulation_utils import simulate_from_theta\n",
+    "\n",
+    "start = time.perf_counter()\n",
+    "# Run simulations in parallel with 10 jobs\n",
+    "with joblib.parallel_config(n_jobs=10):\n",
+    "    thetas, x = simulate_from_theta(simulator, thetas)\n",
+    "end = time.perf_counter()\n",
+    "print(f\"Elapsed time: {end - start} seconds\")\n",
+    "print(x[:10])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Creating a vectorized simulator\n",
+    "\n",
+    "While `simulate_from_theta` is convenient for generating a static dataset, you might sometimes need a **simulator object** that can handle batches of parameters (e.g., to pass into an inference algorithm that runs simulations on-the-fly).\n",
+    "\n",
+    "The `parallelize_simulator` utility wraps your simulator and returns a new function that:\n",
+    "1. Accepts a batch of parameters.\n",
+    "2. Splits them into chunks.\n",
+    "3. Runs the chunks in parallel.\n",
+    "4. Re-assembles the results into a batch.\n",
+    "\n",
+    "This effectively \"vectorizes\" your simulator without rewriting its internal logic."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/local/home/bm284012/Workspace/sbi/sbi/utils/simulation_utils.py:255: UserWarning: Joblib is used for parallelization. It is recommended to use numpy arrays for the simulator input and output to avoid serialization overhead with torch tensors.\n",
+      "  return decorator(simulator)\n",
+      "/tmp/ipykernel_17395/1249179843.py:8: UserWarning: Simulation batch size is greater than 1, but simulator_is_batched is False. Simulations will be run sequentially (batch size 1).\n",
+      "  x = parallel_simulator(thetas)\n"
+     ]
+    },
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "1b32ae874cf147a7b407875b43270b5f",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "  0%|          | 0/100 [00:00<?, ?it/s]"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Elapsed time: 1.0284301520000554 seconds\n",
+      "['/path/to/simulation/output_0-5192.npy', '/path/to/simulation/output_0-0983.npy', '/path/to/simulation/output_0-6823.npy', '/path/to/simulation/output_0-9948.npy', '/path/to/simulation/output_0-3902.npy', '/path/to/simulation/output_0-4715.npy', '/path/to/simulation/output_0-0382.npy', '/path/to/simulation/output_0-9635.npy', '/path/to/simulation/output_0-8807.npy', '/path/to/simulation/output_0-2155.npy']\n"
+     ]
+    }
+   ],
+   "source": [
+    "from sbi.utils.simulation_utils import parallelize_simulator\n",
+    "\n",
+    "# Create a new simulator function that handles batches automatically\n",
+    "batched_simulator = parallelize_simulator(simulator)\n",
+    "\n",
+    "start = time.perf_counter()\n",
+    "# We can now call this new simulator with the full batch of parameters.\n",
+    "# The parallel execution is handled internally, governed by the joblib context.\n",
+    "with joblib.parallel_config(n_jobs=10):\n",
+    "    x = batched_simulator(thetas)\n",
+    "\n",
+    "end = time.perf_counter()\n",
+    "print(f\"Elapsed time: {end - start} seconds\")\n",
+    "print(x[:10])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "sbi",
+   "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.11"
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": false
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

sbi/examples/minimal.py

@@ -46,6 +46,10 @@ def flexible():
     inference = SNPE(prior)
 
     theta, x = simulate_for_sbi(simulator, proposal=prior, num_simulations=500)
+    # The simulator returns a Tensor,
+    # but `simulate_for_sbi` returns `Tensor | List[str]`.
+    # To satisfy type checking, we assert that x is a Tensor.
+    assert isinstance(x, torch.Tensor)
     density_estimator = inference.append_simulations(theta, x).train()
     posterior = inference.build_posterior(density_estimator)
     posterior.sample((100,), x=x_o)