-
Notifications
You must be signed in to change notification settings - Fork 6
/
hw5 really.jl
1235 lines (974 loc) · 39.1 KB
/
hw5 really.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
### A Pluto.jl notebook ###
# v0.19.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)
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
end
# ╔═╡ 2b37ca3a-0970-11eb-3c3d-4f788b411d1a
begin
using Pkg
Pkg.activate(mktempdir())
end
# ╔═╡ 2dcb18d0-0970-11eb-048a-c1734c6db842
begin
Pkg.add(["PlutoUI", "Plots"])
using Plots
gr()
using PlutoUI
end
# ╔═╡ faec52a8-0a60-11eb-082a-f5787b09d88c
using Statistics
# ╔═╡ 19fe1ee8-0970-11eb-2a0d-7d25e7d773c6
md"_homework 5, version 0_"
# ╔═╡ 49567f8e-09a2-11eb-34c1-bb5c0b642fe8
# WARNING FOR OLD PLUTO VERSIONS, DONT DELETE ME
html"""
<script>
const warning = html`
<h2 style="color: #800">Oopsie! You need to update Pluto to the latest version for this homework</h2>
<p>Close Pluto, go to the REPL, and type:
<pre><code>julia> import Pkg
julia> Pkg.update("Pluto")
</code></pre>
`
const super_old = window.version_info == null || window.version_info.pluto == null
if(super_old) {
return warning
}
const version_str = window.version_info.pluto.substring(1)
const numbers = version_str.split(".").map(Number)
console.log(numbers)
if(numbers[0] > 0 || numbers[1] > 12 || numbers[2] > 1) {
} else {
return warning
}
</script>
"""
# ╔═╡ 181e156c-0970-11eb-0b77-49b143cc0fc0
md"""
# **Homework 5**: _Epidemic modeling II_
`18.S191`, fall 2020
This notebook contains _built-in, live answer checks_! In some exercises you will see a coloured box, which runs a test case on your code, and provides feedback based on the result. Simply edit the code, run it, and the check runs again.
_For MIT students:_ there will also be some additional (secret) test cases that will be run as part of the grading process, and we will look at your notebook and write comments.
Feel free to ask questions!
"""
# ╔═╡ 1f299cc6-0970-11eb-195b-3f951f92ceeb
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)
student = (name = "Jazzy Doe", kerberos_id = "jazz")
# you might need to wait until all other cells in this notebook have completed running.
# scroll around the page to see what's up
# ╔═╡ 1bba5552-0970-11eb-1b9a-87eeee0ecc36
md"""
Submission by: **_$(student.name)_** ($(student.kerberos_id)@mit.edu)
"""
# ╔═╡ 2848996c-0970-11eb-19eb-c719d797c322
md"_Let's create a package environment:_"
# ╔═╡ 69d12414-0952-11eb-213d-2f9e13e4b418
md"""
In this problem set, we will look at a simple **spatial** agent-based epidemic model: agents can interact only with other agents that are *nearby*. (In the previous homework any agent could interact with any other, which is not realistic.)
A simple approach is to use **discrete space**: each agent lives
in one cell of a square grid. For simplicity we will allow no more than
one agent in each cell, but this requires some care to
design the rules of the model to respect this.
We will adapt some functionality from the previous homework. You should copy and paste your code from that homework into this notebook.
"""
# ╔═╡ 3e54848a-0954-11eb-3948-f9d7f07f5e23
md"""
## **Exercise 1:** _Wandering at random in 2D_
In this exercise we will implement a **random walk** on a 2D lattice (grid). At each time step, a walker jumps to a neighbouring position at random (i.e. chosen with uniform probability from the available adjacent positions).
"""
# ╔═╡ 3e623454-0954-11eb-03f9-79c873d069a0
md"""
#### Exercise 1.1
We define a struct type `Coordinate` that contains integers `x` and `y`.
"""
# ╔═╡ 0ebd35c8-0972-11eb-2e67-698fd2d311d2
struct Coordinate
x::Int
y::Int
end
# ╔═╡ 027a5f48-0a44-11eb-1fbf-a94d02d0b8e3
md"""
👉 Construct a `Coordinate` located at the origin.
"""
# ╔═╡ 1408f93e-0972-11eb-34b7-fdbac820be10
origin = Coordinate(0,0)
# ╔═╡ b2f90634-0a68-11eb-1618-0b42f956b5a7
# origin = missing
# ╔═╡ 3e858990-0954-11eb-3d10-d10175d8ca1c
md"""
👉 Write a function `make_tuple` that takes an object of type `Coordinate` and returns the corresponding tuple `(x, y)`. Boring, but useful later!
"""
# ╔═╡ 189bafac-0972-11eb-1893-094691b2073c
function make_tuple(c::Coordinate)
return (c.x, c.y)
end
# ╔═╡ d2051d9c-0a68-11eb-07c5-53658fd5b657
# function make_tuple(c::Coordinate)
# return missing
# end
# ╔═╡ 46a602c4-0a28-11eb-156e-3f374e29170d
make_tuple(Coordinate(5,6))
# ╔═╡ 73ed1384-0a29-11eb-06bd-d3c441b8a5fc
md"""
#### Exercise 1.2
In Julia, operations like `+` and `*` are just functions, and they are treated like any other function in the language. The only special property you can use the _infix notation_: you can write
```julia
1 + 2
```
instead of
```julia
+(1, 2)
```
_(There are [lots of special 'infixable' function names](https://github.com/JuliaLang/julia/blob/master/src/julia-parser.scm#L23-L24) that you can use for your own functions!)_
We you call it with the prefix notation, it becomes clear that it really is 'just another function', with lots of predefined methods.
"""
# ╔═╡ 9c9f53b2-09ea-11eb-0cda-639764250cee
md"""
> #### Extending + in the wild
> Because it is a function, we can add our own methods to it! This feature is super useful in general languages like Julia and Python, because it lets you use familiar syntax (`a + b*c`) on objects that are not necessarily numbers!
>
> One example we've see before is the `RGB` type in Homework 1. You are able to do:
> ```julia
> 0.5 * RGB(0.1, 0.7, 0.6)
> ```
> to multiply each color channel by $0.5$. This is possible because `Images.jl` [wrote a method](https://github.com/JuliaGraphics/ColorVectorSpace.jl/blob/master/src/ColorVectorSpace.jl#L131):
> ```julia
> *(::Real, ::AbstractRGB)::AbstractRGB
> ```
👉 Implement addition on two `Coordinate` structs by adding a method to `Base.:+`
"""
# ╔═╡ b5c7eb1a-09ea-11eb-27fe-936c47c19a72
function Base.:+(a::Coordinate, b::Coordinate)
return Coordinate(a.x + b.x, a.y + b.y)
end
# ╔═╡ 96707ef0-0a29-11eb-1a3e-6bcdfb7897eb
+(1, 2)
# ╔═╡ b0337d24-0a29-11eb-1fab-876a87c0973f
+
# ╔═╡ e24d5796-0a68-11eb-23bb-d55d206f3c40
# function Base.:+(a::Coordinate, b::Coordinate)
# return missing
# end
# ╔═╡ ec8e4daa-0a2c-11eb-20e1-c5957e1feba3
Coordinate(3,4) + Coordinate(10,10)
# ╔═╡ e144e9d0-0a2d-11eb-016e-0b79eba4b2bb
md"""
_Pluto has some trouble here, you need to manually re-run the cell above!_
"""
# ╔═╡ 71c358d8-0a2f-11eb-29e1-57ff1915e84a
md"""
#### Exercise 1.3
In our model, agents will be able to walk in 4 directions: up, down, left and right. We can define these directions as `Coordinate`s.
"""
# ╔═╡ 5278e232-0972-11eb-19ff-a1a195127297
possible_moves = [
Coordinate( 1, 0),
Coordinate( 0, 1),
Coordinate(-1, 0),
Coordinate( 0,-1),
]
# ╔═╡ 448f5870-0a29-11eb-395f-7dd102358015
md"""
👉 Use the `rand` function to choose a random move, and use `move` to move `test_move_wanderer` in that direction.
"""
# ╔═╡ 8f829ffc-0a30-11eb-1c1a-bde01253db94
coordinate_after_random_move = Coordinate(3,4) + rand(possible_moves)
# ╔═╡ e8dc6e46-0a68-11eb-0d57-973612b77ae9
# coordinate_after_random_move = missing
# ╔═╡ 3eb46664-0954-11eb-31d8-d9c0b74cf62b
md"""
We are able to make a `Coordinate` perform one random step, by adding a move to it. Great!
👉 Write a function `trajectory` that calculates a trajectory of a `Wanderer` `w` when performing `n` steps., i.e. the sequence of positions that the walker finds itself in.
Possible steps:
- Use `rand(possible_moves, n)` to generate a vector of `n` random moves. Each possible move will be equally likely.
- To compute the trajectory you can use either of the following two approaches:
1. 🆒 Use the function `accumulate` (see the live docs for `accumulate`). Use `move` as the function passed to `accumulate` and the `w` as the starting value (`init` keyword argument).
1. Use a `for` loop calling `move`.
"""
# ╔═╡ 548b9c16-0972-11eb-3ce9-838a63d3524b
function trajectory(w::Coordinate, n::Int)
return accumulate(+, rand(possible_moves, n), init=w)
end
# ╔═╡ edf86a0e-0a68-11eb-2ad3-dbf020037019
# function trajectory(w::Coordinate, n::Int)
# return missing
# end
# ╔═╡ 478309f4-0a31-11eb-08ea-ade1755f53e0
function plot_trajectory!(p::Plots.Plot, trajectory::Vector{Coordinate}; kwargs...)
plot!(p, make_tuple.(trajectory);
label=nothing,
linewidth=2,
linealpha=LinRange(1.0, 0.2, length(trajectory)),
kwargs...)
end
# ╔═╡ 3ebd436c-0954-11eb-170d-1d468e2c7a37
md"""
#### Exercise 1.4
👉 Plot 10 trajectories of length 1000 on a single figure, all starting at the origin. Use the function `plot_trajectory!` as demonstrated above.
Remember from last week that you can compose plots like this:
```julia
let
# Create a new plot with aspect ratio 1:1
p = plot(ratio=1)
plot_trajectory!(p, test_trajectory) # plot one trajectory
plot_trajectory!(p, another_trajectory) # plot the second one
...
p
end
```
"""
# ╔═╡ dcefc6fe-0a3f-11eb-2a96-ddf9c0891873
# ╔═╡ b4d5da4a-09a0-11eb-1949-a5807c11c76c
md"""
#### Exercise 1.5
Agents live in a box of side length $2L$, centered at the origin. We need to decide (i.e. model) what happens when they reach the walls of the box (boundaries), in other words what kind of **boundary conditions** to use.
One relatively simple boundary condition is a **collision boundary**:
> Each wall of the box is a wall, modelled using "collision": if the walker tries to jump beyond the wall, it ends up at the position inside the box that is closest to the goal.
👉 Write a function `collide_boundary` which takes a `Coordinate` `c` and a size $L$, and returns a new coordinate that lies inside the box (i.e. ``[-L,L]\times [-L,L]``), but is closest to `c`. This is similar to `extend_mat` from Homework 1.
"""
# ╔═╡ d07744dc-09a0-11eb-3f65-85fff408b7fb
function collide_boundary(c::Coordinate, L::Number)
return Coordinate(clamp(c.x, -L, L), clamp(c.y, -L, L))
end
# ╔═╡ 0237ebac-0a69-11eb-2272-35ea4e845d84
# function collide_boundary(c::Coordinate, L::Number)
# return missing
# end
# ╔═╡ ad832360-0a40-11eb-2857-e7f0350f3b12
collide_boundary(Coordinate(12,4), 10)
# ╔═╡ b4ed2362-09a0-11eb-0be9-99c91623b28f
md"""
#### Exercise 1.6
👉 Implement a 3-argument method of `trajectory` where the third argument is a size. The trajectory returned should be within the boundary (use `reflect_boundary` from above). You can still use `accumulate` with an anonymous function that makes a move and then reflects the resulting coordinate, or use a for loop.
"""
# ╔═╡ a4e93fe4-099f-11eb-0a29-d10c28a2d9af
function trajectory(c::Coordinate, n::Int, L::Number)
return accumulate(rand(moves, n), init=c) do old, m
new = collide_boundary(old + m, L)
end
end
# ╔═╡ 44107808-096c-11eb-013f-7b79a90aaac8
test_trajectory = trajectory(Coordinate(4,4), 30)
# ╔═╡ 87ea0868-0a35-11eb-0ea8-63e27d8eda6e
try
p = plot(ratio=1, size=(650,200))
plot_trajectory!(p, test_trajectory; color="black", showaxis=false, axis=nothing, linewidth=4)
p
catch
end
# ╔═╡ 51788e8e-0a31-11eb-027e-fd9b0dc716b5
let
long_trajectory = trajectory(Coordinate(4,4), 1000)
p = plot(ratio=1)
plot_trajectory!(p, long_trajectory)
p
end
# ╔═╡ 0665aa3e-0a69-11eb-2b5d-cd718e3c7432
# function trajectory(c::Coordinate, n::Int, L::Number)
# return missing
# end
# ╔═╡ 3ed06c80-0954-11eb-3aee-69e4ccdc4f9d
md"""
## **Exercise 2:** _Wanderering Agents_
In this exercise we will create Agents which have a location as well as some infection state information.
Let's define a type `Agent`. `Agent` contains a `position` (of type `Coordinate`), as well as a `state` of type `InfectionStatus` (as in Homework 4).)
(For simplicity we will not use a `num_infected` field, but feel free to do so!)
"""
# ╔═╡ 35537320-0a47-11eb-12b3-931310f18dec
@enum InfectionStatus S I R
# ╔═╡ cf2f3b98-09a0-11eb-032a-49cc8c15e89c
mutable struct Agent
position::Coordinate
status::InfectionStatus
end
# ╔═╡ 814e888a-0954-11eb-02e5-0964c7410d30
md"""
#### Exercise 2.1
👉 Write a function `initialize` that takes parameters $N$ and $L$, where $N$ is the number of agents abd $2L$ is the side length of the square box where the agents live.
It returns a `Vector` of `N` randomly generated `Agent`s. Their coordinates are randomly sampled in the ``[-L,L] \times [-L,L]`` box, and the agents are all susceptible, except one, chosen at random, which is infectious.
"""
# ╔═╡ 19fdccb6-0a46-11eb-3632-2f94e842414f
function initialize(N, L)
Agent.(Coordinate.(rand(-L:L, N), rand(-L:L, N)), [I, fill(S,N-1)...])
end
# ╔═╡ 0cfae7ba-0a69-11eb-3690-d973d70e47f4
# function initialize(N::Number, L::Number)
# return missing
# end
# ╔═╡ 1d0f8eb4-0a46-11eb-38e7-63ecbadbfa20
initialize(3, 10)
# ╔═╡ b5a88504-0a47-11eb-0eda-f125d419e909
position(a::Agent) = a.position
# ╔═╡ e0b0880c-0a47-11eb-0db2-f760bbbf9c11
color(s::InfectionStatus) = if s == S
"blue"
elseif s == I
"red"
else
"green"
end
# ╔═╡ 87a4cdaa-0a5a-11eb-2a5e-cfaf30e942ca
color(a::Agent) = color(a.status)
# ╔═╡ 49fa8092-0a43-11eb-0ba9-65785ac6a42f
md"""
#### Exercise 2.2
👉 Write a function `visualize` that takes in a collection of agents as argument, and the box size `L`. It should plot a point for each agent at its location, coloured according to its status.
You can use the keyword argument `c=color.(agents)` inside your call to the plotting function make the point colors correspond to the infection statuses. Don't forget to use `ratio=1`.
"""
# ╔═╡ 1ccc961e-0a69-11eb-392b-915be07ef38d
# function visualize(agents::Vector, L)
# return missing
# end
# ╔═╡ f953e06e-099f-11eb-3549-73f59fed8132
md"""
### Exercise 3: Spatial epidemic model -- Dynamics
Last week we wrote a function `interact!` that takes two agents, `agent` and `source`, and an infection of type `InfectionRecovery`, which models the interaction between two agent, and possibly modifies `agent` with a new status.
This week, we define a new infection type, `CollisionInfectionRecovery`, and a new method that is the same as last week, except it **only infects `agent` if `agents.position==source.position`**.
"""
# ╔═╡ d929f7d4-0a4d-11eb-0ff6-5111155c549d
bernoulli(p::Number) = rand() < p
# ╔═╡ e6dd8258-0a4b-11eb-24cb-fd5b3554381b
abstract type AbstractInfection end
# ╔═╡ de88b530-0a4b-11eb-05f7-85171594a8e8
struct CollisionInfectionRecovery <: AbstractInfection
p_infection
p_recovery
end
# ╔═╡ d1bcd5c4-0a4b-11eb-1218-7531e367a7ff
function interact!(agent::Agent, source::Agent, infection::CollisionInfectionRecovery)
if agent.position == source.position && agent.status == S && source.status == I
if bernoulli(infection.p_infection)
agent.status = I
end
elseif agent.status == I
if bernoulli(infection.p_recovery)
agent.status = R
end
end
end
# ╔═╡ 34778744-0a5f-11eb-22b6-abe8b8fc34fd
md"""
#### Exercise 3.1
Your turn!
👉 Write a function `step!` that takes a vector of `Agent`s, a box size `L` and an `infection`. This that does one step of the dynamics on a vector of agents.
- Choose an Agent `source` at random.
- Move the `source` one step, and use `collide_boundary` to ensure that our agent stays within the box.
- For all _other_ agents, call `interact!(other_agent, source, infection)`.
- return the array `agents` again.
"""
# ╔═╡ 24fe0f1a-0a69-11eb-29fe-5fb6cbf281b8
# function step!(agents::Vector, L::Number, infection::AbstractInfection)
# return missing
# end
# ╔═╡ 1fc3271e-0a45-11eb-0e8d-0fd355f5846b
md"""
#### Exercise 3.2
If we call `step!` `N` times, then every agent will have made one step, on average. Let's call this one _sweep_ of the simulation.
👉 Create a before-and-after plot of ``k_{sweeps}=1000`` sweeps.
- Initialize a new vector of agents (`N=50`, `L=40`, `infection` is given as `pandemic` below).
- Plot the state using `visualize`, and save the plot as a variable `plot_before`.
- Run `k_sweeps` sweeps.
- Plot the state again, and store as `plot_after`.
- Combine the two plots into a single figure using
```julia
plot(plot_before, plot_after)
```
"""
# ╔═╡ 18552c36-0a4d-11eb-19a0-d7d26897af36
pandemic = CollisionInfectionRecovery(0.5, 0.00001)
# ╔═╡ 4e7fd58a-0a62-11eb-1596-c717e0845bd5
@bind k_sweeps Slider(1:10000, default=1000)
# ╔═╡ 778c2490-0a62-11eb-2a6c-e7fab01c6822
# let
# N = 50
# L = 40
# plot_before = plot(1:3) # replace with your code
# plot_after = plot(1:3)
# plot(plot_before, plot_after)
# end
# ╔═╡ e964c7f0-0a61-11eb-1782-0b728fab1db0
md"""
#### Exercise 3.3
Every time that you move the slider, a completely new simulation is created an run. This makes it hard to view the progress of a single simulation over time. So in this exercise, we we look at a single simulation, and plot the S, I and R curves.
👉 Plot the SIR curves of a single simulation, with the same parameters as in the previous exercise. Use `k_sweep_max = 10000` as the total number of sweeps.
"""
# ╔═╡ 4d83dbd0-0a63-11eb-0bdc-757f0e721221
k_sweep_max = 10000
# ╔═╡ e4eabe9e-0a63-11eb-0aa4-939536369dbd
# let
# N = 50
# L = 30
# missing
# end
# ╔═╡ 201a3810-0a45-11eb-0ac9-a90419d0b723
md"""
#### Exercise 3.4 (optional)
Let's make our plot come alive! There are two options to make our visualization dynamic:
👉1️⃣ Precompute one simulation run and save its intermediate states using `deepcopy`. You can then write an interactive visualization that shows both the state at time $t$ (using `visualize`) and the history of $S$, $I$ and $R$ from time $0$ up to time $t$. $t$ is controlled by a slider.
👉2️⃣ Use `@gif` from Plots.jl to turn a sequence of plots into an animation. Be careful to skip about 50 sweeps between each animation frame, otherwise the GIF becomes too large.
This an optional exercise, and our solution to 2️⃣ is given below.
"""
# ╔═╡ e5040c9e-0a65-11eb-0f45-270ab8161871
# let
# N = 50
# L = 30
# missing
# end
# ╔═╡ 2031246c-0a45-11eb-18d3-573f336044bf
md"""
#### Exercise 3.5
👉 Using $L=20$ and $N=100$, experiment with the infection and recovery probabilities until you find an epidemic outbreak. (Take the recovery probability quite small.) Modify the two infections below to match your observations.
"""
# ╔═╡ 63dd9478-0a45-11eb-2340-6d3d00f9bb5f
causes_outbreak = CollisionInfectionRecovery(0.5, 0.001)
# ╔═╡ 269955e4-0a46-11eb-02cc-1946dc918bfa
does_not_cause_outbreak = CollisionInfectionRecovery(0.5, 0.001)
# ╔═╡ 4d4548fe-0a66-11eb-375a-9313dc6c423d
# ╔═╡ 20477a78-0a45-11eb-39d7-93918212a8bc
md"""
#### Exercise 3.6
👉 With the parameters of Exercise 3.2, run 50 simulations. Plot $S$, $I$ and $R$ as a function of time for each of them (with transparency!). This should look qualitatively similar to what you saw in the previous homework. You probably need different `p_infection` and `p_recovery` values from last week. Why?
"""
# ╔═╡ 601f4f54-0a45-11eb-3d6c-6b9ec75c6d4a
# ╔═╡ b1b1afda-0a66-11eb-2988-752405815f95
need_different_parameters_because = md"""
i say so
"""
# ╔═╡ 05c80a0c-09a0-11eb-04dc-f97e306f1603
md"""
## **Exercise 4:** _Effect of socialization_
In this exercise we'll modify the simple mixing model. Instead of a constant mixing probability, i.e. a constant probability that any pair of people interact on a given day, we will have a variable probability associated with each agent, modelling the fact that some people are more or less social or contagious than others.
"""
# ╔═╡ b53d5608-0a41-11eb-2325-016636a22f71
md"""
#### Exercise 4.1
We create a new agent type `SocialAgent` with fields `position`, `status`, `num_infected`, and `social_score`. The attribute `social_score` represents an agent's probability of interacting with any other agent in the population.
"""
# ╔═╡ 1b5e72c6-0a42-11eb-3884-a377c72270c7
mutable struct SocialAgent
position::Coordinate
status
num_infected
social_score
end
# ╔═╡ e97e39aa-0a5d-11eb-3d5f-f90a0acfe5a2
position(a::SocialAgent) = a.position
# ╔═╡ b554b654-0a41-11eb-0e0d-e57ff68ced33
md"""
👉 Create a population of 500 agents, with `social_score`s chosen from 10 equally-spaced between 0.1 and 0.5.
"""
# ╔═╡ 1cd76cfc-0a42-11eb-2794-f100233b5ecd
function initialize_social(N, L)
social_agents = SocialAgent.(Coordinate.(rand(-L:L, N), rand(-L:L, N)),
S, 0, rand(LinRange(0.1, 0.5, 10), N))
social_agents[2].status = I
social_agents
end
# ╔═╡ 40c1c1d6-0a69-11eb-3913-59e9b9ec4332
# function initialize_social(N, L)
# return missing
# end
# ╔═╡ b56ba420-0a41-11eb-266c-719d39580fa9
md"""
#### Exercise 4.2
Not all two agents who end up in the same grid point may actually iterect in an infectious way -- they may just be passing by and do not create enough exposure for communicating the disease.
👉 Write a new `interact!` method on `SocialAgent` which adds together the social_scores for two agents and uses that as the probability that they interact in a risky way. Only if they interact in a risky way, the infection is transmitted with the usual probability.
"""
# ╔═╡ 1df4b57c-0a42-11eb-21ba-0d14438912b5
function interact!(agent::SocialAgent, source::SocialAgent, infection::CollisionInfectionRecovery)
p = agent.social_score + source.social_score
if bernoulli(p)
if agent.position == source.position && agent.status == S && source.status == I
if bernoulli(infection.p_infection)
agent.status = I
source.num_infected += 1
end
elseif agent.status == I
if bernoulli(infection.p_recovery)
agent.status = R
end
end
end
end
# ╔═╡ 4b79f25c-0a45-11eb-07fb-530b1628ff3b
function step!(agents::Vector, L::Number, infection::AbstractInfection)
source = rand(agents)
source.position = collide_boundary(source.position + rand(possible_moves), L)
for b in agents
if source !== b
interact!(source, b, infection)
end
end
agents
end
# ╔═╡ ef27de84-0a63-11eb-177f-2197439374c5
let
N = 50
L = 30
agents = initialize(N, L)
SIRs = map(1:k_sweep_max) do t
for _ in 1:N
step!(agents, L, pandemic)
end
sum(a -> a.status == S, agents), sum(a -> a.status == I, agents), sum(a -> a.status == R, agents)
end
p = plot()
plot!(p, 1:k_sweep_max, [a[1] for a in SIRs])
plot!(p, 1:k_sweep_max, [a[2] for a in SIRs])
plot!(p, 1:k_sweep_max, [a[3] for a in SIRs])
p
end
# ╔═╡ 465e918a-0a69-11eb-1b59-01150b4b0f36
# function interact!(agent::SocialAgent, source::SocialAgent, infection::CollisionInfectionRecovery)
# # your code here
# end
# ╔═╡ 76d84954-0a5d-11eb-1223-e976b12486d0
color(s::SocialAgent) = color(s.status)
# ╔═╡ 1e1c1064-0a46-11eb-1e4a-8b89c73b5ed0
function visualize(agents::Vector, L)
coords = make_tuple.(position.(agents))
p = plot(xlim=(-L,L), ylim=(-L,L), ratio=1)
scatter!(p, first.(coords), last.(coords), c=color.(agents), label=nothing)
p
end
# ╔═╡ 1f96c80a-0a46-11eb-0690-f51c60e57c3f
let
N = 20
L = 10
visualize(initialize(N, L), L)
end
# ╔═╡ e822ac2c-0a61-11eb-2abf-4dbca6145ee8
let
N = 50
L = 40
agents = initialize(N, L)
plot_before = visualize(agents, L)
for _ in 1:N*k_sweeps
step!(agents, L, pandemic)
end
plot_after = visualize(agents, L)
plot(plot_before, plot_after)
end
# ╔═╡ a885bf78-0a5c-11eb-2383-9d74c8765847
md"""
Make sure `step!`, `position`, `color`, work on the type `SocialAgent`. If `step!` takes an untyped first argument, it should work for both Agent and SocialAgent types without any changes. We actually only need to specialize `interact!` on SocialAgent.
#### Exercise 4.3
👉 Plot the SIR curves of the resulting simulation.
"""
# ╔═╡ 1f172700-0a42-11eb-353b-87c0039788bd
let
N = 50
L = 40
global social_agents = initialize_social(N, L)
x = social_agents
Ss, Is, Rs = [], [], []
Tmax = 200
@gif for t in 1:Tmax
for i in 1:50N
step!(x, L, pandemic)
end
push!(Ss, sum(a -> a.status == S, x))
push!(Is, sum(a -> a.status == I, x))
push!(Rs, sum(a -> a.status == R, x))
left = visualize(x, L)
right = plot(xlim=(1,Tmax), ylim=(1,N))
plot!(right, 1:t, Ss, color="blue", label="S")
plot!(right, 1:t, Is, color="red", label="I")
plot!(right, 1:t, Rs, color="green", label="R")
plot(left, right, size=(600,300))
end
# Animation(fastplot.(plots), 100)
end
# ╔═╡ 572f3b4c-0a69-11eb-313d-81c46847530c
# ╔═╡ b59de26c-0a41-11eb-2c67-b5f3c7780c91
md"""
#### Exercise 4.4
👉 Make a scatter plot showing each agent's mixing rate on one axis, and the `num_infected` from the simulation in the other axis. Run this simulation several times and comment on the results.
"""
# ╔═╡ 203ee2ee-0a42-11eb-2d10-bdb2210a5969
begin
scores = LinRange(0.1, 0.5, 10)
bar(scores,
legend=nothing,
[mean(x->x.num_infected,
filter(a->a.social_score == x, social_agents)) for x in scores])
end
# ╔═╡ 5924635a-0a69-11eb-05d2-d7bbb1e8117c
# ╔═╡ b5b4d834-0a41-11eb-1b18-1bd626d18934
md"""
Run a simulation for 100 steps, and then apply a "lockdown" where every agent's social score gets multiplied by 0.25, and then run a second simulation which runs on that same population from there. What do you notice? How does changing this factor form 0.25 to other numbers affect things?
"""
# ╔═╡ a83c96e2-0a5a-11eb-0e58-15b5dda7d2d2
# ╔═╡ 05fc5634-09a0-11eb-038e-53d63c3edaf2
md"""
## **Exercise 5:** (Optional) _Effect of distancing_
We can use a variant of the above model to investigate the effect of the
mis-named "social distancing"
(we want people to be *socially* close, but *physically* distant).
In this variant, we separate out the two effects "infection" and
"movement": an infected agent chooses a
neighbouring site, and if it finds a susceptible there then it infects it
with probability $p_I$. For simplicity we can ignore recovery.
Separately, an agent chooses a neighbouring site to move to,
and moves there with probability $p_M$ if the site is vacant. (Otherwise it
stays where it is.)
When $p_M = 0$, the agents cannot move, and hence are
completely quarantined in their original locations.
👉 How does the disease spread in this case?
"""
# ╔═╡ 24c2fb0c-0a42-11eb-1a1a-f1246f3420ff
# ╔═╡ c7649966-0a41-11eb-3a3a-57363cea7b06
md"""
👉 Run the dynamics repeatedly, and plot the sites which become infected.
"""
# ╔═╡ 2635b574-0a42-11eb-1daa-971b2596ce44
# ╔═╡ c77b085e-0a41-11eb-2fcb-534238cd3c49
md"""
👉 How does this change as you increase the *density*
$\rho = N / (L^2)$ of agents? Start with a small density.
This is basically the [**site percolation**](https://en.wikipedia.org/wiki/Percolation_theory) model.
When we increase $p_M$, we allow some local motion via random walks.
"""
# ╔═╡ 274fe006-0a42-11eb-1869-29193bb84957
# ╔═╡ c792374a-0a41-11eb-1e5b-89d9de2cf1f9
md"""
👉 Investigate how this leaky quarantine affects the infection dynamics with
different densities.
"""
# ╔═╡ d147f7f0-0a66-11eb-2877-2bc6680e396d
# ╔═╡ 0e6b60f6-0970-11eb-0485-636624a0f9d7
if student.name == "Jazzy Doe"
md"""
!!! danger "Before you submit"
Remember to fill in your **name** and **Kerberos ID** at the top of this notebook.
"""
end
# ╔═╡ 0a82a274-0970-11eb-20a2-1f590be0e576
md"## Function library
Just some helper functions used in the notebook."
# ╔═╡ 0aa666dc-0970-11eb-2568-99a6340c5ebd
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))
# ╔═╡ 8475baf0-0a63-11eb-1207-23f789d00802
hint(md"""
After every sweep, count the values $S$, $I$ and $R$ and appemd them to an array.
""")
# ╔═╡ f9b9e242-0a53-11eb-0c6a-4d9985ef1687
hint(md"""
```julia
let
N = 50
L = 40
x = initialize(N, L)
# initialize to empty arrays
Ss, Is, Rs = Int[], Int[], Int[]
Tmax = 200
@gif for t in 1:Tmax
for i in 1:50N
step!(x, L, pandemic)
end
# `sum` with two arguments count the number of elements
# of `x` for which the function returns `true`
push!(Ss, sum(a -> a.status == S, x))
push!(Is, sum(a -> a.status == I, x))
push!(Rs, sum(a -> a.status == R, x))
left = visualize(x, L)
right = plot(xlim=(1,Tmax), ylim=(1,N), size=(600,300))
plot!(right, 1:t, Ss, color=color(S), label="S")
plot!(right, 1:t, Is, color=color(I), label="I")
plot!(right, 1:t, Rs, color=color(R), label="R")
plot(left, right)
end
end
```
""")
# ╔═╡ 0acaf3b2-0970-11eb-1d98-bf9a718deaee
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))
# ╔═╡ 0afab53c-0970-11eb-3e43-834513e4632e
still_missing(text=md"Replace `missing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Here we go!", [text]))
# ╔═╡ 0b21c93a-0970-11eb-33b0-550a39ba0843
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]))
# ╔═╡ 0b470eb6-0970-11eb-182f-7dfb4662f827
yays = [md"Fantastic!", md"Splendid!", md"Great!", md"Yay ❤", md"Great! 🎉", md"Well done!", md"Keep it up!", md"Good job!", md"Awesome!", md"You got the right answer!", md"Let's move on to the next section."]
# ╔═╡ 0b6b27ec-0970-11eb-20c2-89515ee3ab88
correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]))
# ╔═╡ ec576da8-0a2c-11eb-1f7b-43dec5f6e4e7
let
# we need to call Base.:+ instead of + to make Pluto understand what's going on
# oops
result = Base.:+(Coordinate(3,4), Coordinate(10,10))
if result isa Missing
still_missing()
elseif !(result isa Coordinate)
keep_working(md"Make sure that your return a `Coordinate`. 🧭")
elseif result.x != 13 || result.y != 14
keep_working()
else
correct()
end
end
# ╔═╡ 0b901714-0970-11eb-0b6a-ebe739db8037
not_defined(variable_name) = Markdown.MD(Markdown.Admonition("danger", "Oopsie!", [md"Make sure that you define a variable called **$(Markdown.Code(string(variable_name)))**"]))
# ╔═╡ 66663fcc-0a58-11eb-3568-c1f990c75bf2
if !@isdefined(origin)
not_defined(:origin)
else
let
if origin isa Missing
still_missing()
elseif !(origin isa Coordinate)
keep_working(md"Make sure that `origin` is a `Coordinate`.")
else
if origin == Coordinate(0,0)
correct()
else
keep_working()
end
end
end
end
# ╔═╡ ad1253f8-0a34-11eb-265e-fffda9b6473f
if !@isdefined(make_tuple)
not_defined(:make_tuple)
else
let
result = make_tuple(Coordinate(2,1))
if result isa Missing
still_missing()
elseif !(result isa Tuple)
keep_working(md"Make sure that you return a `Tuple`, like so: `return (1, 2)`.")
else
if result == (2,1)
correct()
else