Skip to content

Commit

Permalink
Add piecewise functions for Scenario 2
Browse files Browse the repository at this point in the history
  • Loading branch information
@bgyori
bgyori committed Jul 12, 2023
1 parent 48478ed commit f12adec
Show file tree
Hide file tree
Showing 2 changed files with 426 additions and 411 deletions.
293 changes: 154 additions & 139 deletions notebooks/hackathon_2023.07/scenario2.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,128 +26,61 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 1,
"id": "0991209a",
"metadata": {},
"outputs": [],
"source": [
"import sympy\n",
"from copy import deepcopy as _d\n",
"from mira.metamodel import *\n",
"from mira.modeling import Model\n",
"from mira.modeling.askenet.petrinet import AskeNetPetriNetModel"
]
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 2,
"id": "dc86827c",
"metadata": {},
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"times = [0, 82, 198, 282, 327, 370, 429, 468, 527, 548, 562, 609, 670, 731, 1096]\n",
"beta_scales = [1, 0.4, 0.55, 0.6, 0.55, 0.6, 0.55, 0.3, 0.4, 0.15, 0.2, 0.3, 0.35, 0.4, 0.5, 0.4]\n",
"\n",
"\n",
"\n",
"\"\"\"\n",
"% simulate lockdown\n",
"function [beta_scale] = get_beta(t)\n",
"ind = 1;\n",
"\n",
"% first day of simulation\n",
"day0 = '1-jan-2020';\n",
"\n",
"% state of emergency and such [2]\n",
"tlock(ind) = 82; %daysact(day0,'23-mar-2020');\n",
"tunlock(ind) = 198; %daysact(day0,'17-jul-2020');\n",
"scale(ind) = 0.4;\n",
"ind = ind+1;\n",
"\n",
"% stage 3 [3]\n",
"tlock(ind) = 198; %daysact(day0,'17-jul-2020');\n",
"tunlock(ind) = 282; %daysact(day0,'9-oct-2020');\n",
"scale(ind) = 0.55;\n",
"ind = ind+1;\n",
"\n",
"% oops, back to stage 2 [2]\n",
"tlock(ind) = 282; %daysact(day0,'9-oct-2020');\n",
"tunlock(ind) = 327; %daysact(day0,'23-nov-2020');\n",
"scale(ind) = 0.6;\n",
"ind = ind+1;\n",
"\n",
"% strict lockdown, stay at home order [1]\n",
"tlock(ind) = 327; %daysact(day0,'23-nov-2020');\n",
"% tunlock(ind) = daysact(day0,'5-mar-2021');\n",
"tunlock(ind) = 370; %daysact(day0,'5-jan-2021');\n",
"scale(ind) = 0.55;\n",
"ind = ind+1;\n",
"\n",
"% stricter lockdown, stay at home order [1]\n",
"tlock(ind) = 370; %daysact(day0,'5-jan-2021');\n",
"tunlock(ind) = 429; %daysact(day0,'5-mar-2021');\n",
"scale(ind) = 0.3; %0.25;\n",
"ind = ind+1;\n",
"\n",
"% stay at home order lifted [2]\n",
"tlock(ind) = 429; %daysact(day0,'5-mar-2021');\n",
"tunlock(ind) = 468; %daysact(day0,'13-apr-2021');\n",
"scale(ind) = 0.4; %0.4;\n",
"ind = ind+1;\n",
"\n",
"% lockdown again [1]\n",
"tlock(ind) = 468; %daysact(day0,'13-apr-2021');\n",
"tunlock(ind) = 527; %daysact(day0,'11-jun-2021');\n",
"scale(ind) = 0.15;% 0.1;\n",
"ind = ind+1;\n",
"\n",
"% stage 2 reopening [3]\n",
"tlock(ind) = 527; %daysact(day0,'11-jun-2021');\n",
"tunlock(ind) = 548; %daysact(day0,'2-jul-2021');\n",
"scale(ind) = 0.2;%0.15;\n",
"ind = ind+1;\n",
"beta_scales = [1, 0.4, 0.55, 0.6, 0.55, 0.3, 0.4, 0.15, 0.2, 0.3, 0.35, 0.4, 0.5, 0.4]\n",
"\n",
"% stage 3 reopening [4]\n",
"tlock(ind) = 548; %daysact(day0,'2-jul-2021');\n",
"tunlock(ind) = 562; %daysact(day0,'16-jul-2021');\n",
"scale(ind) = 0.3;%0.25;\n",
"ind = ind+1;\n",
"t = sympy.Symbol('t')\n",
"\n",
"% physical distancing [5]\n",
"tlock(ind) = 562; %daysact(day0,'16-jul-2021');\n",
"tunlock(ind) = 609; %daysact(day0,'1-sep-2021');\n",
"scale(ind) = 0.35;\n",
"ind = ind+1;\n",
"\n",
"% physical distancing [5]\n",
"tlock(ind) = 609; %daysact(day0,'1-sep-2021');\n",
"tunlock(ind) = 670; %daysact(day0,'1-nov-2021');\n",
"scale(ind) = 0.4;\n",
"ind = ind+1;\n",
"\n",
"% physical distancing loosening [6]\n",
"tlock(ind) = 670; %daysact(day0,'1-nov-2021');\n",
"tunlock(ind) = 731; %daysact(day0,'1-jan-2022');\n",
"scale(ind) = 0.5;\n",
"ind = ind+1;\n",
"\n",
"% physical distancing tightened [6]\n",
"tlock(ind) = 731; %daysact(day0,'1-jan-2022');\n",
"tunlock(ind) = 1096; %daysact(day0,'1-jan-2023');\n",
"scale(ind) = 0.4;\n",
"ind = ind+1;\n",
"\n",
"ilock = find(tlock<=t);\n",
"iunlock= find(t<tunlock);\n",
"ii = intersect(ilock,iunlock);\n",
"if (length(ii)>1)\n",
" [t,ii]\n",
"end\n",
"if (isempty(ii))\n",
" beta_scale = 1;\n",
"else\n",
" beta_scale = scale(ii);\n",
"end\n",
"\"\"\""
"beta_scale = sympy.Piecewise(\n",
" *[(b, (tl <= t) & (t < tu))\n",
" for (tl, tu), b in zip(zip(times[:-1], times[1:]), beta_scales)],\n",
" (1, True)\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "496d650e",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\begin{cases} 1 & \\text{for}\\: t \\geq 0 \\wedge t < 82 \\\\0.4 & \\text{for}\\: t \\geq 82 \\wedge t < 198 \\\\0.55 & \\text{for}\\: t \\geq 198 \\wedge t < 282 \\\\0.6 & \\text{for}\\: t \\geq 282 \\wedge t < 327 \\\\0.55 & \\text{for}\\: t \\geq 327 \\wedge t < 370 \\\\0.3 & \\text{for}\\: t \\geq 370 \\wedge t < 429 \\\\0.4 & \\text{for}\\: t \\geq 429 \\wedge t < 468 \\\\0.15 & \\text{for}\\: t \\geq 468 \\wedge t < 527 \\\\0.2 & \\text{for}\\: t \\geq 527 \\wedge t < 548 \\\\0.3 & \\text{for}\\: t \\geq 548 \\wedge t < 562 \\\\0.35 & \\text{for}\\: t \\geq 562 \\wedge t < 609 \\\\0.4 & \\text{for}\\: t \\geq 609 \\wedge t < 670 \\\\0.5 & \\text{for}\\: t \\geq 670 \\wedge t < 731 \\\\0.4 & \\text{for}\\: t \\geq 731 \\wedge t < 1096 \\\\1 & \\text{otherwise} \\end{cases}$"
],
"text/plain": [
"Piecewise((1, (t >= 0) & (t < 82)), (0.4, (t >= 82) & (t < 198)), (0.55, (t >= 198) & (t < 282)), (0.6, (t >= 282) & (t < 327)), (0.55, (t >= 327) & (t < 370)), (0.3, (t >= 370) & (t < 429)), (0.4, (t >= 429) & (t < 468)), (0.15, (t >= 468) & (t < 527)), (0.2, (t >= 527) & (t < 548)), (0.3, (t >= 548) & (t < 562)), (0.35, (t >= 562) & (t < 609)), (0.4, (t >= 609) & (t < 670)), (0.5, (t >= 670) & (t < 731)), (0.4, (t >= 731) & (t < 1096)), (1, True))"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"beta_scale"
]
},
{
Expand Down Expand Up @@ -250,36 +183,67 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 12,
"id": "563bef26",
"metadata": {},
"outputs": [],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"40\n"
]
}
],
"source": [
"all_infectious = c['I'], c['IV'], c['IR'], c['A'], c['AR']\n",
"infection_rate_term = I + IV + IR + ai_beta_ratio * (A + AR)\n",
"all_infectious = [c['I'], c['IV'], c['IR'], c['A'], c['AR']]\n",
"inf_terms = [I, IV, IR, ai_beta_ratio*A, ai_beta_ratio*AR]\n",
"\n",
"templates = [\n",
"templates = []\n",
"for inf, inf_term in zip(all_infectious, inf_terms):\n",
" # Infection asym\n",
" GroupedControlledConversion(subject=c['S'], outcome=c['A'],\n",
" controllers=all_infectious,\n",
" rate_law=ai*beta*S*infection_rate_term),\n",
" GroupedControlledConversion(subject=c['V1'], outcome=c['A'],\n",
" controllers=all_infectious,\n",
" rate_law=ai*beta_v1*V1*infection_rate_term),\n",
" GroupedControlledConversion(subject=c['V2'], outcome=c['A'],\n",
" controllers=all_infectious,\n",
" rate_law=ai*beta_v2*V2*infection_rate_term),\n",
" t = ControlledConversion(subject=c['S'], outcome=c['A'],\n",
" controller=_d(inf),\n",
" rate_law=ai*beta*beta_scale*S*inf_term)\n",
" templates.append(t)\n",
" t = ControlledConversion(subject=c['S'], outcome=c['A'],\n",
" controller=_d(inf),\n",
" rate_law=ai*beta*beta_scale*S*inf_term)\n",
" templates.append(t)\n",
" t = ControlledConversion(subject=c['S'], outcome=c['A'],\n",
" controller=_d(inf),\n",
" rate_law=ai*beta*beta_scale*S*inf_term)\n",
" templates.append(t)\n",
"\n",
" # Infection sym\n",
" GroupedControlledConversion(subject=c['S'], outcome=c['I'],\n",
" controllers=all_infectious,\n",
" rate_law=(1-ai)*beta*S*infection_rate_term),\n",
" GroupedControlledConversion(subject=c['V1'], outcome=c['IV'],\n",
" controllers=all_infectious,\n",
" rate_law=ai*beta_v1*V1*infection_rate_term),\n",
" GroupedControlledConversion(subject=c['V2'], outcome=c['IV'],\n",
" controllers=all_infectious,\n",
" rate_law=ai*beta_v2*V2*infection_rate_term),\n",
" t = ControlledConversion(subject=c['S'], outcome=c['I'],\n",
" controller=_d(inf),\n",
" rate_law=(1-ai)*beta*beta_scale*S*inf_term)\n",
" templates.append(t)\n",
" t = ControlledConversion(subject=c['V1'], outcome=c['IV'],\n",
" controller=_d(inf),\n",
" rate_law=ai*beta_v1*beta_scale*V1*inf_term)\n",
" templates.append(t)\n",
" t = ControlledConversion(subject=c['V2'], outcome=c['IV'],\n",
" controller=_d(inf),\n",
" rate_law=ai*beta_v2*beta_scale*V2*inf_term)\n",
" templates.append(t)\n",
" \n",
" \n",
" # Reinfection asym\n",
" t = ControlledConversion(subject=c['R'], outcome=c['AR'],\n",
" controller=_d(inf),\n",
" rate_law=ai_R*beta_R*beta_scale*R*inf_term)\n",
" templates.append(t)\n",
" # Reinfection sym\n",
" t = ControlledConversion(subject=c['R'], outcome=c['IR'],\n",
" controller=_d(inf),\n",
" rate_law=(1-ai_R)*beta_R*beta_scale*R*inf_term)\n",
" templates.append(t)\n",
"\n",
"print(len(templates))\n",
" \n",
"templates += [\n",
" # Vaccination\n",
" # Note: structurally we want to have these in the model but we make the rates 0\n",
" NaturalConversion(subject=c['S'], outcome=c['V1'], rate_law=0*S),\n",
Expand All @@ -292,15 +256,6 @@
" NaturalConversion(subject=c['A'], outcome=c['R'], rate_law=gamma*A),\n",
" NaturalConversion(subject=c['AR'], outcome=c['R'], rate_law=gamma*AR),\n",
"\n",
" # Reinfection asym\n",
" GroupedControlledConversion(subject=c['R'], outcome=c['AR'],\n",
" controllers=all_infectious,\n",
" rate_law=ai_R*beta_R*R*infection_rate_term),\n",
" # Reinfection sym\n",
" GroupedControlledConversion(subject=c['R'], outcome=c['IR'],\n",
" controllers=all_infectious,\n",
" rate_law=(1-ai_R)*beta_R*R*infection_rate_term),\n",
"\n",
" # Second recovery\n",
" NaturalConversion(subject=c['IR'], outcome=c['R2'], rate_law=gamma*IR),\n",
" NaturalConversion(subject=c['AR'], outcome=c['R2'], rate_law=gamma*AR),\n",
Expand All @@ -323,9 +278,18 @@
" NaturalConversion(subject=c['V2'], outcome=c['SVR'], rate_law=nu_v2*V2),\n",
" NaturalConversion(subject=c['R'], outcome=c['SVR'], rate_law=nu_R*R),\n",
" NaturalConversion(subject=c['R2'], outcome=c['SVR'], rate_law=nu_R*R2),\n",
"]"
"]\n",
"print(len(te))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "31e0ac67",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 8,
Expand All @@ -336,15 +300,66 @@
"tm = TemplateModel(templates=templates,\n",
" parameters={p.name: p for p in parameters},\n",
" initial_conditions=initials)\n",
"tm = simplify_rate_laws(tm)\n",
"tm.annotations= Annotations(name=\"Scenario 2a\")\n",
"AskeNetPetriNetModel(Model(tm)).to_json_file('scenario2_a.json')"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "f0d2cd39",
"metadata": {},
"outputs": [],
"source": [
"tm = TemplateModel(templates=[templates[10]])"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "01eaa75c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"ControlledConversion(rate_law=I*S*beta*(1 - ai)*Piecewise((1, (t >= 0) & (t < 82)), (0.4, (t >= 82) & (t < 198)), (0.55, (t >= 198) & (t < 282)), (0.6, (t >= 282) & (t < 327)), (0.55, (t >= 327) & (t < 370)), (0.3, (t >= 370) & (t < 429)), (0.4, (t >= 429) & (t < 468)), (0.15, (t >= 468) & (t < 527)), (0.2, (t >= 527) & (t < 548)), (0.3, (t >= 548) & (t < 562)), (0.35, (t >= 562) & (t < 609)), (0.4, (t >= 609) & (t < 670)), (0.5, (t >= 670) & (t < 731)), (0.4, (t >= 731) & (t < 1096)), (1, True)), type='ControlledConversion', controller=Concept(name='I', display_name=None, description=None, identifiers={'ido': '0000511'}, context={}, units=Unit(expression=1)), subject=Concept(name='S', display_name=None, description=None, identifiers={'ido': '0000514'}, context={'vaccination_status': 'vo:0001377'}, units=Unit(expression=1)), outcome=Concept(name='I', display_name=None, description=None, identifiers={'ido': '0000511'}, context={}, units=Unit(expression=1)), provenance=[])"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"templates[3]"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "ac6816b8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"ControlledConversion(rate_law=I*V1*ai*beta_v1*Piecewise((1, (t >= 0) & (t < 82)), (0.4, (t >= 82) & (t < 198)), (0.55, (t >= 198) & (t < 282)), (0.6, (t >= 282) & (t < 327)), (0.55, (t >= 327) & (t < 370)), (0.3, (t >= 370) & (t < 429)), (0.4, (t >= 429) & (t < 468)), (0.15, (t >= 468) & (t < 527)), (0.2, (t >= 527) & (t < 548)), (0.3, (t >= 548) & (t < 562)), (0.35, (t >= 562) & (t < 609)), (0.4, (t >= 609) & (t < 670)), (0.5, (t >= 670) & (t < 731)), (0.4, (t >= 731) & (t < 1096)), (1, True)), type='ControlledConversion', controller=Concept(name='I', display_name=None, description=None, identifiers={'ido': '0000511'}, context={}, units=Unit(expression=1)), subject=Concept(name='V1', display_name=None, description=None, identifiers={'ido': '0000514'}, context={'vaccination_status': 'askemo:0000018'}, units=Unit(expression=1)), outcome=Concept(name='IV', display_name=None, description=None, identifiers={'ido': '0000511'}, context={'vaccination_status': 'vo:0001376'}, units=Unit(expression=1)), provenance=[])"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"templates[4]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "34c423df",
"id": "3dfe61aa",
"metadata": {},
"outputs": [],
"source": []
Expand Down
Loading

0 comments on commit f12adec

Please sign in to comment.