Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,10 @@ The official documentation can be found here https://memilio.readthedocs.io .

### Marimo notebook setup

We recommend using the editor VSCode and installing the Python and Marimo extensions. Only this this set up was

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
We recommend using the editor VSCode and installing the Python and Marimo extensions. Only this this set up was
We recommend using the editor VSCode and installing the Python and Marimo extensions. Only this setup was

tested by us. If you encounter any issues with the installation process, please open an issue on our
[GitHub](https://github.com/SciCompMod/memilio/issues) page.

First, open a new terminal.
Go to a directory for the new project, then create a virtual environment (here "venv") and install marimo:

Expand Down
60 changes: 37 additions & 23 deletions exercises/exercise07.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
import marimo

__generated_with = "0.20.4"
__generated_with = "0.24.0"
app = marimo.App(width="medium")


Expand Down Expand Up @@ -128,7 +128,8 @@ def _(AgeGroup, model, np, num_age_groups):
model.parameters.CriticalPerSevere[AgeGroup(2)] = 0.25 * 2

# Set contact frequency
model.parameters.ContactPatterns.cont_freq_mat[0].baseline = np.ones((num_age_groups, num_age_groups)) * 10
model.parameters.ContactPatterns.cont_freq_mat[0].baseline = np.ones(
(num_age_groups, num_age_groups)) * 10
return


Expand Down Expand Up @@ -161,8 +162,10 @@ def _(AgeGroup, model, num_age_groups, osecir, total_population_per_region):
# The population is equally distributed among the age groups
for group in range(num_age_groups):
# 1% of the population is initially infected, 0.5% Exposed and 0.5% in the pre- or asymptomatic state
model.populations[AgeGroup(group), osecir.InfectionState.Exposed] = 0.005 * total_population_per_region / num_age_groups
model.populations[AgeGroup(group), osecir.InfectionState.InfectedNoSymptoms] = 0.005 * total_population_per_region / num_age_groups
model.populations[AgeGroup(group), osecir.InfectionState.Exposed] = 0.005 * \
total_population_per_region / num_age_groups
model.populations[AgeGroup(group), osecir.InfectionState.InfectedNoSymptoms] = 0.005 * \
total_population_per_region / num_age_groups
# The rest of the population is Susceptible
model.populations.set_difference_from_group_total_AgeGroup(
(AgeGroup(group), osecir.InfectionState.Susceptible), total_population_per_region / num_age_groups)
Expand All @@ -180,7 +183,7 @@ def _(mo):
@app.cell
def _(graph, model, t0):
# Add node with id 0 and copy beforehand initialized model to it
graph.add_node(id=0, model=model, t0=t0)
graph.add_node(id=0, model=model, t0=t0)
return


Expand All @@ -198,7 +201,8 @@ def _(AgeGroup, model, num_age_groups, osecir, total_population_per_region):
for age in range(num_age_groups):
# No infected individuals
model.populations[AgeGroup(age), osecir.InfectionState.Exposed] = 0
model.populations[AgeGroup(age), osecir.InfectionState.InfectedNoSymptoms] = 0
model.populations[AgeGroup(
age), osecir.InfectionState.InfectedNoSymptoms] = 0
# The total population is Susceptible
model.populations.set_difference_from_group_total_AgeGroup(
(AgeGroup(age), osecir.InfectionState.Susceptible), total_population_per_region / num_age_groups)
Expand Down Expand Up @@ -251,7 +255,8 @@ def _(mo):
### Exercise

Set the mobility coefficients for `Dead` individuals to zero:
Hint: For the first age group this can be done via `mobility_coefficients[0 * int(osecir.InfectionState.Dead) + osecir.InfectionState.Dead] = 0`.

Hint: For the first age group this can be done via `mobility_coefficients[0 * (int(osecir.InfectionState.Dead)+1) + int(osecir.InfectionState.Dead)] = 0`.
""")
return

Expand Down Expand Up @@ -298,21 +303,26 @@ def _(mo):


@app.cell
def _(dt_exchange, graph, osecir, t0):
# Create graph simulation and advance until tmax
sim = osecir.MobilitySimulation(graph, t0, dt=dt_exchange)
def _(dt_exchange, graph, osecir, t0, tmax):
import os
# Silence C++ log output; it can deadlock marimo's output pipe on Windows
devnull = os.open(os.devnull, os.O_WRONLY)
saved_stdout, saved_stderr = os.dup(1), os.dup(2)
os.dup2(devnull, 1)
os.dup2(devnull, 2)
try:
# Create graph simulation and advance until tmax
sim = osecir.MobilitySimulation(graph, t0, dt=dt_exchange)
sim.advance(tmax)
finally:
os.dup2(saved_stdout, 1)
os.dup2(saved_stderr, 2)
os.close(devnull)
os.close(saved_stdout)
os.close(saved_stderr)
return (sim,)


@app.cell
def _(mo):
mo.md(r"""
### Exercise
Please advance the simulation until `tmax`.
""")
return


@app.cell
def _():
# Insert code here

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# Insert code here

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

or maybe even remove this entirely?

Expand Down Expand Up @@ -344,8 +354,10 @@ def _(mo):

@app.cell
def _(osecir, result_region0, result_region1):
result_region0_interpolated = osecir.interpolate_simulation_result(result_region0)
result_region1_interpolated = osecir.interpolate_simulation_result(result_region1)
result_region0_interpolated = osecir.interpolate_simulation_result(
result_region0)
result_region1_interpolated = osecir.interpolate_simulation_result(
result_region1)
return result_region0_interpolated, result_region1_interpolated


Expand All @@ -368,8 +380,10 @@ def _(osecir, plt, result_region0_interpolated, result_region1_interpolated):
# Plot the number of non-symptomatically infected for both regions
fig, ax = plt.subplots()
time = result_array_0[0, :]
InfectedNoSymptoms_0 = result_array_0[1 + int(osecir.InfectionState.InfectedNoSymptoms), :] + result_array_0[1 + int(osecir.InfectionState.InfectedNoSymptoms) + int(osecir.InfectionState.Dead) + 1, :] + result_array_0[1 + int(osecir.InfectionState.InfectedNoSymptoms) + 2 * (int(osecir.InfectionState.Dead) + 1), :]
InfectedNoSymptoms_1 = result_array_1[1 + int(osecir.InfectionState.InfectedNoSymptoms), :] + result_array_1[1 + int(osecir.InfectionState.InfectedNoSymptoms) + int(osecir.InfectionState.Dead) + 1, :] + result_array_1[1 + int(osecir.InfectionState.InfectedNoSymptoms) + 2 * (int(osecir.InfectionState.Dead) + 1), :]
InfectedNoSymptoms_0 = result_array_0[1 + int(osecir.InfectionState.InfectedNoSymptoms), :] + result_array_0[1 + int(osecir.InfectionState.InfectedNoSymptoms) + int(
osecir.InfectionState.Dead) + 1, :] + result_array_0[1 + int(osecir.InfectionState.InfectedNoSymptoms) + 2 * (int(osecir.InfectionState.Dead) + 1), :]
InfectedNoSymptoms_1 = result_array_1[1 + int(osecir.InfectionState.InfectedNoSymptoms), :] + result_array_1[1 + int(osecir.InfectionState.InfectedNoSymptoms) + int(
osecir.InfectionState.Dead) + 1, :] + result_array_1[1 + int(osecir.InfectionState.InfectedNoSymptoms) + 2 * (int(osecir.InfectionState.Dead) + 1), :]
ax.plot(time, InfectedNoSymptoms_0, label='Infected No Symptoms Region 1')
ax.plot(time, InfectedNoSymptoms_1, label='Infected No Symptoms Region 2')
ax.set_xlabel('Time [days]')
Expand Down
98 changes: 60 additions & 38 deletions tutorial07.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
import marimo

__generated_with = "0.19.11"
__generated_with = "0.24.0"
app = marimo.App(width="medium")


Expand Down Expand Up @@ -49,7 +49,7 @@ def _():
@app.cell
def _(mo):
mo.md(r"""
We set the simulation start time `t0`, the end time `tmax` and the initial step size `dt` as:
We set the simulation start time `t0` and the end time `tmax`.:
""")
return

Expand All @@ -58,7 +58,6 @@ def _(mo):
def _():
t0 = 0
tmax = 100
dt = 0.1
return t0, tmax


Expand Down Expand Up @@ -104,20 +103,20 @@ def _(mo):

@app.cell
def _(AgeGroup, model, np, num_age_groups):
for ag in range(num_age_groups):
for i in range(num_age_groups):
# Set infection state stay times (in days)
model.parameters.TimeExposed[AgeGroup(ag)] = 3.2
model.parameters.TimeInfectedNoSymptoms[AgeGroup(ag)] = 2.
model.parameters.TimeInfectedSymptoms[AgeGroup(ag)] = 6.
model.parameters.TimeInfectedSevere[AgeGroup(ag)] = 12.
model.parameters.TimeInfectedCritical[AgeGroup(ag)] = 8.
model.parameters.TimeExposed[AgeGroup(i)] = 3.2
model.parameters.TimeInfectedNoSymptoms[AgeGroup(i)] = 2.
model.parameters.TimeInfectedSymptoms[AgeGroup(i)] = 6.
model.parameters.TimeInfectedSevere[AgeGroup(i)] = 12.
model.parameters.TimeInfectedCritical[AgeGroup(i)] = 8.

# Set infection state transition probabilities
model.parameters.RelativeTransmissionNoSymptoms[AgeGroup(ag)] = 0.67
model.parameters.TransmissionProbabilityOnContact[AgeGroup(ag)] = 0.1
model.parameters.RecoveredPerInfectedNoSymptoms[AgeGroup(ag)] = 0.2
model.parameters.RiskOfInfectionFromSymptomatic[AgeGroup(ag)] = 0.25
model.parameters.DeathsPerCritical[AgeGroup(ag)] = 0.3
model.parameters.RelativeTransmissionNoSymptoms[AgeGroup(i)] = 0.67
model.parameters.TransmissionProbabilityOnContact[AgeGroup(i)] = 0.1
model.parameters.RecoveredPerInfectedNoSymptoms[AgeGroup(i)] = 0.2
model.parameters.RiskOfInfectionFromSymptomatic[AgeGroup(i)] = 0.25
model.parameters.DeathsPerCritical[AgeGroup(i)] = 0.3

# The groups have an increasing risk of severe and critical infections
model.parameters.SeverePerInfectedSymptoms[AgeGroup(0)] = 0.2
Expand All @@ -128,7 +127,8 @@ def _(AgeGroup, model, np, num_age_groups):
model.parameters.CriticalPerSevere[AgeGroup(2)] = 0.25 * 2

# Set contact frequency
model.parameters.ContactPatterns.cont_freq_mat[0].baseline = np.ones((num_age_groups, num_age_groups)) * 10
model.parameters.ContactPatterns.cont_freq_mat[0].baseline = np.ones(
(num_age_groups, num_age_groups)) * 10
return


Expand Down Expand Up @@ -159,13 +159,15 @@ def _(mo):
@app.cell
def _(AgeGroup, model, num_age_groups, osecir, total_population_per_region):
# The population is equally distributed among the age groups
for group in range(num_age_groups):
for j in range(num_age_groups):
# 1% of the population is initially infected, 0.5% Exposed and 0.5% in the pre- or asymptomatic state
model.populations[AgeGroup(group), osecir.InfectionState.Exposed] = 0.005 * total_population_per_region / num_age_groups
model.populations[AgeGroup(group), osecir.InfectionState.InfectedNoSymptoms] = 0.005 * total_population_per_region / num_age_groups
model.populations[AgeGroup(j), osecir.InfectionState.Exposed] = 0.005 * \
total_population_per_region / num_age_groups
model.populations[AgeGroup(j), osecir.InfectionState.InfectedNoSymptoms] = 0.005 * \
total_population_per_region / num_age_groups
# The rest of the population is Susceptible
model.populations.set_difference_from_group_total_AgeGroup(
(AgeGroup(group), osecir.InfectionState.Susceptible), total_population_per_region / num_age_groups)
(AgeGroup(j), osecir.InfectionState.Susceptible), total_population_per_region / num_age_groups)
return


Expand All @@ -180,7 +182,7 @@ def _(mo):
@app.cell
def _(graph, model, t0):
# Add node with id 0 and copy beforehand initialized model to it
graph.add_node(id=0, model=model, t0=t0)
graph.add_node(id=0, model=model, t0=t0)
return


Expand All @@ -195,13 +197,14 @@ def _(mo):
@app.cell
def _(AgeGroup, model, num_age_groups, osecir, total_population_per_region):
# The population is equally distributed among the age groups
for age in range(num_age_groups):
for k in range(num_age_groups):
# No infected individuals
model.populations[AgeGroup(age), osecir.InfectionState.Exposed] = 0
model.populations[AgeGroup(age), osecir.InfectionState.InfectedNoSymptoms] = 0
model.populations[AgeGroup(k), osecir.InfectionState.Exposed] = 0
model.populations[AgeGroup(
k), osecir.InfectionState.InfectedNoSymptoms] = 0
# The total population is Susceptible
model.populations.set_difference_from_group_total_AgeGroup(
(AgeGroup(age), osecir.InfectionState.Susceptible), total_population_per_region / num_age_groups)
(AgeGroup(k), osecir.InfectionState.Susceptible), total_population_per_region / num_age_groups)
return


Expand All @@ -216,7 +219,7 @@ def _(mo):
@app.cell
def _(graph, model, t0):
# Add node with id 0 and copy beforehand initialized model to it
graph.add_node(id=1, model=model, t0=t0)
graph.add_node(id=1, model=model, t0=t0)
return


Expand All @@ -231,35 +234,50 @@ def _(mo):


@app.cell
def _(graph, mio, model, np, osecir):
def _(graph, mio, model, np, num_age_groups, osecir):
# One coefficient per (age group x compartment)
mobility_coefficients = 0.1 * np.ones(model.populations.numel())
# Dead individuals do not commute
mobility_coefficients[osecir.InfectionState.Dead] = 0
for l in range(num_age_groups):
mobility_coefficients[l * (int(osecir.InfectionState.Dead)+1) +
int(osecir.InfectionState.Dead)] = 0
mobility_params = mio.MobilityParameters(mobility_coefficients)
# Add two edges to graph
graph.add_edge(0, 1, mobility_params)
graph.add_edge(1, 0, mobility_params)
# Individuals are exchanged every half day
dt_exchange = 0.5
return
return (dt_exchange,)


@app.cell
def _(mo):
mo.md(r"""
## Model simulation

We now have finished initializing the metapopulation model. The graph-based simulation is created and advanced until `tmax` via:
""")
return


@app.cell
def _(graph, osecir, t0, tmax):
# Create graph simulation and advance until tmax
sim = osecir.MobilitySimulation(graph, t0, dt=dt_exchange)
sim.advance(tmax)
def _(dt_exchange, graph, osecir, t0, tmax):
import os
# Silence C++ log output; it can deadlock marimo's output pipe on Windows
devnull = os.open(os.devnull, os.O_WRONLY)
saved_stdout, saved_stderr = os.dup(1), os.dup(2)
os.dup2(devnull, 1)
os.dup2(devnull, 2)
try:
# Create graph simulation and advance until tmax
sim = osecir.MobilitySimulation(graph, t0, dt=dt_exchange)
sim.advance(tmax)
finally:
os.dup2(saved_stdout, 1)
os.dup2(saved_stderr, 2)
os.close(devnull)
os.close(saved_stdout)
os.close(saved_stderr)
return (sim,)


Expand Down Expand Up @@ -288,16 +306,18 @@ def _(mo):

@app.cell
def _(osecir, result_region0, result_region1):
result_region0_interpolated = osecir.interpolate_simulation_result(result_region0)
result_region1_interpolated = osecir.interpolate_simulation_result(result_region1)
result_region0_interpolated = osecir.interpolate_simulation_result(
result_region0)
result_region1_interpolated = osecir.interpolate_simulation_result(
result_region1)
return result_region0_interpolated, result_region1_interpolated


@app.cell
def _(mo):
mo.md(r"""
## Visualization of model output

Finally, we can compare the trajectories of all infection states for both regions. In the following, we plot the number of `InfectedNoSymptoms` aggregated over all age groups for both regions:
""")
return
Expand All @@ -312,8 +332,10 @@ def _(osecir, plt, result_region0_interpolated, result_region1_interpolated):
# Plot the number of non-symptomatically infected for both regions
fig, ax = plt.subplots()
time = result_array_0[0, :]
InfectedNoSymptoms_0 = result_array_0[1 + int(osecir.InfectionState.InfectedNoSymptoms), :] + result_array_0[1 + int(osecir.InfectionState.InfectedNoSymptoms) + int(osecir.InfectionState.Dead) + 1, :] + result_array_0[1 + int(osecir.InfectionState.InfectedNoSymptoms) + 2 * (int(osecir.InfectionState.Dead) + 1), :]
InfectedNoSymptoms_1 = result_array_1[1 + int(osecir.InfectionState.InfectedNoSymptoms), :] + result_array_1[1 + int(osecir.InfectionState.InfectedNoSymptoms) + int(osecir.InfectionState.Dead) + 1, :] + result_array_1[1 + int(osecir.InfectionState.InfectedNoSymptoms) + 2 * (int(osecir.InfectionState.Dead) + 1), :]
InfectedNoSymptoms_0 = result_array_0[1 + int(osecir.InfectionState.InfectedNoSymptoms), :] + result_array_0[1 + int(osecir.InfectionState.InfectedNoSymptoms) + int(
osecir.InfectionState.Dead) + 1, :] + result_array_0[1 + int(osecir.InfectionState.InfectedNoSymptoms) + 2 * (int(osecir.InfectionState.Dead) + 1), :]
InfectedNoSymptoms_1 = result_array_1[1 + int(osecir.InfectionState.InfectedNoSymptoms), :] + result_array_1[1 + int(osecir.InfectionState.InfectedNoSymptoms) + int(
osecir.InfectionState.Dead) + 1, :] + result_array_1[1 + int(osecir.InfectionState.InfectedNoSymptoms) + 2 * (int(osecir.InfectionState.Dead) + 1), :]
ax.plot(time, InfectedNoSymptoms_0, label='Infected No Symptoms Region 1')
ax.plot(time, InfectedNoSymptoms_1, label='Infected No Symptoms Region 2')
ax.set_xlabel('Time [days]')
Expand Down
Loading