{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "using GraphViz\n", "\n", "using StockFlow\n", "using StockFlow.Syntax\n", "\n", "using Catlab\n", "using Catlab.CategoricalAlgebra\n", "using LabelledArrays\n", "using OrdinaryDiffEq\n", "using Plots\n", "\n", "using Catlab.Graphics\n", "using Catlab.Programs\n", "using Catlab.Theories\n", "using Catlab.WiringDiagrams" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# SEIR model" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "StockAndFlowF with elements S = 1:4, SV = 1:1, LS = 1:4, F = 1:8, I = 1:4, O = 1:7, V = 1:10, LV = 1:8, LSV = 1:2, P = 1:5, LVV = 1:2, LPV = 1:8\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Ssname
1S
2E
3I
4R
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
SVsvname
1N
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
LSlsslssv
111
221
331
441
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Ffvfname
11fbirth
24fincid
35fdeathS
46finf
57fdeathE
68frec
79fdeathI
810fdeathR
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Iifnis
111
222
343
464
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Oofnos
121
231
342
452
563
673
784
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Vvnamevop
1##v_fbirth#292*
2##v_fincid#293*
3##v_fincid#294*
4##v_fincid#295/
5##v_fdeathS#296*
6##v_finf#297/
7##v_fdeathE#298*
8##v_frec#299/
9##v_fdeathI#300*
10##v_fdeathR#301*
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
LVlvslvvlvsposition
1122
2332
3151
4261
5271
6381
7391
84101
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
LSVlsvsvlsvvlsvsvposition
1112
2142
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Ppname
1μ
2β
3tlatent
4trecovery
5δ
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
LVVlvsrclvtgtlvsrcposition
1231
2341
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
LPVlpvplpvvlpvpposition
1111
2221
3552
4362
5572
6482
7592
85102
\n", "
\n" ], "text/plain": [ "StockAndFlowF with elements S = 1:4, SV = 1:1, LS = 1:4, F = 1:8, I = 1:4, O = 1:7, V = 1:10, LV = 1:8, LSV = 1:2, P = 1:5, LVV = 1:2, LPV = 1:8\n", "┌───┬───────┐\n", "│\u001b[1m S \u001b[0m│\u001b[1m sname \u001b[0m│\n", "├───┼───────┤\n", "│ 1 │ S │\n", "│ 2 │ E │\n", "│ 3 │ I │\n", "│ 4 │ R │\n", "└───┴───────┘\n", "┌────┬────────┐\n", "│\u001b[1m SV \u001b[0m│\u001b[1m svname \u001b[0m│\n", "├────┼────────┤\n", "│ 1 │ N │\n", "└────┴────────┘\n", "┌────┬─────┬──────┐\n", "│\u001b[1m LS \u001b[0m│\u001b[1m lss \u001b[0m│\u001b[1m lssv \u001b[0m│\n", "├────┼─────┼──────┤\n", "│ 1 │ 1 │ 1 │\n", "│ 2 │ 2 │ 1 │\n", "│ 3 │ 3 │ 1 │\n", "│ 4 │ 4 │ 1 │\n", "└────┴─────┴──────┘\n", "┌───┬────┬─────────┐\n", "│\u001b[1m F \u001b[0m│\u001b[1m fv \u001b[0m│\u001b[1m fname \u001b[0m│\n", "├───┼────┼─────────┤\n", "│ 1 │ 1 │ fbirth │\n", "│ 2 │ 4 │ fincid │\n", "│ 3 │ 5 │ fdeathS │\n", "│ 4 │ 6 │ finf │\n", "│ 5 │ 7 │ fdeathE │\n", "│ 6 │ 8 │ frec │\n", "│ 7 │ 9 │ fdeathI │\n", "│ 8 │ 10 │ fdeathR │\n", "└───┴────┴─────────┘\n", "┌───┬─────┬────┐\n", "│\u001b[1m I \u001b[0m│\u001b[1m ifn \u001b[0m│\u001b[1m is \u001b[0m│\n", "├───┼─────┼────┤\n", "│ 1 │ 1 │ 1 │\n", "│ 2 │ 2 │ 2 │\n", "│ 3 │ 4 │ 3 │\n", "│ 4 │ 6 │ 4 │\n", "└───┴─────┴────┘\n", "┌───┬─────┬────┐\n", "│\u001b[1m O \u001b[0m│\u001b[1m ofn \u001b[0m│\u001b[1m os \u001b[0m│\n", "├───┼─────┼────┤\n", "│ 1 │ 2 │ 1 │\n", "│ 2 │ 3 │ 1 │\n", "│ 3 │ 4 │ 2 │\n", "│ 4 │ 5 │ 2 │\n", "│ 5 │ 6 │ 3 │\n", "│ 6 │ 7 │ 3 │\n", "│ 7 │ 8 │ 4 │\n", "└───┴─────┴────┘\n", "┌────┬─────────────────┬─────┐\n", "│\u001b[1m V \u001b[0m│\u001b[1m vname \u001b[0m│\u001b[1m vop \u001b[0m│\n", "├────┼─────────────────┼─────┤\n", "│ 1 │ ##v_fbirth#292 │ * │\n", "│ 2 │ ##v_fincid#293 │ * │\n", "│ 3 │ ##v_fincid#294 │ * │\n", "│ 4 │ ##v_fincid#295 │ / │\n", "│ 5 │ ##v_fdeathS#296 │ * │\n", "│ 6 │ ##v_finf#297 │ / │\n", "│ 7 │ ##v_fdeathE#298 │ * │\n", "│ 8 │ ##v_frec#299 │ / │\n", "│ 9 │ ##v_fdeathI#300 │ * │\n", "│ 10 │ ##v_fdeathR#301 │ * │\n", "└────┴─────────────────┴─────┘\n", "┌────┬─────┬─────┬─────────────┐\n", "│\u001b[1m LV \u001b[0m│\u001b[1m lvs \u001b[0m│\u001b[1m lvv \u001b[0m│\u001b[1m lvsposition \u001b[0m│\n", "├────┼─────┼─────┼─────────────┤\n", "│ 1 │ 1 │ 2 │ 2 │\n", "│ 2 │ 3 │ 3 │ 2 │\n", "│ 3 │ 1 │ 5 │ 1 │\n", "│ 4 │ 2 │ 6 │ 1 │\n", "│ 5 │ 2 │ 7 │ 1 │\n", "│ 6 │ 3 │ 8 │ 1 │\n", "│ 7 │ 3 │ 9 │ 1 │\n", "│ 8 │ 4 │ 10 │ 1 │\n", "└────┴─────┴─────┴─────────────┘\n", "┌─────┬───────┬──────┬───────────────┐\n", "│\u001b[1m LSV \u001b[0m│\u001b[1m lsvsv \u001b[0m│\u001b[1m lsvv \u001b[0m│\u001b[1m lsvsvposition \u001b[0m│\n", "├─────┼───────┼──────┼───────────────┤\n", "│ 1 │ 1 │ 1 │ 2 │\n", "│ 2 │ 1 │ 4 │ 2 │\n", "└─────┴───────┴──────┴───────────────┘\n", "┌───┬───────────┐\n", "│\u001b[1m P \u001b[0m│\u001b[1m pname \u001b[0m│\n", "├───┼───────────┤\n", "│ 1 │ μ │\n", "│ 2 │ β │\n", "│ 3 │ tlatent │\n", "│ 4 │ trecovery │\n", "│ 5 │ δ │\n", "└───┴───────────┘\n", "┌─────┬───────┬───────┬───────────────┐\n", "│\u001b[1m LVV \u001b[0m│\u001b[1m lvsrc \u001b[0m│\u001b[1m lvtgt \u001b[0m│\u001b[1m lvsrcposition \u001b[0m│\n", "├─────┼───────┼───────┼───────────────┤\n", "│ 1 │ 2 │ 3 │ 1 │\n", "│ 2 │ 3 │ 4 │ 1 │\n", "└─────┴───────┴───────┴───────────────┘\n", "┌─────┬──────┬──────┬──────────────┐\n", "│\u001b[1m LPV \u001b[0m│\u001b[1m lpvp \u001b[0m│\u001b[1m lpvv \u001b[0m│\u001b[1m lpvpposition \u001b[0m│\n", "├─────┼──────┼──────┼──────────────┤\n", "│ 1 │ 1 │ 1 │ 1 │\n", "│ 2 │ 2 │ 2 │ 1 │\n", "│ 3 │ 5 │ 5 │ 2 │\n", "│ 4 │ 3 │ 6 │ 2 │\n", "│ 5 │ 5 │ 7 │ 2 │\n", "│ 6 │ 4 │ 8 │ 2 │\n", "│ 7 │ 5 │ 9 │ 2 │\n", "│ 8 │ 5 │ 10 │ 2 │\n", "└─────┴──────┴──────┴──────────────┘\n" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "seir = @stock_and_flow begin\n", " :stocks\n", " S\n", " E\n", " I\n", " R\n", "\n", " :parameters\n", " μ\n", " β\n", " tlatent\n", " trecovery\n", " δ\n", "\n", " :flows\n", " CLOUD => fbirth(μ * N) => S\n", " S => fincid(β * S * I / N) => E\n", " S => fdeathS(S * δ) => CLOUD\n", " E => finf(E / tlatent) => I\n", " E => fdeathE(E * δ) => CLOUD\n", " I => frec(I / trecovery) => R\n", " I => fdeathI(I * δ) => CLOUD\n", " R => fdeathR(R * δ) => CLOUD\n", "\n", " :sums\n", " N = [S, E, I, R]\n", "\n", "end\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", "\n", "\n", "G\n", "\n", "\n", "\n", "s1\n", "\n", "S\n", "\n", "\n", "\n", "v2\n", "β * S\n", "\n", "\n", "\n", "s1->v2\n", "\n", "\n", "\n", "\n", "\n", "v4\n", "((β * S) * I) / N\n", "\n", "\n", "\n", "s1->v4\n", "\n", "\n", "\n", "\n", "\n", "\n", "v5\n", "S * δ\n", "\n", "\n", "\n", "s1->v5\n", "\n", "\n", "\n", "\n", "\n", "\n", "s1->v5\n", "\n", "\n", "\n", "\n", "\n", "sv1\n", "\n", "N\n", "\n", "\n", "\n", "s1->sv1\n", "\n", "\n", "\n", "\n", "\n", "s2\n", "\n", "E\n", "\n", "\n", "\n", "v6\n", "E / tlatent\n", "\n", "\n", "\n", "s2->v6\n", "\n", "\n", "\n", "\n", "\n", "\n", "s2->v6\n", "\n", "\n", "\n", "\n", "\n", "v7\n", "E * δ\n", "\n", "\n", "\n", "s2->v7\n", "\n", "\n", "\n", "\n", "\n", "\n", "s2->v7\n", "\n", "\n", "\n", "\n", "\n", "s2->sv1\n", "\n", "\n", "\n", "\n", "\n", "s3\n", "\n", "I\n", "\n", "\n", "\n", "v3\n", "(β * S) * I\n", "\n", "\n", "\n", "s3->v3\n", "\n", "\n", "\n", "\n", "\n", "v8\n", "I / trecovery\n", "\n", "\n", "\n", "s3->v8\n", "\n", "\n", "\n", "\n", "\n", "\n", "s3->v8\n", "\n", "\n", "\n", "\n", "\n", "v9\n", "I * δ\n", "\n", "\n", "\n", "s3->v9\n", "\n", "\n", "\n", "\n", "\n", "\n", "s3->v9\n", "\n", "\n", "\n", "\n", "\n", "s3->sv1\n", "\n", "\n", "\n", "\n", "\n", "s4\n", "\n", "R\n", "\n", "\n", "\n", "v10\n", "R * δ\n", "\n", "\n", "\n", "s4->v10\n", "\n", "\n", "\n", "\n", "\n", "\n", "s4->v10\n", "\n", "\n", "\n", "\n", "\n", "s4->sv1\n", "\n", "\n", "\n", "\n", "\n", "p1\n", "\n", "μ\n", "\n", "\n", "\n", "v1\n", "μ * N\n", "\n", "\n", "\n", "p1->v1\n", "\n", "\n", "\n", "\n", "\n", "p2\n", "\n", "β\n", "\n", "\n", "\n", "p2->v2\n", "\n", "\n", "\n", "\n", "\n", "p3\n", "\n", "tlatent\n", "\n", "\n", "\n", "p3->v6\n", "\n", "\n", "\n", "\n", "\n", "p4\n", "\n", "trecovery\n", "\n", "\n", "\n", "p4->v8\n", "\n", "\n", "\n", "\n", "\n", "p5\n", "\n", "δ\n", "\n", "\n", "\n", "p5->v5\n", "\n", "\n", "\n", "\n", "\n", "p5->v7\n", "\n", "\n", "\n", "\n", "\n", "p5->v9\n", "\n", "\n", "\n", "\n", "\n", "p5->v10\n", "\n", "\n", "\n", "\n", "\n", "fs_1u\n", "\n", "\n", "\n", "\n", "fs_1u->v1\n", "\n", "\n", "\n", "\n", "\n", "\n", "fs_3d\n", "\n", "\n", "\n", "\n", "fs_5d\n", "\n", "\n", "\n", "\n", "fs_7d\n", "\n", "\n", "\n", "\n", "fs_8d\n", "\n", "\n", "\n", "\n", "v1->s1\n", "\n", "\n", "\n", "\n", "fbirth\n", "\n", "\n", "\n", "v2->v3\n", "\n", "\n", "\n", "\n", "\n", "v3->v4\n", "\n", "\n", "\n", "\n", "\n", "v4->s2\n", "\n", "\n", "\n", "\n", "fincid\n", "\n", "\n", "\n", "v5->fs_3d\n", "\n", "\n", "\n", "\n", "fdeathS\n", "\n", "\n", "\n", "v6->s3\n", "\n", "\n", "\n", "\n", "finf\n", "\n", "\n", "\n", "v7->fs_5d\n", "\n", "\n", "\n", "\n", "fdeathE\n", "\n", "\n", "\n", "v8->s4\n", "\n", "\n", "\n", "\n", "frec\n", "\n", "\n", "\n", "v9->fs_7d\n", "\n", "\n", "\n", "\n", "fdeathI\n", "\n", "\n", "\n", "v10->fs_8d\n", "\n", "\n", "\n", "\n", "fdeathR\n", "\n", "\n", "\n", "sv1->v1\n", "\n", "\n", "\n", "\n", "\n", "sv1->v4\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "Graph(\"G\", true, \"dot\", Catlab.Graphics.Graphviz.Statement[Catlab.Graphics.Graphviz.Node(\"s1\", OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(:label => \"S\", :shape => \"square\", :color => \"black\", :style => \"filled\", :fillcolor => \"#9ACEEB\")), Catlab.Graphics.Graphviz.Node(\"s2\", OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(:label => \"E\", :shape => \"square\", :color => \"black\", :style => \"filled\", :fillcolor => \"#9ACEEB\")), Catlab.Graphics.Graphviz.Node(\"s3\", OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(:label => \"I\", :shape => \"square\", :color => \"black\", :style => \"filled\", :fillcolor => \"#9ACEEB\")), Catlab.Graphics.Graphviz.Node(\"s4\", OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(:label => \"R\", :shape => \"square\", :color => \"black\", :style => \"filled\", :fillcolor => \"#9ACEEB\")), Catlab.Graphics.Graphviz.Node(\"p1\", OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(:label => \"μ\", :shape => \"circle\", :color => \"black\")), Catlab.Graphics.Graphviz.Node(\"p2\", OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(:label => \"β\", :shape => \"circle\", :color => \"black\")), Catlab.Graphics.Graphviz.Node(\"p3\", OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(:label => \"tlatent\", :shape => \"circle\", :color => \"black\")), Catlab.Graphics.Graphviz.Node(\"p4\", OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(:label => \"trecovery\", :shape => \"circle\", :color => \"black\")), Catlab.Graphics.Graphviz.Node(\"p5\", OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(:label => \"δ\", :shape => \"circle\", :color => \"black\")), Catlab.Graphics.Graphviz.Node(\"fs_1u\", OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(:label => \"\", :shape => \"point\", :color => \"white\")) … Catlab.Graphics.Graphviz.Edge(Catlab.Graphics.Graphviz.NodeID[Catlab.Graphics.Graphviz.NodeID(\"v3\", \"\", \"\"), Catlab.Graphics.Graphviz.NodeID(\"v4\", \"\", \"\")], OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}()), Catlab.Graphics.Graphviz.Edge(Catlab.Graphics.Graphviz.NodeID[Catlab.Graphics.Graphviz.NodeID(\"v2\", \"\", \"\"), Catlab.Graphics.Graphviz.NodeID(\"v3\", \"\", \"\")], OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}()), Catlab.Graphics.Graphviz.Edge(Catlab.Graphics.Graphviz.NodeID[Catlab.Graphics.Graphviz.NodeID(\"p5\", \"\", \"\"), Catlab.Graphics.Graphviz.NodeID(\"v10\", \"\", \"\")], OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}()), Catlab.Graphics.Graphviz.Edge(Catlab.Graphics.Graphviz.NodeID[Catlab.Graphics.Graphviz.NodeID(\"p5\", \"\", \"\"), Catlab.Graphics.Graphviz.NodeID(\"v9\", \"\", \"\")], OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}()), Catlab.Graphics.Graphviz.Edge(Catlab.Graphics.Graphviz.NodeID[Catlab.Graphics.Graphviz.NodeID(\"p4\", \"\", \"\"), Catlab.Graphics.Graphviz.NodeID(\"v8\", \"\", \"\")], OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}()), Catlab.Graphics.Graphviz.Edge(Catlab.Graphics.Graphviz.NodeID[Catlab.Graphics.Graphviz.NodeID(\"p5\", \"\", \"\"), Catlab.Graphics.Graphviz.NodeID(\"v7\", \"\", \"\")], OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}()), Catlab.Graphics.Graphviz.Edge(Catlab.Graphics.Graphviz.NodeID[Catlab.Graphics.Graphviz.NodeID(\"p3\", \"\", \"\"), Catlab.Graphics.Graphviz.NodeID(\"v6\", \"\", \"\")], OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}()), Catlab.Graphics.Graphviz.Edge(Catlab.Graphics.Graphviz.NodeID[Catlab.Graphics.Graphviz.NodeID(\"p5\", \"\", \"\"), Catlab.Graphics.Graphviz.NodeID(\"v5\", \"\", \"\")], OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}()), Catlab.Graphics.Graphviz.Edge(Catlab.Graphics.Graphviz.NodeID[Catlab.Graphics.Graphviz.NodeID(\"p2\", \"\", \"\"), Catlab.Graphics.Graphviz.NodeID(\"v2\", \"\", \"\")], OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}()), Catlab.Graphics.Graphviz.Edge(Catlab.Graphics.Graphviz.NodeID[Catlab.Graphics.Graphviz.NodeID(\"p1\", \"\", \"\"), Catlab.Graphics.Graphviz.NodeID(\"v1\", \"\", \"\")], OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}())], OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(:rankdir => \"LR\"), OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(), OrderedCollections.OrderedDict{Symbol, Union{String, Catlab.Graphics.Graphviz.Html}}(:splines => \"splines\"))" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "GraphF(seir)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 1. Measles model \n", "## Time unit: Month" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 89070.0\n", " :E => 0.0\n", " :I => 930.0\n", " :R => 773545.0" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# define parameter values and initial values of stocks\n", "# define constant parameters\n", "p_measles = LVector(\n", " β=49.598, μ=0.03/12, δ=0.03/12, tlatent=8.0/30, trecovery=5.0/30\n", ")\n", "# define initial values for stocks\n", "u0_measles = LVector(\n", " S=90000.0-930.0, E=0.0, I=930.0, R=773545.0\n", ")" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deZgU1b038N+prfdZevaFYdgXNxAQXJDFhUXFBVBxSUASMTEaE2OSq14To8+9Ro2aqImvGp9cI2rEGFEUAQmCS0BFYhAB2WZh9n16q67tvH/U0Iysg8LUzNT38/DwdNVUV52uPl3fqnNOdTPOOQEAALiV4HQBAAAAnIQgBAAAV0MQAgCAqyEIAQDA1RCEAADgaghCAABwNQQhAAC4GoIQAABcDUEIAACuhiAEAABX61NB+MQTT+zatcvpUjjMNE2ni+Awy7LwxYGoBpxzy7KcLoXDUA2IqCtHgz4VhG+99db27dudLoXD4vG400VwmKZp+PyjGhiGoWma06VwGKqBZVmGYRx1sT4VhAAAAMcKQQgAAK6GIAQAAFdDEAIAgKshCAEAwNUQhAAA4GoIQgAAcDUEIQAAuJrkdAGOp0k500Jri9Z/uJWImMhEjyD5RNEryAFRDklyQFLSZSVN8mTInkxZ8olOlxcAAJzXp4Lw05aPTrti+AXTLiAibnBTs4yEaSQsI2boUUOPmtGKuNZuqM261qZzi3syFW+W7M1SvFmKL1vxZivesMJE5vTrcAVucT1qaq26HjOMhGUmTUvnlrH/O7Ekn0iMST5BVATRJ0o+UfaLUkBkAt4ggONA1/X//Oc/Tpfi2IwcOdLn8x331fapIIwZEa5YqUs9+YgLG6qZbNbVJk1t1hINyZatkUSjprXqSrrsy1F8uR5frseXo/hyPJ6MI68JuiTZqkfKE9G9iXiNGq9Paq265BeVNFkOSZJXEL2iIDNB2t9WrzbpxLmdkYZqGXHTiJt6zJADkhwUPRmyHJI8GbKSJnsyZVzln2iWwY24aSRMM2mZqkmcDHX/99iJHpEJJEhMUATRI4geUfQKooKelx7t3XffnTdv3uDBg50uSFeVlZU99thj11133XFfc58KwmMieUWpUAwUejvP5CZXm/VEQzLRkIxVqY3/bkvUJU3N8uV6/Lkef77Hl+Px53u8Wbhw7BIjYbZsi7Zui7bujFoGT+vvDxZ788Zn+vO+6T7kpEcNLWpobYYeMZItWqxGbdkaUZv1ZKvGiHnCspIhecKyP8drX+t7sxRBwpt1DJIterxWjdcn1QZNbdaSrbrWZpiaJQdEySeKHkH0imRfsu9jJk1ukWVwS7PMpGmqlpG0uMllvyj5RSkgyUFRSZPkoKSEJCVNVtIkJV1WQhLhnfnWLIPrUSN1mmLplqla3OJEFI8k2iSVGJO8HSclokdgIpO8IpNYrC5xzjkTly1709HiH4MFCxZ05YtDvwH3BuEhMZH5chRfjkIUSs00VDNRryXqkvG6ZP2nrfG6ZLJV94Zlf57Hl+vxF3j9uR5fngdH2xQjYTZ+3t7477ZIRTx9UCBzeKj4vGxfruc4rJqRHJLkkBQoOPR21WY9WhvT28xEQ7JlW1Rt1JItmhyUvNl267fHbgP3ZSsCrlf2MeJm+554pDweqUhEKxOCzPwFXn+ex5/vyRwZ8oYVJU2S/Md8td1xERkz9LipRw09YuhRM1qZ0CIRrc3Q2nU9Zir2ZX2G7EmXPJmKJ7Pj+l4O4tB0IFO1EvXJRGNSbdTUZj3ZqmttutZmWIYlBSQ5IIpeUfQIoiKIXkGQGREzTFPwWZxztbGj08FMmtwkI2FaBq/+d5Oluf0HOmyobUcnecVQiS9Usr9hmps8UZ+M1yfjtcmmze1765OJBk1Jk/z5Hn+e15/v8ed7fLketzUNcYu3bIvWf9LS+lUsY2gw/6zwyBtKujNvJJ8YLBKlLJIkSZL21W1OaoumNmlqo5Zo1CLl8USDpjYmO9Ixx+PLVrw5HQ/cc6FvJq22nbHWHdG2HbFkix4q9YX6+4smZQX7+Y5XCAkSU9IkJe2wa+Mm1yJGslXXWvVkq642aW07Y8kWLdmqmzr3ZsqesOLJlL1h2RtWPJmyJ6woIbccsrjFE3XJaLUar1aj1WqiLmkkTF+Oxz6NC5X6c8KyPfrvCD0CkUgkFAod7q8VywrE/+euY9ThuKVWHV9MZP4Cr7/AS6d1zOEWV5u0eG0yXpds2RatXtsYb9CUoOTL8wTyPb48jz/f68/ziJ6+We20Nr12fUvdhhYlXc4bnzn4qiLJ22O66xh5w4o3rNCQTjM5JVv1RGNSbdASDVrbrliiUUs2a3YPsTfb489VvNkeX47iyZT70vCcWLXasi3Ssi0arUyESvwZwwKDryoKFnsdeY1MZJ4M+ZB98KZmJZt1tVlLtujJZq2xsl1t0ZLNupm0vGHZE1Y60rHjf0UO9Jj69o1xitcno5WJaGUiUpmIV6tKhhwo9AaLvEXnZvnzvZ5MDFY4URCExwcTmC/H48vxZJ2ybxYntVmL1ybjtWrbzljNh82JuqQUEH25djR6/Xkef57nG7Q49SCcWrZHa//V3LYrlnN6+sjv9w8UeI/+rJ6Akd0K1zkducWTHT3EWrwu2fRFRG3UtHbdE1Z8OYov2+PLUbw5ii/b48mQe1Hnlpm0Wr+KtmyNtGyLMpFljggVT8lOHxToyY3DoiLYLSsHzE8FpNqkJZv19vJ4sklTm/cNArczMlP2hGUpJEhpAvXs+qg2aanki+5NyEEpVOIL9vNlnZoWLPb11fPmHghBeMIwskdqhE/a1zTBSW3REnXJeG0yUh6v29CSqE8yifnzPHYvo93p2CsGqWoRo/7jltr1LZJPzD8rPPTa4j7QDswE5s1WvNlK5oj9My2Dq01aoj6ZaNSiexMN/25TGzU9anizlFS/Y8eonLDcc1pWucUj5fHWHbHW7dHo3kSo1B8eESqaku3LOR49tc45bECqln35qDZryWatvTxuT1oat0cUe9JlT1hW0js6I79Zr+e3ZCatRH0yVqPGa9RolRqrUkWPEOznC/bz9TsvJ1jiw7BnpyAIu9G+NrrMEftb7bV2I16nJuq0eJ3a/EV7vC5pJjsGqfryPL4cj30vRw8ZiWMZvHlLpP7Tlvbd8ezT0oZf3y9Ycvzv6elRBPtMJe9rR17L4GpDMtGoqY1arEZt2tyuNmlamy6HJG+W4g0rqeY7T6aipEvd0/BoqlakMhEpj7fvibfvjnmzlIwhgeKpPf3i77gQvUKg0Nt5ELiu66ZpKqKiNnf0QSZb9Uh5vKnNHmZiWAa3uzCVkKSky1JAVEKSHJQkvygHJMkvSj5BkL/JfjNUU2839KiRbOno+1Sb9ERjUo+Z9uBzf4G334hQoNCLMUHHxNKsihX1m+N7BFkQFUGQmSALxOiQHTHGvpt8QgN9RWdnH3nNeBscpqRJSlowo1MD3f5BqvXJhk2tifqk2qgpabI3W/Hl7mug6957/w3VbN0ebdocadkaCfbz5Y7NGHZ9vz5wCfiNCdK+TuJOuMVTR71ks9a6PWpflGgRQwlJnkxZSZOVdMkeEmkff6WgJH+zrwjgpEWMZLOmNmvxei1Rq0arVT1iBIp8of6+/AmZhbMyMvPSj9sL7rUEWTj4PMZm6ZY9eFWLGFq7oUeN6F5Vjxp6zDRiHXcjEJF9u4joEZjAUheR9iQ3ualZRB33jRiqaaqWkTDtUUJySPZkyJ4MKVjsyx6V7s1WvJlKL2pRd0RFRUVdXV1WVtbAgQMP/qsgC3mnZZRcmGvplqlZlsYtwzrgllYbE5k3LBMjQWbefOWo20UQ9jiHGKRq8WSLnmjQ7O6rlm377v1Pk71ZsieseLPswXWKJ0PmAj8uxdCjRqQiESmLt+2KxarVtIGBrJNDAy/Lxzns4TCBHWJUDhG3uNZuJFt0rU23H0SrVL1d1yKmETP0mCkqgn39IXpFQRFERWAiiR6x8xpM1bIMy0xahn0rQsyU/aInU/ZmKb5cT/ao9JIZeb4cJZWpkUikO197byTIgt0SfoRlLIObqmkmLTNpcYsb8Y4Drj0pSPZdCh1f6Ch6RMkrSH6x57SQ9yLRaPTyyy/funXr8OHDy8vLGWNfffXVgQsx8mQq6YMDXV+tZVmmeWBMHgwHtV6ACczuhcocHkzNtNNRbdLUZl1t0lq2RpMtHdcfckC0x1Xb9y9LAVH2S6JXEL1Cx/eW7bu71j6ltTTLSFh63NBajWSLpjZp8bqkpfNgiS+tv79kem5aqf+btREBETHhsGMjbYZqGjHTUC1TNS3NMjWLm2Qm9396BUkQFCZIgv31uXJQlIMSjrbdQJCYEJTk4NGXhG/p8ccfV1W1rKzMvvFp+/bt3bl1BGFvlUrHA+a3t7V7yZdsM/SIoUXsr1/RY1WqqVpm0jISJufcVC0iEmRGnARZEGQm+UQ5ICnpHc04vWXMTt8gecUedLcJgBNqampKSkpSt/8OGzasO7eOIOxrmMCUkKykI8YA4NjUJ+g/zcenb+XIPCJNzP9ak8bVV199/vnnX3jhhVOnTj3//PPHjBnDWPe1eSAIAQCAiOijeuvJL7vjS9c8Ao3L+Vo7yFlnnbVly5aXXnrpvffe+81vfjN58uSlS5fKcjed0CMIAQCAiOiy/sJl/R0bDTBgwIA777zzzjvvrKioOO20015//fW5c+d2z6YxAgIAAHqQkpKSoqKi9vb2btsirggBAMBh//Vf/6Wq6sSJE9PS0t58882amppp06Z129YRhAAA4LD58+cvWbJkyZIliURi6NChH3/8cXFxcbdtHUEIAAAOGzZs2N133+3U1tFHCAAAroYgBAAAV0PTKAAAOKypqSkajaYmFUUpKCjotq0jCAEAwGG333770qVLs7Ky7MkhQ4YsX76827aOIAQAAOd973vfe+ihhxzZNPoIAQDA1XBFCAAARERmW5O258tu2BDzeL3Dx9LXv1b73Xffvfnmm+3HY8eOXbBgQTeUxIYgBAAAIiK9anf83+u6YUNMlDyDT2Py135FLicn59RTT7Ufl5SUdEMxUhCEAABAROQdOc47cpxTWz/ttNMWLVrkyKbRRwgAAK6GK0IAAHDexx9//Nvf/tZ+7PF4brvttm7bNIIQAAAcNmvWrM8++ywSidiTuq5359YRhAAA4LArrrjiiiuucGrr6CMEAABXQxACAICrIQgBAMBh1dXV9fX1Tm0dQQgAAA678847nfqiUUIQAgCAyyEIAQDA1XD7BAAAEBHVxeq3Nu3ohg1JgnR28RmM2NEX7RYIQgAAICIqb9+7pvyDbtiQX/aNLzhdFuVu2FZXIAgBAICI6IyC088oON3pUjgAfYQAAOBqCEIAAHA1NI0CAIDD5syZ4/f7ndo6ghAAABx28cUXO7h1NI0CAICrIQgBAMDV0DQKAAAOKysra2pqIqK0tLT+/fsritKdW0cQAgCAw+699941a9YMGTKkurq6qalp8eLF5513XrdtHU2jAADgvCuvvHLVqlVbtmxZuHDhz372s+7cNIIQAAB6kFGjRtXU1HTnFtE0CgAARERqo9a6I9oNGxI9Ys7o9AO+c7u2tnbjxo1tbW0PPfTQpZde2g3FSEEQAgAAEVGyVY/uVbthQ0xkWaemCdLXknDlypVffvlldXV1OBz+3e9+1w3FSEEQAgAAEVH64ED64IBTW//Od77z4IMPGoYxa9ase+6555FHHum2TaOPEAAAegpJkh5//PE//elPO3Z0xy8j2hCEAADQgwwaNGjevHn3339/t20RTaMAAOCwuXPn+ny+1OSvf/3rV155xTAMSeqOkEIQAgCAw2bOnNl5sqSkpDtvJUTTKAAAuBqCEAAAXA1BCAAAroYgBAAAV0MQAgCAq/WOUaMbNmx477337MdXXXVVaWmpk6UBAIA+pHcE4dq1a3ft2nXBBRcQUed7TQAA4BsrLy9/+umnnS5FV23fvn3SpEknYs29IwiJ6NRTT73ooov8fr/TBQEA6AtGjRp19tlnb9y40emCdNWoUaPOPPPME7Hm3hGEkiS9+OKLL7zwAmPslVde6devn9MlAgDo3YqLi//0pz85XYoeoXcMlvnJT37y0UcfrV+/ft68eb/61a+cLg4AAPQdxxaEGzZs+L//+78VK1bE4/GD/2qa5sqVK59//vnKysqur7OpqWnt2rV79uw5YFWrVq16/vnnKyoqiIixjp+tOuOMM45p5QAAAEfW1aZRwzCuueaaTz/9dOLEibW1tdXV1QsWLOi8gGVZF110UVNT0ymnnPKTn/zktdde60qv5ty5c998801Zln/xi1/cfffd9kzO+SWXXFJfX3/qqaf+9Kc/feWVV9ra2gYNGmRZ1q9+9asZM2Yc64sEAAA4nK4G4RNPPLFr167NmzcHAof+2caVK1du3779yy+/9Pl8Tz755N133/3+++8TUSQSWbJkyQ033JBa8tlnn503b569ngcffPDFF1+85pprOq/q3Xff3bJly9atW/1+/1NPPXX33Xf/8Ic/fOihh0RRvO6666699tpv+FoBAAAOwjjnXVnujDPOuPnmmydMmFBbW3v66aeHQqEDFrj55psFQXj88ceJqK6uLj8/v7m5OTMzs7Gx8bzzzps9e/Y999xDRHfdddc777zz7rvvZmZmpp47d+7c0047LXVFeOuttxqG8cc//pGIGhsbc3JyGhoasrOzj1rIkpKSzg2n/fv337x5c1deXV8SjUaDwaDTpXCSqqqSJHXPr7f0WKgGuq6bpun1ep0uiJNQDSzLkmX5qLcbdPVgsXv37ueff/6ZZ55JT0//7LPPli1bNmbMmM4LVFVVjR8/3n6cl5enKEpVVVVmZmZ2dva777573nnnEVEymXznnXdWrVrVOQUPVlVVNXr0aPtxdna21+utqqrqShCedNJJTz311AE/5+FCB5+muIosywhCcn01QBDaXF4NLMsyTfOoi3X1YJFIJLKzs//2t78R0a9+9as77rjjn//8Z+cFTNMURTE1KQiCYRj245ycnNWrV48cOVJRlM2bN4fD4SNv64BViaKYWhUAAMDx1dVRowUFBZMnT7YfT5069eAmx4KCgvr6evtxW1ubqqqFhYWpvz722GP9+/fPzs5+8sknu7Kt1Kqi0WgsFuu8KgAAgOOoq0E4derUXbt22Y937NhRXFxsP25pabEvPCdNmvTuu+/aPY4rV64cOXJkbm6uvcw999yzfPnyFStWvPfee0uXLv3Nb35z5G3Zq7Isy17VsGHD8vPzj/2lAQAAHF1Xm0Zvv/32iRMnBgKBtLS0Bx980L6w0zQtHA5/+umnY8aMmTNnzv3333/ttdeOGzfuwQcffPjhh+0n7tmzZ+3atatXr7b7Bd95553LLrvshhtusKP0pZdeeu+99zZu3FheXl5ZWXn99defc845V1xxxX333XfttdeOHz/+wQcffOCBB1L3EQIAABxfXQ3CYcOGrV+/fvHixbFY7K233ho7diwRSZL01FNP9e/fn4g8Hs9HH3303HPP1dfXv/TSS6l21AEDBqxduza1nuzs7A8++CA1WVJSMmbMmNS4G3tEjKIoH3744XPPPVdXV7d48eIpU6YchxcKAABwKF29faJXmDFjxi233OLyUaORSMTl48Rw+wShGmDUKBGhGuwbNSrL8pEX6x3fNQoAAHCCIAgBAMDVEIQAAOBqCEIAAHA1BCEAALgaghAAAFwNQQgAAK6GIAQAAFdDEAIAgKshCAEAwNUQhAAA4GoIQgAAcDUEIQAAuBqCEAAAXA1BCAAAroYgBAAAV0MQAgCAqyEIAQDA1RCEAADgaghCAABwNQQhAAC4GoIQAABcDUEIAACuhiAEAABXQxACAICrIQgBAMDVEIQAAOBqCEIAAHA1BCEAALgaghAAAFwNQQgAAK6GIAQAAFdDEAIAgKshCAEAwNUQhAAA4GoIQgAAcDUEIQAAuBqCEAAAXA1BCAAAroYgBAAAV0MQAgCAqyEIAQDA1RCEAADgaghCAABwNQQhAAC4GoIQAABcDUEIAACuhiAEAABXQxACAICrIQgBAMDVEIQAAOBqCEIAAHA1BCEAALgaghAAAFwNQQgAAK6GIAQAAFdDEAIAgKshCAEAwNUQhAAA4GoIQgAAcDUEIQAAuBqCEAAAXA1BCAAAroYgBAAAV0MQAgCAqyEIAQDA1RCEAADgaghCAABwNcnpAnTJihUrXnvtNfvxbbfdNmLECGfLAwAAfUbvCMLPP/88Kytr9uzZRFRUVOR0cQAAoO/oHUFIRIqiKIoydOhQj8fjdFkAAKDv6B19hFlZWTt27HjggQdOOeWUL774wuniAABA39E7rggXLly4cOFCInrxxRfvv//+l19+2ekSAQBAH3HMQbh06VLDMOzuugO0t7cvXry4rq5u+vTpEyZM6Mraksnk559/vmXLllGjRo0ePTo1PxKJLF68uLa2dtq0aWeeeWZqfnFxcWtr67GWGQAA4HCOLQjff//96667bvDgwQcHoa7rEydOHDRo0JgxY2bNmvXHP/5xzpw5R13hjBkz9u7d297e/qMf/SgVhIZhnHvuuf379x83btxll132hz/8oaGhoaioyLKs+++//6c//ekxlRkAAOAIjiEIE4nEj370o5///OepOxk6+8c//mGa5quvvioIwsCBA++77z47CGtqap544on77rtPEAQiMk3z7rvv/ulPf5qTk0NEy5cv93g8c+fO7byqpUuXJpPJv//976IoDhky5De/+c3jjz++bt06QRCeeOKJs88++3Al1HW9trZ29+7d9qTX6y0sLOz6CwQAABc6hiD87//+7+uuu+5wdy+sXr162rRpdtrNmDHjmmuuaWhoyMnJCYfDn3/++cKFC//85z9zzufPn9/Y2BgKhexnHXII6OrVqy+88EJRFO1VXXXVVSNGjJgyZcpRS7hnz54777zT5/PZk3l5eatWrer6C+wbotGo00VwmKqqkiRJUu/o/z5BUA10XTdNU9d1pwviJFQDy7JkWZZl+ciLdfVgsWHDhjVr1qxfv37JkiWHXKCmpiZ1rZaRkeHxeKqrq3Nycjwez5IlSy699NLvfe97yWSyra3t9ddfP/ItEDU1NePGjbMfh0Ihv99fXV2dn59/1EIOHTr0lltumTlzZhdfVF+VOs9wJ1mWEYTk+mpgB6HX63W6IA5zeTWwLMs0zaMu1qWDRTKZXLRo0Z///Ocj5KogCJZlpSY55/YlHRH5fL6///3vAwcOVBRl586dR70R8IBVWZaVWhUAAMDx1aUg/Ne//lVWVvbLX/6SiGpra8vLyy+44II333yz89lWYWFhTU2N/bihoUHTtFT/nGmaN9100/jx4zVNu/nmm5999lm7BfVwOq+qublZVVV09QEAwAnSpRvqx4wZs3r16gceeOCBBx64+uqrCwoKHnjgAUVROOeffvppPB4nopkzZy5fvlzTNCJ6/fXXJ0yYEA6Hicg0zfnz5zc3N7/66qtvvvlmfX39woULO1/wHcxeVTKZtFc1btw4e2QNAADAcdelK8JQKDRmzBj78fbt2/1+vz2padq4ceM+/fTTMWPGzJgx4+GHH546deqoUaNeeumll156yV6+qqpKkqRUv+CSJUtuuumm2tpa+yLv8ccff+ONN/7zn/9s2rRp7dq1P/nJT2bOnDlt2rTS0tKpU6eefvrpL7744uLFi0/ISwcAACBinPNjeoLdNDp+/Hgi4pyvXr16/Pjxdn+srutvvfVWU1PT5MmTBw0a1JW1bdu2be/evanJESNG2KNSdV1/++23GxsbJ02aNHjw4C6WbcaMGRgsE4lEXN49jlGjhGqAwTJEhGqwb7DMcRs1mpKfn58awMkYO//881N/kmX5sssuO6a1DR8+fPjw4QfPl2X50ksvPdayAQAAHKve8aXbAAAAJwiCEAAAXA1BCAAAroYgBAAAV0MQAgCAqyEIAQDA1RCEAADgaghCAABwNQQhAAC4GoIQAABcDUEIAACuhiAEAABXQxACAICrIQgBAMDVEIQAAOBqCEIAAHA1BCEAALgaghAAAFwNQQgAAK6GIAQAAFdDEAIAgKshCAEAwNUQhAAA4GoIQgAAcDUEIQAAuBqCEAAAXA1BCAAAroYgBAAAV0MQAgCAqyEIAQDA1RCEAADgaghCAABwNQQhAAC4GoIQAABcDUEIAACuhiAEAABXQxACAICrIQgBAMDVEIQAAOBqCEIAAHA1BCEAALgaghAAAFwNQQgAAK6GIAQAAFdDEAIAgKshCAEAwNUQhAAA4GoIQgAAcDUEIQAAuBqCEAAAXA1BCAAAroYgBAAAV0MQAgCAqyEIAQDA1RCEAADgaghCAABwNQQhAAC4GoIQAABcDUEIAACuhiAEAABXQxACAICrIQgBAMDVEIQAAOBqCEIAAHA1BCEAALgaghAAAFwNQQgAAK6GIAQAAFfrHUH4yiuvXLDPZ5995nRxAACg75CcLkCXlJWVzZgxY8GCBUQUDAadLg4AAPQdveOKkIh27969bt26RCIhy7LTZQEAgL6jdwThoEGDwuHwxo0bJ06cuHbtWqeLAwAAfUfvaBqdPXv27NmzieiMM8549NFHJ02a5HSJAACgj+hqEK5fv37ZsmUVFRXFxcWLFi3q37//wcuUl5c//vjjdXV1M2bMuOaaa7qy2j179mzYsOGLL76YOnXq1KlTU/MrKyv/8Ic/1NbWTp8+/ZprrmGMdRRXkizL6mKZAQAAjqqrQXjrrbdOnz592rRpH3zwwejRoz///PN+/fp1XiAWi51zzjlz58696KKL7r777paWlptvvvmoq/35z38eiUR27Njh9XpTQRiPx88+++wrrrjioosu+u///u+mpqba2tpwOGxZ1jPPPPPEE08c64sEAAA4nGO4IhQEgYiuvfbajRs3Llu27Ac/+EHnBV566aXi4uJHHnmEiEKh0K233vqDH/xAEITy8vJf/vKXzz33nM/nI6JYLLZgwYJHH320qKiIiJYsWUJEc+fO7byqv/3tb/n5+Y899hgRZWRk3HTTTcuXL9+4cSNjbPXq1SUlJYcrZDQaff/997OcHDUAACAASURBVKPRqD2ZlZU1ZcqULu+KPsKyLJdfNFv7OF0QJ2EPoBoQqgFRF19+V4PQTkF7vY2Njbm5uQcs8K9//SvVdTd58uTdu3fX1dUVFBSUlJRkZGRMnz797bffZoxdcsklgwYNKigoOMK2Oq9q0qRJ5eXloVDouuuuO2ohm5qa1q1bt337dnsyGAyOHz++iy+wz0gkEqIoOl0KJ6mqKkmSJPWO/u8TBNVA13XTNF0eA6gGlmXJsnzUew2O+WDxP//zP+np6ZdeeukB82tra4cPH24/DgQCPp+vpqamoKCAMfbkk08uXLjw8ssv1zRtxIgRf/zjH1N9fodUW1s7cOBA+7HP5wsGgzU1NcXFxUctW//+/W+55ZaZM2ce64vqSzjnLr/VUtrH6YI4CdXADkKv1+t0QZyEamBZlmmaR13s2A4Wzz333NNPP71u3bqDjzI+n0/TNPsx51zXdb/fb08KgvD73/9+0KBBiqLY14VH3krnVRGRpmmpVQEAABxfx3Af4eLFi++555533323tLT04L8WFRVVVFTYj6uqqizLKiwstCdjsdisWbPmzJkzffr0yy67LJFIHHlDnVdVU1Oj63pqVQAAAMdXV4Pw1VdfveOOO1asWDF06NDUTMuy/va3v7W0tBDR7Nmzly1b1traSkQvvPDC+eefn5aWRkTxePySSy4ZOnTok08++eyzzw4aNGj69OmxWOwI25o9e/Zbb71lr/aFF16YOnVqZmbmN36FAAAAR8K7JhQKZWRkDNznoYce4pwnk0ki+vTTT+1lrr/++oEDB06fPj0vL2/jxo32zLq6uvvvv9+yLHvSNM177723sbHRnrznnnsGDhwYCATC4fDAgQNffvlle/78+fMHDBgwY8aMvLy8Tz75pIuFnD59+ltvvdXFhfuq9vZ2p4vgsEQioeu606VwGKqBpmmJRMLpUjgM1cA0TU3TjroY45x3JS/37NnTecmMjIxwOGzPLyoqUhTFnr958+bGxsaxY8eGQqGurLaxsbG9vT01mZOTk3riF1980dDQMGbMGPvKsitmzJiBwTKRSKSLO7+vwqhRQjXAYBkiQjXYN1jmuI0aHTBgQFfmn3LKKV1coS07Ozs7O/uQfzr55JOPaVUAAADfQO/40m0AAIATBEEIAACuhiAEAABXQxACAICrIQgBAMDVEIQAAOBqCEIAAHA1BCEAALgaghAAAFwNQQgAAK6GIAQAAFdDEAIAgKshCAEAwNUQhAAA4GoIQgAAcDUEIQAAuBqCEAAAXA1BCAAAroYgBAAAV0MQAgCAqyEIAQDA1RCEAADgaghCAABwNQQhAAC4GoIQAABcTXK6AHBCqCbVJ3ijSk1JatV4a5LiBqkmcaJWjRNRusIEoqBMaQpleVihn4oDLNvrdLkBALodgrAXS5pUFuVlEaqI8soYL49SdZxXRjy1qq6alONl2V7K8lDYw9IV8kvkk0hklKkwImrTuEVUEaN2jZqSVnWMKmPc5DQsnZ0SZqdnsQm57LQsJjKnX6T7mJzatY7HGR7COwBwoiEIe4eITjva+I52vrOddrXzXe18V4QaVd4vwEpD1D/I+gXY1EIq9AsZpA7JCWYo32QrLUna1sb/08w/beB/3GpVxfikAmFGP3ZxP1YUwAH5W2lOUnWcV8epNs7rVaqL88YkNanUnORtGkV0atd53KCkSZJAIbnjWa1J4kQBiYIyZXpYlofyfKw4QP2CbGCIhqazoelMQf8GwLeDIOxxDIvKonxbK21v49vb+Fdt/Ks23q7R4DQ2JJ0NTqOz8tj1Q4RBISoOMOGgeIpEeOgbpSARZXrozFx2Zi5bNJyIqD5Bq6uttyr5XZ+YwzPYVQOFqwcJOWg+PTzNooooL4tQeZSXR3lFlCqifG+c9sa4V6RCPyvwU4GP5fgoz8dGZlKWhzI9QoZCIZnSFeaTyCseYrVRnaIGtSR5k0p1CV4Vp7II/7CWtrdZZVE+IMhGZ7Nx2Wx8Ljs9G7l4bKI6Vcd5g0oNKm9SqVGl5iRv0ag1SW0ajxqkmtSaJIuoTeMHPz1TYaJAaTKlK5ShsGwv5fkoz8f6Bag4wEpDLOzp/tfkCvZZIxGpJiXMjreGEWUojOhrZ5NdgSB0WMyg7a18Wxv/soVva6NtrXxXOy/0s2EZNDydjclm8wYJQ9Ko2IkLslwfzRskzBtEuiW+W8Vf2mXds1GfVCAsHMZm9hNc3mranNx/ab6rne+O8D0RqkvwQj8bEKL+QdY/yCYXUElQKA5QvwDzfYuPWlCmoEz5PnuPf22/axZtbeWbGvknjfz/dlg72vnYbDalUDgzQ5gSIAmhSKRbVBPnFVHaG+NVcdob49Vxqo7zmjjVxLnAKN/Hcn2U7WXZHsryUpaHDUqjTA+lyUJQJq9ImZ79R1hRIM7J2peJLRo3LIro1KpRa5I3qFSv0petfMVeqohaZVEuMhqcxoZnsBEZ7ORMdnImlYbc/ck5jKhOdYmOM5LmJDWp1JTkzUlqSVKrxtu0jtPBNo0nDFJNss8aPQKJAvPsq+ep8xXdoqhOPokCEi0aRvePO8rWGeeHOM3ppWbMmHHLLbfMnDnT6YIcVkuStrbyL1v51taO5KtP8KHpHZ+TERk0LJ0NS2eeQ10WdFEkEgmFQsevyF8T1WnJHuuZbVZljL4/XPjeMFbo73GfalVVJUmSpON2ktec7GiX/qqN72qnne18Zzu3OA1KY4PS2KAQDUpjA0JsYIiKA8zZ7Ino9EEtX1Njraw0ymPClEJhejGb0Y/1c0HLdlTv6DIvj/KKKC+LWBVRXhkXGhI8z8dKglQcYEUB6hdgBX4q8rN8PxX4WeAEXws0qrSznW9t5Vtb+eZmvrmFYjofncVOz2ZjstnYbDYk/QS+NSf0aHBMLE71KtXEeVWMahO8Kkb1Kq+KUX2C1yaoNsEFohwfy/NRtoeyvCzLQ1leFvZQhkKZHpauUFCioEwZCvOK1MXTyrhBEc2SyQz7j3J5iCA8gRpV2tLCt7Zy+/8vW3ncoOEZ7KQMNjyDjcxkIzKoNHiI5s1vo3uq/uZm/tQ26+Vd1tRC4aYRwtRC1nMOtN8mCFN9sV+10VdtfEc739HGLU5D0tmQNDY4jYaks8FpbHBaTx9hG4lEElJoVZW1fC9fsdcq9LOZ/diMfsJZuQ5H9bcXN2hPhJdFaU+El0V4WZTKIrw8yuMGlQY7usxLgqzYZxb5rCFhpcDfs8Z8Nar0WRP/rJFvbOSfNPB2nY/LZuNyOv4d3zPLbg7ChEGVMV4Tp8oYr45TVYzvjVFNnFfGqC7BszyU72dF/o7/c32sKEC5XpbvpzzfiTojsSzLNE1ZRhB2l7rE12JvSws3OI3MYCdlshEZ7KRMNjyDuuHEvDurfkSnxTutP221VJO+P1z4zmAh19c9Wz6SLgZh0qSd7fyrNr6jff8FX7vWkXlD0mlIGhuWwYb0+Mw7pM7VwOL0cQN/q9JaXsl3R/h5hcKFReyCItbD2+hUk1I5Z1/q7YnwsiiP6FQaZANCVBpi/YOsNEgDQqw0xA7ovdZ13TRNr7env3kNKn3cwD9psD5p4J80cFlgY7LZmGw2OotOy2L9g9/qPTruRwNOVJeg2nhHwlXF+d4YVcd5ZZSq4zxhUnGAFfioJNhx2V0coAI/6xekfB+TnTgDQxCeWHtjfGsrfbnvUu/LFs4YpTLP/j/fiVRwpDHkX/X8mW3W6+XWufnC9YPZRSXCIQd9dI+DgzBp0p4I39lO9hXezna+o51q47w0xIamsyFpNDSdDUlnTvXFngiHqwZ1CVqx11pZxVdXWX6JTS5g5xawM3PZ0BPZQHdkLUmqjPHyKC+PUsW+/8sivEUjO+dKQ6w01PFgQKirH6veEoQHKI/yjY18YyPf1MQ/b6K4wU/KZCdlsmHpbFgGGxiiASHW9Q/XNzga2H119WpHo6Ude9VxqkvwvTGqT/BMD+X7OhLObm0u9LN+ASoMsKyeNzIIQXjcmJzKInxrK21t5dta+ZZWvq2V+8SO2Bu5L/Z6yHBKB3sFIjr9fY+1eJf1aQOf3k+4tIRdWCx086i5Vo22NyXLY0JZXLTHsOxsp9o47xdkg9NoSBpLXfD1D/asFrPjqyvVYEsLf6+Gf1jHP6rjEZ2PyWajs9jJYTYigw1NZ2nHMujuyFST6hK8Ok71Cb43RnUJXhGlqjivilFFlMsC9Quy/kEqCbKSACsJ2pd6VOD/Vu3tvTQID9CUpM3NfFsr39bGt7fy3REqj/JMhYoDrDDAcryU5+u4UThdoZDMJEbpCqV6W+LxuN/vp31jLJMWxQ0e0SmmU7tO7Rpv0aglSc1J3qRSU5IaVS4yyvOxXB/leFmBn/J9lOtjxYGO//N8vWxYMoLwG2rXyb5jYWsr395K29v4jnae72PDM2ik3beXwUZksMyed+5j6wnd4/UJWlpuvVnB19ZYwzLYufnsrDx2etZxa4vjRLVx2hvje2O8IkrlUV5m/x/hukUDgjTQHsaSxgalscFpVOL0GJbud6zVoC5BGxv5583cPuzuaOdekUqCrMjPcn2U66UMD8tQKHWDhyyQyEg1O57eppFuUbtGbRpv16lVo5Ykb0pSo0p1CW5YlONlRZ0OpiUBKgwwezztMQ1z77q+EYQHS1X+mjivS1D9vltR2zVq17nJqU3bP6jVNE1RFCWBFIE8IikCBSQWkikgU1CidIVleijTQ2GPPTiFsj3fanhzD9TFIOxbL/oYRfSOEfA722lnO9/Rzre38ohOwzLY0HQ2PJ3NHkDD0oVh6X2tcpxouT76/nDh+8NJs8QN9XxdLf/LV/zWJiuq8+EZbHBax/l+jpeyPCxNIZFRmkyiQJQ6dTUpblCLxts1ak5SU5I3JKhe5fUJqopTfYJneexBgB1XEmflUWlIKA2ybG+qadS5xtleKM9HM/uxmf32n6nUJ+whD7w2Tg0qNal8V/v+m7c0iyy+/67HNIUUgdIUSpNZYYDSFcpUhCwvZXmowH+ios6dGFGBnwr8jLrwpUORSMLx0+Jeoe8f4JMmNah8b4xqE7wiSpVRXh6l8ijfE+Fxgwam2YMAaVwOu2awY3fs9VWKQBPz2cT8jl3anKTtbXxnOy+P0LZWvi5BLZpln8C268SIBEayQIpAikABmTIUliZ3nLEOSaMcL8vzscKAYx3vrpLro1wfG4OveAMX6FNBuHvo7HubRz2xwkia1GZfSahcNSnXxwr9VODvuJdobA71DwoDQiyvBwxxdJXwvm+ucbogAAD79akgLKx8b/Lk08afVCAL9m2YlOVh6d/0+8YAAMAN+lQQehNN4zwN04txwQEAAF2FnhYAAHA1BCEAALgaghAAAFwNQQgAAK6GIAQAAFdDEAIAgKshCAEAwNUQhAAA4GoIQgAAcDUEIQAAuBqCEAAAXA1BCAAAroYgBAAAV0MQAgCAqyEIAQDA1RCEAADgaghCAABwNQQhAAC4GoIQAABcDUEIAACuhiAEAABXQxACAICrIQgBAMDVEIQAAOBqCEIAAHA1BCEAALgaghAAAFwNQQgAAK6GIAQAAFdDEAIAgKtJThegSyoqKlavXh0Oh2fMmKEoitPFAQCAvqMXXBF+8MEHo0aN+te//vXwww9PnTpV13WnSwQAAH1HLwjCX//613fdddfTTz/9z3/+s7m5eenSpU6XCAAA+o6e3jSaSCT++c9/PvPMM0Qky/KsWbPeeuutOXPmOF0uN+LJBDdN+7GlxohbBy8j+IJEjIhIYII30J3FAwD4Znp6ENbU1HDOCwsL7cnCwsJNmzYdbmHfycY/vvp/q3f+xWMxIpIkKRwOM0EgQbQXkJjgZV9/yYyRJJ+o0ncPQyfOU1OmaYiiRERJburctGfxA0LL0HUyk2Tt+ysnIo0bmr28ZXLiMWaRZRIRcZ4QOGeMGCNiFqNEx+5khygMt4hIIPKZxC2LiIgxmTHFJGICMUZEQSbbz2WCQEzwkyAQI3H/+yIzySMe2BPMBJGE/Q0YAmN+dug3zrIsJghMdqgv2TTIOsQpQmcqN4zUO8Kt1OnF1xi6ya0EM/ctZaVWazIrYaY6CDiZJhFFBYuIOt6yfSxuJQ7T6MMYk+z3JSX1MSHmJZGYsH9RQVSYIFv73jXhaysNCAoTDn0kOeBdS60vcJj37jiyLIuIhIO2fgCNmx3VvhNuHL7/hXNu6JxRnIyvzez0JnJGca7zA6oB73gHOaMYO0oNOYDfEgQiEsUDPnT2J0gi5uEdL5MxRqKYWsCyLEEQGLGAoKTe38NhTOj83EM6wufuG4hxnXc6dh2CaRxlAZtlcNMyGVfpoLfSNE/LP2XWjNuPvIKeHoSmaVKn2iyKomEYh1s496t4QPIraVJS5ERkcSGUlkuWmTo66NyK8q8/3bK4qp6gwh93jAmcOlULzu2DFAmCyAQiMrnFiYkC54wpTAkKMhEnJuw/qAmCh0mWwAQmeQSRM4GJoiRIoiCTICiyh5ggS14SBK+oMI+PiBTZKys+gQlE5BEUi1seyXPkchqWqZr79qplaWosaWlk6Nw0iCiSaLXzkhsGWUbUVIlz0rXU05OmltQPelM6vY9EZBFFLe3AZTr+xBmZpEaOtjuPHWNEtP+0g7GOx/YDxohzJor2EUdkgskt+3+BMYtzgYgzxjn3Ch6vnRz2s76eIpIoc+KCIHIm+AWZiEiUiDiJksRExgQi5vX4iTrOLTyK37CMgBIgJthndQIJmm54vB4SRb9y6Otybui6piY77UOuqfbLMSwjocX2723L4oammbpGJhFxw2Bfj9uIHjt0ltOB71rHhjiP8uPf02/v6v1bIeLEZZGZ+07O7J0vkH0CSIwxzrnC5CDzkiAQt59ERET+QxzrBVEUuWCSSaLMGOsnKEREgkCMcYszWSYihUSLMS4wn+BhkpR6JhETmcgkybJMxphf8pEdY9zeIdb+unQocSNhESct+bWP/75PkMFN1dq3Py2LOqW4puuKLFvEY3r8sO/R/tVZZB726Nqx+sN/7o7K3uGd5wSYj7F90S4IZFmUmuSc7PNjdqgTbiIikgXJsvcaYyRKIjGfsO+NYwIJjCzOBaG44OSjlq2nB2FBQQER1dfXFxUVEVFdXV3q6vBg5VraLRMXzZw5s/vK1/NEIpFQKOR0KZykqqokSZLU0+v2CYVqoOu6aZper9fpgjgJ1cCyLPOoZwA9f7BMMBgcN27cihUriIhzvmLFiqlTpzpdKAAA6Dt6wVnzXXfdtXDhwrq6us2bN7e2tl511VVOlwgAAPqOnn5FSESXXnrp22+/rarq+PHj169f7/f7D7dkS0tLLBbrzrL1NJqmrV692ulSOGzz5s0VFRVOl8Jh77zzjtNFcFhZWdmWLVucLoXDVq1adYRBFW7Q1NS0fv36oy7WC4KQiM4444x77733xz/+cWZm5hEW271799atW7utVD1QWVnZHXfc4XQpHPbcc8+9/fbbTpfCSaZpzps3z+lSOGzp0qXPP/+806Vw2G233VZVVeV0KZz04YcfPvzww0ddrHcEIcAx6dKQa+jTUAeAulwNEIQAAOBqCEIAAHC1A+9w7NVCoVBeXt6AAQOcLohj4vH4v//977POOsvpgjhp27ZtgUCgX79+ThfEMZzzNWvWuPxGo4qKClVVhw4d6nRBnPTBBx+MHTvWzTdTNjQ0ENG///3vIy/Wp4LwL3/5S0ZGRjAYdLogjuGcl5eXl5aWOl0QJzU0NPh8PjdXAyLas2ePm88Iiai9vV3X9aysLKcL4qSysrL+/fuzw385S5+XTCb9fv+UKVOOvFifCkIAAIBjhT5CAABwNQQhAAC4GoIQAABcDUEIAACuJv761792ugzHR2Vl5dKlS6uqqkpLS8Wj/bxkn1FTU7Nq1apNmzb5fL5wOJyaH4vFli1btmnTpsLCwiN8O2tfsnfv3o8//rigoECW5dSc119/3SVVQtO0VatWvf/++7FYrKioyP4JT8MwVqxY8dFHH2VlZaWnpztdxhMrkUisWLFi/fr1jLG8vLzU/Kqqqtdff72ysnLAgAF9shq0trZ+8skniUQiOzs7NZNzvnbt2jVr1vh8vs7zGxoali5dunPnztLS0tQnpbfjnO/YsWPTpk15eXmK0vGL3JqmrVu37v33329ubj5g6OyGDRtWrVpF+37mr2MVfcD777+fmZl5ww03TJgwYcqUKYZhOF2i7rBq1arMzMxLL7302muvTU9P//3vf2/Pb2lpGT58+LRp0+bNm5eXl7dz505ny9kNNE0bO3YsEe3Zs8ee88EHH4TD4RtuuOHMM8+cPHmyruuOFvDEqq6uHjly5IQJExYsWDBhwoStW7dyzg3DOO+8884444wbbrghHA6vXbvW6WKeQLW1taWlpTNnzvzxj39cVFR0zz332PPXr18fDocXLFhw9tlnn3POOZqmOVvO4+5HP/qRoijp6ek/+MEPOs9fsGDByJEjFy1alJOT8+KLL9ozt23blp2dfe21115wwQWnnHJKe3u7E0U+zpqbm9PT0+2wt2u+beDAgWeeeeb8+fNHjhx51llnxeNxe/7dd99dWlq6aNGioqKiRx55xJ7ZR4JwypQp9ktKJpNDhw5dunSp0yXqDnV1damqvGzZskAgYJ8BPPzww+eff75lWZzzH/7whzfeeKOTpewW99577y9+8YvOQTh16tSHH36Yc55MJocPH/6Pf/zDyfKdYJdddtlNN910wMxly5YNGjQokUhwzv/whz+ce+65ThStmzzxxBMTJkywH69bty4YDNr1f/r06f/7v//LOdc07eSTT37llVecLOUJUFlZqarqj3/8485B+MUXX6SlpTU2NnLO33zzzdLSUtM0Oeff/e53b7vtNs65ZVmTJk16/PHHnSr2caTrellZGef8gCBMXQCoqlpaWvrCCy9wzuvq6rxer/2nzz77LC0tLRKJcM77Qh9hNBp97733Zs+eTUSKolxyySXLli1zulDdITc3N/Xz0wUFBaZpWpZFRMuWLZs9e7bdFDBnzpw333zTyVKeeFu3bn3ttdfsILTFYrE1a9a4pErEYrE33njj9ttvX79+/ccff6zruj1/2bJlF198sf2tInPmzFm3bl17e7ujJT2BsrKyEomEXf9jsVg4HGaMaZq2cuVKuxrIsjxr1qy+Vw2Ki4s9Hs8BM5ctWzZ58mT7ywSmT5/e0NCwefNm2ndkICLG2BVXXNE39oYkSf379z94/qBBg+wHHo8nHA5rmkZEq1atOumkk+w/jR49Ojs7e+3atdQrfpj3qKqrq4mosLDQniwqKnLhjzHdd9993/nOd+xG/6qqqqKiInt+UVFRXV2dYRiS1Bfe64OZprlgwYInn3yy8+GgpqaGc955J9gHgj5pz549oih+//vfD4fDlZWVuq6vWbMmIyOjqqpq4sSJ9jL5+fmSJFVVVaWlpTlb2hNk7ty5GzduHDduXGlp6VdfffXKK68QUW1trWVZxcXF9jJFRUUbNmxwtJjdpKqqKvWqJUnKy8urqqoaNmxYU1NT573hkl9oWrZs2d69e2fNmkVf3zPUaSf0hStC0zQZY6m+UFEU3fZblHfddVd5eflDDz1kT5qmaY+VICJRFDnnpmk6V7oT63e/+924cePOPvvszjNdVSVUVdV1/corr/z73/++YcOGrKysRx99lL5eDYiIMdaHd8K2bduWLFly5ZVXXnnllUVFRU899RQR2dXeJdWgM7v+pyYlSTIMw517Y9OmTQsXLvzrX/9qXx8fcs9Q3wjC/Px8y7IaGxvtybq6uv1jgVzgvvvue+ONN1asWJE62S8oKKivr7cf19XVZWVlHdx40mc89NBDTU1NixYtuuWWW4jozjvv3LBhQ35+PufcJVXCbguZNGkSETHGJk2aZP8ye+dq0NLSout6qtWk73n00UcvvvjiX/ziF1ddddWrr776wgsvbN++PT8/n/Z97TIR1dXV9eE90Fnnt55zXl9fX1hYGAgE0tLSXLU3Nm/ePHPmzKeeeurCCy+053TeM9RpJ/SFIMzMzBw1atTKlSuJiHO+atWqo37Fap/xyCOPvPDCCytXrszJyUnNnDx58ooVK+zHK1eu7Nt745lnnrn88svPP/98+2Wec845+fn56enpo0ePdkmVKCgoGDFixM6dO+3JHTt22I0/kydPXrlypT2IYOXKlaeeemof/gZqURTtTiAisoeGiqLo8/nGjx/f+bMwefJkx4rYjSZPnrx27dpkMklEGzZskCTplFNOoYOODH17b2zfvn3GjBmPPPLI5Zdfnpp57rnnfvbZZ01NTURUUVGxe/fujt/q6cbRPSfQkiVLcnJyHnzwwWuvvXbo0KGpkbJ92/Lly4nokksuuXEfe5zY3r17s7Ozb7vttnvvvTctLe3jjz92uqTdIRaLUadRo6+++mp2dvaDDz543XXXDRkyJBaLOVq6E+vll18uLi5+9NFHb7/99uzs7F27dnHOE4nEsGHDrr766oceeig3N/fll192upgn0Pr16wOBwB133PHkk0+eeeaZF154oT1q9I033giHw7/97W/nz58/cOBAe4hgX/L222/feOONJ5988siRI2+88cY33njDnj9lypRp06Y98sgjgwcP/u1vf2vP/PDDD9PS0u67775bbrklLy/P7krvA372s5/deOONRHTllVfeeOON9oe9uLh4yJAhqWPj66+/bi983XXXTZgw4bHHHhs9evStt95qz+w7vz7xwQcfvPPOO+FweP78+Z1vLe/DduzYsWbNms5z5s2bZ48jraioeOGFFzRNmzNnzsknn+xQAbuVYRjPPfdcag8Q0Ycffrh8+fLMzMwFCxb0+Srx0Ucf2S/26quvTjV5tbS0/OUvf2lqapo2bVpq4ExftWvXrtdee629vX3EiBFz585N3S2+fv36ZcuWZWRkzJ8/v/OtWwsNAwAAAbdJREFU5X3Dpk2bPvnkk9TkmDFjxowZQ0SJROIvf/lLZWXlWWeddfHFF6cW+Pzzz1977TWfz3f99denRpP1dn/9618TiURq8rvf/a7H4/nzn//ceWzE6NGjx40bR0SGYSxevHjr1q2jRo266qqr7C7DvhOEAAAA30Bf6CMEAAD4xhCEAADgaghCAABwNQQhAAC4GoIQAABcDUEIAACuhiAEAABXQxAC9Cn19fVPP/20/ZMsANAVCEKAPmX37t2LFi366quvnC4IQK+BIAQAAFdDEAL0HStXrrR/cWbWrFnhcDgcDts/wA0AR4DvGgXoOxobG1966aVbb7314YcfPu2004jo9NNP7/NfOA7wLUlOFwAAjpvs7Gz7K/bHjBnTt39tDuA4QtMoAAC4GoIQAABcDUEIAACuhiAE6FOCwSARdf7BbgA4MgyWAehTSktLA4HAM888EwwG/X7/kCFD0tLSnC4UQI+GK0KAPiUYDD777LNbtmw5//zzx44d+/HHHztdIoCeDvcRAgCAq+GKEAAAXA1BCAAAroYgBAAAV0MQAgCAqyEIAQDA1RCEAADgaghCAABwtf8P6rswD5750MkAAAAASUVORK5CYII=", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# solve the ODEs\n", "# The model results are compared with the same model built by Anylogic, and the resules are the same!\n", "prob_measles = ODEProblem(vectorfield(seir),u0_measles,(0.0,120.0),p_measles);\n", "sol_measles = solve(prob_measles,Tsit5(),abstol=1e-8);\n", "plot(sol_measles)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 2. Chickenpox model \n", "## Time unit: Month" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 295354.0\n", " :E => 0.0\n", " :I => 1000.0\n", " :R => 567191.0" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# define parameter values and initial values of stocks\n", "# define constant parameters\n", "p_chickenpox = LVector(\n", " β=18.0, μ=0.03/12.0, δ=0.03/12.0, tlatent=14.0/30.0, trecovery=5.0/30.0\n", ")\n", "# define initial values for stocks\n", "u0_chickenpox = LVector(\n", " S=295354.0, E=0.0, I=1000.0, R=567191.0\n", ")" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dd3wUdf4/8PeUbekVUkACikgHIx2lKQGkCQKifE+acCootvNx4s8D0bvvCXJiOT08PL6AcoJ0kCpIkaMKyAkovYQQkpCwKbs77fP747NZV0CyKGGSndfze3ffmdnJ7HuXz8xr5jOf3RUYYwQAAGBVotkFAAAAmAlBCAAAloYgBAAAS0MQAgCApSEIAQDA0hCEAABgaQhCAACwNAQhAABYGoIQAAAsDUEIAACWFlZB+P777x8/ftzsKkym67rZJZjMMAx8cSCaAWPMMAyzqzAZmgERhXI0CKsgXLVq1Q8//GB2FSYrKyszuwSTKYqC/R/NQNM0RVHMrsJkaAaGYWiaVuFqYRWEAAAANwpBCAAAloYgBAAAS0MQAgCApSEIAQDA0hCEAABgaQhCAACwNAQhAABYmmx2AVDpmMF0n6F5dN1nGArTfbruMwyV6T5dVwymM63MIMY0r8EMpnsNItK8OjFiOtN9RmAL5ZsjzXPDH1eXXRIJ/mlRFkSbSEQkkOySiEgQSXJIRCTZRUEWBEmQ7KIgkuQUBVGQHKIoi6JdEG2iKAuSUxIkQXaKol0UZYFvAcAidJ/BDEZEhmIYmv87UwyVGdrPvkbHUJmhGmVlHj1CuMZWOEGQXVdeC0l2UZD8fyLaRNHmn5adP+3C4cfqQVhWVnb48GGzq7gBgiA0a9aMVEEr07VSXfPoapmueXStTNc8uu7RvcUKqQWaV9c9hubVda9haIbslGSXJNpFyS5ITklyiKJNkBySZBcESZQjRCLBmSQGdgzJIQkiCZIgOUQSBEEkyS4SEd9DJMfPdh5BFIiI75y/RPPoVP64oTFDNYh+ylSmM10xiEhXDKYxPst00soMPm1ohqEwQzUMjele3dCZ7jX4gUDz6DxZZZfkD0uHIDslySFJdlF2iaJdlByiZBflCElySJJD5P+RIyTJIfLiASqP7jN0n2Eo/v1RVwxDNTSPYaj+CaYxvpCfnjLDv7/w81HdZzCdBfaRQKMVbf4TSlEWSBBE2d+SJadoqIYgCqJN1HVdkn75y2UY0zxXfgsdPznm07wkPh04/fWfvMoCPybITpEEgZ+V+k9hJUFyiIJEkl0SZEEKnMI6REESZJckiILkFPlC2SkJsnDFIeXWs3oQzp8//49//ONtt91mdiFBGDHG/5//f/h/GSNidPzksTd6vdPqjra2CFmOlGSXJEdIcoQkuyRblOxKdjhIik6Ikp2S5BQllyQ7RNFufgd4pV638WTVynRDMQyVlbk9pAmCLuiKoXl0XWFKkar7DM1j8Kth3WfoXv8lMt9pJafIzxUkpyg5RckhyU5RckqyS5RdkuSUJKcoO0TJJfE3tvJeC1RxmoefcRqaV9fKdN1raF5d9+iaz9A9ulZ+9ql7/S1N8+j8PEx0iLzxiDbRf4pmE0WbKDlFMVKQHHxW4OegvAeFX4Tx/OAZc6PVFhcXR0dH3/R3gPhOpxiB2eC05rlefgrLtDLDUDVDY7pPZzppHp33PPHzWs2j8/Ulh7+/R7SJol2QnZJoE0RZlCMkyS4Isii7RMnuf8ckhyTaBMnBT39F0SbIzt90hLF6EGqaNmDAgI8++sjsQkLV/YHud42u1a57o19aobi4ODo66laWZDpRFkRZCmSt7BVkWZblkNq2oRiaz9C9hu4rP675+IHM0Eo1b4Ghlem6T9e8hu7RtfIDnOwKpGZ5ZPLTDpf/YOePTx6xLim8u5WqNUNjukfn8aaV6ZrH0Mp0HnK8i0XzlC/06JpHl52S5BIDDSAwYYuSXUl2f6+DU5SD/vXNfok3WSW9It3n7+/RFYOpTPPqvHdX8+i6Yhgq08oMX5FqKEz3GbpP1xVmKIbm0Q3F0FVD9xr8ZEJyiP77JryXyCbE3BmRck9CBS+qMl4SVCKBBAHH1JtGtIt2u0g3dNLMSPPqmkf3X1Z6Dd3rz0itzPAWqIFZ/+WCt/yywJ+aouQoz1GHyPNSdvrPbf3dti7JfxlRBa7mqxdDMfib7y1R1DK1RPdpHsOfZ179p+nyeGMG83erlMcbn7VHy64aDtklyRFBy104oakskkOUHGSL/PUpq/sMQ/VfhRsq433Rmld3xNsq/FsEIcANEsh/TLwR/h4zn/+a0p+jXl33GWqJ5s33d9XqPsPgF6B8WjP4tYVo/1l/mhwh8VtEsksUZVG0BY0hkgXRLvpU1aYovFetet0KNRSD3wNmevBlQeCOmv9SIDD+y399UH5SIkgCvyATnYLkEO2RNtnFe7klV7TMbyLIrvIr+AjpV/Q3QtXETyJtP+8OMwwjlN+iCfMg5MMd+ehHPrCK903rin+i8Ehx4M4wQOWRndKvuY3B6OcjLHQ+aMh/zqsaP92ACYwh0pihGKpXE1i+7tWZ4R9nyOMwMMKCBzlPTQoeyksku0Qq73UQBJJ+1d0XPtDD/yI0/0APKr+lFBifzGvjgzL4TaPyscGiIAq8a0u0iTIfW2H3D4myx8iB+0OBS23eKR0Y8aiqqq7rTqfzVxQPVhNWQdg0NjNqZ8quXUcMlfEBFHzwPT8p5ru65OTjl/xn2YbK0NcBVdevuvqka42SKI8cZqiGf0Ri0JCHn4byEmkeg8p/y5QZpPtu7NMyoiwyxpwJNr5nCZIgCIJoLx+F7//AjMDHHJXHsyjahOAwBriVwioIi9RLnrsK2/ZqKtpFURJCub+SmBdzZh+SEMKff4S6y+w6AKqesArCs2Un9VifPbbiW6MAAABcWAVhGGCMnTx5sqCgICUlpXbt2maXAwAQ/hCEVciFCxf69Olz6dKlOnXqnDhxom7dups2bTK7KACAMIcgrELeeOONunXr7tq1i39S8IcffjC7IgCA8IcxWlVITk5ORkZG4PPyDRo0MLceAAArwBXhlc6Wsh+KbsUTxdipdfLPBqwOGzZs6NChP/zwQ+fOnR944IEmTZrcijoAAKwNQXilNWfZgpNXfiN7ZUhwCJ93/dnnwx566KG9e/d+/vnny5cvf/nll4cNG/bJJ5/cgkoAAKwMQXilJ+4Sn7jLtB7jxo0bv/7660R08ODBzMzMJ554ol27dmYVAwBgBbhHWEU1adIkKirK7XabXQgAQJjDFWEVMnr06Jo1a7Zp08bhcMybNy8mJgaXgwAAlQ1XhFXI008/LYri3LlzP/7447p1627fvj0mJsbsogAAwhyuCKuQli1btmzZ0uwqAACsBVeEAABgaQhCAACwNHSNViEXL170eDyBWZfLVaNGDRPrAQCwAgRhFfLoo48eOHAgNjaWz3bs2HH27NmmVgQAEP4QhFXLxIkTJ0yYYHYVAAAWgnuEAABgabgivJKWn6OeO3YLnkiMjHHUb37FwkWLFh09epRPd+vWbcCAAbegEgAAK0MQXkk5ecjz/Y5b8ESiM/LqIExLS2vWrFlg+haUAQBgcQjCK0W06hbRqptZz96uXbuxY8ea9ewAABaEe4QAAGBpuCKsWjZu3Ojz+fh0UlLSqFGjzK0HACDsIQirkGHDhh07dqy4uJjPulwuc+sBALACBGEVMnz4cLNLAACwHNwjBAAAS0MQAgCApSEIq5AzZ85cunTJ7CoAAKyletwj3Llz59dff82nhwwZkpGRYWY1lWbkyJG9e/fGd40CANxK1eOKcPPmzcePH69Xr169evUwlhIAAG6i6nFFSEQtW7bs27evw+EwuxAAAAgr1SMIRVGcNWvWrFmzIiIi/v3vf1fql3CecZ87UXS68rYfEGFztU69+xY8EQAAXEf1CMLnnnvuxRdfJKK33377tdde++c//1l5z/XjpeNbz96KL91OcMUjCAEATBdqEJ47dy43Nzcw27JlS1G88v6iYRhbtmwpKCjo0KFDSkpKiFsuLi4+duxYSkpKampq8Ka2bt2an5/fvn371NRUSZL48o4dO65evTrELf8692d0uj+jU6U+BQAAVB2hBuE777zz2WefBfokt23b5nQ6g1dgjPXv3//06dONGzceO3bssmXLOnToUOFmH3300UWLFkmS9Morr7z66quBTQ0YMOD48ePNmjUbO3bs4sWL3W53/fr1GWNTpkzJysq6kRcIAABwPTfQNTpy5Mg33njjlx7dsGHDgQMHDh06FBkZOWPGjIkTJ/IPPJSWli5btuzRRx8NrDlv3ryBAwfywZ+TJk2aNWvW7373u+BNbdq0ae/evYcPH46Kinr//fcnTpw4YsSIuXPnCoLQr1+/kSNH3uiLBAAA+CU3EIR5eXmbNm3KyMioW7fu1Y8uW7asd+/ekZGRRDRkyJDnnnuuqKgoLi6utLT0zTffPHfu3B/+8Aciev311xctWtSzZ08ehHfeeec1N/Xggw9GRUXxTY0fP37p0qWh5N/333//4IMPBmZr16596NCh6/9J4KceqoJRo0bdcccd11+HMebxeEpKSn5phdLSUkEQbnZp1YnX65VlWZarx/3vSoJmoKqqruuappldiJnQDAzDsNlsNpvt+quFerAQBGHnzp0//vjjd99916VLl/nz51+x6ezs7NatW/PpmjVr2my27OzsuLi4GjVqbNy4sUuXLkTk9Xq/+OKLr776KjEx8TrPlZ2d3aJFCz6dnJzsdDqzs7Ov/ydc48aNP/roo169eoX4ooioSn0eY+jQoRWuIwiCy+XiZwnXxBi7zqNWIJczuxAzoRnwILziDo7VoBkYhqHreoWrhXqweOONN6ZOnUpEhYWFbdq0+cc//jFu3LjgFVRVDQxpEQRBFEVFUfhszZo1N27c2LRpU5vNduDAgeTk5Os/V/CmiEiSpMCmAAAAbq5Qv1kmcOUUHx/fr1+/PXv2XLFCampqfn4+n3a73V6vN/jTfjNnzkxPT4+Pj//Xv/5V4XMFb6q0tLS0tLRSPzgIAABW9mu6jw4ePJiZmcmny8rKXC6XIAj33nvvjBkzGGOCIHz11VcNGjSoUaMGX+ett95avHjxhg0bDMPo1q2bqqoTJ068zvbvvffeqVOnBjZVv3794E9WhLGjR4+63W4iio+Pv+222yzeuQcAcGuEeqjt06dPmzZt4uLivvrqq3379s2aNYuIFEWJjIzcs2dPZmbm4MGDp0yZMmrUqFatWr355ptvvvkmv0l78uTJJUuWbNiwISkpiYg2bNjQp0+fxx9/vFatWkS0bNmy//znPwcPHszNzS0pKRkwYEDr1q0ffvjh119/fcSIEW3atPnzn/88efJki9zvfeqpp86cOXPbbbedOnVK07TFixe3bNnS7KIAAMJcqF2jo0aN8ng8P/74Y9euXY8cOZKenk5Esiy/8847tWvXJiKn07l9+/Y77rjjyJEjs2bNevzxx/kf1q1bd/v27TwFiahmzZo7d+7kKUhEERER8fHxI0aMePDBB+Pj43kHrMPh+Oabbxo0aHDkyJGZM2da6vMSTz755Pr1648ePdqpU6fXXnvN7HIAAMJfqFeE/fv379+//xULRVF89tlnA7PJycmvvPLK1X97xfVc8OwDDzzwwAMPXP0nSUlJf/zjH0OsLSy1aNFi7ty5ZlcBABD+cBfqSqXnvcWny27BE9mi5MSmMVcsPHfu3N69e/Py8j788MPHHnvsFpQBAGBxCMIrefOVknPeW/BEslO8OggXLFjw9ddfnzp1qnHjxte8vAYAgJsLQXilxGYxic2uzKdb5vnnn58wYUJZWVmnTp1mzJjxwgsvmFUJAIBFVI9fqLeaiIiIGTNmTJ48uaCgwOxaAADCHIKwimrfvn27du2mT59udiEAAGEOQViFjB49umPHjoHZ6dOn84+pAABA5cE9wipkyJAhwbONGzdu3LixWcUAAFgErggBAMDSEIQAAGBpCEIAALA0BCEAAFgaghAAACwNQQgAAJaGj0/Q4cOHZ86caXYVoTp37pzZJQAAhBWrB2HHjh2//fbbvXv3ml1IqLp27dqkSROzqwAACB9WD8LGjRv/4x//MLsKAAAwDe4RAgCApSEIAQDA0hCEAABgaQhCAACwNAQhAABYGoIQAAAsDUEIAACWhiAEAABLQxACAIClIQgBAMDSEIQAAGBpCEIAALA0BCEAAFgaghAAACwNQQgAAJaGIAQAAEtDEAIAgKUhCAEAwNIQhAAAYGkIQgAAsDQEIQAAWBqCEAAALA1BCAAAloYgBAAAS5PNLiAkmzdvXr16NZ8ePXr0HXfcYW49AAAQNqpHEO7cudPr9fbu3ZuIEhISzC4HAADCR/UIQiKqVatWo0aN0tLSzC4EAADCSvW4RxgREbFhw4aRI0c2bdr0+PHjZpcDAADho3oE4bhx49asWbNmzZpnn3120qRJZpcDAADh44aDcM+ePTt37rzmQz6fb+HChe+///6RI0dC3+CJEyc2bNhw8uTJ4IWKovBNHT58OHh5o0aNcnNzb7RmAACAX3Jj9wj37t3bqVOnO++8c9++fVc8pOt6t27dbDZbixYt/vSnP3366ac9evSocIO9evXasWMHET3//POvvvoqX2gYxgMPPEBEd99996RJk+bMmVNQUFCnTh3DMP7f//t/Q4cOvaGaAQAAruMGglDTtKeeemrcuHHr1q27+tFVq1ZdvHjx+++/t9lsTZo0+dOf/sSDsKCgYM6cOc899xxfjTE2ffr0UaNGxcXFEdHs2bNr1KgxaNCg4E2tXr06Ozv70KFDdru9efPmr7322quvvrpw4UJRFF966aW+ffv+UoWGYZSUlBQWFvJZm80WFRUV+gsEAAALEhhjIa46adIkxliDBg2mTp169RXhk08+abfbZ8yYQUT5+fnJycn5+fmJiYlutzsrK6tdu3bTp09njD377LPffvvt6tWro6OjA387aNCg5s2bB64Ix40bR0Tvv/8+ERUWFiYkJOTm5taoUaPCCjMyMvLz8202G5+tVavW9u3bQ3x1YaOkpMTi8e/1emVZluVqMyK6MqAZqKqq67rT6TS7EDOhGRiGYbPZIiIirr9aqAeLgwcPLlq0aPfu3YsXL77mCufPn2/Xrh2fTkpKstvt58+fT0xMjImJWbNmTVZW1osvvuj1eg8cOLBmzZrr/9ucP38+MzOTT8fHx7tcrvPnz4cShA0bNhw/fnyvXr1CfFHhKvgkw4JsNhuCkCzfDBCEnMWbgWEYuq5XuFpIBwtN00aPHv3hhx9ep1UxxgRBCMwKwk/XmrGxsWvWrKlfv77dbj9y5EiFZyhXbIovCaVOAACAGxXSqNHt27cfPXr03XffHTx48Lvvvnvq1KnBgwf7fL7gdVJTUwPjOQsLC30+X2pqKp9ljL366qsNGzasXbv25MmTK3y64E253W6PxxPYFAAAwM0VUhDedddd//jHPwYNGjRo0KA2bdrExcUNGjRIkiQiOnXqlKIoRHT//fevXbuWX4R++eWXLVq0SE5OJiJ+X3D//v2rVq1at27dN9988/zzz1//6e6///5169YFNtW0adOUlJTf+DoBAACuKaSu0eCBnaqqbtmyhc8qilK3bt09e/ZkZmb279//L3/5S//+/TMzMz/44IOZM2fy9U+fPn3y5MnAfcE1a9Y8+uij586dq1WrFhHNmjVr7dq1O3bsOHLkyHfffTd27Nhu3br17dv3z3/+c79+/Vq1avXBBx/8/e9/r5SXDgAA8Cu+a7Rjx478Uo+IZFmeP39+vXr1iMhms23dunX+/Pl5eXlr1qwJjHbJyMhYsWJF4M9jY2NXrVoVmL377rtjYmICKVu3bl2+2S1btnz22Wd5eXlffvnlPffc82tfHQAAQAVu4OMTVV/Pnj0xarS4uNji48Tw8QlCM8CoUSJCMygfNRr4TN0vqR7fNQoAAFBJEIQAAGBpCEIAALA0BCEAAFgaghAAACwNQQgAAJaGIAQAAEtDEAIAgKUhCAEAwNIQhAAAYGkIQgAAsDQEIQAAWBqCEAAALA1BCAAAloYgBAAAS0MQAgCApSEIAQDA0hCEAABgaQhCAACwNAQhAABYGoIQAAAsDUEIAACWhiAEAABLQxACAIClIQgBAMDSEIQAAGBpCEIAALA0BCEAAFgaghAAACwNQQgAAJaGIAQAAEtDEAIAgKUhCAEAwNIQhAAAYGkIQgAAsDQEIQAAWBqCEAAALA1BCAAAloYgBAAAS0MQAgCApclmFxCStWvXLl68mE9PmDChYcOG5tYDAABho3oE4YEDBxITEwcOHEhE6enpZpcDAADho3oEIRHZbDZJkho0aOByucyuBQAAwkf1uEeYlJR04sSJ6dOnN2vW7LvvvjO7HAAACB/V44pw5MiRI0eOJKL58+e/+eabn3/+udkVAQBAmAg1CA8dOrRp06acnJyaNWsOGDDgmjfq3G733Llzc3Nze/To0b59+1A26/F4Dhw4cOjQoZYtW7Zs2fLqTWVlZXXo0CGwPD093e12h1gzAABAhULtGv30008PHz4cHR29Z8+ehg0bfv/991esoChKhw4dNm7c6HK5HnrooQULFoSy2QcffHD48OETJ05ctWpVYKGqqvfee+9XX33lcrkGDhw4f/789957b9GiRQsXLnzmmWeGDRsWYs0AAAAVCvWK8M033wxMFxUVLVy4sHHjxsEr8I83LFy4UBTFjIyMN954Y/DgwUSUk5Pz7rvvvvnmm6IoEpGu66+88soLL7xQo0YNIlq7dq3NZhs0aFDwppYuXaqq6hdffCGK4u233z558uQPPvhg69atoih++OGH7dq1+6UiFUU5duzY3r17+Wx8fHy9evVCfIEAAGBNN3yPsLi4+MSJE4888sgVyzdt2tS9e3eedj179nz00Ufz8vKSk5MTExP/+9//PvbYY/PmzSOi4cOHFxQUxMbG8r+y2WxXP8UVmxoyZMhdd93VuXPnCms7ffr0jBkz/vWvf/HZuLi45cuX3+gLrO5KSkrMLsFkXq9XlmVZrh73vysJmoGqqrquq6pqdiFmQjMwDMNms10zaILdwMFi9uzZU6ZMycnJGTt27NChQ694NCcnJ3AzLy4uzul0nj9/Pjk52W63L1y4sF+/fqNGjfL5fCUlJUuXLrXb7dd5opycnFatWvHp6OjoiIiI8+fPp6SkVFhh/fr1x48f36tXr9BfVFiKjo42uwQz2Ww2BCFZvhnwIHQ6nWYXYjKLNwPDMHRdr3C1GzhYPPTQQ/fee++hQ4eefvrptm3bDhkyJPhRURQNwwh+ekmS+LTT6Vy0aFHdunXtdvuxY8eun4LX39T1XajVYVNZSswF5pIp2kaySDE2wSZSbAVPCAAA1nUDQRgbGxsbG3v77bcfOnRozpw5VwRhWlpaTk4On87Ly1MUJS0tjc/quv7EE0+0a9dOVdUnn3zyk08+4d2evyR4U5cuXfJ6vYFNXZ9NKd3hq/mf3bpXJ7dKukGXFaYxuqyQQ6IImWJsgkumSJni7BRlE6JsFG2jWDvF2YUYG8XYKdZO8XYhzkFxdop3CJFWv6iAasark0cjj84KSgVSmWKQRyOvTopBpRojoiIfMSLFoFKViH5aTkSXFTL8k2Qwuqz8mgLsEgXvNfx8lIicErkkIbDEJlKUjUSiWLsgChRrJ0mgGDs5RIqQhcBfAdwaoR7pFUUJXMn997//5cnEGNu1a1fTpk0jIiJ69eo1YcKEqVOnOhyOJUuWtGvXLiEhgYh0XR8+fPjly5eXLFkiCMLAgQP5/cLrXOT16tXr6aef9nq9TqdzyZIlrVu3Tk5ODqXIxIv7x8cfuGbXKD9AXFaYR6cyjQp9VKKxEpWKVXIrVKiws6V0WaEiHytSiP+n0Mc0g+IdlOAQEhyU6BQSHJTooASHkOSkJCclOYUkJyU6hCQnSUKIbyTAL7qsUGl5s+TTZRpdVqhEpTKNSjTmVsijU6lKbpV5df9yn05FClMNKlb9J3wuSbAL9ii77pDIKZFLJptIUbJARHEOEsifQ0RkFyneLggCyQLVifI3Y4dEjFHkz++q2ERijHRGApFxdenlFJ1KtZ9m+fkoEXl1KlQYEZ0pJc3wJ7FBdFkxdIPcKmkGFavk1cmjM7dCjCjWThGy4JIo1k6RMkXIFG0T4hwUKVOUjaJtQpydYmwUYxei+VmsjWLtAn+BADck1CBs3LhxZmZmfHz8gQMHzp8///XXXxORqqpt27bds2dPZmZmz549p0+f3qVLl+bNmy9YsCDwmffz5887nc7AfcGFCxc+9dRTFy9eTE1NJaL33ntv+fLl33333b59+zZv3vzcc8/16tUrKyurXr16Xbp0adGixYIFC+bPn//bX6dTIqdE8Y7gfaTi/cWnU6FCBV52yUeXfPx/qcDLThZTgY8KvEaely75WJ6XEhyU7BQSnZTkEGq4qIaLkp1CkpNqOAU+nezESa6FFPrIrbLLCrkVcqvkVphb9S/k517Fqn+FYpWKVSpR2WWFYmwUaaMomxBjo2gbRdkoQhbi7BQpU6SN4u3CbZHkkinKRjE20SFRtI0iZHJIFGsX7CJFB0VXcXFxtb45pDNyK1SqMY9GbpVKVCrVqFRjhT4q1ahEpXwvO+6mywoVq4ZbJbdClxVyq6zQRzE2inMIMTaKsYkJTi3OLsTZKc5BcXYh3kFx9p8m4h1CHO6bAJHAGKt4LaJz58795z//cbvdtWrV6tKlC081xtimTZtatWrFdzlVVdesWZOfn9+5c+e6deuGstkjR46cO3cuMNuwYUP+UX1N09asWZOXl9epU6fQPwLRs2dPswbL5Hspz8sKvJTvYxfK/LN5Xsr1sDwv5XtZnpfi7OXp6BJquijZSUlOIdlJKS4hcH35268sq/sR8Le7uaNG+RG5UGGX+aFWYZdVPkGXFX//gVthlxVy+5czt0pxdoq1CzH+SxaKsQmxdoqzU7TNvzDa5l8h2p95N/mIbOVmUKTQZYUVlGmFXqOMbEUKFfl4Nw8r5BM+VqhQoY8KfaxYpfjyjIx3/Cwj4+0UZ6c4PuEgHqiOkMyFl4AAAB/6SURBVIYrVBXh2gzKNCrVqFhlboXKNCrTqFBhfMKtkFv1Txf6qFQ1uqXS880rGDUaahBWCyYGYSh4OuZ76aKH5XqIB+RFHpYeKvCxfC/FO/yJmOQQkpz+sEx0UqK/e9bfT3udvAzXph+66wThZYWKg67JChVWfnFGlxXGryqKFBa4wuB96fzOMb+FzCMt1k4xNoq1C/EO/8IYO8XZf8o806EZhDhq1GBUWB6NRT665GNFChX6glIzaKLIR5JIsTaKcwix5Z2xwW2AjzngbSDKRlE2ipKFeMetecXXUAWbAT+z9OjMq9NlhXw6lajkVplP510j5NXJrbJS1X/pX6wyj04lKl1WqEwjj84KfRQh+8d88E6RSBvF2YUImSJlirFTtM0/HWcnl0R3xej142/exyfgN+KXfUT0S72yjKjAS/lelu+lfB+PTMouY99dogKfUeClAh9d8rECL8XYKdHhP4cNjO7xn7EaUo0Yg++WfJ+MspGzWp3GXkeJSj6DLivMo5FHpyIflWnMq1NR0I20gjLRY1CZrherrMhHxSqVaD/1PUbzW0o2iikfJBVtpxgbpUYId8VSnINi7WJM+T0nfv8JwpgoUKKDEn+6aVJBnwy/a1uksKLyHoJAr8BFd6DT2yjy8VZHJSorUvyd21GyEGsnV/lBnHdu85u40TZBFijWTmL5//IlLtm/8/KFRBQlCzaRiEgQ6Cb3Iqikld/+LdGYahAR+XQq0342UawyjflvBjOiIh8R0WWFGURuhXTmH3VVpDCDUZHiv/vLky9w95cPnoq1k13inSKCQ6IYG0XZyCFSvF2oFcnvClO0TXRKxPfZSBu5pBu7DWwYRgifnkAQViVCCGHJFfroUnn3Dj9j5RPnSim/VPScN8pvTdFlhZVopBn+EyWX5D9RtYkUZxdsIkXb/AP2iCjeLhD5bzsFTxD5d8tfoVRjSvnexYcsEZFqUIlKRP7dgzEqUojKd0X+J3z6skI6oyKF8R0vUiaHRHF2gR8+4h3kkgSXTHHl4yni7UK6g0U7WKxDjLaJ/MS8MvoewZr4tUhqRKjByQXO0op85NWpTPtpuBMf1svT5WQxMaIihRgjt2rozP8oBY3p5WsS/bTXBPCdIhhjDkG49lcKBLZMRDaRBCKH5B/HIBJFlMdtYBgwH4dF5YcCPt5KIIpzEBHF2ASpfLxVjI0kkWLtoiRQnJ1ksTzvJSHGXkXHFSIIq6V4R/DAn5+1rOLisujoK7uD+EmZW2Uef986aQYVKkwrz5tilRj5x/Vll5FPJyL/iET/ZlVDu8FOdIdIRGQT/XsUHzTPT28Doxadkn/U4u0xROVj6yNl0S760zrGRnJ5ZkeFdnHm9TJZFmSMTYIqg5991iCBfuqnvPmBEDjRDCgpKYmKigpeYhNJFEhnP+2MQAhCi5DFK7KTq5LnZgDwq/Cx8cFkhUWbd4eyGsFZMwAAWBqCEAAALA1BCAAAloYgBAAAS0MQAgCApSEIAQDA0hCEAABgaQhCAACwNAQhAABYGoIQAAAsDUEIAACWhiAEAABLQxACAIClIQgBAMDSEIQAAGBpCEIAALA0BCEAAFgaghAAACwNQQgAAJaGIAQAAEtDEAIAgKUhCAEAwNIQhAAAYGkIQgAAsDQEIQAAWBqCEAAALA1BCAAAloYgBAAAS0MQAgCApSEIAQDA0hCEAABgaQhCAACwNAQhAABYGoIQAAAsDUEIAACWhiAEAABLqx5BuGDBggfKffvtt2aXAwAA4UM2u4CQnDp1qmfPniNGjCCiqKgos8sBAIDwUT2uCIno3LlzO3fu1HXdZrOZXQsAAISP6hGEderUkSRp48aNbdq02b59u9nlAABA+KgeXaNDhgwZMmQIEXXo0GHq1KlLliwxuyIAAAgToQbhgQMHVq1adebMmdq1a48cOTI1NfXqdXJycv7+978XFBR07969f//+oWw2Jydn9+7dP/zwQ4cOHdq3bx9YfuHChQ8++CA/P7979+4PPfRQYHlERISqqiHWDAAAUKFQg/B3v/vd/fffn5mZuW3btubNm+/fvz8tLS14BY/H0759+6ysrNatWz/77LP5+fmjR4+ucLNjxoxxu92nT5/2+XyBIPR6vR06dOjWrVvbtm2fe+65ixcvXrhwISkpiTH23nvvTZ069UZfJAAAwC8JNQj37t0ryzIRPfHEE5mZmStWrBg7dmzwCp9//nlSUtJHH31ERImJiS+99NKoUaMEQTh79uykSZM+/PBDu91ORF6v9/e///2f//xnnqMrVqwgokGDBgVvauHChbGxsTNnziSiGjVqTJgwYeHChbt37xYEYcWKFXfeeecvFVlWVrZ//37+REQUHR3dpk2bkN8KAACwolCDkKcgETHG3G53fHz8FSts27ata9eufLpr165Hjx7Nzc1NSUlJS0tTFKVfv35LliwRBOHhhx9OTk5OSUm5znNdsanjx4/XqFFjzJgxFRZ54cKFzz//fMOGDXw2Li5uzpw5Ib7AsFFaWioIgtlVmMnr9cqyHGix1oRmoKqqruuappldiJnQDAzDsNlsFX7W4IYPFtOmTXM4HFffArxw4UKDBg34dHR0tMvlysnJSUlJkSRp9uzZw4YNGzRokKIoqamps2bNEsXrjVa9cOFC3bp1+XRERERUVFROTk56enqFtdWrV2/8+PG9evW60RcVThhjFv+opVzO7ELMhGbAg9DpdJpdiJnQDAzD0HW9wtVu7GDx2WefvfPOO19//XWg+zHA4XAoisKnGWOqqgaaoCRJH3/8cd26de12+9KlS6+fgldsiogURbF4awYAgMpzA58jXLRo0Ysvvrh27dr69etf/Wh6enp2djafzsnJ0XU9MJrG6/UOHDiwd+/enTp1Gjx4cHDIXVPwpnJzc1VVvWJgDgAAwM0SahCuXr366aefXrFiRZMmTQILGWMrVqy4fPkyEfXv33/FihXFxcVENH/+/C5dusTGxhKRz+d7+OGH09LSZs2aNXfu3Li4uH79+nm93us8V//+/VeuXOl2u/mm7rvvvoSEhF/9CgEAAK4j1K7RRx55JCIiIjBSdMSIEU8//bSqqn379t2zZ09mZmaXLl3atm3bpk2bhg0bbt26dfny5XzNS5cuNW/efMqUKbxHdPbs2RMnTiwuLua9nX/9618XLlx44sSJ7du3L1269NVXX+3fv/99993XsWPHNm3aNGrUaOvWrUuXLq2EFw4AAEBEJDDGQllv3759hmEEZlNSUtLT0xljhw4duv3223mqMcZ2796dm5vbvn37xMTEUDZ75syZvLy8wGydOnWSkpL4pvbs2ZOTk9OhQ4cQN0VEPXv2xGCZ4uLi6Ohos6swE0aNEpoBBssQEZpB+WCZmzZqtGXLllcvFAShcePGwbOtW7cOvUQiuu2222677bZrbrlVq1Y3tCkAAIBfoXp86TYAAEAlQRACAIClIQgBAMDSEIQAAGBpCEIAALA0BCEAAFgaghAAACwNQQgAAJaGIAQAAEtDEAIAgKUhCAEAwNIQhAAAYGkIQgAAsDQEIQAAWBqCEAAALA1BCAAAloYgBAAAS0MQAgCApSEIAQDA0hCEAABgaQhCAACwNAQhAABYGoIQAAAsDUEIAACWhiAEAABLQxACAIClIQgBAMDSEIQAAGBpCEIAALA0BCEAAFgaghAAACwNQQgAAJaGIAQAAEtDEAIAgKUhCAEAwNIQhAAAYGkIQgAAsDQEIQAAWBqCEAAALA1BCAAAllY9gnDOnDn3lNu1a5fZ5QAAQPiQzS4gJBcuXBg9evTvf/97swsBAIBwUz2uCIloz549c+bMOXnypNmFAABAWKkeQdiwYcMGDRrk5uZ279597dq1ZpcDAADho3p0jfbp06dPnz5E1KxZs/fffz8rK8vsigAAIEyEGoTnz5/fvHnz999/f8cddwwfPvya6/z4449vv/12Xl5ejx49nnjiCUEQKtzsoUOHvvnmm+PHj3fv3r1r166B5ceOHZs2bdrFixezsrLGjBkT2JSqqpIkhVgzAABAhUINwk8//XTTpk0lJSX79++/ZhC63e777rtvzJgxffr0+cMf/lBSUvL8889XuNn//d//9Xg8Bw8ejIqKCgRhSUnJvffeO2rUqN69e7/88stut/v8+fMxMTGGYcyfP3/WrFkhvzoAAIAKhBqEL7300ksvvfT2229v2rTpmit8+umn9evXf/3114nIbrePHj16woQJoiieOnVq/Pjx8+fPj4qKIqLLly8PHTp05syZtWrVIqI5c+YQ0aBBg4I39dlnn9WrV++NN94gIpfL9fjjj2/btm3//v2SJD3zzDPJycm/VKTb7d6yZUtxcTGfjY6OtmAnqq7ruq6bXYWZdF0XBCGUDokwhmaglzO7EDPhHTAMgzFW4Wo37R7h7t27O3bsyKc7dux49uzZCxcupKWl1alTp27dullZWWvWrDEMIysrq3379jwFr7OpDh06BDaVnZ0tSVL//v0rrKGoqGjr1q1Hjx7ls0lJSZ07d/5Nr6oaUhTF5/OZXYWZfD4f9n80A1VV+SmR2YWYCc3AMIxQ7qbdtCDMzc1t1KgRn46IiIiIiOBBKAjCjBkzxo8f36tXL1VVO3bsOG3atAo3dccdd/Bph8MRHR194cKF2rVrV1jDbbfdxp/oN76Wak3X9YiICLOrMJMoirIsy3L1GAhWSdAMeBA6nU6zCzETmoFhGKGcE9+0g0VERITX6w08t8/ni4yM5LOCILzxxhv169e32+1/+tOfbmhTjDGv1xvYFAAAwM1104Kwdu3ap0+f5tNnz54lorS0ND57+fLlHj16jBgxwuv19urVa/Xq1fx+YSibys7O1nU9PT39ZtUJAAAQ7Dd9oN4wjH/+858FBQVENGjQoBUrVuTn5xPRv/71rx49ekRHRxPR5cuXs7KyOnTo8NZbb82YMePuu+/OysoKjGe5pkGDBq1cuTIvL49vqnv37rGxsb+lTgAAgF/EQvP555/Hx8e7XC6bzRYfHz9u3DjGGL8Nu2fPHr7Ok08+mZaW1qFDh1q1ah08eJAvLCgoeOeddwLbMQxj+vTpRUVFfPbFF1+Mj4+32+0ulys+Pv7//u//+PJx48alpqZ27NgxPT39wIEDIRbZo0ePVatWhbhyuHK73WaXYDKPx6OqqtlVmAzNQFEUj8djdhUmQzPQdV1RlApXE1gIQ0uJSFGU0tLSwKzD4eD3YPPz8+Pj4wPDck6dOnXx4sXmzZs7HI5QNltWVhY8qCkiIiLwh3xTzZo1C/12d8+ePTFYpri4mF+LW5bX68VgGTQDDJYhNIPywTI2m+36q4V6sLDb7Xa7/erlSUlJwbMZGRkZGRkhbpPKx5de86Eb3RQAAMCvUD2+dBsAAKCSIAgBAMDSEIQAAGBpCEIAALA0BCEAAFgaghAAACwNQQgAAJaGIAQAAEtDEAIAgKUhCAEAwNIQhAAAYGkIQgAAsDQEIQAAWBqCEAAALA1BCAAAloYgBAAAS0MQAgCApSEIAQDA0hCEAABgaQhCAACwNAQhAABYGoIQAAAsDUEIAACWhiAEAABLQxACAIClIQgBAMDSEIQAAGBpCEIAALA0BCEAAFgaghAAACwNQQgAAJaGIAQAAEtDEAIAgKUhCAEAwNIQhAAAYGkIQgAAsDQEIQAAWBqCEAAALA1BCAAAloYgBAAAS0MQAgCApSEIAQDA0hCEAABgabLZBYTkzJkzX331VUJCQs+ePe12u9nlAABA+KgGV4Tbtm1r0aLFf/7zn2nTpnXt2lVVVbMrAgCA8FENgnDSpEkTJ06cOXPmxo0bL126tGzZMrMrAgCA8FHVg9Dj8WzcuHHAgAFEZLPZ+vbtu2rVKrOLAgCA8FHV7xHm5OQwxtLS0vhsWlravn37fmllRyt94emPvpnxD4GIiAQSbHbbNdbTNCL222tjjJGuhb6+Jgq+q048HAbTBcGhh1qPIF/rFQVXRSSEuC1RJFEKcd1r1iLIFbcfiQlOVkFFgij9tkp+hjFGAgkhvw3XUNGb/Otpt6hj3zAMUQztNNfQmaGHsiIjKhONile7wf2iYr/qTeN71G9oBLeKKFGI/1I3sk1BFImIGUwQK34PHEyQK9pJfyIIglRpO8jNxog1SW/Zvccz11+tqgehrutEFNilJUnStF/cx9L2FAsNIi7VjOf7gMPhuLfVfbxBEOP/x4hItDuus4NIJOrk39tFQTT4XwmCQGQwg0gQiBgxkQQmChTUIISK0lUWZZfsvGKhR/PKouTVfNf90yCqcv3HvV6v03nls1wTM3T2W47LjJFWQTFEpDHdq1ewGtM1YhUfYUOka5ogiqHGwDWqYZUYV7KNhFtxcFZV1WazEZFMoiYQMYMEgbf/4BZORKIos6D3SiCBBTXk4FlBECMlx8+eRhCpfAchop8m5Bse0SaIIjMMQRSJEQuqlohEe0jtubxgEgTRYIau64Zh8DchQBRE4+a1tJuC6SrTQzoRuQG6RoZORD6f4nBU/G/h1VWNhXzuwgx2c090KhVj0SkNK1yrqgdhamoqEV28eDE9PZ2IcnNzA1eHVzvJ4sY/MKZXr163rr6qp7i4ODo62uwqzOT1emVZlkO4Wg1jaAaqquq6HuJJYbhCMzAMQw/hPKOq3yOMiopq1arV2rVriYgxtnbt2q5du5pdFAAAhI9qcNY8ceLEUaNG5ebmHjx4sKioaMiQIWZXBAAA4aOqXxESUb9+/b788kuv19umTZsdO3ZERET80pqFhYWlpaW3sraqRlGUr776yuwqTHbw4MEzZ86YXYXJ1qxZY3YJJjt16tT3339vdhUmW79+/XUGVVhBQUHBjh07KlytGgQhEbVu3Xry5MnPPvtsfHz8dVY7ceLE4cOHb1lVVdCpU6deeukls6sw2SeffPLll1+aXYWZdF0fOnSo2VWYbNmyZXPmzDG7CpNNmDAhOzvb7CrM9M0330ybNq3C1apHEALcEMZuwsdjoFpDGwAKuRkgCAEAwNIQhAAAYGlCOHUgREdH16xZs27dumYXYpqysrL9+/e3b9/e7ELMdOTIkcjIyNq1a5tdiGkYY5s2bbL4B43OnDnj9XrvvPNOswsx07Zt2+655x4rf5gyLy+PiPbv33/91cIqCGfPnh0XFxcVFWV2IaZhjJ0+fTojI8PsQsyUl5fncrms3AyI6OTJk1Y+IyQit9utqmpiYqLZhZjp1KlTderUEW7JlxlVTT6fLyIiokuXLtdfLayCEAAA4EbhHiEAAFgaghAAACwNQQgAAJaGIAQAAEuTJk2aZHYNN8fZs2eXLVuWnZ2dkZEhSTftV16ruJycnPXr1+/bt8/lciUkJASWl5aWrly5ct++fWlpadf5dtZwcu7cuV27dqWmpgZ+gu7cuXNLly61SJNQFGX9+vVbt24tLS1NT0/nP8eoadratWu3b9+emJgYGxtrdo2Vy+PxrF27dseOHYIg1KxZM7A8Ozt76dKlZ8+erVu3blg2g6Kiot27d3s8nqSkpMBCxtjmzZs3bdrkcrmCl+fl5S1btuzYsWMZGRlX/Fhj9cUYO3r06L59+2rWrGm3+39/UVGULVu2bN269dKlS1cMnd25c+f69eup/Gf+/JsIA1u3bo2Pjx85cmTbtm27dOmiaZrZFd0K69evj4+P79ev32OPPRYbGztjxgy+vLCw8K677srKyho6dGjNmjWPHTtmbp23gKIo99xzDxGdPHmSL9m2bVtCQsLIkSPbtWvXuXNnVVVNLbBynT9/vlGjRm3bth0xYkTbtm0PHz7MGNM0rVu3bq1btx45cmRCQsLmzZvNLrMSXbhwISMjo1evXs8++2x6evprr73Gl+/YsSMhIWHEiBEdOnTo2LGjoijm1nnTjRs3zm63x8bGPvnkk8HLR4wY0ahRo7FjxyYnJ3/22Wd84ZEjR5KSkh577LEHHnigadOmbrfbjJJvskuXLsXGxvKw5y2fq1evXrt27YYPH96oUaP27duXlZXx5a+++mpGRsbYsWPT09OnT5/OF4ZJEHbp0oW/JJ/Pd+eddy5btszsim6F3NzcQFNeuXJlZGQkPwOYNm3a/fffbxgGY+ypp54aM2aMmVXeEpMnT3755ZeDg7Br167Tpk1jjPl8vrvuumvJkiVm1lfJ+vfv//vf//6KhStXrrz99ts9Hg9j7N13373vvvvMKO0Wef/999u2bcunt2zZEhUVxdt/jx49/vKXvzDGFEVp0qTJggULzKyyEpw9e9br9T777LPBQfjf//43JiYmPz+fMbZixYqMjAxd1xljjz/++IQJExhjhmF06tTpvffeM6vsm0hV1VOnTjHGrgjCwAWA1+vNyMiYN28eYyw3N9fpdPKHvv3225iYmOLiYsZYONwjLCkp+frrrwcOHEhEdru9T58+K1euNLuoW6FGjRqBn59OTU3Vdd0wDCJauXLlwIEDeVfAww8/vGLFCjOrrHyHDx9evHgxD0KutLR006ZNFmkSpaWly5cvf+GFF3bs2LFr1y5VVfnylStX9u7dm3+ryMMPP7xlyxa3221qpZUoMTHR4/Hw9l9aWpqQkCAIgqIo69at483AZrP17ds3/JpBrVq1HA7HFQtXrlzZuXNn/mUCPXr0yMvLO3jwIJUfGYhIEIQBAwaEx7shy3KdOnWuXn777bfzCYfDkZCQoCgKEa1fv75x48b8oZYtWyYlJW3evJmqxQ/zVuj8+fNElJaWxmfT09Mt+GNMU6ZM+d3vfsc7/bOzs9PT0/ny9PT03NxcTdNkORz+ra+m6/qIESM++OCD4MNBTk4OYyz4TeAHgrB08uRJSZKeeOKJhISEs2fPqqq6adOmuLi47Ozse++9l6+TkpIiy3J2dnZMTIy51VaSQYMG7d27t1WrVhkZGT/++OOCBQuI6MKFC4Zh1KpVi6+Tnp6+c+dOU8u8RbKzswOvWpblmjVrZmdnN2jQoKCgIPjdsMgvNK1cufLcuXN9+/aln78zFPQmhMMVoa7rgiAE7oVKkmS136KcOHHi6dOnp06dymd1XedjJYhIkiTGmK7r5lVXud5+++1WrVp16NAheKGlmoTX61VVdfDgwYsWLdq5c2diYuLf/vY3+nkzICJBEML4TThy5MjChQsHDx48ePDg9PT0jz76iIh4s7dIMwjG239gVpZlTdOs+W7s27dv1KhRc+fO5dfH13xnKDyCMCUlxTCM/Px8Ppubm/vTWCALmDJlyvLly9euXRs42U9NTb148SKfzs3NTUxMvLrzJGxMnTq1oKBg7Nix48ePJ6JXXnll586dKSkpjDGLNAneF9KpUyciEgShU6dO/JfZg5tBYWGhqqqBXpPw87e//a13794vv/zykCFDvvjii3nz5v3www8pKSlU/rXLRJSbmxvG70Cw4H96xtjFixfT0tIiIyNjYmIs9W4cPHiwV69eH330Uffu3fmS4HeGgt6EcAjC+Pj4Fi1arFu3jogYY+vXr6/wK1bDxvTp0+fNm7du3brk5OTAws6dO69du5ZPr1u3LrzfjY8//vihhx66//77+cvs2LFjSkpKbGxsy5YtLdIkUlNTGzZseOzYMT579OhR3vnTuXPndevW8UEE69ata9asWRh/A7UkSfwmEBHxoaGSJLlcrjZt2gTvC507dzatxFuoc+fOmzdv9vl8RLRz505Zlps2bUpXHRnC+9344YcfevbsOX369Iceeiiw8L777vv2228LCgqI6MyZMydOnPD/Vs8tHN1TiRYuXJicnPzWW2899thjd955Z2CkbHhbvXo1EfXp02dMOT5O7Ny5c0lJSRMmTJg8eXJMTMyuXbvMrvRWKC0tpaBRo1988UVSUtJbb701bNiw+vXrl5aWmlpd5fr3v/9dq1atv/3tby+88EJSUtLx48cZYx6Pp0GDBo888sjUqVNr1Kjx73//2+wyK9GOHTsiIyNfeumlDz74oF27dt27d+ejRpcvX56QkPDXv/51+PDh9erV40MEw8mXX345ZsyYJk2aNGrUaMyYMcuXL+fLu3TpkpWVNX369DvuuOOvf/0rX/jNN9/ExMRMmTJl/PjxNWvW5LfSw8CLL744ZswYIho8ePCYMWP4zl6rVq369esHjo1Lly7lKw8bNqxt27bvvPNOy5Ytn3nmGb4wfH59Ytu2bWvWrElISBg+fHjwR8vD2NGjRzdt2hS8ZOjQoXwc6ZkzZ+bNm6coysMPP9ykSROTCrylNE375JNPAu8AEX3zzTerV6+Oj48fMWJE2DeJ7du38xf7yCOPBLq8CgsLZ8+eXVBQkJWVFRg4E66OHz++ePFit9vdsGHDQYMGBT4tvmPHjpUrV8bFxQ0fPjz4o+XhYd++fbt37w7MZmZmZmZmEpHH45k9e/bZs2fbt2/fu3fvwAoHDhxYvHixy+X6n//5n8Bosupu7ty5Ho8nMPv44487HI5Zs2YFj41o2bJlq1atiEjTtE8//fTw4cMtWrQYMmQIv2UYPkEIAADwK4TDPUIAAIBfDUEIAACWhiAEAABLQxACAIClIQgBAMDSEIQAAGBpCEIAALA0BCFAWLl48eLMmTP5T7IAQCgQhABh5cSJE2PHjv3xxx/NLgSg2kAQAgCApSEIAcLHunXr+C/O9O3bNyEhISEhgf8ANwBcB75rFCB85Ofnz58//5lnnpk2bVrz5s2J6O677w77LxwH+I1kswsAgJsmKSmJf8V+ZmZmeP/aHMBNhK5RAACwNAQhAABYGoIQAAAsDUEIEFaioqKIKPgHuwHg+jBYBiCsZGRkREZGfvzxx1FRUREREfXr14+JiTG7KIAqDVeEAGElKirqn//85/fff3///fffc889u3btMrsigKoOnyMEAABLwxUhAABYGoIQAAAsDUEIAACWhiAEAABLQxACAIClIQgBAMDSEIQAAGBp/x9Qr4v7SY6PVQAAAABJRU5ErkJggg==", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# solve the ODEs\n", "# The model results are compared with the same model built by Anylogic, and the resules are the same!\n", "prob_chickenpox = ODEProblem(vectorfield(seir),u0_chickenpox,(0.0,120.0),p_chickenpox);\n", "sol_chickenpox = solve(prob_chickenpox,Tsit5(),abstol=1e-8);\n", "plot(sol_chickenpox)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "retcode: Success\n", "Interpolation: specialized 4th order \"free\" interpolation\n", "t: 298-element Vector{Float64}:\n", " 0.0\n", " 2.8133919327141074e-9\n", " 3.0947311259855185e-8\n", " 3.122865045312659e-7\n", " 3.125678437245373e-6\n", " 3.1259597764386445e-5\n", " 0.00031259879103579717\n", " 0.003125990723749904\n", " 0.015721220942788605\n", " 0.04026863845950178\n", " 0.07423904882138671\n", " 0.11791864254410875\n", " 0.1723420300368352\n", " ⋮\n", " 115.65914474054375\n", " 116.08964142668084\n", " 116.52011842425044\n", " 116.95057564863274\n", " 117.38101312924142\n", " 117.81143080997383\n", " 118.24182866323792\n", " 118.67220671844858\n", " 119.10256500501667\n", " 119.53290360934692\n", " 119.96322253233947\n", " 120.0\n", "u: 298-element Vector{LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}}:\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 295354.0\n", " :E => 0.0\n", " :I => 1000.0\n", " :R => 567191.0\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 295353.99998667586\n", " :E => 1.7320507818745795e-5\n", " :I => 999.9999831126152\n", " :R => 567191.000012891\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 295353.9998534344\n", " :E => 0.00019052556412603002\n", " :I => 999.9998142387877\n", " :R => 567191.0001418012\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 295353.9985210214\n", " :E => 0.001922573939183367\n", " :I => 999.9981255026566\n", " :R => 567191.0014309019\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 295353.9851970526\n", " :E => 0.01924283889007708\n", " :I => 999.9812383556265\n", " :R => 567191.0143217528\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 295353.8519734983\n", " :E => 0.19242361030166427\n", " :I => 999.8123883115076\n", " :R => 567191.1432145798\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 295352.5213499726\n", " :E => 1.9220453840944718\n", " :I => 998.126028662952\n", " :R => 567192.4305759802\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 295339.3749437633\n", " :E => 19.00152652208452\n", " :I => 981.4746974066799\n", " :R => 567205.1488323078\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 295283.93378345115\n", " :E => 90.83469714947545\n", " :I => 911.4633089972851\n", " :R => 567258.768210402\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 295190.2782579937\n", " :E => 211.30502569745482\n", " :I => 794.147911380377\n", " :R => 567349.2688049284\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 295086.5711864948\n", " :E => 342.8704085749534\n", " :I => 666.2506099590504\n", " :R => 567449.3077949712\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 294986.7359150886\n", " :E => 466.5215703090174\n", " :I => 546.4334327500894\n", " :R => 567545.3090818522\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 294899.2517716976\n", " :E => 570.3926060558011\n", " :I => 446.38889145874026\n", " :R => 567628.9667307878\n", " ⋮\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 286890.0357270133\n", " :E => 505.8635280900825\n", " :I => 180.9533416981066\n", " :R => 575968.147403198\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 287045.6077336715\n", " :E => 504.28784443079377\n", " :I => 180.36400939159776\n", " :R => 575814.7404125056\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 287202.1814522802\n", " :E => 502.9027884934137\n", " :I => 179.84284695065017\n", " :R => 575660.0729122752\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 287359.58047511696\n", " :E => 501.7078205713716\n", " :I => 179.38966267303041\n", " :R => 575504.3220416381\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 287517.6291542\n", " :E => 500.7023967219029\n", " :I => 179.00426283902502\n", " :R => 575347.6641862384\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 287676.1525249288\n", " :E => 499.8859736702451\n", " :I => 178.68645323720042\n", " :R => 575190.2750481631\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 287834.9763419593\n", " :E => 499.25801156993975\n", " :I => 178.43604045982232\n", " :R => 575032.3296060103\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 287993.92708219093\n", " :E => 498.817976441694\n", " :I => 178.2528329333104\n", " :R => 574874.0021084335\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 288152.8319148042\n", " :E => 498.5653420410812\n", " :I => 178.13664159720813\n", " :R => 574715.466101557\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 288311.51871292066\n", " :E => 498.4995906622028\n", " :I => 178.08728039535794\n", " :R => 574556.8944160213\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 288469.81599174725\n", " :E => 498.6202134875822\n", " :I => 178.10456626442502\n", " :R => 574398.4592285003\n", " 4-element LArray{Float64, 1, Vector{Float64}, (:S, :E, :I, :R)}:\n", " :S => 288483.304863531\n", " :E => 498.662023366161\n", " :I => 178.08630510848926\n", " :R => 574384.9468079938" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_chickenpox" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n" ], "text/plain": [ "HTML{String}(\"\\n\")" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# to have the figures plotted fix to the wider of the cells\n", "HTML(\"\"\"\n", "\n", "\"\"\")\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.9.1", "language": "julia", "name": "julia-1.9" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.9.1" } }, "nbformat": 4, "nbformat_minor": 2 }