Models/Predator-Prey/predprey_client.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ca79c985",
   "metadata": {},
   "source": [
    "# Monte Carlo Simulation - Client for Predator-Prey Model\n",
    "\n",
    "## Overview\n",
    "\n",
    "The Uncertainty Quantification (UQ) problem in this context is to evaluate the uncertainty in the predator-prey model outcomes due to **uncertainty in the initial conditions or parameters**, specifically we want to estimate the expected value of the prey population $\\mathbb{E}[u_1(T)]$ at a future time T given some uncertainty in the initial conditions.\n",
    "\n",
    "As explained in the file [predprey_UQ.ipynb](https://github.com/MathSEE-Modeling-Week/Modeling-Week/blob/main/Models/Predator-Prey/predprey_UQ.ipynb), the uncertainty in the predator-prey model comes from the initial conditions. The true initial condition $\\bar{u}_0 = [0.5, 2]$ is perturbed with a uniform distribution $\\mathcal{U}(\\Gamma)$ over a square region $\\Gamma$ of radius $\\delta = 0.2$. This perturbed initial condition is then used to run a Monte Carlo simulation.\n",
    "\n",
    "The key objective is to calculate an estimate of $\\mathbb{E}[u_1(T)]$, the expected prey population at time $T = 6$, using the **Monte Carlo (MC) method**.\n",
    "\n",
    "## Step-by-Step Process\n",
    "\n",
    "### 1. Sampling Initial Conditions\n",
    "\n",
    "The third file defines that the uncertain initial condition follows a uniform distribution $\\mathcal{U}(\\Gamma)$, where $\\Gamma$ is a square region centered at $\\bar{u}_0 = [0.5, 2]$ with side lengths $2\\delta = 0.4$. We generate $N = 1000$ Monte Carlo samples of initial conditions within this region.\n",
    "\n",
    "### 2. Solving the Lotka-Volterra Equations\n",
    "\n",
    "For each sampled initial condition, we solve the **Lotka-Volterra equations** over the time interval $[0, T]$, where $T = 6$. The equations describe the dynamics between prey and predator populations, with fixed parameter values $\\theta_1 = \\theta_2 = \\theta_{12} = \\theta_{21} = 1$, as set in [predprey_UQ.ipynb](https://github.com/MathSEE-Modeling-Week/Modeling-Week/blob/main/Models/Predator-Prey/predprey_UQ.ipynb).\n",
    "\n",
    "We use the **umbridge model** to solve the system of ODEs for each initial condition.\n",
    "\n",
    "### 3. Monte Carlo Estimator for Expected Prey Population\n",
    "\n",
    "Once the prey population $u_1(T)$ is computed for all $N = 1000$ samples, we estimate the expected value $\\mathbb{E}[u_1(T)]$ by averaging the results:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbb{E}[u_1(T)] \\approx \\frac{1}{N} \\sum_{i=1}^N u_1^{(i)}(T)\n",
    "\\end{equation}\n",
    "\n",
    "### 4. Comparison with Unperturbed Initial Condition\n",
    "\n",
    "We also compute the prey population $u_1(T)$ for the **unperturbed initial condition** $\\bar{u}_0 = [0.5, 2]$ to provide a baseline for comparison with the perturbed trajectories.\n",
    "\n",
    "### 5. Visualization of Results\n",
    "\n",
    "The results are visualized using a histogram of the Monte Carlo samples of $u_1(T)$. The mean (MC estimator) and the prey population for the unperturbed trajectory are marked on the plot for comparison.\n",
    "\n",
    "\n",
    "## Conclusion\n",
    "\n",
    "This file demonstrates the use of the Monte Carlo method for **Uncertainty Quantification** in the Predator-Prey model. By generating multiple samples of the uncertain initial conditions and solving the Lotka-Volterra equations for each, we obtain an estimate of the expected prey population at a future time T. The results are compared to the unperturbed scenario to highlight the effects of uncertainty on the system dynamics."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d492a9e",
   "metadata": {},
   "source": [
    "### Task:\n",
    "Try implementing a MC simulation using the umbridge framework, where the predator-prey model is running as a server see [here](https://github.com/MathSEE-Modeling-Week/Modeling-Week/blob/main/Models/Predator-Prey/predprey_server.ipynb) and the MC simulation as an umbridge client.\n",
    "\n",
    "**In the next section you cand find an example code.**\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ced2f8b9",
   "metadata": {},
   "source": [
    "## Client for MC simulation\n",
    "In the following you can see the client code corresponding to the MC simulation. Be sure that your model is running [predprey_server.ipynb](https://github.com/MathSEE-Modeling-Week/Modeling-Week/blob/main/Models/Predator-Prey/predprey_server.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cd80bf85",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MC estimator for expected value of u1(T) at T=6: 1.5181587659330416\n",
      "Value of u1(T) for the unperturbed trajectory for T=6: 1.4098300518932472\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "import matplotlib.pyplot as plt\n",
    "import umbridge\n",
    "\n",
    "# Connect to the umbridge model\n",
    "model = umbridge.HTTPModel(\"http://0.0.0.0:4242\", \"forward\")\n",
    "\n",
    "# Define parameters\n",
    "theta1, theta2, theta12, theta21 = 1.0, 1.0, 1.0, 1.0\n",
    "T = 6\n",
    "u0_bar = np.array([0.5, 2])  # Example unperturbed initial condition\n",
    "\n",
    "# Generate perturbed initial conditions\n",
    "delta = 0.2\n",
    "N = 1000  # Number of MC samples\n",
    "\n",
    "np.random.seed(42)\n",
    "\n",
    "initial_conditions = u0_bar + delta * (2 * np.random.rand(N, 2) - 1)\n",
    "\n",
    "# Solve ODE for each initial condition using the umbridge model\n",
    "u1_T_values = []  # List for storing values of u1(T), i.e., prey, at T=6\n",
    "for u0 in initial_conditions:\n",
    "    # Create the input for the model\n",
    "    parameters = [T, u0[0], u0[1], theta1, theta2, theta12, theta21]\n",
    "    result = model([parameters])\n",
    "    u1_T_values.append(result[0][0])  # Extract u1(T) from the result\n",
    "\n",
    "# Calculate MC estimator\n",
    "u1_T_mean = np.mean(u1_T_values)\n",
    "print(f\"MC estimator for expected value of u1(T) at T={T}: {u1_T_mean}\")\n",
    "\n",
    "# Solve ODE for the unperturbed initial condition using the umbridge model\n",
    "parameters_unperturbed = [T, u0_bar[0], u0_bar[1], theta1, theta2, theta12, theta21]\n",
    "result_unperturbed = model([parameters_unperturbed])\n",
    "u1_T_unperturbed = result_unperturbed[0][0]\n",
    "print(f\"Value of u1(T) for the unperturbed trajectory for T=6: {u1_T_unperturbed}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5f6d81a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}