From 8a07acd46db507ff90d976b85988bdb92adeb37b Mon Sep 17 00:00:00 2001 From: gvhaelew <162558233+gvhaelew@users.noreply.github.com> Date: Wed, 12 Feb 2025 16:54:21 +0100 Subject: [PATCH] Add files via upload --- exercises/ode_model_catalyst_intro.jl | 112 ++++---- exercises/ssa_model_bike_sharing.jl | 351 ++++++++++++++++++++++++++ 2 files changed, 417 insertions(+), 46 deletions(-) create mode 100644 exercises/ssa_model_bike_sharing.jl diff --git a/exercises/ode_model_catalyst_intro.jl b/exercises/ode_model_catalyst_intro.jl index 008548e3..ddb0cc1d 100644 --- a/exercises/ode_model_catalyst_intro.jl +++ b/exercises/ode_model_catalyst_intro.jl @@ -168,8 +168,8 @@ The substrates and the products may contain one or more reactants, separated by # ╔═╡ bb5578b8-40d4-48c9-be87-cfaeadb8c5f9 md""" -!!! hint "Hint" -The Greek letters can be visualized by typing a **backslash** followed by the **name of the Greek letter** and then the **TAB** key. For example `\alpha` followed by the TAB key results in in a list where you can choose `α`. +!!! tip "Tip" + The Greek letters can be visualized by typing a **backslash** followed by the **name of the Greek letter** and then the **TAB** key. For example `\alpha` followed by the TAB key results in in a list where you can choose `α`. """ # ╔═╡ af6b169d-55c3-4384-9f57-ad629e02bad6 @@ -196,7 +196,7 @@ parameters(infection_model) # ╔═╡ 2bd02577-7321-46b5-9f42-1484ed86d1d3 md""" !!! important "Important" -You can also get the different species and parameters using `@unpack` followed by comma separated species and/or parameter names followed by the equal sign and the name of the reaction network model. For example: + You can also get the different species and parameters using `@unpack` followed by comma separated species and/or parameter names followed by the equal sign and the name of the reaction network model. For example: """ # ╔═╡ 64455b26-4374-435f-b71b-8b9600cab263 @@ -212,7 +212,7 @@ osys = convert(ODESystem, infection_model) # ╔═╡ eeb25e50-165b-4e93-8397-5b09fe8e7242 md""" -Note that the model equations are essencially: +Note that the model equations are essentially: $$\cfrac{dS(t)}{dt} = -\alpha \beta S(t) I(t)$$ $$\cfrac{dI(t)}{dt} = \alpha \beta S(t) I(t) - r I(t)$$ @@ -248,18 +248,18 @@ parameters(osys) md""" ### Simulating the system as an ODE-problem -We first need to load the Differential and Plot package, which is required for simulating the system and plotting the results. +We first need to load the `DifferentialEquations` and `Plots` packages, which are required for simulating the system and plotting the results. """ # ╔═╡ 3197244f-655b-4dca-80f3-794b30722551 md""" -Now we wish to simulate our model. To do this, we need to provide some the following information: +Now we wish to simulate our model. To do this, we need to provide the following information: - Initial conditions for the state variables $S$, $I$, $D$ and $R$. - The parameter values for $\alpha$, $\beta$, $r$ and $m$. - The timespan, which is the timeframe over which we wish to run the simulation. -Assume in this example that there are $10\,000\,000$ people in the country, and that initially $1\,000$ person are infected. Hence, $I_0 = 1\,000$, $S_0 = 10\,000\,000-I_0 = 9\,999\,000$, $D_0 = 0$ and $R_0 = 0$.\ +Assume in this example that there are $10\,000\,000$ people in the country, and that initially $1\,000$ persons are infected. Hence, $I_0 = 1\,000$, $S_0 = 10\,000\,000-I_0 = 9\,999\,000$, $D_0 = 0$ and $R_0 = 0$.\ Furthermore, we take the following values for the parameters: $\alpha = 0.08\;person/contact$, $\beta = 10^{-6}\;contact/(person^2\,day)$, $r = 0.2\;day^{-1}$ (i.e. a person is contagious for an average of $5\;days$) and $m=0.4$. The following table summarizes the above values: |Initial conditions |Parameters | @@ -279,7 +279,7 @@ md""" # ╔═╡ 1ba859fa-46a5-434f-a99c-e710ba85caf8 md""" -The initial conditions are given as a *Vector*. This is a type which collects several different values. To declare a vector, the values are specific within brackets, `[]`, and separated by `,`. Since we have four species, the vector holds four elements. E.g., we set the value of $I$ using the `:I => 1` syntax. Here, we first denote the name of the species (with a colon `:` pre-appended), next follows a `=>` and then the value of `I`.\ +The initial conditions are given as a *Vector*. This is a type which collects several different values. To declare a vector, the values are specified within brackets, `[]`, and separated by `,`. Since we have four species, the vector holds four elements. E.g., we set the value of $I$ using the `:I => 1` syntax. Here, we first denote the name of the species (with a colon `:` pre-appended), next follows a `=>` and then the value of `I`.\ The vector holding the initial conditions for $S$, $I$, $D$ and $R$ can be created in the following way: """ @@ -288,7 +288,7 @@ u0 = [:S => 9_999_000.0, :I => 1_000.0, :D => 0.0, :R => 0.0] # ╔═╡ 95c1c0ea-51d0-47ee-8948-a5b87c42a70d md" -Note that the order of the vector elements doesn't matter, because the initial values of each of the species is indicated using its variable name. +Note that the order of the vector elements doesn't matter here, since the initial values of each of the species is indicated using its variable name. " # ╔═╡ 57036f49-2f1f-4327-89ed-2c96098a1c22 @@ -298,7 +298,7 @@ md""" # ╔═╡ 56ebb72d-2351-4e67-b268-1f48bbb77cb3 md""" -The timespan sets the time point at which we start the simulation (typically `0.0` is used) and the final time point of the simulation. These are combined into a two-valued Tuple. Tuples are similar to vectors, but are enclosed by `()` and not `[]`. Again, we will let both time points be decimal valued. +The timespan sets the time point at which we start the simulation (typically `0.0` is used) and the final time point of the simulation. These are combined into a two-valued *Tuple*. Tuples are similar to vectors, but are enclosed by `()` and not `[]`. Again, we will let both time points be decimal valued. """ # ╔═╡ a25d5925-a254-488b-b782-d29cff4470a2 @@ -311,7 +311,7 @@ md""" # ╔═╡ a235c7ce-f14a-4c7c-86a1-08aa5f2d9c85 md""" -Similarly, the parameters values are also given as a vector. We have four parameters, hence, the parameter vector will also contain four elements. We use a similar notation for setting the parameter values as the initial condition (first the colon, then the parameter name, then an arrow, then the value). +Similarly, the parameter values are also given as a vector. We have four parameters, hence, the parameter vector will also contain four elements. We use a similar notation for setting the parameter values as the initial condition (first the colon, then the parameter name, then an arrow, then the value). """ # ╔═╡ 9a9440fa-d8a3-44bc-8037-4bf1f8af40b0 @@ -358,7 +358,7 @@ md""" # ╔═╡ 0af3c166-46b6-455d-af6b-a72c4d2a5ce4 md""" -Finally, we can plot the solution through the plot function. +Finally, we can plot the solution through the `plot` function. """ # ╔═╡ 513c037b-c54c-47fa-b97a-06f69a983386 @@ -374,7 +374,7 @@ plot(osol, idxs=[:S, :I]) # brackets [ ] # ╔═╡ f3cf01a1-3c65-4a37-bf9a-cd2233cb470b md""" -If you want a fase plot of for example just $I$ versus $S$, you can specify this with the option `idxs=(:S, :I)` (*notice the parentheses*) in the plot function. You can indicate the $S$ and $I$ axes with the additional options `xlab="S"` and `ylab="I"`. +If you want a phase plot of, for example, just $I$ versus $S$, you can specify this with the option `idxs=(:S, :I)` (*notice the parentheses*) in the plot function. You can indicate the $S$ and $I$ axes with the additional options `xlab="S"` and `ylab="I"`. """ # ╔═╡ ff8f4c23-5695-48aa-9965-02677103f2c9 @@ -437,7 +437,7 @@ You may have noticed that while using the Pluto notebooks, when you change the v md""" ### Example 1 - Influence of $r$ -Influence of the duration of infection $1/r$ for average infection periods of between $10$, days and $1$ day contagious ($r$ between $0.1$ and $1.0$, step $0.1$, default value $0.1$). +Influence of the duration of infection $1/r$ for average infection periods of between $10$ days and $1$ day of being contagious ($r$ between $0.1$ and $1.0$, step $0.1$, default value $0.1$). """ # ╔═╡ b6baafc2-6d5e-43c3-8ef9-845961cdd20b @@ -450,14 +450,14 @@ We will create a slider for the $r$-values between $0.1$ and $1.0$, stepsize $0. # ╔═╡ ae38c663-0ee4-409e-bfca-5f13ed88b67d md""" -We will create een new parameter value vector, ODE problem and solution object by putting `1` at the end of the corresponding variable names. In that way, the previous simulation results will be unaffected! The model, the initial conditions and the timespan are identical as before. In there we also use the variable `r` coupled to the slider. +We will create a new parameter value vector, ODE problem and solution object by putting `1` at the end of the corresponding variable names. In that way, the previous simulation results will be unaffected! The model, the initial conditions and the timespan are identical as before. In there we also use the variable `r` coupled to the slider. """ # ╔═╡ d44da6c6-c93d-4c61-8125-9eee464c897e params1 = [:α => 0.08, :β => 1.0e-6, :r => r, :m => 0.4] # ╔═╡ 00a72697-d36a-41cc-9eec-8e821829ce0e -# put semi-colon at end of instruction to avoid seeing its output. +# type a semi-colon at end of an instruction to avoid seeing its output oprob1 = ODEProblem(infection_model, u0, tspan, params1); # ╔═╡ cf39b4cf-9cd0-4755-80db-ca4aea7c1084 @@ -476,28 +476,36 @@ md""" Try to interpret the results yourself. Ask yourself the following questions: -1. What are the trends in the results obtained? +!!! question + 1. What are the trends in the results obtained? """ -# ╔═╡ ff9ce0c5-931f-44f4-aa0d-4aead0284160 -md"- Answer: missing" +# ╔═╡ 4242ec68-b36d-4f7f-9001-2c90c66d5b9b +md""" +- Answer: +""" # ╔═╡ f5f98168-119c-4b49-a8e3-a3cc775faeb0 -md"""2. How can this be explained from the model structure?""" +md""" +!!! question + 2. How can this be explained from the model structure? +""" -# ╔═╡ a308de53-80a2-4b37-a752-cb2c383cf7a4 -md"- Answer: missing" +# ╔═╡ 30b0e2b3-2b8d-43a5-9587-9fcb0b89b258 +md""" +- Answer: +""" # ╔═╡ ade413c2-d7d9-4250-8490-75534900a389 md""" ### Example 2 - Discrete Event -Suppose that regulations are such that on day 14, people need to reduce their contacts by 50%. Hence, this means that the parameter value $\beta$ needs to be divided by a factor of 2 at timepoint 14. In order to realize that we need to now the order of the parameters in the model because we will need to address the value of $\beta$ by means of an index. +Suppose that regulations are such that on day 14, people need to reduce their contacts by 50%. Hence, this means that the parameter value $\beta$ needs to be divided by a factor of 2 at timepoint 14. """ # ╔═╡ 58730ac6-d83b-420d-a000-2f50545f0d39 md""" -We need to state that the parameter $\beta$ needs to be reduce by $50\%$ at time $t=14\,days$. We put this in a condition named `condition2`. +We need to state that the parameter $\beta$ needs to be reduced by $50\%$ at time $t=14$. We include this condition in a variable named `condition2` in the following way: """ # ╔═╡ 7411474c-fac7-4b7c-8ded-4c2df5956fb0 @@ -505,7 +513,7 @@ condition2 = [14.0] => [infection_model.β ~ infection_model.β/2] # ╔═╡ e6949052-6d12-4414-ac67-cb9435b46290 md""" -The discrete time event needs to be included in our model. +The discrete time event needs to be now included in our model. """ # ╔═╡ 8752373f-a602-437b-9b3a-2602c8babf87 @@ -513,7 +521,7 @@ The discrete time event needs to be included in our model. # ╔═╡ ee596a6c-29bf-4d18-b687-be943c128aa7 md""" -After that, we need to *complete* our *reaction network model*. +After that, we need to *complete* our *reaction network model*, so that the model can be simulated. """ # ╔═╡ 40f849cb-20f9-4bbe-806e-512abd6f3210 @@ -555,19 +563,25 @@ osol2.u[end] md""" Try to interpret the results yourself. Ask yourself the following questions: -1. What are the trends in the results obtained? +!!! question + 1. What are the trends in the results obtained? """ -# ╔═╡ a736a17b-ee8e-491b-840a-6319619a8dab -md"- Answer: missing" +# ╔═╡ 6eb368f5-8d43-4fb8-b6eb-c783675d1dd9 +md""" +- Answer: +""" # ╔═╡ 2ba54805-e901-4281-9b07-cc7137b8c809 md""" -2. How much less casualties are there compared to not altering the contact rate? +!!! question + 2. How much less casualties are there compared to not altering the contact rate? """ -# ╔═╡ 20d56d2f-2532-4286-88cc-b42ccb22d246 -md"- Answer: missing" +# ╔═╡ fec6e165-bd1a-4e57-b7f3-ae712e03276c +md""" +- Answer: +""" # ╔═╡ cb8a6f77-f08c-4fc8-9445-bd1c17521fcc md""" @@ -578,7 +592,7 @@ Suppose that when the number of infected individuals reaches $1\,000\,000$, then # ╔═╡ 5005a09d-a844-4f9d-a058-0d93583d5bab md""" -Normally in a continuous event the value of one or more species can be changed when a certain condition is met. In our specific case we want the change in the species happening only once! So, if you want that the continuous event: ''when $I$ reaches $10^6$ then $999000$ is subtrated from $I$'' happens only once, then we need to include a ficticious new *species* in our *reaction network model*. We will call this ficticious *species* `pwc` (a short for _**p**roceed **w**ith **c**ondition_) and we set it default to `true`. +Normally in a continuous event the value of one or more species can be changed when a certain condition is met. In our specific case we want the change in the species happening only once! So, if you want that the continuous event "when $I$ reaches $10^6$ then $999000$ is subtracted from $I$" happens only once, then we need to include a ficticious new *species* in our *reaction network model*. We will call this ficticious species `pwc` (a short for _**p**roceed **w**ith **c**ondition_) and we set its default value to `true`. """ # ╔═╡ f77d2d75-468f-4746-b81e-0bbbf33fc8d7 @@ -638,7 +652,7 @@ Now we can plot the results. """ # ╔═╡ 02173eaa-1178-4383-b57a-03e53ce38af8 -plot(osol3) +plot(osol3, idxs=[:S, :I, :D, :R]) # ╔═╡ 6d1aa79b-8614-4a77-a327-f7e6962d1944 md""" @@ -652,19 +666,25 @@ osol3.u[end] md""" Try to interpret the results yourself. Ask yourself the following questions: -1. What are the trends in the results obtained? +!!! question + 1. What are the trends in the results obtained? """ -# ╔═╡ 7872701b-8506-4303-b615-32d3c1afa50d -md"- Answer: missing" +# ╔═╡ 33ff901b-6c7e-4b8c-a163-663d9443a213 +md""" +- Answer: +""" # ╔═╡ 0d6c2233-35c5-4257-9dfd-398704080ba3 md""" -2. How much less casualties are there compared to not putting $999\,000$ individuals into isolation? (Hint: you also need to take into account the casualties in the $999\,000$ individuals that had been put into isolation.) +!!! question + 2. How much less casualties are there compared to not putting $999\,000$ individuals into isolation? (Hint: you also need to take into account the casualties in the $999\,000$ individuals that had been put into isolation.) """ -# ╔═╡ 029a2024-1ec7-4af9-9e73-717f13483468 -md"- Answer: missing" +# ╔═╡ f5d23dba-02a9-4d27-8eb1-ff89d2902cdf +md""" +- Answer: +""" # ╔═╡ Cell order: # ╠═be326550-25ea-4c5f-ac4a-dd74d12bc89a @@ -755,9 +775,9 @@ md"- Answer: missing" # ╠═52901bbf-e47b-4da1-95c4-f0869812398c # ╟─102b4fbc-23b1-46ed-bb72-124eb88517ce # ╟─f6cbafbf-98b4-4286-86d2-fe95821a5ff4 -# ╟─ff9ce0c5-931f-44f4-aa0d-4aead0284160 +# ╟─4242ec68-b36d-4f7f-9001-2c90c66d5b9b # ╟─f5f98168-119c-4b49-a8e3-a3cc775faeb0 -# ╟─a308de53-80a2-4b37-a752-cb2c383cf7a4 +# ╟─30b0e2b3-2b8d-43a5-9587-9fcb0b89b258 # ╟─ade413c2-d7d9-4250-8490-75534900a389 # ╟─58730ac6-d83b-420d-a000-2f50545f0d39 # ╠═7411474c-fac7-4b7c-8ded-4c2df5956fb0 @@ -774,9 +794,9 @@ md"- Answer: missing" # ╟─18151ab1-1e4e-48d8-be70-4fa4c8a43af1 # ╠═7d7ed974-7c9c-4c30-ab2a-e3028ce702dd # ╟─7375edbc-24a4-4300-bdef-2686cb377cfc -# ╟─a736a17b-ee8e-491b-840a-6319619a8dab +# ╟─6eb368f5-8d43-4fb8-b6eb-c783675d1dd9 # ╟─2ba54805-e901-4281-9b07-cc7137b8c809 -# ╟─20d56d2f-2532-4286-88cc-b42ccb22d246 +# ╟─fec6e165-bd1a-4e57-b7f3-ae712e03276c # ╟─cb8a6f77-f08c-4fc8-9445-bd1c17521fcc # ╟─5005a09d-a844-4f9d-a058-0d93583d5bab # ╠═f77d2d75-468f-4746-b81e-0bbbf33fc8d7 @@ -796,6 +816,6 @@ md"- Answer: missing" # ╟─6d1aa79b-8614-4a77-a327-f7e6962d1944 # ╠═d04b8d97-e9d3-4428-a695-bfde4b44a291 # ╟─70fd136e-5ded-4c30-818e-a5de61fbbf86 -# ╟─7872701b-8506-4303-b615-32d3c1afa50d +# ╟─33ff901b-6c7e-4b8c-a163-663d9443a213 # ╟─0d6c2233-35c5-4257-9dfd-398704080ba3 -# ╟─029a2024-1ec7-4af9-9e73-717f13483468 +# ╟─f5d23dba-02a9-4d27-8eb1-ff89d2902cdf diff --git a/exercises/ssa_model_bike_sharing.jl b/exercises/ssa_model_bike_sharing.jl new file mode 100644 index 00000000..25d1754c --- /dev/null +++ b/exercises/ssa_model_bike_sharing.jl @@ -0,0 +1,351 @@ +### A Pluto.jl notebook ### +# v0.20.4 + +using Markdown +using InteractiveUtils + +# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error). +macro bind(def, element) + #! format: off + quote + local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end + #! format: on +end + +# ╔═╡ 309035dd-5653-48a6-a53d-817e743279fa +begin + # add this cell if you want the notebook to use the environment from where the Pluto server is launched + using Pkg + Pkg.activate(".") +end + +# ╔═╡ 6b342f14-e7d5-11ef-1ea0-77ceb0d78f32 +using Markdown + +# ╔═╡ 7bfcc024-7d7f-4c0b-b918-8f7626b10974 +using InteractiveUtils + +# ╔═╡ 71140c81-af29-4857-8020-4f94c8bd64b3 +using Catalyst, DifferentialEquations, Plots, Distributions + +# ╔═╡ 284f5847-9c15-41f3-a595-1e12a22df69f +using PlutoUI; TableOfContents() + +# ╔═╡ 1f975552-b0b8-4830-8dcc-214574d4fc38 +md""" +# Exercise: Modeling a simple Bike Share System +""" + +# ╔═╡ d2f32eab-0b35-4794-9219-5bcbb4c069c5 +md""" +Imagine a bike share system for students traveling between Olin College and Wellesley College, which are about three miles apart in eastern Massachusetts. Suppose the system contains 12 bikes and two bike racks, one at Olin and one at Wellesley, each with the capacity to hold 12 bikes. As students arrive, check out a bike, and ride to the other campus, the number of bikes in each location changes. +Initially there are 10 bikes at Olin and, hence, 2 bikes at Wellesley. For this simple model, we will aslo assume that the changes in the number of bikes at both locations is instantenuously. The rate at which a bike is moved from Olin to Wellesley is denoted as $p_1$ ($min^{-1}$); the rate at which a bike is moved from Wellesley to Olin is denoted as $p_2$ ($min^{-1}$). Both processes are zeroth-order and we want to see the evolution of bikes during $1\,h = 60\,min$. + +This is a discreet and stochastic problem and you need to solve it with SSA. +""" + +# ╔═╡ 016842c9-9479-4061-a27e-9dc006121f23 +md""" +Create a *reaction network object* model for the aforementioned problem in order to simulate the evolution of the number of bikes at Olin ($O$) and Wellesley ($W$) with time. Name it `bike_sharing`. +""" + +# ╔═╡ 6c97bf81-ef32-45a4-aa7c-c8c26ba2d2c3 +# bike_sharing = @reaction_network begin +# missing +# ... +# end +bike_sharing = @reaction_network begin + @species O(t)=10 W(t)=2 + p₁, O => W + p₂, W => O +end + +# ╔═╡ 9c7ab7fb-7380-41a3-85ea-714478ade218 +md""" +Convert the system to a symbolic differential equation model and verify, by analyzing the differential equation, that your model is correctly implemented. +""" + +# ╔═╡ 1536fe23-0f8d-4b86-98d2-076248b35954 +# osys = missing +osys = convert(ODESystem, bike_sharing) + +# ╔═╡ e6e2ff5c-38eb-4ba3-b430-c9031483a0a5 +md""" +Initialize a vector `u0` with the initial conditions: +""" + +# ╔═╡ ab6af765-1cde-4da8-bbc1-a5fab391db54 +# u0 = missing +u0 = [:O => 10, :W => 2] + +# ╔═╡ 378878a0-5c09-4eb0-ac43-1031014ff12a +md""" +Set the timespan for the simulation: +""" + +# ╔═╡ 3ae98e83-7beb-4597-89be-80c813d4349b +tspan = (0.0, 60.0) + +# ╔═╡ 988f79c0-9c7b-4752-a7f2-d4473ad73ce6 +md""" +Create a slider for the variable `p₁` in the range of $0.0$ and $1.0$ with a step of $0.1$. Take a default value of $0.0$. +""" + +# ╔═╡ 0d8f53f8-0a14-4ac6-bd0c-2190d4db0909 +# @bind p₁ missing +@bind p₁ Slider(0.0:0.1:1, default=0.0, show_value=true) + +# ╔═╡ 08d43ac8-a973-4d7b-baf7-4c37e54cfe24 +md" +Initialize vector `parms` with parameter values, `p₁` is the slider value and assign a constant value of `0.3` to `p₂`. +" + +# ╔═╡ e20e4dd8-bdbb-4005-af68-6bf7e4ec130e +# parms = missing +parms = [:p₁=>p₁, :p₂=>0.30] + +# ╔═╡ 238e1120-34af-4d57-8efa-aa80ab28a874 +md""" +Create a DiscreteProblem and store it in `dprob`: +""" + +# ╔═╡ d4c45709-70c9-4ba0-8fb8-6b600473723d +# dprob = missing +dprob = DiscreteProblem(bike_sharing, u0, tspan, parms) + +# ╔═╡ d06fb076-76e4-4248-a940-96804ea68833 +md""" +Create a JumpProblem and store it in `jdprob`. Use the simulation method `Direct()` and an additional option `save_positions=(false, false)`. The latter prohibits to save the values just before and after the jump event (later, when solving the problem we will namely use `saveat=1.0`). +""" + +# ╔═╡ 7644adf4-d992-48b1-b40a-12fdf30f6cb5 +jdprob = JumpProblem(bike_sharing, dprob, Direct(), save_positions=(false, false)); +# https://docs.sciml.ai/JumpProcesses/dev/jump_solve/#JumpProcesses.jl + +# ╔═╡ 59d2d3e1-354b-4444-b8a1-16ad8ea2ba94 +md""" +If we would now solve the problem, you might for example encounter negative values for the number of bikes at Olin and, hence, a value larger than 12 at Wellesley. Of course, this makes no sense. Therefore, we will create the so-called `condition` function that needs to invoke the so-called `affect!` function at each jump event in order to check on $O$ and $W$ and setting them to valid values if necessary. +""" + +# ╔═╡ 53e15d08-7d1d-4bc2-9fd3-4bc6fbb9de84 +md""" +Create the `condition` function. +""" + +# ╔═╡ 1511d269-9706-41b1-b8e2-f85ebcedc2d8 +# function condition(u, t, integrator) +# missing +# end +function condition(u, t, integrator) + true +end + +# ╔═╡ fd9b4521-30d2-4af3-a080-40e9dbf27aa4 +md""" +Create the `affect` function. + +Hints: +- If the number of bikes at Olin is larger than 12 (this implies that the number of bikes at Wellesley is negative), then the number of bikes at Olin should be set to 12 and the number of bikes at Wellesley should be set to 0. +- If the number of bikes at Olin is negative (this implies that the number of bikes at Wellesley is larger than 12), then the number of bikes at Olin should be set to 0 and the number of bikes at Wellesley should be set to 12. +- The number of bikes at Olin is accessed through `integrator.u[1]` and the number of bikes at Wellesley is accessed through `integrator.u[2]`. +""" + +# ╔═╡ 3cac94d6-14ae-42f2-b0b6-f47d87cdb518 +# function affect!(integrator) +# if missing +# missing +# missing +# end +# if missing +# missing +# missing +# end +# end +function affect!(integrator) + if integrator.u[1] > 12 + integrator.u[1] = 12 + integrator.u[2] = 0 + end + if integrator.u[1] < 0 + integrator.u[1] = 0 + integrator.u[2] = 12 + end +end + +# ╔═╡ b2797fd2-b6e1-4bb0-af23-26931fe8be69 +md""" +Create the discrete callback function. Store it in `cb`. Again use the option `save_positions=(false,false)`. +""" + +# ╔═╡ f947c2d9-9123-422d-8972-157717c85b3c +# cb = missing +cb = DiscreteCallback(condition, affect!, save_positions=(false,false)); + +# ╔═╡ 74708270-b1ec-48c7-af32-3b970b92c706 +md""" +Solve the problem and store it in `jdsol`. Use the `SSAStepper()` stepping algorithm, the option `saveat=1.0` and the callback function. +""" + +# ╔═╡ 2b00df5d-994e-47a1-8068-c93ce3f1a618 +jdsol = solve(jdprob, SSAStepper(), saveat=1.0, callback=cb); + +# ╔═╡ 9d06c31e-3525-4889-a1de-3fe02413c7d8 +md""" +Plot the solution. +""" + +# ╔═╡ 9a90f800-3669-4831-b50b-c5405bbb9a03 +plot(jdsol, ylim=(0, 12)) + +# ╔═╡ a554fd16-aa3d-48ca-8de6-5582725c27d8 +md""" +Analyse the results. See what happens when you: +- run the notebook cell with the `solve` function repeatly +- change the value of `p₁` using the slider +""" + +# ╔═╡ d6872046-b5ef-4c2d-a9bb-2418f57f715d +md""" +!!! question "Question" + From what value of `p₁` do you start to get empty bike racks at Olin? +""" + +# ╔═╡ 747e20c4-b06b-4e78-a09a-55053cf42bf4 +md"- Answer: missing" + +# ╔═╡ 92181028-60fc-4830-afba-2380ac91455d +md""" +You can inspect the actual number of bike values at Olin by using `jdsol[:O]`: +""" + +# ╔═╡ f8942b10-773a-4b22-baad-8004fba8bd34 +# missing +jdsol[:O] + +# ╔═╡ 43d41284-053d-4dfe-8d5b-96be70c0495c +md""" +If you want to have a `true` boolean value on positions where the vector value is zero (and `false` on non-zero values), then you would compare `jdsol[:O]` element wise with `0`. In Julia, if you want to do element wise operations with/on vectors, you always need to place a dot (`.`) in front of the operator, like for example `.==`. + +Compare in that way `jdsol[:O]` with `0`: +""" + +# ╔═╡ a999ae2a-7567-41e7-9c0c-e94fad6f5d46 +# missing +jdsol[:O] .== 0 + +# ╔═╡ 9ebb5b44-04d7-4b89-acdb-e40a245703d2 +md""" +Furthermore, if you want the count the number of `true` values in the latter (hence, the zero element values), you can simply use the function `count(...)`. Count the number of zeros: +""" + +# ╔═╡ 049de8d5-b221-452b-b2c4-9bc1e0c17f48 +# missing +count(jdsol[:O] .== 0) + +# ╔═╡ a73a2853-1f48-4179-9771-083794d3f137 +md""" +Using the aforementioned way to count zeros in a vector, we will now count the zeros for a range of $p$ values. Because of the stochastic behaviour of the system, for each $p$ values we will count the zeros for a $1000$ simulations and then storing only the average value. + +To introduce a new value for $p_1$ you can take a deepcopy of the problem and remake the problem like this: +- `jdprob_re = remake(deepcopy(jdprob); p=[:p₁=>p_val])` +and then solving the problem and store it in `jdsol_re`. + +In the layout below, `mean_zero_counts` while contain the final mean values of the averaged numbers of zeros from a 1000 simulations using a specific $p$ value, `zero_counts_p_val` will contain the actual number of zeros for a 1000 simulations using a specific $p$ value. + +Use the layout below to fill in `mean_zero_counts`. +""" + +# ╔═╡ 682e9120-0e1c-4dfa-9ec6-66bb0a3f4374 +# begin +# p_values = 0.0:0.1:1.0 # different p-values +# mean_zero_counts = [] # vector to store the corresponding mean zero values +# for p_val = p_values # p_val will be each of the p_values +# zero_counts_p_val = [] # vector to store the zeros for the 1000 simulations +# for i = 1:1000 # do a 1000 simulations +# # take a deepcopy and remake the problem for the specific p-value +# jdprob_re = missing +# # solve the problem +# jdsol_re = missing +# # append the number of zeros to zero_counts_p_val +# append!(..., ...) +# end +# # append the mean number of zeros to mean_zero_counts +# append!(..., ...) +# end +# end +begin + p_values = 0.0:0.1:1.0 # different p-values + mean_zero_counts = [] # vector to store the corresponding mean zero values + for p_val = p_values # p_val will be each of the p_values + zero_counts_p_val = [] # vector to store the zeros for the 1000 simulations + for i = 1:1000 # do a 1000 simulations + # take a deepcopy and remake the problem for the specific p-value + jdprob_re = remake(deepcopy(jdprob); p=[:p₁=>p_val]) + # solve the problem + jdsol_re = solve(jdprob_re, SSAStepper(), saveat=1.0, callback=cb); + # append the number of zeros to zero_counts_p_val + append!(zero_counts_p_val, count(jdsol_re[:O] .== 0)) + end + # append the mean number of zeros to mean_zero_counts + append!(mean_zero_counts, mean(zero_counts_p_val)) + end +end + +# ╔═╡ 5968317a-6c07-4655-8137-6702656bb3b4 +mean_zero_counts + +# ╔═╡ 48be49d0-0b60-44f3-8152-1ca917a4232e +plot(p_values, mean_zero_counts) + +# ╔═╡ Cell order: +# ╠═6b342f14-e7d5-11ef-1ea0-77ceb0d78f32 +# ╠═7bfcc024-7d7f-4c0b-b918-8f7626b10974 +# ╠═309035dd-5653-48a6-a53d-817e743279fa +# ╠═71140c81-af29-4857-8020-4f94c8bd64b3 +# ╠═284f5847-9c15-41f3-a595-1e12a22df69f +# ╟─1f975552-b0b8-4830-8dcc-214574d4fc38 +# ╟─d2f32eab-0b35-4794-9219-5bcbb4c069c5 +# ╟─016842c9-9479-4061-a27e-9dc006121f23 +# ╠═6c97bf81-ef32-45a4-aa7c-c8c26ba2d2c3 +# ╟─9c7ab7fb-7380-41a3-85ea-714478ade218 +# ╠═1536fe23-0f8d-4b86-98d2-076248b35954 +# ╟─e6e2ff5c-38eb-4ba3-b430-c9031483a0a5 +# ╠═ab6af765-1cde-4da8-bbc1-a5fab391db54 +# ╟─378878a0-5c09-4eb0-ac43-1031014ff12a +# ╠═3ae98e83-7beb-4597-89be-80c813d4349b +# ╟─988f79c0-9c7b-4752-a7f2-d4473ad73ce6 +# ╠═0d8f53f8-0a14-4ac6-bd0c-2190d4db0909 +# ╟─08d43ac8-a973-4d7b-baf7-4c37e54cfe24 +# ╠═e20e4dd8-bdbb-4005-af68-6bf7e4ec130e +# ╟─238e1120-34af-4d57-8efa-aa80ab28a874 +# ╠═d4c45709-70c9-4ba0-8fb8-6b600473723d +# ╟─d06fb076-76e4-4248-a940-96804ea68833 +# ╠═7644adf4-d992-48b1-b40a-12fdf30f6cb5 +# ╟─59d2d3e1-354b-4444-b8a1-16ad8ea2ba94 +# ╟─53e15d08-7d1d-4bc2-9fd3-4bc6fbb9de84 +# ╠═1511d269-9706-41b1-b8e2-f85ebcedc2d8 +# ╟─fd9b4521-30d2-4af3-a080-40e9dbf27aa4 +# ╠═3cac94d6-14ae-42f2-b0b6-f47d87cdb518 +# ╟─b2797fd2-b6e1-4bb0-af23-26931fe8be69 +# ╠═f947c2d9-9123-422d-8972-157717c85b3c +# ╟─74708270-b1ec-48c7-af32-3b970b92c706 +# ╠═2b00df5d-994e-47a1-8068-c93ce3f1a618 +# ╟─9d06c31e-3525-4889-a1de-3fe02413c7d8 +# ╠═9a90f800-3669-4831-b50b-c5405bbb9a03 +# ╟─a554fd16-aa3d-48ca-8de6-5582725c27d8 +# ╟─d6872046-b5ef-4c2d-a9bb-2418f57f715d +# ╟─747e20c4-b06b-4e78-a09a-55053cf42bf4 +# ╟─92181028-60fc-4830-afba-2380ac91455d +# ╠═f8942b10-773a-4b22-baad-8004fba8bd34 +# ╟─43d41284-053d-4dfe-8d5b-96be70c0495c +# ╠═a999ae2a-7567-41e7-9c0c-e94fad6f5d46 +# ╟─9ebb5b44-04d7-4b89-acdb-e40a245703d2 +# ╠═049de8d5-b221-452b-b2c4-9bc1e0c17f48 +# ╟─a73a2853-1f48-4179-9771-083794d3f137 +# ╠═682e9120-0e1c-4dfa-9ec6-66bb0a3f4374 +# ╠═5968317a-6c07-4655-8137-6702656bb3b4 +# ╠═48be49d0-0b60-44f3-8152-1ca917a4232e