---
author:
- "[Eugeniy E. Mikhailov](http://physics.wm.edu/~evmik/) and [Greg
  Bentsen](https://physics.wm.edu/~gbentsen/)"
title: From simulation to optimization
---

## Agenda

Today

-   Review 2D NumPy arrays and indexing
-   Simulate a virus outbreak on a lattice
-   Optimize a vaccination policy
-   Simulate a forest fire and discover percolation
-   Optimize tree removal to prevent fire spread

\

**Main ideas:**

-   A **simulation** predicts what happens for a choice of parameters.
-   **Optimization** selects parameters to achieve a desired outcome.

# From Simulation to Optimization

## Simulation and Optimization

A **simulation** takes a model and parameters and predicts an outcome.

$$\text{parameters }\theta\ \longrightarrow\ \text{simulation}\ \longrightarrow\ \text{outcome}.$$

An **objective function** assigns a numerical cost to that outcome:

$$J(\theta)=\text{cost associated with the outcome of parameters }\theta.$$

**Optimization** means finding the parameters with the lowest cost:

$$\boxed{\theta^*=\arg\min_{\theta}J(\theta).}$$

\

Examples:

-   Simulating the path of a rocket from earth to mars
-   Simulating the behavior of the stock market
-   Simulating the flow of traffic and traffic lights

## Simulation and Optimization

This builds on our iterated maps from last week:

-   Simple rules
-   Iteration
-   Complex behavior

\

But this time we'll be adding a few new ingredients

-   randomness
-   objective functions
-   modification of parameters to optimize the objective function

\

Today: two spatial simulations, two objective functions.

## Review: Indexing a 2D array

A 2D NumPy array uses **row, column** indexing:

``` {.python .numberLines}
A = np.array([
    [0.3, 0.8, 0.1, 0.6, 0.9],
    [0.7, 0.2, 0.5, 0.4, 0.1],
    [0.9, 0.6, 0.3, 0.8, 0.2],
    [0.4, 0.1, 0.7, 0.5, 0.6]
])

print(A[2, 0])       # one element

print(A[1, :])       # one row
print(A[:, 2])       # one column

print(A[1:3, 2:4])   # subarray (what size?)

print(A[4,4])
print(A[:,5])
```

**Predict before running:** What shape and value(s) does each expression
return?

\

Key distinction: `A[i, j]` selects a cell; `A[i, :]` and `A[:, j]`
select a row or column.

# Simulation 1: Virus Spreading

## A population on a 2D grid

Divide a region into an $L\times L$ lattice. Each site contains **either
zero or one person**.

Let $n_{ij}\in\{0,1\}$ indicate whether site $(i,j)$ is occupied:

$$n_{ij}\sim\operatorname{Bernoulli}(p),$$

where $p$ is the **probability a site is occupied** (and the average
number of people per site).

\

``` {.python .numberLines}
import numpy as np

rng = np.random.default_rng(4)

L = 50
p = 0.65
occupied = rng.random((L, L)) < p
```

## A spatial SIR model

::: columns
::: {.column width="45%"}
Use one integer-valued array `state` to describe each site:

-   `0`: empty
-   `1`: susceptible
-   `2`: infected
-   `3`: recovered
-   `4`: vaccinated

An occupied site holds **exactly one person**, whose state changes over
time; an empty site remains empty.
:::

::: {.column width="45%"}
![Source: generated with
ChatGPT](./pics/virus_population.png){width="70%"}
:::
:::

Initially, vaccinate a fraction $v$ of occupied sites, then seed a few
infections among the remaining susceptible people.

## A spatial SIR model

**Control parameter:**

-   Fraction of population $v$ that is vaccinated

\

**Other parameters:** (not in our control)

-   Population density $p$
-   Number of initial infections (5)
-   Infection probability $\beta$
-   Recovery probability $\gamma$

## Transmission and recovery

Let $k_{ij}$ count infected people on the **four nearest-neighbor
sites**. Since there is at most one person per site,
$0\leq k_{ij}\leq4$.

A susceptible person becomes infected with probability

$$P_{ij}=1-(1-\beta)^{k_{ij}},$$

where $\beta$ is the transmission probability per infected neighbor per
step.

Each infected person recovers with probability $\gamma$ per step.

\

**Update rules:**

-   susceptible $\to$ infected with probability $P_{ij}$
-   infected $\to$ recovered with probability $\gamma$
-   Empty, recovered, and vaccinated sites do not change.

\

Determine all transitions from the **old state**, then apply them
simultaneously.

## Review: Count nearest neighbors with slicing

How can we add the four neighboring infected counts without writing a
loop over every grid cell?

``` {.python .numberLines}
infected = (state == 2)         # Boolean array

k = np.zeros_like(state, dtype=int)

k[1:, :]  += infected[:-1, :]  # neighbor above
k[:-1, :] += infected[1:, :]   # neighbor below
k[:, 1:]  += infected[:, :-1]  # neighbor to the left
k[:, :-1] += infected[:, 1:]   # neighbor to the right
```

\

**Check:** Why do the two slices on each line have the same shape?

**Boundary condition:** No wraparound at the edges.

## Boolean masks and vectorized updates

Rather than updating sites one at a time, use **Boolean masks** and
vectorized operations.

``` {.python .numberLines}
susceptible = (state == 1)
infected = (state == 2)
P = 1 - (1 - beta)**k

new_infections = susceptible & (rng.random(state.shape) < P)
new_recoveries = infected & (rng.random(state.shape) < gamma)

state[new_infections] = 2
state[new_recoveries] = 3
```

## Demo: Simulate an outbreak

Hold $\beta$ and $\gamma$ fixed. Compare different site-occupancy
probabilities $p$ and vaccination fractions $v$.

-   Plot infected people versus time:
    $I_{\mathrm{tot}}(t)=\sum_{ij} I_{ij}(t)$.
-   Animate the spatial distribution of infection.
-   Compare the final number ever infected.

\

**Predict:** What happens when we increase $p$? What happens when we
increase $v$?

## Optimization 1: Vaccination costs versus infection costs

Vaccination prevents infections, but vaccination itself has a cost.

Let $C(v)$ be the cumulative number of people infected over the
outbreak, **including initial infections**. Let $N$ be the total
population and $N_V(v)$ the number vaccinated.

Define the objective

$$\boxed{J_{\mathrm{virus}}(v)=\mathbb{E}\!\left[\frac{C(v)}{N}\right]+\lambda_v\,\mathbb{E}\!\left[\frac{N_V(v)}{N}\right].}$$

-   First term: expected infected fraction
-   Second term: vaccination cost, weighted by $\lambda_v$
-   Choose $0\leq v\leq1$ to **minimize** $J_{\mathrm{virus}}(v)$

\

Why not just minimize the number of infections?

## Optimize by scanning parameters

Evaluate the objective for a set of candidate vaccination fractions.

``` {.python .numberLines}
v_values = np.linspace(0, 1, 21)
J_values = np.array([virus_objective(v)[2] for v in v_values])

best_index = np.argmin(J_values)
```

Here `virus_objective(v)` estimates the expected cost by averaging over
many independent outbreaks.

\

**Discussion:** How would the optimum change if vaccination became more
expensive ($\lambda_v$ increased)?

# Simulation 2: Forest Fire and Percolation

## A simpler spreading model

Represent a forest by a Boolean occupancy array:

$$A_{ij}=\begin{cases}1 & \text{tree},\\0 & \text{empty} .\end{cases}$$

Each site independently contains a tree with probability $p$.

``` {.python .numberLines}
p = 0.65
forest = rng.random((L, L)) < p
```

Fire begins on occupied sites along the **left edge**.

At every step, burning trees ignite **all** adjacent unburned trees, and
then become burned.

**Important:** Fire spread is deterministic once the forest is
generated.

## Demo: Vary the tree density

Try $p=0.45$, $p=0.59$, and $p=0.75$.

-   Does the fire reach the **right edge**?
-   How many trees burn?
-   Do repeated random forests give the same outcome?

\

Even though fire propagation is deterministic, the **forest itself is
random**.

## Percolation: A connectivity transition

A **connected cluster** consists of trees joined by nearest-neighbor
steps.

In this model, the trees that eventually burn are the connected clusters
that are **connected to the ignition boundary**.

Therefore,

$$\Pr(\text{fire reaches right edge})
=\Pr(\text{a left-right connected cluster exists}).$$

For independent site percolation on a large square lattice,

$$\boxed{p_c\simeq0.5927.}$$

Below $p_c$, spanning is unlikely on a large lattice; above $p_c$, it is
likely.

![](./pics/fire_percolation.png){width="60%"}

## Estimating the spanning probability

A single random forest answers **yes or no**. Many trials estimate a
probability.

$$P_{\text{span}}(p)\approx\frac{\text{number of trials that span}}{\text{total number of trials}}.$$

``` {.python .numberLines}
p_values = np.linspace(0.35, 0.85, 21)
probabilities = [
    np.mean([fire_spans(p) for _ in range(100)])
    for p in p_values
]
```

\

**Question:** Why is the transition smooth for a finite grid rather than
perfectly sharp at $p_c$?

## Optimization 2: Preventing forest fires

Suppose we start with a dense forest of tree density $p_0$.

Remove a fraction $q$ of the existing trees at random. The remaining
density is

$$p=p_0(1-q).$$

Removing trees costs money; letting the fire reach the opposite side
causes damage.

Define

$$\boxed{J_{\text{fire}}(q)=P_{\text{span}}\big(p_0(1-q)\big)+\lambda_f q.}$$

-   First term: risk of a crossing fire
-   Second term: cost of clearing a fraction $q$ of trees
-   Choose $0\leq q\leq1$ to **minimize** $J_{\text{fire}}(q)$

## Predict the optimal removal fraction $q$

From the percolation threshold,

$$p_0(1-q_c)=p_c.$$

Thus,

$$\boxed{q_c=1-\frac{p_c}{p_0}.}$$

For $p_0=0.80$,

$$q_c\approx1-\frac{0.5927}{0.80}\approx0.26.$$

\

**Prediction:** Removing roughly $26\%$ of the existing trees moves the
remaining forest through the percolation threshold.

Does that guarantee the minimum of $J_{\text{fire}}$ occurs at $q_c$?
**No**---the answer depends on $\lambda_f$ and finite-size effects.

## Demo: Minimize fire risk and removal cost

For each value of $q$:

1.  Generate many forests with density $p=p_0(1-q)$.
2.  Compute the fraction of fires that reach the far edge.
3.  Add the removal cost $\lambda_f q$.
4.  Plot $J_{\text{fire}}(q)$ and locate its minimum.

``` {.python .numberLines}
q_values = np.linspace(0, 0.6, 31)
J_values = np.array([fire_objective(q) for q in q_values])
q_best = q_values[np.argmin(J_values)]
```

# Review

## Review

**Virus outbreak**

-   Random population: $n_{ij}\sim\operatorname{Bernoulli}(p)$
-   Stochastic infection and recovery
-   Choose vaccination fraction $v$ to balance health and intervention
    costs

\

**Forest fire**

-   Random occupancy: $A_{ij}\sim\operatorname{Bernoulli}(p)$
-   Deterministic spreading on a fixed forest
-   Choose tree-removal fraction $q$ to balance fire risk and clearing
    costs

\

## Looking ahead: How do we find the best parameters?

Today, we evaluated the objective on a grid of candidate parameter
values:

$$J(\theta_1),J(\theta_2),\ldots,J(\theta_m).$$

But what if there are **many parameters**, or each simulation is
expensive?

-   Can we search more efficiently than testing every possibility?
-   How do noise, local minima, and constraints affect optimization?
-   How can we use derivatives to guide the search?

\

**Next:** Numerical optimization algorithms.