{
"cells": [
{
"cell_type": "markdown",
"id": "a40f6857",
"metadata": {},
"source": [
"## Carry out a few single indivual steps of diffusion, \n",
"### and directly verify that the values satisfy the diffusion equation\n",
"\n",
"In this \"PART 2\", we quickly repeat all the steps discussed at length in part1,\n",
"but using a much-finer resolution.\n",
"This time, we'll start with slightly more complex initial concentrations built from 2 superposed sine waves.\n",
"\n",
"**We'll also explore the effects of:** \n",
"-Spatial resolution (\"delta x\") \n",
"-Temporal resolution (\"delta t\") \n",
"-Alternate methods of estimating numerical derivatives\n",
"\n",
"LAST REVISED: June 4, 2023"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "d142d450-f2af-42fa-9026-2b44919c5b40",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Added 'D:\\Docs\\- MY CODE\\BioSimulations\\life123-Win7' to sys.path\n"
]
}
],
"source": [
"import set_path # Importing this module will add the project's home directory to sys.path"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "9897632a-191c-4cdc-9bc2-a1143c2f0d36",
"metadata": {},
"outputs": [],
"source": [
"from experiments.get_notebook_info import get_notebook_basename\n",
"\n",
"from src.life_1D.bio_sim_1d import BioSim1D\n",
"from src.modules.reactions.reaction_data import ChemData as chem\n",
"from src.modules.movies.movies import MovieArray\n",
"from src.modules.numerical.numerical import Numerical as num\n",
"\n",
"import numpy as np\n",
"\n",
"import plotly.express as px\n",
"import plotly.graph_objects as go"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "0c291bcf-fd3a-443d-9df1-1b8ce9fefcae",
"metadata": {},
"outputs": [],
"source": [
"# We'll be considering just 1 chemical species, \"A\"\n",
"diffusion_rate = 10.\n",
"\n",
"chem_data = chem(diffusion_rates=[diffusion_rate], names=[\"A\"])"
]
},
{
"cell_type": "markdown",
"id": "82504138-da32-486a-930d-bc9419cb1dc2",
"metadata": {},
"source": [
"# BASELINE\n",
"This will be our initial system, whose adherence to the diffusion equation we'll test.\n",
"Afterwards, we'll tweak individual parameters - and observed their effect on the closeness of the approximation"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "bca9a060-f62d-433a-9924-1fb2314870f3",
"metadata": {},
"outputs": [],
"source": [
"# Parameters of the simulation run (the diffusion rate got set earlier, and will never vary)\n",
"delta_t = 0.01\n",
"n_bins = 300\n",
"delta_x = 2 # Note that the number of bins also define the fraction of the sine wave cycle in each bin\n",
"algorithm = None # \"Explicit, with 3+1 stencil\""
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "696122b0-c8c4-4d43-a4fc-319f95221602",
"metadata": {},
"outputs": [],
"source": [
"# Initialize the system\n",
"bio = BioSim1D(n_bins=n_bins, chem_data=chem_data)\n",
"\n",
"# Initialize the concentrations to 2 superposed sine waves\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=1, amplitude=10, bias=50)\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=2, amplitude=8)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "338bd51c-6b12-4cf0-8526-797389fbc55c",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/vnd.plotly.v1+json": {
"config": {
"plotlyServerURL": "https://plot.ly"
},
"data": [
{
"hovertemplate": "Chemical=A
Bin number=%{x}
concentration=%{y}",
"legendgroup": "A",
"line": {
"color": "red",
"dash": "solid"
},
"marker": {
"symbol": "circle"
},
"mode": "lines",
"name": "A",
"orientation": "v",
"showlegend": true,
"type": "scatter",
"x": [
0,
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
],
"xaxis": "x",
"y": [
50,
50.544429428667165,
51.088179283950524,
51.630571063807565,
52.170928407051974,
52.70857815921861,
53.242851432961885,
53.77308466117989,
54.29862064106969,
54.81880956733467,
55.333010052784,
55.84059013458597,
56.34092826446227,
56.83341428113769,
57.317450363391465,
57.79245196208926,
58.25784870961231,
58.713085305138996,
59.157622374276286,
59.59093730158345,
60.01252503457715,
60.42189885785704,
60.81859113604284,
61.20215402426825,
61.57216014503327,
61.92820323027551,
62.26989872758102,
62.596884369518115,
62.90882070514144,
63.20539159277945,
63.48630465328596,
63.751291683005014,
64.00010902576838,
64.2325379033164,
64.44838470360496,
64.64748122653477,
64.82968488671305,
64.99487887293242,
65.1429722641268,
65.27390010164028,
65.38762341772012,
65.48412922022146,
65.5634304335874,
65.62556579624368,
65.67059971462322,
65.6986220741107,
65.70974800727217,
65.704117619809,
65.6818956747483,
65.64327123545478,
65.5884572681199,
65.51769020445475,
65.43122946538178,
65.3293569465877,
65.21237646686652,
65.08061318024517,
64.93441295294744,
64.77414170631202,
64.60018472683976,
64.41294594460129,
64.21284718129132,
64.00032736926742,
63.775841742961624,
63.53986100410003,
63.29287046221103,
63.03536915194446,
62.767868928764315,
62.49089354461371,
62.204977705185655,
61.910666110463616,
61.60851248022481,
61.29907856622386,
60.98293315279715,
60.66065104764711,
60.33281206458206,
60.00000000000001,
59.66280160491487,
59.32180555433006,
58.977601415768284,
58.63077861876623,
58.28192542714066,
57.93162791582594,
57.580468954074114,
57.229027196796395,
56.87787608580947,
56.527582862731656,
56.178707595252185,
55.83180221847259,
55.48740959299154,
55.146062581373926,
54.80828314461176,
54.47458146014833,
54.14545506299808,
53.82138801145301,
53.50285007882258,
53.19029597260686,
52.884164582453884,
52.58487825820056,
52.29284211924252,
52.00844339642252,
51.73205080756888,
51.46401396775535,
51.204662835292,
50.954307194392996,
50.71323617540191,
50.481717813388244,
50.2599986458607,
50.04830335027312,
49.846834421928385,
49.65577189281366,
49.47527309182776,
49.305472446787945,
49.146481328529276,
48.9983879373349,
48.86125723186071,
48.735130900642396,
48.62002737619715,
48.51594189165739,
48.42284657979783,
48.34069061424248,
48.269400392563504,
48.20887976090919,
48.159010279725436,
48.11965153006181,
48.09064145988159,
48.071796769724486,
48.062913337001035,
48.063766678128985,
48.07411244765571,
48.093686973444406,
48.12220782693885,
48.1593744274586,
48.20486867941727,
48.25835564129731,
48.31948422515935,
48.387887925409686,
48.46318557549796,
48.544982131167615,
48.63286947883482,
48.72642726762714,
48.82522376357119,
48.92881672437982,
49.03675429325235,
49.148575910068345,
49.263813238324204,
49.38199110613446,
49.502628459594334,
49.625239326778704,
49.74933379063347,
49.874418968999976,
50,
50.125581031000024,
50.25066620936653,
50.374760673221296,
50.497371540405666,
50.618008893865536,
50.736186761675796,
50.851424089931655,
50.96324570674764,
51.07118327562018,
51.17477623642881,
51.27357273237286,
51.36713052116517,
51.455017868832385,
51.53681442450203,
51.61211207459031,
51.68051577484064,
51.74164435870269,
51.79513132058273,
51.8406255725414,
51.87779217306115,
51.90631302655559,
51.92588755234429,
51.936233321871015,
51.937086662998965,
51.92820323027551,
51.90935854011841,
51.88034846993819,
51.840989720274564,
51.791120239090816,
51.7305996074365,
51.659309385757524,
51.57715342020218,
51.484058108342616,
51.37997262380286,
51.26486909935761,
51.13874276813929,
51.00161206266512,
50.85351867147074,
50.69452755321206,
50.52472690817225,
50.344228107186346,
50.15316557807162,
49.951696649726884,
49.74000135413931,
49.518282186611756,
49.2867638245981,
49.04569280560701,
48.795337164708016,
48.53598603224466,
48.26794919243113,
47.9915566035775,
47.70715788075749,
47.41512174179946,
47.11583541754612,
46.80970402739315,
46.497149921177424,
46.17861198854701,
45.85454493700193,
45.525418539851685,
45.19171685538825,
44.85393741862609,
44.51259040700847,
44.16819778152742,
43.82129240474783,
43.47241713726835,
43.122123914190546,
42.77097280320363,
42.41953104592591,
42.068372084174065,
41.71807457285935,
41.36922138123378,
41.02239858423173,
40.678194445669945,
40.337198395085146,
40,
39.66718793541796,
39.3393489523529,
39.017066847202855,
38.70092143377615,
38.3914875197752,
38.0893338895364,
37.79502229481435,
37.509106455386295,
37.23213107123569,
36.964630848055556,
36.70712953778898,
36.460138995899975,
36.22415825703838,
35.99967263073258,
35.787152818708684,
35.58705405539872,
35.39981527316025,
35.22585829368798,
35.065587047052574,
34.919386819754834,
34.7876235331335,
34.670643053412306,
34.568770534618224,
34.48230979554525,
34.41154273188011,
34.35672876454522,
34.31810432525169,
34.295882380191,
34.29025199272783,
34.3013779258893,
34.329400285376764,
34.374434203756316,
34.436569566412594,
34.515870779778524,
34.612376582279865,
34.726099898359706,
34.8570277358732,
35.00512112706758,
35.17031511328693,
35.35251877346523,
35.551615296395035,
35.76746209668359,
35.99989097423161,
36.24870831699498,
36.51369534671404,
36.794608407220544,
37.09117929485854,
37.40311563048187,
37.73010127241896,
38.07179676972448,
38.427839854966706,
38.79784597573175,
39.18140886395714,
39.57810114214295,
39.987474965422834,
40.409062698416534,
40.84237762572369,
41.286914694861004,
41.74215129038767,
42.20754803791071,
42.68254963660851,
43.16658571886228,
43.659071735537715,
44.159409865414,
44.666989947215995,
45.1811904326653,
45.701379358930296,
46.226915338820085,
46.75714856703811,
47.29142184078137,
47.82907159294801,
48.36942893619241,
48.911820716049476,
49.45557057133282
],
"yaxis": "y"
}
],
"layout": {
"autosize": true,
"legend": {
"title": {
"text": "Chemical"
},
"tracegroupgap": 0
},
"template": {
"data": {
"bar": [
{
"error_x": {
"color": "#2a3f5f"
},
"error_y": {
"color": "#2a3f5f"
},
"marker": {
"line": {
"color": "#E5ECF6",
"width": 0.5
},
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "bar"
}
],
"barpolar": [
{
"marker": {
"line": {
"color": "#E5ECF6",
"width": 0.5
},
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "barpolar"
}
],
"carpet": [
{
"aaxis": {
"endlinecolor": "#2a3f5f",
"gridcolor": "white",
"linecolor": "white",
"minorgridcolor": "white",
"startlinecolor": "#2a3f5f"
},
"baxis": {
"endlinecolor": "#2a3f5f",
"gridcolor": "white",
"linecolor": "white",
"minorgridcolor": "white",
"startlinecolor": "#2a3f5f"
},
"type": "carpet"
}
],
"choropleth": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "choropleth"
}
],
"contour": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"type": "contour"
}
],
"contourcarpet": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "contourcarpet"
}
],
"heatmap": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"type": "heatmap"
}
],
"heatmapgl": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"type": "heatmapgl"
}
],
"histogram": [
{
"marker": {
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "histogram"
}
],
"histogram2d": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"type": "histogram2d"
}
],
"histogram2dcontour": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"type": "histogram2dcontour"
}
],
"mesh3d": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "mesh3d"
}
],
"parcoords": [
{
"line": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "parcoords"
}
],
"pie": [
{
"automargin": true,
"type": "pie"
}
],
"scatter": [
{
"fillpattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
},
"type": "scatter"
}
],
"scatter3d": [
{
"line": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatter3d"
}
],
"scattercarpet": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattercarpet"
}
],
"scattergeo": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattergeo"
}
],
"scattergl": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattergl"
}
],
"scattermapbox": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattermapbox"
}
],
"scatterpolar": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterpolar"
}
],
"scatterpolargl": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterpolargl"
}
],
"scatterternary": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterternary"
}
],
"surface": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"type": "surface"
}
],
"table": [
{
"cells": {
"fill": {
"color": "#EBF0F8"
},
"line": {
"color": "white"
}
},
"header": {
"fill": {
"color": "#C8D4E3"
},
"line": {
"color": "white"
}
},
"type": "table"
}
]
},
"layout": {
"annotationdefaults": {
"arrowcolor": "#2a3f5f",
"arrowhead": 0,
"arrowwidth": 1
},
"autotypenumbers": "strict",
"coloraxis": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"colorscale": {
"diverging": [
[
0,
"#8e0152"
],
[
0.1,
"#c51b7d"
],
[
0.2,
"#de77ae"
],
[
0.3,
"#f1b6da"
],
[
0.4,
"#fde0ef"
],
[
0.5,
"#f7f7f7"
],
[
0.6,
"#e6f5d0"
],
[
0.7,
"#b8e186"
],
[
0.8,
"#7fbc41"
],
[
0.9,
"#4d9221"
],
[
1,
"#276419"
]
],
"sequential": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"sequentialminus": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
]
},
"colorway": [
"#636efa",
"#EF553B",
"#00cc96",
"#ab63fa",
"#FFA15A",
"#19d3f3",
"#FF6692",
"#B6E880",
"#FF97FF",
"#FECB52"
],
"font": {
"color": "#2a3f5f"
},
"geo": {
"bgcolor": "white",
"lakecolor": "white",
"landcolor": "#E5ECF6",
"showlakes": true,
"showland": true,
"subunitcolor": "white"
},
"hoverlabel": {
"align": "left"
},
"hovermode": "closest",
"mapbox": {
"style": "light"
},
"paper_bgcolor": "white",
"plot_bgcolor": "#E5ECF6",
"polar": {
"angularaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"bgcolor": "#E5ECF6",
"radialaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
}
},
"scene": {
"xaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
},
"yaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
},
"zaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
}
},
"shapedefaults": {
"line": {
"color": "#2a3f5f"
}
},
"ternary": {
"aaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"baxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"bgcolor": "#E5ECF6",
"caxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
}
},
"title": {
"x": 0.05
},
"xaxis": {
"automargin": true,
"gridcolor": "white",
"linecolor": "white",
"ticks": "",
"title": {
"standoff": 15
},
"zerolinecolor": "white",
"zerolinewidth": 2
},
"yaxis": {
"automargin": true,
"gridcolor": "white",
"linecolor": "white",
"ticks": "",
"title": {
"standoff": 15
},
"zerolinecolor": "white",
"zerolinewidth": 2
}
}
},
"title": {
"text": "Initial System State"
},
"xaxis": {
"anchor": "y",
"autorange": true,
"domain": [
0,
1
],
"range": [
0,
299
],
"title": {
"text": "Bin number"
},
"type": "linear"
},
"yaxis": {
"anchor": "x",
"autorange": true,
"domain": [
0,
1
],
"range": [
32.544724436364255,
67.45527556363574
],
"title": {
"text": "concentration"
},
"type": "linear"
}
}
},
"image/png": "iVBORw0KGgoAAAANSUhEUgAABrsAAAFoCAYAAADn+EAcAAAAAXNSR0IArs4c6QAAIABJREFUeF7s3QmYVNWd//9PLd1d1XSD7LtsbriFCAhGRUSNCiKKYsg/mYkxMUazuQQHTIwxJsJP4vIkMxrGiZNxYmJAEUUwjqIILijR4BLAyCo7Cgi9VHV1Lf/n3LKaXunq7ltV91a97zM83XTfe+73vL6XMPLpc64nkUgkxIEAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIICACwU8hF0u7BolI4AAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIWAKEXTwICCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACrhUg7HJt6ygcAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEECAsItnAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAwLUChF2ubR2FI4AAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIEHbxDCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCLhWgLDLta2jcAQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAcIungEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAHXChB2ubZ1FI4AAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIEDYxTOAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCDgWgHCLte2jsIRQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQIu3gGEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEXCtA2OXa1lE4AggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAYRfPAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAgGsFCLtc2zoKRwABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQIOziGUAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEHCtAGGXa1tH4QgggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAoRdPAMIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAKuFSDscm3rKBwBBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQICwi2cAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEDAtQKEXa5tHYUjgAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggQdvEMIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIuFaAsMu1raNwBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABwi6eAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAdcKEHa5tnUUjgACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAgggQNjFM4AAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIOBaAcIu17aOwhFAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBAi7eAYQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQRcK0DY5drWUTgCCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggABhF88AAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIICAawUIu1zbOgpHAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBAg7OIZQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQcK0AYZdrW0fhCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAAChF08AwgggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAq4VIOxybesoHAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAgLCLZwABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQMC1AoRdrm2d/YVv3LpT1916r274xhRNnTjO/hswIgIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCBgswBhl82g9YdbuHSFbr/nkQZ3+MMDMzV6xAltumtqnPrX3jdvvt5as14PzblJXbuU14134GCFrp95v47u10t3zrhGwUBx2vdqa9gVCkd0x9xHtGTZqgb3OGX40CZ1pV1EKye2tUa77tuWcVavWa+rb5zT5JK7br2mQYhox1zMs/Hg/zyteffcomGD+rWlTM5FAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQCAvBAi7MtzGjoRPqdKcGHalgprTTj62QaiWmu+n+w9mJICxIyDKZMtNCPn7Py9V41Az1cNvfXWibr7uKqsEO+ZC2JXJbjI2AggggAACCCCAAAIIIIAAAggggAACCCCAgBsECLsy3CU7wq7mSmxpZVdHptOW8KW1+z+//C0dM2SA7auN2lJjRyzac21rtZlnYemyVfra1AsIu9oDzDUIIIAAAggggAACCCCAAAIIIIAAAggggAACCDQjQNiV4ceipbCrfjBiSqi/3WFzq4Lqb1WXWj1Uv/S+vbtbK6m6HVVubWN4+ogT6lYQtbStXuP7tBbWpO6X2r7Q/L61rRJT869fT2qc5r7XXK2pOlP17dqzr0HX6q+Wam5rxcbbB5rVUPMXL9fds67VbbMf1vvrNlnjpbZf3P9ZhfXustR90t2WMVV7OltVtjaXdHrW3DaZZh71758yTs2x8fcz/PgzPAIIIIAAAggggAACCCCAAAIIIIAAAggggAACGRcg7MowcWthlwlU6ocxzW1L19zXWntnV+Owa+Wb79WFX2bKzQUz6YZd5vpU4NY4SGqOs6Wt9kwNs2Y/XLfdYXM1Nb72SDWmvjdxwpi6uTYXqKVCovohVv1QqLmvNxfWNZ5r6v49unVJ651lR5qLsUinZ0faxjDlWb9Hbelxhv9oMDwCCCCAAAIIIIAAAggggAACCCCAAAIIIIAAArYIEHbZwtjyIK2FXTd8Y4qmThxXN0Dq/Ksmj6/7ekfDruaqS62A6tOrW7veIdXcCipzn0nnjW2y2qu5OTW3OqylAO+DDzcrGCixtkQ8Ulhjrt+9d3+T+zf2S63semjOTerapbyOp61fb841nRVZqevaGjw117OWwq7mzk3dt7UtKDP8R4LhEUAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBCwVYCwy1bOpoO1N+yqv5LIjrCrue3sTLX1twBsa/iSmm1z2+k1Dr0aByype82eda1GjzjBGio1TnOBWWsBUXOBWuNrUvdqa6jV0vlHenRa8q6/xWBr3un0rKWwqznfVL2NV9Rl+I8AwyOAAAIIIIAAAggggAACCCCAAAIIIIAAAgggkFEBwq6M8kpOCLtSIVJz77Zq78quI7Gl7tfc9nmpwOlIq7Dqv7/M3CedbfhaegdW/TpTQVM2wq7GPqnnwHw9taLsSGFXuj070haRV984p8U2pd7xZlbLcSCAAAIIIIAAAggggAACCCCAAAIIIIAAAggg4GYBwq4Mdy/XYVdz76wyU+7oNoZHYmvpnqmA66bvTNNNP/8P3XLdVXWrupobr/5WiamgqqWA6EgruxqPnYuwy9SQ7vvH2tKz9qzsyvAjz/AIIIAAAggggAACCCCAAAIIIIAAAggggAACCGRVgLArw9yZCrtaCmwaByUthUMdCbvMtU8ueUVXTDpHwUBxE8GW7pnaPm/c2C9o7T+31K1wSg3w/PK3NG7siAZjNh6rpVCruXeAtdTaTIVdZn7mSG3L2Pj+jbdybGkubelZS1sSthSYZfhxZ3gEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBLIuQNiVYfJMhV0m5DDb1NV/B5SZSuOQo7kQqP6Kqfa8s6v+9fW3GDT3TwU1p518rO6ccU2D4OpI15lrG4dB5muNVy4dKdRq6d7mmrkPPa6vTT1fZtu+TIZdpienDB/aJMgzc/v9n5c26FdLc7GrZ6lnpH6PU8/I7N88plk//Jq6dinP8J8AhkcAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAILMChF0Z9E29d6n+LdLdju/0ESfo5uuusi5taau6+uOn3sHU7ahyXT/zftW/vn7IZMYz5z7wi+/r0fnPqyPv7EoFOI0JGwdw9b/fUtBkzmlcp/lac8FRKtB7f90ma+jm3kW2ZNmqBmXVPydTYZe5YUvvDpt03tgm4Z85v6W5pNszM0Yq1EpNuL5/4/FT5zQOKTP4x4ChEUAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBDIqABhV0Z5Gby+AFvr8TwggAACCCCAAAIIIIAAAggggAACCCCAAAIIIICA3QKEXXaLMl6LAi2tUIMMAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEGivAGFXe+W4rk0CrOpqExcnI4AAAggggAACCCCAAAIIIIAAAggggAACCCCAQJoChF1pQnEaAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIICA8wQIu5zXEypCAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBIU4CwK00oTkMAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEHCeAGGX83pCRQgggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAmkKEHalCcVpCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACzhMg7HJeT6gIAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAgTQHCrjShOA0BBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQMB5AoRdzusJFSGAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCKQpQNiVJhSnIYAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIOE+AsMt5PaEiBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQACBNAUIu9KE4jQEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAHnCRB2Oa8nVIQAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIJCmAGFXmlCchgACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggg4DwBwi7n9YSKEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEE0hQg7EoTitMQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQScJ0DY5byeUBECCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggECaAoRdaUJxGgIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAgPMECLuc1xMqQgABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQSFOAsCtNKE5DAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBwngBhl/N6QkUIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAJpChB2pQnFaQgggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAs4TIOxyXk+oCAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAIE0Bwq40oTgNAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEDAeQKEXc7rCRUhgAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAgikKUDYlSYUpyGAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCDhPgLDLeT2hIgQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAgTQFCLvShOI0BBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAAB5wkQdjmvJ1SEAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCQpgBhV5pQnIYAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIOA8AcIu5/WEihBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBNIUIOxKE4rTEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEnCdA2OW8nlARAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIBAmgKEXWlCcRoCCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggIDzBAi7nNcTKkIAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEhTgLArTShOQwABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQcJ4AYZfzekJFCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACaQoQdqUJxWkIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAALOEyDscl5PqAgBBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQCBNAcKuNKE4DQEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAwHkChF3O6wkVIYAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIpClA2JUmFKchgAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAgg4T4Cwy3k9oSIEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAIE0BQi70oTiNAQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAecJEHY5rydUhAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAgggkKYAYVeaUJyGAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCDgPAHCLuf1hIoQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQTSFCDsShOK0xBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBJwnQNhlQ0927gvZMApDIICAmwQCxT6Vlvi0vyLiprKpFQEEbBDweT3q0aVEew6EbRiNIRBAwG0CfbsFtXt/SAm3FU69CCDQYYEenUt0KFSrSG28w2MxAAIIuEugS6ciRWMJVYWj7iqcahFAoM0C/boH23wNFzhDgLDLhj4QdtmAyBAIuEyAsMtlDaNcBGwUIOyyEZOhEHChAGGXC5tGyQjYJEDYZRMkwyDgQgHCLhc2jZIRaKcAYVc74RxwGWGXDU0g7LIBkSEQcJkAYZfLGka5CNgoQNhlIyZDIeBCAcIuFzaNkhGwSYCwyyZIhkHAhQKEXS5sGiUj0E4Bwq52wjngMsIuG5pA2GUDIkMg4DIBwi6XNYxyEbBRgLDLRkyGQsCFAoRdLmwaJSNgkwBhl02QDIOACwUIu1zYNEpGoJ0ChF3thHPAZYRdNjSBsMsGRIZAwGUChF0uaxjlImCjAGGXjZgMhYALBQi7XNg0SkbAJgHCLpsgGQYBFwoQdrmwaZSMQDsFCLvaCeeAywi7bGgCYZcNiAyBgMsECLtc1jDKRcBGAcIuGzEZCgEXChB2ubBplIyATQKEXTZBMgwCLhQg7HJh0ygZgXYKEHa1E84BlxF22dAEwi4bEBkCAZcJEHa5rGGUi4CNAoRdNmIyFAIuFCDscmHTKBkBmwQIu2yCZBgEXChA2OXCplEyAu0UIOxqJ5wDLiPssqEJhF02IDIEAi4TIOxyWcMoFwEbBQi7bMRkKARcKEDY5cKmUTICNgkQdtkEyTAIuFCAsMuFTaNkBNop4Kawa/Wa9bp33nw9NOcmde1S3s4Z5/6yAwcrdP3M+3XLdVdp9IgT2l0QYVe76Q5fSNhlAyJDINBIwPvpJ/Ie/EzeQ4fkqayU98C+5OcHD8pjfT318ZA8FRXyfnYg+bVDB63zWjrinTsrESxVIhBUIhio+1yBgOKlpbI+Bs3Hz79f2qnZc/2dy1TSpUwHPSVKlJUp0alM8U6drI8cCCCQ3wKEXfndX2aHQGsChF2tCfF9BPJXgLArf3vLzBBoTYCwqzUhvo9A/gg4KewKhSO6Y+4jWrJsVQPgu269RlMnjlMuw66FS1do/uLltgRthF0O+vND2OWgZlCKYwW8+/fLv3WzfDt3yLtnt7y7d8lnPu7dI+/+fclgywRZJtCqCTt2HukUFus/QLH+AxXr20+xfv0V79/f+hjrNyD5tf4D0hmGcxBAwKEChF0ObQxlIZAlAcKuLEFzGwQcKEDY5cCmUBICWRIg7MoSNLdBwAECTgm7Nm7dqetuvVcTJ4zRzdddVSdjgqFZdz+sGTdM1/4Dh1jZVe+ZYWWXDX+ACLtsQGSIvBDwb9kk39at8m/eKJ8JtrZsln/rFvk3bbBWZ7XliHftqnjno5To0kVx86tzFyU6f/6xa1clyssV73JU8uvWx85KdO5sXWOube4wIZqnOiRPOCRPqFqecFiekPm9+Vid/Nz6/ee/qkOS+bzanPv59z7/vq8mLPMrduAzeauq5KmskMd8TDOoi/Xpq9iAgbI+9jWB2ADF+iWDMCsQO3pwW7g4FwEEsihA2JVFbG6FgAMFCLsc2BRKQiBLAoRdWYLmNgg4UICwy4FNoSQEMiTghLArtaKrT69uDYKuxlNOrewy2//Nmv2wdu3ZZ53yhwdmNtgOMBWcpb7/ra9OrBs3NcY3v3KRbv75g9b1fXt317x7btG7/9ig2+95xPraKcOHNljFZVZ2rXp7re6ccY2CgWLrHPO11Pnm96n7NL6/+V5qdZr5nJVdGXqY2zMsYVd71LjGjQJme0Dfpg3WCi3/1q3ybt6Y/HzLZvk+3nrEKZltA2ODBik6cJDivfso1ruP4n36Jj/v1UuJ8lSo1dkVWwEe6Z1dJtzz7dol347t8u3aKe+ObdbnZhWb/+Ot8u77NK32x3v0VPToQVYAFj3mOMWGDlP0mGMVHTxU8Z690hqDkxBAwH4Bwi77TRkRATcJEHa5qVvUioC9AoRd9noyGgJuEiDsclO3qBWBjgk4IexKhUOzZ117xHdYmaDq6hvnaNJ5Y+tCp8bbC5qxfjL7Yf1q1rUaNqifGgdpqTHqB2D3zZuv3/95aV1YZUTN18yRWmXWOOxqfF9znyeXvKIrJp2jnXs+1bKVb+s7X59sjdF4foRdHXtmbb2asMtWTgZzgID/w/Uq+sd78q/9h/ybNsr/8Rb5Nm+0thk80hHr1VuxwUMVHTRIsaHHWKuTokOHWh9NuJVPx5HCrnTmaUJC766d8u3YYW3taIVhu3bIZ7Z33PaxfHv3HHEYs5Kt9tQRqv3CFxUdfrJqTzzJ+j0HAghkXoCwK/PG3AEBJwsQdjm5O9SGQGYFCLsy68voCDhZgLDLyd2hNgTsFXBC2GUCKLNSy6yuMgFVS0dz7+xqHG6ZkGrwwD7WO75SR/3rNmze0WQrxObGbfy1+mFXuKZG18+8X2aF2egRJ6TVkPp1EXalRZadkwi7suPMXTIj4Nu2VcV/e0vFb65S0XtrVPT+u9aWfi0dqZVFscFDFBs0WNFhxyo2cJCiQ4bIrN4qlKOjYVc6TlbotWuHFTj6Nm6Qf8NHyZV0H65vcbvE2lO+YAVf0RNPqQvDTDDGgQAC9gkQdtlnyUgIuFGAsMuNXaNmBOwRIOyyx5FREHCjAGGXG7tGzQi0TyCfwq5+vXvojrmPaMmyVU0wUtsS2hF2mZVbcx98XLNvu1Zdu5Q3C59aQVb/m6nVZIRd7XtWM3IVYVdGWBk0QwJF7/5dJW++rqI3XlfJqtfk/WRvkzuZd17VjjhNkVNHKGa2zxs8RNFBQxTr1z9DVblv2GyEXUdSMe9F869bK9NPK6Rc+4G1Oqy5w7wDrPakU1R78qmKfnGkak8YboWUHAgg0D4Bwq72uXEVAvkiQNiVL51kHgi0XYCwq+1mXIFAvggQduVLJ5kHAq0LOCHsass2hvfOm9/gXVr1V3alwq6xI09ssLKrvkI6q7jM+Uda2dVa2GVWcS196c0GK9Xqb4tI2NX6c5m1Mwi7skbNjdoo4KmqVPFbb6r4zddV/NYb1ueecKjBKInyckVOG6XIqDGqPW2UFYiYd0RxHFkg12FXc9V5KipUbIVff5d/7Qcqev89Ff3j/WYnkigJHN4G8fgTrL7XnnJqQa3O4xlHoL0ChF3tleM6BPJDgLArP/rILBBojwBhV3vUuAaB/BAg7MqPPjILBNIRcELY1fi9Wo3rfn75WzpmyADtP3CoyRaEzW1jaK5PvWur8Vh2hF1H2sYwNZdpk8c32OKQsCudpzEH5xB25QCdWzYrYN73ZAVbr61U8eo3rS0JGx+1w09U7WmjVTt6jCIjR6t2+ElotkPAiWFXS9OwVn6tWyv/+++qaO37KlrzjryHDjV7enTIsOQKsFNOtd4HFhk9VmyD2I4HhEvyWoCwK6/by+QQaFWAsKtVIk5AIG8FCLvytrVMDIFWBQi7WiXiBATyRsAJYZfBTK3umjhhTIOgqv4qqXTCrtT2gXfdek3d6i6zkuq/H39O13/jMn2wflOH39kVDBTL1PXWmvV1q8xMyPXkklc08bwzNOe3f1SfXt3q5pGqiW0MHfjHhrDLgU0pkJJMiFH81ioVv/mGFXL5du5oMPP4UUclV2yZYGtUMtxKlJUViE5mp+mmsKs5CfOsFK37h7UFov+9d60VYP5NG5pFM1sgRsacocjZ41XzpbMU794js7iMjoDDBQi7HN4gykMgwwKEXRkGZngEHCxA2OXg5lAaAhkWIOzKMDDDI+AgAaeEXYYktSqq/ju3Uu/aMu/Gam5VVuOVXfWDs1179tVJp8IvO1Z2mbDLHCbw+v2fl9bdo3GY9f66Tdb3zNdTh1lxxjaGDvoDQNjloGbkeSklb7yq4jdeU/HrK5NbElZXNZhx7akjFBl1umpHnq7IqNG8lymDz4Pbw67maMwWl9YqsH98IP+H61X09lsq/vvbTU6NHne8asZNUOTsc1Rz9nhWfmXwOWNoZwoQdjmzL1SFQLYECLuyJc19EHCeAGGX83pCRQhkS4CwK1vS3AeB3As4KezKvYa7KvCVH6QDAAAgAElEQVQkEomEu0p2XrWEXc7rSb5UVPzO31Ty1yUqef1VFa96rcG04j17qeb0sao1K7dGjVbkiyN531IWG5+PYVezAVioWsWrzNaYK1Ty+mvWu98aH+adbzVnn6PIOeepZtz4LHaBWyGQGwHCrty4c1cEnCJA2OWUTlAHAtkXIOzKvjl3RMApAoRdTukEdSCQeQHCrswbZ+oOhF02yBJ22YDIEHUCJS+/qMDzSxVc9KS8n37SQCYy9kzVfPlihS+4kHdt5fiZKZSwqzGzJxV+vfqKSla+IhPINj6s4Gv8+ckA7LRROe4Ut0fAfgHCLvtNGREBNwkQdrmpW9SKgL0ChF32ejIaAm4SIOxyU7eoFYGOCRB2dcwvl1cTdtmgT9hlA2IBD2G2IjThVuCvSxR4bkmDrQnN+7XC512o8GVTVXPOeWwX56DnpFDDribhV1WlSszWmib8WrHc2gax/pEoL7e2OqwZd671y2yByIGA2wUIu9zeQepHoGMChF0d8+NqBNwsQNjl5u5ROwIdEyDs6pgfVyPgJgHCLjd1q2GthF029I6wywbEAhvCu+9TBZ99WiVLnlHJq6/IE4nUCcR691H4okmqmXSpFRIkiooKTMcd0yXsar5P3kOHku+UW/mK9WwXrf1AqrdbbqxvP2urw8i4Cao59zzFevV2R8OpEoF6AoRdPA4IFLYAYVdh95/ZF7YAYVdh95/ZF7YAYVdh95/ZF5YAYZd7+03YZUPvCLtsQCyAIfwf/VPBpYtVsvSZ5NZv9QKA6DHHKjRxsmomXqrIyNGSx1MAIu6eImFXev3zHjigktdWWCu/zK+iD9c3fPaPPU4150ywVn1FzjqH1YvpsXJWjgUIu3LcAG6PQI4FCLty3ABuj0AOBQi7cojPrRHIsQBhV44bwO0RyKIAYVcWsW2+FWGXDaCEXTYg5ukQxatXKbDkGWt7Qv/GjxrMMjJ6jGouukShSZfKhF0c7hIg7Gpfv8yqRiv8WrFcJa+vlP+fHzb8czFytGrOGmet/IqMPUOJkkD7bsRVCGRQgLArg7gMjYALBAi7XNAkSkQgQwKEXRmCZVgEXCBA2OWCJlEiAjYJEHbZBJmDYQi7bEAn7LIBMY+GCLz4vEr+ukTBZxbKu39/g5mZbdvCF05S6PIrFe/eI49mXXhTIeyyp+e+vXusLQ+td369tlL+TRsahl+jTlfNmWer9oyzVDP2TJn32HEgkGsBwq5cd4D7I5BbAcKu3PpzdwRyKUDYlUt97o1AbgUIu3Lrz90RyKYAYVc2te29F2GXDZ6EXTYgungIT0WFgs8uUsnzSxVY9oI8oeq62SRKOyk0eYpqLpqk8IQLlOjEP9S7uNUNSifsykwnfTt3WKFX8fIXFfjrEnkPHmxwI7Pqy7zPLjT5csX69M1MEYyKQCsChF08IggUtgBhV2H3n9kXtgBhV2H3n9kXtgBhV2H3n9kXlgBhl3v7TdhlQ+8Iu2xAdNkQnqpKlT71hAIL56tkxfIG1Zt/gA9PulThiy6RWcnFkZ8ChF3Z6avZCtQKkl94XkX/eL/BTSOnjVL40qkKT56i6KAh2SmIuyAgibCLxwCBwhYg7Crs/jP7whYg7Crs/jP7whYg7Crs/jP7whIg7HJvvwm7bOgdYZcNiC4ZovjN11X6p/9V8Mn58oRDdVVHjz1OoUlTVHPxJYqMHO2S2VBmRwQIuzqi175rfdu3Kbj4KQWeWSQTgtU/aoefpPCUqQpPnKzaE09u3w24CoE0BQi70oTiNATyVICwK08by7QQSEOAsCsNJE5BIE8FCLvytLFMC4FmBAi7Wn8sDhys0PUz79fR/XrpzhnXKBgobv2iLJxB2GUDMmGXDYgOHsJ74IBK//gHlf7pf+T/6J91lZp/UA9PnabQpZcrOvQYB8+A0jIhQNiVCdX0x/R++omCS55R4JmFKnnl5QYXRgcPVejSy1Qz8VJFRp2e/qCciUCaAoRdaUJxGgJ5KkDYlaeNZVoIpCFA2JUGEqcgkKcChF152limhUAzAoRdrT8Wq9es14LFy3WoslozbpiuYYP6tX5RFs4g7LIBmbDLBkQHDlGy8hWVPvqIgk8tqKsuesyxCl05XaHLrpT5nKNwBQi7nNN772efWe/3CphVX88vbVBYrF9/hSdeam0tWnP2Oc4pmkpcLUDY5er2UTwCHRYg7OowIQMg4FoBwi7Xto7CEeiwAGFXhwkZAAHXCBB2td6q++bN19ljTtXKN9/T4IF9NHXiuNYvysIZhF02IBN22YDokCHqVnE9+oj8mzdaVSU6lan6q19X9f/3r6o9dYRDKqWMXAsQduW6A83f31NdpcCLzytgVn09t0Tm96kjftRRyeBr8mUKX3CRMydAVa4QIOxyRZsoEoGMCRB2ZYyWgRFwvABhl+NbRIEIZEyAsCtjtAyMgOMEHBV27dkjrVuXfaPevaXhw5u9r9nCcPZvHtOsH35NGzbvsFZ4OWUrQ8IuGx4Vwi4bEHM8hLWK639+r8DSZ+SJRKxqIqeNUvXV31Zo6jQlAsEcV8jtnSZA2OW0jjStx/xZLnnlJQUWL1Lgr8/Ku39/3UmJ8nKFL7hY4clTFD7/QiWCpc6fEBU6RoCwyzGtoBAEciJA2JUTdm6KgCMECLsc0QaKQCAnAoRdOWHnpgjkRMBRYdf//q/0r/+afYd/+Rfp0Uebva/ZwtCs6Lr5uquUenfXLdddpdEjTsh+nY3uSNhlQwsIu2xAzMEQ1iqux/7H2qrQv2mDVYH5B/DqK6er6tvXK3p87v+A5oCFW6YpQNiVJpRTTovFVPL6qwo8u0iBpYvl27WzrjITZtece57Cl1ym8MWXKN65s1Oqpg6HChB2ObQxlIVAlgQIu7IEzW0QcKAAYZcDm0JJCGRJgLArS9DcBgEHCDgq7HrhBenuu7OvcsEF0m23NblvKBzRHXMf0bTJ4+vCLbOloTlM+JXrg7Dr8w6YRPLqG+dYvztl+FA9NOcmde1Sbv1+4dIVuv2eR6zPJ503tsmyPMKuXD/Gbbs/q7ja5sXZzQsQdrn4yUgkVPz26mTw9ewz8m/ZVDeZRFGRImeebQVfoUumKN6jp4snSumZEiDsypQs4yLgDgHCLnf0iSoRyIQAYVcmVBkTAXcIEHa5o09UiYAdAo4Ku+yYkI1jbNy6U9fdeq927dnXYNTGeYqNt2zTUIRdkkzQde+8+Q0CrpRi4+81l1QSdrXpmcvJyc2+i4tVXDnpRb7clLArXzopFa39ILnV4bOLVLRu7eGJeb2KjB5rbXUYuuQyxQYMzJ9JM5MOCRB2dYiPixFwvQBhl+tbyAQQaLcAYVe76bgQAdcLEHa5voVMAIG0BQi7WqYyi4JWvb22wWKg5lZ7pY1t84kFH3aZfSVn3f2wZtwwXcMG9WvCa8KtwQP7aOrEcdb3mgvGCLtsfirtGi6RUN0qrucW8y4uu1wZxxIg7MrPB8Gs8go+s0glzy5S8d/flhKJuonWnjpC4cmXKTRpiqLHHZ+fAMwqLQHCrrSYOAmBvBUg7Mrb1jIxBFoVIOxqlYgTEMhbAcKuvG0tE0OgiQBhV/MPRSrUGjvyxLqsJHWmCcG2bNud860MCz7sam7p3be+OtFqTHMNNOf/ZPbD+tWsa+vCsd37w/zPgoMEvPv2qfSP/63gI/8p3/ZtVmWJTmWq/urXFbrmO4oeP9xB1VKKWwVKir0qLfbpQGWtW6dA3a0I+HbvUuCZp1SyeJGKX1/Z4OzokGEKX3aFai67QrUnn4plgQl4vR5171ysTz6rKbCZM10EEDACvbsGtPdAWId/HAIXBBAoFIFu5cWqCNeqtpb/BSiUnjNPBFICnTv5FY0lVB2OgYIAAnku0KdbIM9nmL/TK/iwy6zUWrB4ed3SO7PS6/qZ9+uqyeN18YSxTV641lzYFa/3k//5+6g4fGaJhDwvvyz97nfyPP20VJsMIBKnny595ztKfPWrUjDo8ElQnpsEPPLI/F+CP/9ualv7a92/X55Fi6QnnpDnpZekSOTwWMcfr8RXvqLE174mHXts++/Bla4R8EjyeDzi73/XtIxCEbBVwMuff1s9GQwBNwmYP//m//1PEHe7qW3UioAtAta/AZh/Z+LPvy2eDIKAkwXM3/cc7hQg7GoUdpk2pvaenPmDr2vOb/+o+kvzmgu72MYwdw+/99NP1OmxRxV89BH5t25O/j8eZhXXVV9V9TXfUe3wk3JXHHfOawG2Mczr9h5xcp6KCgWfX6rAoicU+OuSBueaVV7hSy9XaOo0RQcPLVykPJ852xjmeYOZHgKtCLCNIY8IAoUrwDaGhdt7Zo4A2xjyDCBQOAJsY+jeXhd82GXCq7kPPq7Zt12rrl3KrU7W32OSd3Y58+EO/N9zCv7pUQWffbquwMgXRyr0jW+p+sqvKBFgFZczO5c/VRF25U8vOzITT1Wlgk8vVMD8WvZ/DYMv846vy660wvdYn74duQ3XOkyAsMthDaEcBLIsQNiVZXBuh4CDBAi7HNQMSkEgywKEXVkG53YI5FCAsCuH+B28dcGHXan3cvXp1c16T1dqG8NbrrtKo0ecILPN4b3z5uuhOTdZYZgJv8xhzk0drOzq4FOY5uW+HdtV+of/Uumf/1fmXTqpo+rrV6v6uu+xiitNR06zR4Cwyx7HfBrF+9lnCix5WsEn/6KSFcsbTC0yeozCU6YqdPk0xXr3yadpF+RcCLsKsu1MGoE6AcIuHgYECleAsKtwe8/MESDs4hlAoHAECLvc2+uCD7tM61IB1/vrNlmdvOvWazR14ri6rpqVXrff84j1+0nnja17v1fqBMKuzP4B8G/aoLJ756j0L3+qu1GsV29Vf+s6VX3zWsW7dc9sAYyOQDMChF08FkcS8B44oOCzixQwwderKxoGX2PPVPVXv67wxEsV79oVSBcKEHa5sGmUjICNAoRdNmIyFAIuEyDsclnDKBcBGwUIu2zEZCgEHC5A2OXwBh2hvJyGXY1Dpvp1njJ8aN1qKqfzEnZlpkNFaz9Q2X33KLjoibobmHdwVf3gZmtbMA4EcilA2JVLfXfd27dnt4JPLVDgmUUqfuuNBsXXjJ+g0BVfUXjSFMU7d3bXxAq4WsKuAm4+U0dAEmEXjwEChStA2FW4vWfmCBB28QwgUDgChF3u7XVOw67mtgR0IyVhl71dK1rzjsrv+ZXMe7lSR+T0M1R50wyFL7jI3psxGgLtFCDsaidcgV/m27tHwYXzFXjqCRW/vbpOI1FcrJrx5yl8+TSFJk1WorRTgUs5e/qEXc7uD9UhkGkBwq5MCzM+As4VIOxybm+oDIFMCxB2ZVqY8RFwjgBhl3N60dZKchZ2mVVds+5+WDNumK5hg/q1tW5HnU/YZU87zDtuyu6b02DLr/D5F6ryxh8rMvZMe27CKAjYJEDYZRNkAQ/j27lDwSfnW6u+it5bczj4KgkofMGFCk+dpvCFE5UoCRSwkjOnTtjlzL5QFQLZEiDsypY090HAeQKEXc7rCRUhkC0Bwq5sSXMfBHIvQNiV+x60twLCrvbK1buOsKtjiIG/LlHZ/fccXuXg8yl06VRV3nyrzLaFHAg4UYCwy4ldcW9Nvm1bVfrkfAUWLpDZwjV1mBVe4YsmKXT5lao578syK8A4ci9A2JX7HlABArkUIOzKpT73RiC3AoRdufXn7gjkUoCwK5f63BuB7AoQdmXX28675SzsMpMw2xgOHthHUyeOs3NOWR+LsKsd5PG49S6usvvnqmjdP6wBEsUlqv7q11X5o1sUO3pwOwblEgSyJ0DYlT3rQruTf8smBRc8bv1vpP/D9XXTN+/0Ck+cbG11GB5/nuTzFRqNY+ZL2OWYVlAIAjkRIOzKCTs3RcARAoRdjmgDRSCQEwHCrpywc1MEciJA2JUTdltumtOwa+PWnXps4Yuacf10BQPu/Wl1wq70n0VPba2C8/+ksgd+Lf/mjcmQq6xMVVd/W5Xfu1Hxnr3SH4wzEcihAGFXDvEL6NZF69Yq8NQCBRc9Kf+mDYeDr27dFJo0ReGpV6nmzLMlr7eAVHI/VcKu3PeAChDIpQBhVy71uTcCuRUg7MqtP3dHIJcChF251OfeCGRXgLAru9523i1nYZd5Z9f1M+/X++s2NTufU4YP1UNzblLXLuV2zjcjYxF2tc7qqQmr06OPqNNv75d5T4054t17qOq731fVt74rs2KBAwE3CRB2ualb+VFr0fvvKrDoCZU+9aR8H2+pm5T5IYHQlKnWVoeR08+QPJ78mLCDZ0HY5eDmUBoCWRAg7MoCMrdAwKEChF0ObQxlIZAFAcKuLCBzCwQcIkDY5ZBGtKOMnIVd7ajVsZcQdrXcGk9lpcr+63fq9LvfyvvpJ9aJsQEDVfn9G1X9L99UoiTg2L5SGAJHEiDs4vnIpUDx399W4KknFHxmoXzbt9WVEuvXX6EpVyhsgq/TRuWyxLy+N2FXXreXySHQqgBhV6tEnIBA3goQduVta5kYAq0KEHa1SsQJCOSNAGGXe1tJ2GVD7wi7miJ6DxxQp4d+o06//528Bw9aJ0SPP0GVP/yxqq+4SvL7bZBnCARyJ0DYlTt77txQoHj1m4eDr927DgdfRw9W9eVXWFsd1p50Cmw2ChB22YjJUAi4UICwy4VNo2QEbBIg7LIJkmEQcKEAYZcLm0bJCLRTgLCrnXAOuCznYdfqNet19Y1zGlD84YGZGj3iBAfwpFcCYddhJ9+e3dZWhWbLQk91lfWNyMjRqvzRjxW++BK210rvkeIsFwgQdrmgSYVWYiKhklWvJYOvxYvk/WRvnUD0mGMVunyaQtOmKzr0mEKTsX2+hF22kzIgAq4SIOxyVbsoFgFbBQi7bOVkMARcJUDY5ap2USwCHRIg7OoQX04vzmnYZYKue+fNb/Buro1bd+q6W+/VDd+YoqkTx+UUJ92bE3bJ2kar7L7/p9LHH5MnUmPR1YyfoMobb1XNWe7oY7r95jwEjABhF8+BowXicZW8tlKBpxYo+OzT8u7fV1du7fCTFL7iKlVfOd3aVpaj7QKEXW034woE8kmAsCufuslcEGibAGFX27w4G4F8EiDsyqduMhcEjixA2OXeJyRnYVcoHNEdcx/RtMnjm6ziMiHYgsXLdeeMaxQMFDtet5DDLv+Gj1R272yVLlwgxWKS16vQpCmqvPlW1Z7yBcf3jgIRaK8AYVd75bgu6wKxmEpWvKzgU08osPQZeT/77HDwNeK05IqvK65SrE/frJfm1hsSdrm1c9SNgD0ChF32ODIKAm4UIOxyY9eoGQF7BAi77HFkFATcIEDY5YYuNV9jzsKuAwcrNOvuhzXjhukaNqhfg+rM6q65Dz6u2bddq65dyh2vW4hhV9H776rs17MVXLpYSiSUKC5W6MrpqrxphqJDhjm+ZxSIQEcFCLs6Ksj1ORGIRhV4+cXkiq/nnpWnoiJZhsejyKgxCl0xTaHLrlS8R8+clOeWmxJ2uaVT1IlAZgQIuzLjyqgIuEGAsMsNXaJGBDIjQNiVGVdGRcCJAoRdTuxKejXlLOxiZVd6DXLaWcV/e0tl985R4IW/WqUlgqWquvrbqvr+jYr17uO0cqkHgYwJEHZljJaBsygQWLpYgcVPKbj4aXnCobo713zpbIWnTFXoiq8oftRRWazIHbci7HJHn6gSgUwJEHZlSpZxEXC+AGGX83tEhQhkSoCwK1OyjIuA8wQIu5zXk3QrylnYZQpcuHSF5i9ezju70u1WDs8reXmZyu+do+JVr1lVxLt2VdU116nquz+wPudAoNAECLsKreP5PV9PTViBF55XYOECBZ9Z2GCy5v2LZuWu2aI2Ue781dbZ6BRhVzaUuQcCzhUg7HJub6gMgUwLEHZlWpjxEXCuAGGXc3tDZQjYLUDYZbdo9sbLadhlpmnez3X1jXMazPgPD8xs8h6v7JG0/U75vI1h4LlnVfbAXBW/vToZcnXrrsof/VhV3/y2EqWd2o7FFQjkiQBhV540kmk0EfCEqmX+t996x9dzzzb4fvj8CxW+7AqFLr28oP8OIOziDw4ChS1A2FXY/Wf2hS1A2FXY/Wf2hS1A2FXY/Wf2hSVA2OXefuc87HIv3eHK8zHsKnl9pTr/bJaK1rxjTTTWt5+qfnCTqr7xLSVKAvnQNuaAQIcECLs6xMfFLhHwVFUquGSx9Y6vkldekicSsSpPBIIKX3ChwlOnKXzhRCWKS1wyI3vKJOyyx5FREHCrAGGXWztH3Qh0XICwq+OGjICAWwUIu9zaOepGoO0ChF1tN3PKFYRdNnQin8KuovffVedf3K6Sl1+sC7kqb5lpvZeLAwEEDgsQdvE0FJqAp6JCwWcXKfjkX1Sy/KW66Zt3N4YnXqLQ5dMUvmhSQbAQdhVEm5kkAi0KEHbxcCBQuAKEXYXbe2aOAGEXzwAChSNA2OXeXhN22dC7fAi7/Fs2qfyXP1fw6SelRELxo45S5Y9mqOo717OSy4ZnhCHyT4CwK/96yozSF/B+9pn1bq/AoidU8tpKKRazLo537qzwpEsVNsHXORMkny/9QV10JmGXi5pFqQhkQICwKwOoDImASwQIu1zSKMpEIAMChF0ZQGVIBBwqQNjl0MakUVbWw64DByt0/cz79c2vXKT//stf9f66Tc2WecrwoXpozk3q2qU8jWnk9hQ3h13eT/aqfM5d6vTY/0jRqLU1VdW131XFzTOVKHe+fW47z90LWYCwq5C7z9zrC3j371Nw0ZMKPPWESt58XYrHk8FXt+4KXTJF4alXqeZLZ0leb97AEXblTSuZCALtEiDsahcbFyGQFwKEXXnRRiaBQLsECLvaxcZFCLhSgLDLlW2zis562JWiMqHXrLsf1owbpmvYoH4NBFevWa8Fi5frzhnXKBgodryuG8Musx1V2QP3qGzeg/KEQ9ZP31dP/7oqZv1MsT59HW9OgQjkWoCwK9cd4P5OFPB++olKn1qgwMInVPy3N62VwuaI9eqt8KVTFZp6pSKnn+HE0ttUE2FXm7g4GYG8EyDsyruWMiEE0hYg7EqbihMRyDsBwq68aykTQqBFAcIu9z4cjgy7Nm7dqbkPPq7Zt13Lyi6bny0TbJX950Pq9Jtfy2xDZY7QpEtV8bO7FB12rM13YzgE8leAsCt/e8vM7BHw7d6l4ML51oqv4r+/XTdorP8AhaZcofDlVyryxZH23CzLoxB2ZRmc2yHgMAHCLoc1hHIQyKIAYVcWsbkVAg4TIOxyWEMoB4EMChB2ZRA3w0M7MuxauHSFVr29lpVddjY/GlWnPz2qsnt+JfMPkOYwP11/6K45iowcbeedGAuBghAg7CqINjNJmwR827epdOECBZ5aoKL3360bNTpoiEKXX2FtdVh74sk23S3zwxB2Zd6YOyDgZAHCLid3h9oQyKwAYVdmfRkdAScLEHY5uTvUhoC9AoRd9npmc7Ssh11m1dZ1t96rXXv2tTjPvr27a949tzTZ3jCbMG25l6O3MUwkFHxmocp/daf8mzZY06odfpK1kit8wUVtmSbnIoBAPQHCLh4HBNon4Nu2VaULHrdWfBWt+0fdINFjj1No6lUKTZuu6OCh7Rs8S1cRdmUJmtsg4FABwi6HNoayEMiCAGFXFpC5BQIOFSDscmhjKAuBDAgQdmUANUtDZj3sSs3rSO/sytLcbbuNU8OukpeXqfNdt6vovTXWXGMDj1bFrDtUPW265PHYNn8GQqAQBQi7CrHrzNluAfNDGMEFjyu46An5P/pn3fC1J52i8OXTVH3lVxQbMNDu23Z4PMKuDhMyAAKuFiDscnX7KB6BDgkQdnWIj4sRcLUAYZer20fxCLRJgLCrTVyOOjlnYZejFDpYjNPCrqI176jzz2ap5PWV1szi3bqr8pZ/U9U3v6NEcXEHZ8vlCCBgBAi7eA4QsFegaO0H1mqv4KIn5d+8sW5w814v834vs+or1qevvTdt52iEXe2E4zIE8kSAsCtPGsk0EGiHAGFXO9C4BIE8ESDsypNGMg0E0hAg7EoDyaGnEHbZ0BinhF3+jR+p/K6fKfjs09asEqWdVPnd76vyh7coUVZmw0wZAgEEUgKEXTwLCGROwKxIDj61QMGnF8r38dbkjTweRUaPVWjqNOuX+UGOXB2EXbmS574IOEOAsMsZfaAKBHIhQNiVC3XuiYAzBAi7nNEHqkAgGwKEXdlQzsw9chp2Hen9XacMH6qH5tykrl3KMzNzG0fNddjl271L5bN/odLH/yjFYpLfr6p/+aYqbv2J4j172ThThkIAgZQAYRfPAgLZESh+52/JFV9PPynfzh3Jm3q9qjnzbGurw9ClUxU/6qjsFPP5XQi7ssrNzRBwnABhl+NaQkEIZE2AsCtr1NwIAccJEHY5riUUhEDGBAi7Mkab8YFzFnaFwhHdMfcRjR15or5w0jF6bOGLmnH9dAUDxbpv3nydPeZUjR5xQsYB7LhBrsIu72efqey+/6dOv58nT03Y+ql3849+FbffqejgoXZMjTEQQKAFAcIuHg0Esi9Q/NYbCi58QoHFT8m3Z3eyAL9f4XHnKmxWfE2aokR55n9IhrAr+73njgg4SYCwy0ndoBYEsitA2JVdb+6GgJMECLuc1A1qQSCzAoRdmfXN5Og5C7sOHKzQrLsf1owbplvzm/vg45p927XWSq7Va9ZrweLlunPGNVb45fQj22GXp7pKZb/7d5X99j55Kiosnpovna1Dd89V7cmnOp2L+hDICwHCrrxoI5Nwq0AioZI3Xk2u+Fq8SN5PP7FmYt5LWTPhAmubw/DFlygRLM3IDAm7MsLKoAi4RoCwyzWtolAEbBcg7LKdlAERcI0AYZdrWkWhCHRYgLCrw4Q5G8ARYVe3o8o1+zePadYPv2aFXWZ7w/rhV8500rxxtsIuT22tSh99ROW/ni3vJ3ut6mpPHaFDP/ulasZPSLNaTkMAAX+F/HoAACAASURBVDsECLvsUGQMBGwQiMVU8uoKBRY9oeCSp+Xdv98aNBEIKvzlixW+/EqFL7xYieISG26WHIKwyzZKBkLAlQKEXa5sG0UjYIsAYZctjAyCgCsFCLtc2TaKRqBdAoRd7WJzxEU5C7vqb2M4deI4a+vCwQP7yHy+cOkKrXp7LSu76j0ipU/8ReV33ynfx1usr5ptCit+codCl09zxINEEQgUmgBhV6F1nPm6RSCw7P+SwdfiRfJUVlplmxVeVd/6jqInnarw+Rcq3rVrh6ZD2NUhPi5GwPUChF2ubyETQKDdAoRd7abjQgRcL0DY5foWMgEE0hYg7EqbynEn5izsaixhtjW8fub9en/dJvXt3V3z7rlFwwb1cxxYcwVlcmVX4MXnVf6L21W09gPr1vEePVUx4zZVfes6V9hQJAL5KkDYla+dZV75JBB44a8KLnhcgeeelSdUXTe1yOixqvnyRar65ncUP+qoNk+ZsKvNZFyAQF4JEHblVTuZDAJtEiDsahMXJyOQVwKEXXnVTiaDwBEFCLvc+4A4JuxyL6GUibCrePWb6nznT1W86jWLJtGpTJU/vFmV1/9AidJObuaidgTyQoCwKy/ayCQKSMAEXsGnnlBg6WJ5wqG6mZtVXjVnn6uaCy5UzbnnK9anb6sqhF2tEnECAnktQNiV1+1lcggcUYCwiwcEgcIVIOwq3N4z88ITIOxyb89zFnaZlVyz7n5YM26Y7poVXC212c6wy//henX+xU8VeH5p3e0qr71BlTNmKd6tu3ufNCpHIM8ECLvyrKFMp2AEPDVh6+9Ys+Kr5JWX5amuajD36HHHq+bs8aqZcIFqzhpn/bBJ44Owq2AeFyaKQLMChF08GAgUrgBhV+H2npkjQNjFM4BA4QgQdrm314RdNvTOjrDLt2O7yu/+uUoXPC7F45LXq+orrlLFbXcoNnCQDVUyBAII2ClA2GWnJmMhkCOBWEzF76y2Qq+Sl5ep6J3V8tTWHi7G71fki6NUM36C9Ssyaozk84mwK0f94rYIOESAsMshjaAMBHIgQNiVA3RuiYBDBAi7HNIIykAgCwKEXVlAztAtchZ2mfncN2++zh5zqkaPOCFD08vOsB0Ju7z7PlXZr+eo7OEH64oNn3+hKu74pWqHn5SdCXAXBBBoswBhV5vJuAABxwuYVV4lr7+q4mUvKLDiJZnV1vUPs42wWe1Ve+556nTJhdrdd5jj50SBCCBgvwBhl/2mrY1o3rvoqayUt6pSnqoqeSor5P38o/m9t7LC+roqDjX9enVInkMHk18311dWyqzyNUeiJCD5/Ur4fZLPfPRLRUWHP/f7lPD5Pz/HfCxSwvq+N/m11PfqrvFZX68/VvIcX3Js83kwqESnToqXlVsfE2XlinfubL2bOdaztxJlTVcUt+bD97MnQNiVPWvuhIDTBAi7nNYR6kEgcwKEXZmzzfTIOQ27Nm7dqccWvqgZ109XMFCc6blmbPz2hF3mP7LK/uMBlf37AzL/8WaOyBdHquIXs1VzxlkZq5WBEUDAHgHCLnscGQUBJwv49u5RyUsvquSlF1SyfJm8+/c1KDfes5fC5385ueXhuHMV797DydOhNgQQsEmAsKt1SO9nn9UFS1ZAVV0tT6MgygRWVvBUP7iqrpK3wgRXyVDKuraiovUb5tEZJoCL9+plvUMy3rVbMgTr3UeJHj0U69lL8V59FO/ePfl5j555NHN3TIWwyx19okoEMiFA2JUJVcZEwJkChF3O7Es6VeUs7DLv7Lp+5v16f92mZus8ZfhQPTTnJnXtUp7OPHJ6TlvDrrL//A9rNVfqH82iw45Vxe2/UOiSKTmdBzdHAIH0BQi70rfiTATyRaDoH++r5JWXkr9ef1UKhRpMrfbEk1VzzgTVTDhfkTPOVCIQzJepMw8EEKgnUOhhl2/nDlm/tn8s3/btyY87tsu7e5f8WzfLu39/Rp6XeLduyVVQZjVU56OUCAYaro7q1Ekq79xk1VRG/rc4GpUnFpXMR7P9berzaMz63BONSrW1decoGjv8eahaXhPy1Q/0Kivl27VT3k8+qVt1li5ickVYT8V79Va8u/m8lxI9eynWKxmGmR/MiJlzBh6d7pCcdwQBwi4eDwQKV4Cwq3B7z8wLT4Cwy709z1nY5V6yppWnG3aVPv6Yyuf8Qr7t26xBzE/rVc68XVVfvzqfOJgLAgUhQNhVEG1mkgg0K5B6Z9dnS15QsVn5teJlFb/ztybn1pw5ThHzvq9x5yoycjSaCCCQJwL5HHZ5Dx6Ub8c2+Xbvkm/rFnl3bJd/21Z5d+6U34Ra2z5Oq4tmK754pzIrnLK27Ps8hEp0Ml9Pbt2XKDffKzu8lZ8VYnVRIhhsGGIddVRa98yXk8xKNt+ne+X99BP59u6V95O9Mlvfe/fstsIw87lvr/l8r0y/2nIYc7NKLNa7r+LWxz6K9++v2uOHKzZ4iKJD2J63NU/CrtaE+D4C+StA2JW/vWVmCDQWIOxy7zORs7DLrOyadffDmnHDdA0b1K+B4Oo167Vg8XLdOeMaV2xv2FrYFVi6WJ1/dUfduz/iXbqo6oc/VuV3v5fcJ54DAQRcJ0DY5bqWUTACtgmkwq49B5LvfDGH99AhlaxcrmKz6mvFy/Jv+KjB/cw/MIZN6GVWfp1zrsyqbg4EEHCngJvDLrPyyqzC8m3bJt/O7fJuS67KsoKsHdutrQNbO8wKq+igIYr36atY/4HWiqFY//6KDTAfByjWt+F/27U2Ht/vmIBZZWeCMbOizoSUVghm/dojn/m6CcxMcLZnd1o3Mn8/RYcMVfT4ExQbMsz6+8r8PjZgYFrX5/tJhF353mHmh0DLAoRdPB0IFI4AYZd7e+3IsMu8y2vug49r9m3Xunobw5I3XlX5nT9V8d/eqntCKr/3I1Xe9G+KF9hPKLr3jwiVI9C8AGEXTwYChSvQXNjVWMO3basCy15Q8Ssvq2TlyzLvr6l/1HzpbOsfiiNnjbPe1Rk99rjCBWXmCLhMwKlhlwk4ksHVNmt1lnfXTvm2bpVv1w75Pt4q8x7C1g7zg3gm1IgNGKBYvwGKHT0oGWCZIKtff0WPIahvzdAN3y9at1a+rZtV9OF6+TZ+JN/mTSra8E8rJDvSUXvyqYqaFWDHHKvY0GMUHXaMokOPsbZKLJSDsKtQOs08EWgqQNjFU4FA4QgQdrm3144MuxYuXaFVb6/N2squ++bN1+//vLRBF++69RpNnTjO+pqp5/Z7HrE+n3Te2CZ1NV7ZZf7jofzOnyjw4vN1Y5qtCiv/7af8pKN7/6xQOQINBAi7eCAQKFyBdMKuxjpF7/5dgZeXqfiVZSpZ+UoTvHj3Hqr612uU6N5dkdNGKXL6GYULzMwRcLhALsIuT004GVyZMMsEWdbqrK3Jd2dt2yb/5o1pqZlt1K0wywRZgwYpbgIsszrLCrQGyvxvEUfhCnhC1Spav06+LZvk37xJvo8+lH/jBvk3bWjyQxv1lRKlnVRrVoINHno4CBs6TLXHnWBtV5lPB2FXPnWTuSDQNgHCrrZ5cTYCbhYg7HJv97IedplVW9fdeq927dnXolrf3t01755bmmxvmClmE3aZ4+brrmpyC7Ol4r3z5uuhOTdZq8yaOzcVdvm3bFLZnLtU+sRf6sYJT5ysQz/9haLHHZ+p8hkXAQRyIEDYlQN0bomAQwTaE3Y1Lr3k9ZUqfv1VFb/xqorfelPmHxgbH7UjTlNk1BhZH0eOZvWXQ/pPGQhkIuyythY078qyVmXtkM9sK2g+N7/MdoMHDrQKb0KFqAmyPg+v4gOTq7KiA0241V+xowe3OgYnINCSgHf/Pvk3bbS26fVt3ij/R/+0QjDzNU91VYtw1raXxxyv2pNPUfTkUxU5dYT195pbD8Iut3aOuhHouABhV8cNGQEBtwgQdrmlU03rzHrYlSrhSO/syjbnkcIu873BA/vUrfJqHH6ZWves26Kye36lTn/4r7rSI2PP1KE7fqnI6DHZng73QwCBLAgQdmUBmVsg4FABO8KuxlMzWx4Xr16lotVvyawCM+/VaXwkOpUpMmp0MgAbPUaRkacr3rWrQ5UoC4H8FWhr2GW2MbWCrN275NuyWd5dO+T/+OPPV2UlV2ulc5gtBaPWe7GS78eKmyBr0GDF+/RTdNAgmdU1HAjkQsA822Z1oRWEmQDsw/UyPwjq/+eHLZZjtkQ0/60cPXWEIiNOU+0pX8hF6W2+J2FXm8m4AIG8ESDsyptWMhEEWhUg7GqVyLEn5CzscpJI420MU1sYhsIR3TH3EY0deWJd2GVWpv1k9sP61axrkyvPfvpTJX59r8zWIuYwL/I9dMevFP7yxU6aIrUggIDNAoRdNoMyHAIuEshE2NV4+uYfx4v/9qaKTAj2tgnC3pSnsrKJklmpEfnCF1V7+hgrBOOHbFz0IFGqawUah13mH/l9O3fKt80EWNvl/Xir/HUrtbYfcdVLCsFsH5gMrvpa2wrGzWqs1LuyzMc+fV3rReGFLeD7eIuKNnyU/EGONe+o+P13rXfINT4SgaC1itn8UEft53+fxXv0dBweYZfjWkJBCGRNgLAra9TcCIGcCxB25bwF7S6AsKsRXWqbxdmzrtXJJwy1wq5pk8dr9IgTrDObhF0ej/X1eP/+ivz8LtV+7V/a3QwuRAAB9wj4fR4V+bwKRWLuKZpKEUDAFgGvRyot8asyHLVlvHQH8a5bK9/f35H3nbfle+N16/PmjtjIUYqNGau4+TjqdMWPPTbdW3AeAgg0EvDs2S3vtm3ybt8uz7aP5dm5U8Xbtii2fYe8W7bIs3dP62bBoOLmfVgDBioxcKAS5l1Z5vOjj1a83wD+jLYuyBl5JuA5dFA+83fZm6vkXf2WfKtXy/PJ3iazTPTrJ+vvtLFnWH+fxUaNloLBnGqUlvhUUxtXLJ7IaR3cHAEEsi8QKPLK/NGPROPZvzl3RACBrAqUB/1ZvR83s08gp2GX2crw+pn36/11m5rM6JThQ+vek2XfdNMbKbV14cUTxra+suuyy1Rz5jhFvveD9AbnLAQQyAsBv8+rIr9HoRrCrrxoKJNAoA0CHo9HpQGfqkLZDbuaK9H8Q6HvrTflXfWGfGvesf7xvfGRKCtXbMwYxU4fo9jESxQ7bWQbZsupCOSxQCgk79Yt8u7YYa3G0o7t8m7dKs/2bcmAa3PT/0ZpTsP8g3z86MGKm4+DB0v9Byg+eIgSffpYK7QS3XvkMSJTQ8AeAfNnz/p7zPx99tZb8q16vdmBY18YodjoMYp/8YvJH+o45VR7CkhzFPPDLjXRmGIxwq40yTgNgbwRKCn2KR5PqJawK296ykQQaEmgvLQIHJcK5DTsOtK7snLpWf89Xem8s2vnvlAuy+XeCCCQAwG2McwBOrdEwCEC2djGsL1Trdv+cPWb1jaIxW+vbn77wz59rfej1I44TZFx41VzxlntvSXXIeBYAd/2bda2gr5t26xtBb3bP5b1NfNrxzZ5Dx5stfZ4586K9R+o2IDkr/iAo1V+/BDt69pH0X4DFBt4dKtjcAICCLRPoPjvb1t/l5ntD4veW6OidWubDJQIlioycpS1BWLt6LHWdr5mW9BMHWxjmClZxkXA+QJsY+j8HlEhAnYJsI2hXZLZHydnYZdZ1TXr7oc144bpyXdf5egwdSxdtkpfm3qBVUHjbQpXr1mve+fNr1tl1lxAR9iVo+ZxWwRyKEDYlUN8bo1AjgWcHHY1R1O0fp2KzD8Uml/m/V9/f7tFwejQYxQdMlTR44crZj4ec5z1e/OP/BwIOEnAu39/MsjavUu+LZvl3bVD/o8/ls+8K8uEWjt3pFWuee9d9OijFevXX7GBgxTv11/Rowcp3re/ooMHy7xHqPHR+J1dad2IkxBAoMMC5t2VxeZdlql3Wv7tTZkf8mh8mD/P1g9zWOHX6YqMPbPD904NQNhlGyUDIeA6AcIu17WMghFok4CnJqxO/zVPJS88p5KVr7TpWk52jkDBh12hcMTaqnDJslV1XfnDAzPr3tFlvrhw6Qrdfs8j1vcnnTdWd864RsFAcd35hF3OeaCpBIFsCRB2ZUua+yDgPAG3hV3NCRat/cD6KXn/u2tU9MF7KtrwT3mbeV9K6tpESUDRoUNlhWHHHqfYkGGKHnOs9THWq7fzmkRFrhYwAZZ3z2759u6xwiyv+fjxVms1VjLM2i5PqLrVOcZ79FR00GDF+vS1VmeZLQVj/QccXqnVu0+rYzR3AmFXu9i4CIGMCPi3blbR26tV9Pd3kkHY6sP/XV//hibwqjnzbEW+dLYiY8Y2G2SnUyBhVzpKnINAfgoQduVnX5kVAkag9PE/qvyun8m3Z3cSJMF2xW59MnIWdhmwxlsEuhWRsMutnaNuBNovQNjVfjuuRMDtAvkQdjXXA091lYo+XC/flk3yb94k30cfyr9ls4rWr5WnoqLFtiVKO6n2+BMUG2xWgh2rmAnEhg5T7XEnKFFe7vZ2U79NAuY/HE2g6t23T75dO5Of790j7yd75Pv8c99e8/1P07qj2bosNmCAYgM+X5VlQi2zOsv6fT8rmM3UQdiVKVnGRcAegeJ3/pbc/tBsg/j23+TftKHJwJExX1LN+AmqOedcRU4/I+0bE3alTcWJCOSdAGFX3rWUCSGgwHPPqvNdt8v/zw8tDbPV/6Gf363uUyei41KBnIZdZsvAxxa+qBnXT2+wUsptloRdbusY9SLQcQHCro4bMgICbhXI17DrSP3w7t8n/6aN8m/cIN+mDfJv+Eh+E4p9uF6ecMvvLo1366boMccrOmSIai64SPGevRUvL1f0+BNkVotxuFvAWnWVCrDMdoJ798hjVmF9skfeTz9Jrsjas0fm+WnLEe/WXbHevRU3K7J69bZWZpltBWODhyRXafUbIPNs5eog7MqVPPdFoH0C5n+bil9boeLXX1XJ6yutv7vqH+aHNmrOGqfI+AkKn3Oe9XdUSwdhV/t6wFUI5IMAYVc+dJE5IJAUKH7rDXW+4zYVr37T+n108FBV/OQOhS6fZv2ed3a590nJWdhl3pV1/cz79f66Tc3qnTJ8aN17spzOS9jl9A5RHwL2CxB22W/KiAi4RaAQw64j9cYEHv7NySDMu2mDtTosFYQd6bpY335SPC5PPK5Yj56K9+ols+1crGdvJXr2UqxXr2TI0bVb8uu8Nywrf0TqVl2Zd2Lt3imz2spjthQ0odanydVY1gqsNgdY3RS3+tpb8R7JjwkTZFl97meFV3Hze/NcOPwg7HJ4gygPgVYETPhV8tKLKnnpBZUsX9bkf89ivfuoZsL5qplwgWrGnat49x51IxJ28XghULgChF2F23tmnj8C5gdeOv/yZ9aKLnOY/86smHm7qq7+doNJEna5t+c5C7vcS9a0csKufOomc0EgPQHCrvScOAuBfBQg7Eq/q75tW+XfvFn+jz6Ub/vH8prVPyY4MVvamV+ffZb+YGbr9PJyxXr0Urx3H8W6d7fCk8O/eiaDlO49k1/r0qVNYxfCyb6Pt1jvtDH+nkOHrCmb8Kpk1WvJFVr797eJwVqBZULKIwVYJuByQYDVlokTdrVFi3MRcL5A0T/eV8krL6l4+Usqef3VJiuWzZZGZrvDmnMmqPzLE3RIfkVq486fGBUigICtAoRdtnIyGAJZFTDv/S2f8wuV/vmP1n0TncpU+aNbVPnd78us8G58EHZltT223oywywZOwi4bEBkCAZcJEHa5rGGUi4CNAoRdNmJKMv/hYW15Z1YS7THb4iW3vzOriHzWVnh7rSDG/CR+Ww/zDifz03rhCy6SvN4Gl9eOOl2JQLBtQ3o9SviLlCgqkoqKrM9V5E/+PvV166Pf+g+olg5PTVie6pDMe9I8odTHautr3sZfC4Wk6mp5qqo+P7/auiZ5Xlieqkp5QsmvWeOZcWvCbZtXvbPNCobDAVZypZ0JF61VWWal3echo/m8UA/CrkLtPPMuFIGSV1eoePkylax4Web9X42P2rPGKXyOed/XBEVOG1UoLMwTgYIXIOwq+EcAABcKeA8cUNmvZ6ts3r/XVV/5ne+p8sczZX5wr6WDsMuFzf685JyGXaFwRHfMfURLlq1S397dNe+eW9Svdw/ra2NHnqipE8e5QpawyxVtokgEbBUg7LKVk8EQcJUAYVfu2mW2zjNhmBWCfWLeDbVXnk8/TW6t9/nKJOs9Up9+YgU/hXyYIC9RWpr8FTSfd1K8U5m1HZcVXpl3X/X+/H1Yn6+QM1t3cbQuQNjVuhFnIJAvAt5Dh1Sycnky/HrlZfk3bWgwtXjnzqo5e7wiVvh1rqLDjs2XqTMPBBBoJEDYxSOBgHsEzHuly373H+r0m1/L/F1ujupp01XxkzvT2h6fsMs9vW5caU7DrvvmzdfggX108YSxmvvQ4/ra1PM1bFA/rV6zXgsWL9edM65RMFDseF3CLse3iAIRsF2AsMt2UgZEwDUChF2uaZXMtn3W+6UOHUyugKqstD56KyulmhrJrLCqjcpTWytFa+t9jEq1tfJEa+t9bO488/1og/OssczXwqEjQpltFhNBE0R9HkjVBVP1vhYsVbxTqVTaqS6wanx+vPEYnwdb7umS+yol7HJfz6gYAbsEeu7fqciS5+R96SWVrHy5yXa80UFDrPd9Rc4+RzXnnMeWunbBMw4CDhAg7HJAEygBgTQEOv33w9ZqLrN1uznMLh8VP7tLtcNPSuPq5CmEXWlTOe7EnIVdBw5WaNbdD2vGDdOt1Vz1w66NW3dq7oOPa/Zt16prl3LHoTUuiLDL8S2iQARsFyDssp2UARFwjQBhl2taRaEIZESAsCsjrAyKgCsEenQu0aFQbfKdXfG4it5bo8Dyl1T8yjIVv7lKnkjN4Xl4vao9dYTC4ycoMv58RU4fq0Sx83+Y1xWNoEgEciBA2JUDdG6JQBsEgoufUvkvfy7/xo+sq2pHnKaDd89V5PQz2jBK8lTCrjaTOeYCR4ZdrOxyzPNBIQgg0IIAYRePBgKFK0DYVbi9Z+YIGAHCLp4DBApXoEHY1YjBrOgtfuM1lbzykkpeelFFaz9ocEaiJKCas8apZvJlqp5yhRLlzv/B3sLtNDNHoKkAYRdPBQLOFDDv2uz889tUtOYdq8Do0GOslVyhS6a0u2DCrnbT5fzCnIVdZuYLl67QqrfXatYPv6bfPvKUtY1ht6PKdf3M+3XV5PG8syvnjwcFIIBASwKEXTwbCBSuAGFX4faemSNA2MUzgEBhCxwp7Gos492/XyWvLLOCr5KXXqjbTil1XvjCiQpNvUrhiydZ71XkQAABZwsQdjm7P1RXeAJFH7ynznf+VCUvv2hNPtanrypn3q6qr1/dYQzCrg4T5myAnIZdZtZmFdfVN85pAPCHB2Zq9IgTcobS1huzjWFbxTgfAfcLEHa5v4fMAIH2ChB2tVeO6xDIDwFWduVHH5kFAu0RaEvY1Xh8/4frFVj2vEpe+KtKVr7S4Nsm+ApfdqXMx3jnzu0pjWsQQCDDAoRdGQZmeATSFDDvZe78q58r+OR86wrzLuTKm2eq6tvXyayituMg7LJDMTdj5Dzsys207b0rYZe9noyGgBsECLvc0CVqRCAzAoRdmXFlVATcIkDY5ZZOUScC9gt0JOyqX42nqlIBs+Lr/5Yq8Nyz8n72Wd23a8ZPUHjy5QpNulTxHj3tnwQjIoBAuwQIu9rFxkUI2Cbg3b9P5ffcrU7/9VDdmJXfv1GVN96q+FFH2XYfMxBhl62cWR0sp2HXffPma/fe/bpzxjUKBpIvag2FI7pj7iMaO/JEtjHM6qPAzRBAoC0ChF1t0eJcBPJLgLArv/rJbBBoqwBhV1vFOB+B/BGwK+xqIBKLqeS1lQo8s1DBJc/I+8ne5Le9XkVGjVF48hSFJl+u2ICB+QPJTBBwoQBhlwubRsl5IeCprlLZg79R2b/fL09lpeTzqXr611Ux62fW1oWZOAi7MqGanTFzFnalQq1pk8c32bLQbG24YPHyBiFYdjjadxdWdrXPjasQcLMAYZebu0ftCHRMgLCrY35cjYDbBQi73N5B6keg/QIZCbsalWOCr5Jnn1bpoicOB1+Sak8+VeEpUxWaeKmix7vntQ/t1+ZKBJwlQNjlrH5QTWEIlD38oMrmzpZZ1WWO8MTJOnTHLxUddmxGAQi7Msqb0cFzFnYdOFihWXc/rBk3TNewQf0aTHLj1p2a++Djmn3bterapTyjAHYMTthlhyJjIOAuAcIud/WLahH4/9u7FzCnqnvv4/9MZibJXETEK3dBrWipCFI5WikFryAiw0Vs9Wi1lmrtqZcHC55aq54KB18v7+vbeiivaFuryACDIlRFFO8oF7FYoFXuyk0FcS5JZibJ++wdEjNhYJLMSvZayTfPwwNM9l77vz7/nTDkN3ttlQKEXSo1GQsB8wQIu8zrGRUjoEogF2FXYq2l770j3heeF9+CueLeuSP+VPOJvcU/qkoCl46Spn79VU2PcRBA4DAChF2cHgjkTqBs7rNSMfVeKd662T5o8N++J7W//Z00DhiYkyIIu3LCnJWDOBZ2cWVXVvrJoAggkCMBwq4cQXMYBDQUIOzSsCmUhEAOBQi7cojNoRDQTCDXYVeL4Gv1yuhSh88vEPe2LfGnQt2621d7WcsdNp59jojLpZka5SCQHwKEXfnRR2ahsUAkIr4XnpOK/75PSjastwttOu3bUvub+yRw/kU5LZywK6fcSg/mWNhlzcJarnDK1JkyY/rt8au7rKu6Jt7xoNx0zSju2aW01QyGAAIqBQi7VGoyFgJmCRB2mdUvqkVAtQBhl2pRxkPAHAEnw65EpZK/r4kGXzXz4j/1bj0f7nS0+EdcJsFLR0lg6AXmwFIpAgYIEHYZ0CRKNFbAt7DGXq6wZN1H9hxC3XtI7Z2/lYaxVzgyJ8IuR9iVHNTRsMuaQSzcSJgRwgAAIABJREFU2rk7uvam9XjykckH3cdLyWyzNAjLGGYJlmER0FiAsEvj5lAaAlkWIOzKMjDDI6C5AGGX5g2iPASyKKBL2JU4ReuDQe9z88X3wgIp/ueG+FPhDh0kcOElErz0cgmcf6FEPN4syjA0AvkvQNiV/z1mhrkX8L7yklROvVdKPvzAPnioazepnXSnNFzxI5Hi4twXdOCIhF2O0bf7wI6HXe2egQYDEHZp0ARKQCDHAoRdOQbncAhoJEDYpVEzKAUBBwQIuxxA55AIaCKgY9iVSFO86RPx1cy1w6/YT8dbz0d8ZXbgFbCCrwsvkUil/vdG16TllIFAXICwi5MBAXUCnteWSuW0e6V01Qp70NBxx0vdbb+Shqt/LJHSUnUHynAkwq4M4TTYjbBLQRMIuxQgMgQChgkQdhnWMMpFQKEAYZdCTIZCwEABwi4Dm0bJCCgS0D3sahF8bdkkvgXz7eUOrWUPYw/rQ8TgeUMkOPJy8Q8fKeGjOinSYRgE8luAsCu/+8vsciPgefN1+0qu0vfftQ9oLb9b98vbpf66n0rE68tNESkchbArBSRNN3E07Nq3v1ZunPywrF2/6SCevn16yWPTbpWOHfT/iSPCLk3PbspCIIsChF1ZxGVoBDQXIOzSvEGUh0CWBQi7sgzM8AhoLGBS2JXI6P50u/gWzBPvwhopXb1SJBKJPu12S/DscyQ4cpT4R4ySUOcuGutTGgLOChB2OevP0c0WKH3vHTni3t+I9bv1sJbarb/5Vqn72c321ce6PQi7dOtI6vU4GnY9NGOOXeltE8enXrGGWxJ2adgUSkIgywKEXVkGZngENBYg7NK4OZSGQA4ECLtygMwhENBUwNSwK5HTvWun+Gqqxfv8AildsbyFdFO//vZShw1jr7Dvm8IDAQS+ESDs4mxAIH2B0pXvS+V//5d4XnvF3jlSVi51N/5C6m6+VesldQm70u+1Lns4FnZZV3VNuX+mTLppgvTu0VkXj4zqIOzKiI2dEDBagLDL6PZRPALtEiDsahcfOyNgvABhl/EtZAIIZCyQD2FXi+Brz27xPl9j//Isf1skHI4/3XjmAAlUjRP/mCskdOxxGZuxIwL5IkDYlS+dZB65EChZ+6FU/tfd4l36cjTk8vqk/vqfSt0td0i4Y8dclNCuYxB2tYvP0Z0JuxTwE3YpQGQIBAwTIOwyrGGUi4BCAcIuhZgMhYCBAoRdBjaNkhFQJJBvYVciS9HeL8W3cIF4a6rF89YbLcQazz5H/KPHiv/yMRI++hhFmgyDgFkChF1m9YtqnREoWf+PaMj10uJ4AfU/uVHqbrvDqB+cIOxy5vxRcVTHwi6reGsZw57djpeq4YNVzMWxMQi7HKPnwAg4JkDY5Rg9B0bAcQHCLsdbQAEIOCpA2OUoPwdHwFGBfA67Dgq+5s8R37xqKV35Xst7fJ17ngSqxov/0sslfOSRjvaDgyOQSwHCrlxqcyzTBIr/uUEqp94jvkXPR//NKC6W+glXSd2vfi2hE8xb0Y2wy7Qz8Jt6HQ27Nm7dIX+d/4pMunGC+LylxioSdhnbOgpHIGMBwq6M6dgRAeMFCLuMbyETQKBdAoRd7eJjZwSMFiiUsCuxSe5Pt4tvwTz7Pl8lH37Qon+BCy6OBl/DL5VIeYXRvaV4BNoSIOxqS4jnC1Gg+JOPpXLqveJ7bl58+g1X/FDqJt0pzT17GUtC2GVs68SxsMu6Z9eNkx+Wtes3tarXt08veWzardKxQ6X2uoRd2reIAhFQLkDYpZyUAREwRoCwy5hWUSgCWREg7MoKK4MiYIRAIYZdycFX2dzZ4p1fLSXrPoo/Zd2LJRp8jZPAhRdLxOM1op8UiUA6AoRd6Wixbb4LFG/eKJXT7hNfzdzo/R5dLvGPHC21d94tzSedbPz0CbvMbaFjYZe5ZAdXTtiVT91kLgikJkDYlZoTWyGQjwKEXfnYVeaEQOoChF2pW7ElAvkmUOhhV2I/izd9Ir7q2fYVX9ZP9scekYoK8V9yqX3FV3DIMImUlOTbacB8ClSAsKtAG8+0Wwi4t22Ryum/k7Lq2SKhkP1c4KLhUvvre6Spz+l5o0XYZW4rCbsU9I6wSwEiQyBgmABhl2ENo1wEFAoQdinEZCgEDBQg7DKwaZSMgCIBwq7WIUv+sVa8NdXiq5knxVs3xzcKd+xo39vLDr7OPU+kqEhRJxgGgdwLEHbl3pwj6iNgLWlb+cD9UvbsX0Wam+3CgkOGyte/vlea+vXXp1BFlRB2KYJ0YBjHw64VazbItbdMazH1Jx+ZLAP7neoAR2aHJOzKzI29EDBZgLDL5O5ROwLtEyDsap8feyNgugBhl+kdpH4EMhcg7GrbrvSDVfYyh9b9W9w7PovvEDr2OAlcViX+MeOk8ayz7SWveCBgkgBhl0ndolZVAu6dO6TiwWlS9tc/iaupyR628azvytf33C+NZ5+j6jDajUPYpV1LUi7I0bDLCroenDGnxb25Nm7dIRPveFBuumaUVA0fnPJEnNyQsMtJfY6NgDMChF3OuHNUBHQQIOzSoQvUgIBzAoRdztlzZAScFiDsSq8Dpe+9Y9/Pxft8jbj37I7vHOrWXfyXjxF/1Xhp6ntGeoOyNQIOCRB2OQTPYR0RcO/eJRUPT5eyPz8hrsagXUPjmQOk9s7fSvAHwxypKZcHJezKpbbaYzkWdvkDjXL3A7Nk3MghB13FZYVg1QuXyT2TrhOft1TtjLMwGmFXFlAZEgHNBQi7NG8Q5SGQRQHCriziMjQCBggQdhnQJEpEIEsChF0ZwobD4nn7zehShy8skKK9e+MDNfc+Wfyjx9rBV/Mp38rwAOyGQPYFCLuyb8wRnBco+uJzqXjkf0n5EzPFFQzYBTWd3ldqp/xGAhePcL7AHFVA2JUj6CwcxrGwa9/+Wply/0yZdNME6d2jc4upWVd3PfCH2TL1zhukY4fKLExb7ZCEXWo9GQ0BEwQIu0zoEjUikB0Bwq7suDIqAqYIEHaZ0inqREC9AGGXGlPv0pfFW/2M+F5cJK66uvig1geqgdFjpWHMeAl166HmYIyCgCIBwi5FkAyjrcARv5ki5Y/PiIdczSedLLV33i3+y6q0rTlbhRF2ZUs2++M6FnZxZVf2m8sREEAgewKEXdmzZWQEdBcg7NK9Q9SHQHYFCLuy68voCOgsQNilvjveFxeJr6ZavIsWiivgjx+gsf9ZdvDlHz1OQsefoP7AjIhAmgKEXWmCsbkRAvaVXI8+JOWzZorL32DXHOreU2p/9WtpuOKHRswhG0USdmVDNTdjOhZ2WdObv/gNmbNwGffsyk2vOQoCCCgUIOxSiMlQCBgmQNhlWMMoFwHFAoRdikEZDgGDBAi7stcsa7ks75KXold8LXq+xYEaB50r/vFXin/kaAl37Ji9IhgZgcMIEHZxeuSTQOyeXOX/73/i0wp17SZ1t0+W+qt/nE9TzWguhF0ZsWmxk6NhlyVg3Z/r2lumtcB48pHJB93HSwutQxTBMoY6d4faEMiOAGFXdlwZFQETBAi7TOgSNSKQPQHCruzZMjICugsQduWmQ66GevEtfkG88+eI9+W/tThocMhQ8Y+5wg6+IhUVuSmIoyAgIoRdnAb5IFC8dbP4nvqTVD48PT6dxoGDpP6mX9jvqzyiAoRd5p4Jjodd5tJ9UzlhVz50kTkgkJ4AYVd6XmyNQD4JEHblUzeZCwLpCxB2pW/GHgjkiwBhV+476aqtFd8LC8Q3f454XlvaooDARcPFXzVeAiNGSsTry31xHLGgBAi7CqrdeTdZzztvSvn//F/xLl4Yn5u1TGH9xJul6Tv98m6+7Z0QYVd7BZ3b39Gw66EZc2TXnr1yz6TrxOcttRVi9/IaNOA0qRo+2DmZNI5M2JUGFpsikCcChF150kimgUAGAoRdGaCxCwJ5JEDYlUfNZCoIpClA2JUmmOLNi776SnzPzRPv/GrxvPuWSDhsHyFSVi528DVmvASHXiCR0ujnSzwQUClA2KVSk7FyIhAKie/5+VLx+/8tJWtW24cMd+gg9ddcLw0//Tn3QzxMEwi7cnKGZuUgjoVdsVBr3MghBy1ZaC1tWL1wWYsQLCuzVzQoYZciSIZBwCABwi6DmkWpCCgWIOxSDMpwCBgmQNhlWMMoFwGFAoRdCjHbOVTRF59LWU21eOfPldKV74lEIvEPcgPDL7Ov+AoOHiLidrfzSOyOQFSAsIszwRQBV32dlP95lpTP+L24P91ul93c40Spv/FmafjRNRLxlZkyFcfqJOxyjL7dB3Ys7Nq3v1am3D9TJt00QXr36NxiIhu37pAH/jBbpt55g3TsUNnuSWZ7AMKubAszPgL6CRB26dcTKkIgVwKEXbmS5jgI6ClA2KVnX6gKgVwIEHblQjn9Y7h37RTfvDniq6mOX71gjRLudLT4LxstgapxEhx0rojLlf7g7IHAAQHCLk4F3QXcOz6zlyos/8sssZaAtR6NA8+W+pt+Kf4Rl4kUFek+BW3qI+zSphVpF+JY2MWVXWn3ih0QQEAjAcIujZpBKQjkWICwK8fgHA4BzQQIuzRrCOUgkEMBwq4cYmd4KOsqhrK5s+2lDkvWfRQfJXRCZ/FfPtYOvhrPHJDh6OxWyAKEXYXcfb3nXrL2Q6l49GF7yUJpbravaPUPv0zqfnm7NPXrr3fxmlZH2KVpY1Ioy7Gwy6rNWq5wytSZMmP67fGru6yruibe8aDcdM0o7tmVQgPZBAEEnBEg7HLGnaMioIMAYZcOXaAGBJwTIOxyzp4jI+C0AGGX0x1I7/jFmz4RX/Vs+4qv4k8+ju9sLeflHz1WAmPGS1Of09MblK0LVoCwq2Bbr+fEIxHxvrRYyv/wf8Tzzpt2jZHyCqn/4b9L/U3/IaFu3fWs25CqCLsMaVQrZToadln1xMKtnbu/jJf35COTD7qPl87ELGOoc3eoDYHsCBB2ZceVUREwQYCwy4QuUSMC2RMg7MqeLSMjoLsAYZfuHTp0fSX/WCve+XPEt2C+FG/dHN+w+Vun2vf38leNk+YTe5s7QSrPugBhV9aJOUAKAq6AX8pm/1XKH3tUijdGQ3zrytX6n/5cGq65XsJHHJHCKGzSlgBhV1tC+j7veNilE01saUWrpnsmXSc+b6ld3vzFb8hd02fZfx4xbFCL56yvEXbp1EVqQSA3AoRduXHmKAjoKEDYpWNXqAmB3AkQduXOmiMhoJsAYZduHcmsntIPVtnLHPqemyfWPW5ij6YzzowGX5ePkVCXrpkNzl55K0DYlbetNWJiRV9+IeV//IOUPzFTivZGLxhp6nuGfT+uhtFjRYqLjZiHKUUSdpnSqYPrJOw6YBILuhYtXd4i0LKWWnxwxhx5bNqt0rFDpTw0Y469x20Tx8c1CbvMfQFQOQKZChB2ZSrHfgiYL0DYZX4PmQEC7REg7GqPHvsiYLYAYZfZ/Wut+tL33hFfzVzxPl8j7j27o5u4XNI4cJB9tZcVfIWPPib/Js6M0hYg7EqbjB0UCFhLsFY8+pC9JKurMWi/PwXOv0jqf36LBL83WMERGKI1AcIuc88Lwq4DvbNCrJ7djrf/tnzVuvjVW7GvVw2PvoEkh1/W1wi7zH0BUDkCmQoQdmUqx34ImC9A2GV+D5kBAu0RIOxqjx77ImC2AGGX2f07bPXhsHjeflO8NdXie2GBFO3dG93c7bY/ULau+AqMGCXhI4/MYwSmdjgBwi7Oj1wKeN5YJuW/f0S8ry4RiUQk4vFKw/grpf4/bmPJ1Rw0grArB8hZOgRhl0iLq7WsJQtjYZdlfvcDs2TQgNMkFnZZ9xj7z6kz5XdTbpDePTrbbSHsytLZybAIaCxA2KVxcygNgSwLEHZlGZjhEdBcgLBL8wZRHgJZFCDsyiKuZkN7l74s3upnxPfiInHV1cWrC1xwsfjHTZDARcMlUl6hWdWUk00Bwq5s6jJ2TMC75EWpvO83UrLuI/tL4U5HS/11P5X6n/zM/jOP3AgQduXGORtHKfiwywq3tmzfFV+WsLWwa9zIITKw36m2f2thVzgSyUZvGBMBBDQWcGlcG6XlQoAzIBfKOh/D5bJ/wI4HAggUoACv/wJsOlNG4ICA/fq3/sz3AAV0TkTEtXChRObNE1d1tYjfH597ZOxYcU2YIJHhw0W83gIyYaoIFIZAzt7q160T1x//KK6//EVk374o7imnSGTSJIlcf31hYGs2yyLrH3weRgoUfNhlLVP4+DOLD2reiGGDZPIvrpJpjz7V5pVdO7/85psdI88CikYAgbQFPKVuKfO4ZW9tY9r7soP5AnzfY34P2zMD68quTkd4ZM9XgfYMw74IIGCowPEdfbJ7n5/Pug3tH2Uj0B6BTpUeqfU3SWNzuD3DsK9BAok/3OQKBsS38DnxLKwR36Ln47OIlJWLf+QoCVhLHQ670KDZUWo6AlzZlY5Wfmyb7bjDt2CulD32qJSuWhEHs5ZNbZh4swQuuTQ/EA2dxQmdfIZWTtkFH3YlnwKJV3b5vKX2EofWvby4ZxcvFgQQSBRgGUPOBwQKV4BlDAu398wcAUuAZQw5DxAoXAGWMSzc3ifP3NVQL77FL4h37mzxvvJS/Olwhw7iv6zKDr6C530fsDwSIOzKo2Y6OBXPm6+L79mnxPdcjbj8DfFKGn50jdT/7BfS1Oc0B6vj0DEBljE091wg7ErqXXLYtWLNBnlwxhx5bNqt0rFDZYv7e8V25Z5d5r4AqByBTAUIuzKVYz8EzBcg7DK/h8wAgfYIEHa1R499ETBbgLDL7P5lq3pXba34Fi8UX/XT4ln2avwwoWOPk8BlVeIfM04aBw7K1uEZN0cChF05gs7DwxRv3ii+p/8iZdXPiPvT7d+8R3TuIg3X/kTqr/2JhI/qlIczN3dKhF3m9o6wq42wy3raCsDumj7L3tJa3vCeSdeJddVX7EHYZe4LgMoRyFSAsCtTOfZDwHwBwi7ze8gMEGiPAGFXe/TYFwGzBQi7zO5fLqov+uor8T03T7zzq8Xz7lsi4eiSl6Fu3cV/+RjxV42Xpr5n5KIUjqFYgLBLMWieD+eqq5Oyec+Kb/ZTUrrivfhsI74y8V96mfgnXB29+rOoKM8lzJweYZeZfbOqJuxS0DvCLgWIDIGAYQKEXYY1jHIRUChA2KUQk6EQMFCAsMvAplEyAooECLsUQRbIMEVffC5lNdXinT9XSle+J3LgBmDNvU8W/+ix4h83Qaw/8zBDgLDLjD45WmUoJN7XXhHfM0+J98VFYt3nz364XBI85zzxX3mV+C8bLdZ9/njoLUDYpXd/DlcdYZeC3hF2KUBkCAQMEyDsMqxhlIuAQgHCLoWYDIWAgQKEXQY2jZIRUCRA2KUIsgCHce/aKb55c8RXUy0la1bHBZpO72vf36uhapx99RcPfQUIu/TtjdOVFf/rn1L21JPimztb3Ht2x8tpPrG3+Cf8SBomXCWhLl2dLpPjpyFA2JUGlmabEnYpaAhhlwJEhkDAMAHCLsMaRrkIKBQg7FKIyVAIGChA2GVg0ygZAUUChF2KIAt8GOuePWVzZ9tLHZas+yiu0ThgoASqxol/9Dix7vfFQy8Bwi69+uF0NUX79klZ9dP2VVwlaz+MlxPu0EH8o8bYV3Fxrz6nu5T58Qm7Mrdzek/CLgUdIOxSgMgQCBgmQNhlWMMoFwGFAoRdCjEZCgEDBQi7DGwaJSOgSICwSxEkw8QFij/5WHzVz4hvwTwp3vix/fXwkUdK0xlnSuO5gyVwwcXc40uT84WwS5NGOFyGb2GN+J7+i3iXvNiiksDQC8T/w6vFf/lYhyvk8CoECLtUKDozBmGXAnfCLgWIDIGAYQKEXYY1jHIRUChA2KUQk6EQMFCAsMvAplEyAooECLsUQTJMqwIlH/1dvAvmiW/hgnjwZYdfRx8jDWMnSHD4pdL07TMkfMQRCDogQNjlALomhyxdtUJ8c562lyks2r8/XlVTn9PEP+FqabjyKgkf1UmTailDhQBhlwpFZ8Yg7FLgTtilAJEhEDBMgLDLsIZRLgIKBQi7FGIyFAIGChB2Gdg0SkZAkQBhlyJIhmlToHjLJvG+tFg8S14Uz7JXW2zf2P8sCX5/qDR+f6gEvze4zbHYQI0AYZcaR1NGse69Vfb0n8U3+ymxrsCMPawrLxsmXC3+8VdK03f6mTId6kxTgLArTTCNNifsUtAMwi4FiAyBgGEChF2GNYxyEVAoQNilEJOhEDBQgLDLwKZRMgKKBAi7FEEyTFoCLn+DeN5YJp7Xlopn2SstPniPeLwSPPc8aRwy1A7Amk7vm9bYbJy6AGFX6lambumqrRXfouek9M1lUvbs0y2m4R852r4PV+DCS0ydHnWnIUDYlQaWZpsSdiloCGGXAkSGQMAwAcIuwxpGuQgoFCDsUojJUAgYKEDYZWDTKBkBRQKEXYogGaZdAu5tW8S75CXxvP6qeN54TVx1dfHxrKXUgoN/IMGh50vw/IskdOxx7ToWO38jQNiVn2eDq6FefM/XiNf69fLfWkzSuorS/8N/F3/VeJYPzc/2H3JWhF3mNpywS0HvCLsUIDIEAoYJEHYZ1jDKRUChAGGXQkyGQsBAAcIuA5tGyQgoEiDsUgTJMEoFSpe/LZ5XXxHPsqVSunpli7Gbv3WqBAYPlcZhF0jwnO9JpKxc6bELaTDCrvzptnv7tugyoa+9IiXr/yHubVvjk2s+6WTxj50g/qpx0tzrpPyZNDNJS4CwKy0urTYm7FLQDsIuBYgMgYBhAoRdhjWMchFQKEDYpRCToRAwUICwy8CmUTICigQIuxRBMkzWBIr27xfPa0vEs3SJHX65d+6IHytSUiJN/QdK8AfDovf86n+WiNudtVrybWDCLoM7GgqJ5/13xfPy3+z74JVsWN9iMlao5b98jARGj5WmPqcbPFFKVyVA2KVKMvfjEHYpMCfsUoDIEAgYJkDYZVjDKBcBhQKEXQoxGQoBAwUIuwxsGiUjoEiAsEsRJMPkTMD6UN8Ov159RUrffVtcAX/82OEjjpDGcwdLcMgwCQ4ZKs29T85ZXSYeiLDLrK4V7dsn3iUvRgOu15aIFQTHHta97hrP+Z4Ez79QAj+4QJpP+ZZZk6ParAsQdmWdOGsHIOxSQEvYpQCRIRAwTICwy7CGUS4CCgUIuxRiMhQCBgoQdhnYNEpGQJEAYZciSIZxRMDVGJTSt98Sz6vWVV/W8m3rWtQR6toter8v+8qvYRI+6ihH6tT1oIRdunbmm7pK1n4YD7jsJT3D4fiTzSefIoGhF9oBlxV0WYEXDwQOJUDYZe65QdiloHeEXQoQGQIBwwQIuwxrGOUioFCAsEshJkMhYKAAYZeBTaNkBBQJEHYpgmQYLQTce3aL57Wl4nlxkXheXypFX3/doq6mvmdIcPAQ8U+4iqXdRISwS4vTtkURroZ68Sx7Vbx/WyjeJS9J0Refx5+37k9nXbEYHHahBIZeIKFu3fWbABVpK0DYpW1r2iyMsKtNorY3IOxq24gtEMg3AcKufOso80EgdQHCrtSt2BKBfBQg7MrHrjInBFITIOxKzYmtzBQoXbXCvs+XveThe+8cNInGgWdL8Nzz7KUPG787SCLlFWZONMOqCbsyhFO8W/HmjeJ9abF4XnnJDroSH41nDojej27o+RI85zzFR2a4QhIg7DK324RdCnpH2KUAkSEQMEyAsMuwhlEuAgoFCLsUYjIUAgYKEHYZ2DRKRkCRAGGXIkiG0V7AVV8nnjdfjy55+NbrUvyvf7aoOXz0MdJ0yqkS6tFDms44U5r6D5TG/mdpP6/2FEjY1R69zPct/vhf4nn3LSl5/11xb9smnnfejA8WPvJICQ45X4IXXiKB8y9i6c3MmdkzSYCwy9xTgrBLQe8IuxQgMgQChgkQdhnWMMpFQKEAYZdCTIZCwEABwi4Dm0bJCCgSIOxSBMkwxgkUffmFeN55S0rfflNK33lTStZ91Oocmvr1t0Ov5u/0k0YrBOt7hnFzPVTBhF3Zb6XL3yCla1ZLyepV4tr/lZT/6XGxzr3ER+N3/02CF1wUvYIrzwPW7ItzhEMJEHaZe24QdinoHWGXAkSGQMAwAcIuwxpGuQgoFCDsUojJUAgYKEDYZWDTKBkBRQKEXYogGcZ4AVddnZR+sEpKV6+U4jWr7IDCvX1bq/NqHDjomwDszAHSfMq3jJw/YZf6tlnnUIn9a6WUfLim1RA1fNRRErhouH3vLesqLutqLh4IZFuAsCvbwtkbn7BLgS1hlwJEhkDAMAHCLsMaRrkIKBQg7FKIyVAIGChA2GVg0ygZAUUChF2KIBkmLwWK9u2TUiu0sK7MWbXC/t29e9dBc42UlYt1b6Wm/gOiSyD26y/NPXtpb0LY1b4WlaxfJyVrVkVDrVXv22Fpa4/m3idL45n9pdm6StAKSgcMbN+B2RuBDAQIuzJA02QXwi4FjSDsUoDIEAgYJkDYZVjDKBcBhQKEXQoxGQoBAwUIuwxsGiUjoEiAsEsRJMMUjEDRF59L6cr3o1fvrFltBxxFe788aP7hjh3t0KvxjP7S1P8s+8+hzl20ciLsSr0d7m1b7F6XrF4pJR+sltIPPxDrPnDJj1DXbnbw2XzmWfbv1q9IRUXqB2JLBLIkQNiVJdgcDEvYpQCZsEsBIkMgYJgAYZdhDaNcBBQKEHYpxGQoBAwUIOwysGmUjIAiAcIuRZAMU9AC7u1bpfTDNVK8eoW9DKK1BKK1LGLyI3xUJ2kYO0EiHTpI+OhjpPmkU6S5Vy8JdevhiB9h18HsVt9K/rleXIGAFG/8WLwLF0jJ6hVStH//wf088khpHPBdabKu6hsw0P6ztUQhDwR0FCDs0rErqdVE2JWa02G3IuxSgMgQCBgmQNhlWMMoFwGFAoRdCjEZCgEDBQi7DGwaJSOgSICwSxEkwyCQJGAFJSUf/V1KVq20AzAFp2yIAAAYQ0lEQVTrKjBXwN+qU9OpfSR81NEipSUS6txVQj1PjP7eo4c0d+2WtTCsUMOu4k8+Fvdnn0rxti1SZC1LGQpJ8Yb14nn7jVav0rOaFluqsnHAWdLcf6A0ntEva33hxYRANgQIu7KhmpsxCbsUOBN2KUBkCAQMEyDsMqxhlIuAQgHCLoWYDIWAgQKEXQY2jZIRUCRA2KUIkmEQSEGgZN1HUvyvDVK8aaO4P/mX/Xvxpk+kaO/eNvcOdekqzd17SqhrVwl17ynhrt2kuWt3CXXvLtY9oTJ55GPYZQWKxZs3iXv7NvtX0Y5PpXjzZnHv/EyKt2yWos/3tEnVfNLJ0tzrJGk++RRpPq1v9D5sp/Zpcz82QEBnAcIunbtz+NoIuxT0jrBLASJDIGCYAGGXYQ2jXAQUChB2KcRkKAQMFCDsMrBplIyAIgHCLkWQDINAOwRc/oZoQLNzh7g3b4oGNFu2iNv+fbNY9wlr6xE+5lhpPnBFWHPf70jjwEH2LpGSYgkfe5yEjj9BIl5fi2FMDrt8Lzwn1vKRRdu2SfGn2w7YbWx1ucFku3hw2KWLhHqcKGErSOxxooR69rR/54FAPgoQdpnbVcIuBb0j7FKAyBAIGCZA2GVYwygXAYUChF0KMRkKAQMFCLsMbBolI6BIgLBLESTDIJBlgfjSe1s3S9Fnn4p762Zxf/aZvRSftSRfKo/gOee12KzY7ZJISYk0+colUlEhkfIKCVdUiFRU2n8Pl1e0+Lr1fOJ21tJ+rT1cDfX2Pa9cDQ3R3wN+sQI9VyAY/d3+FTjwNb+9jTQ2tj5WwC+epUvEvXuXFH35RSrTlOYTe9tXvIW6dGu5JGSX6FVxPBAoRAHCLnO7TtiloHeEXQoQGQIBwwQIuwxrGOUioFCAsEshJkMhYKAAYZeBTaNkBBQJEHYpgmQYBBwWcG/bIsV2CLbVviKsaMdn0YDo8z3Rq55SDMQcnkabhw8ffYyErCvVOneRULfuEu7WQ5q7Re9rFrLub3b8CW2OwQYIFKIAYZe5XSfsUtA7wi4FiAyBgGEChF2GNYxyEVAoQNilEJOhEDBQgLDLwKZRMgKKBAi7FEEyDAIGCNhXVNXXS1Fdrf17ZbNfwl/XSXDfVy2+Lg0NUvT1fnEd2M7e3roS6+v9UlRXJ676OnFZvwcDrc464vFKpMxnL5to//Il/Ir93euViK8s+nyZTyT259i2sefLyiV8xBESPu54CZ3Q2QBlSkRATwHCLj37kkpVhF2pKLWxDWGXAkSGQMAwAcIuwxpGuQgoFCDsUojJUAgYKEDYZWDTKBkBRQKEXYogGQYBAwVMvmeXgdyUjICjAoRdjvK36+CEXe3ii+5M2KUAkSEQMEyAsMuwhlEuAgoFCLsUYjIUAgYKEHYZ2DRKRkCRAGGXIkiGQcBAAcIuA5tGyQhkKEDYlSGcBrsRdiloAmGXAkSGQMAwAcIuwxpGuQgoFCDsUojJUAgYKEDYZWDTKBkBRQKEXYogGQYBAwUIuwxsGiUjkKEAYVeGcBrsRtiloAmEXQoQGQIBwwQIuwxrGOUioFCAsEshJkMhYKAAYZeBTaNkBBQJEHYpgmQYBAwUIOwysGmUjECGAoRdGcJpsBthl4ImEHYpQGQIBAwTIOwyrGGUi4BCAcIuhZgMhYCBAoRdBjaNkhFQJEDYpQiSYRAwUICwy8CmUTICGQoQdmUIp8FuhF0KmkDYpQCRIRAwTICwy7CGUS4CCgUIuxRiMhQCBgoQdhnYNEpGQJEAYZciSIZBwEABwi4Dm0bJCGQoQNiVIZwGuxF2KWgCYZcCRIZAwDABwi7DGka5CCgUIOxSiMlQCBgoQNhlYNMoGQFFAoRdiiAZBgEDBQi7DGwaJSOQoQBhV4ZwGuxG2KWgCYRdChAZAgHDBAi7DGsY5SKgUICwSyEmQyFgoABhl4FNo2QEFAkQdimCZBgEDBQg7DKwaZSMQIYChF0ZwmmwG2GXBk2gBAQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAgcwECLsyc2MvBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABDQQIuzRoAiUggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAghkJkDYlZkbeyGAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCGggQNiVYRPmL35D7po+y957xLBBcs+k68TnLc1wNHZDAAFdBfbtr5UbJz8sa9dvipd4wnGdZMb026V3j87213g/0LV71IVA5gLW63rL9l1y28TxLQZZsWaDXHvLNPtrffv0ksem3SodO1TGt+H9IHNz9kRAF4FDvf4fmjFHHn9mcYsy77vjOqkaPtj+WlvvD7rMjzoQQOBggcTXb2v/x/cHGuXuB2bJoqXL7Z0TX/u8/jmjEDBbYOPWHTLxjgdl5+4vW/0eP/n55P8HJH9m8OQjk2Vgv1PNRqF6BBBAwFABwq4MGmd9I/zgjDnxD7is//haj+QPxDIYml0QQEAzgdg3rrdPHN/qN6y8H2jWMMpBoJ0CiR92XX/l8Bb/tlv/0f3PqTPld1NusMNu6wPx5avWxX/ghfeDduKzOwIOCxzu9W+Vdrjv+dt6f3B4ahweAQTaELD+Te/W+Vj7+/1YsHX8sUfFvw9IfP0n//+A1z+nFwJmC1j//m/fsSf+wyvW633Xnr3x7/GTX+OJs429XwwacJq9/+G2NVuJ6hFAAAEzBAi7MuiT9Q9fz27Ht/gpzsTwK4Mh2QUBBDQVaCvs4v1A08ZRFgLtFGjtyo7kryX/Z5b3g3aiszsCmggc7souq8TWfsCtrfcHTaZGGQggkKJA4g+0BIJBmXL/TJl004T4yg6J4Rev/xRR2QwBQwSSf4DtcAGW9dwDf5gtU++8wV7tITn8MmTKlIkAAgjkjQBhV5qtbO0fLn5yI01ENkfAIIHkJQkSlzDk/cCgRlIqAmkKtPZhd/JVHYlh+LdP7WUvbxT7qU7rcHx/kCY6myOgiUCqyxgmLmN2uPcHljLSpLGUgUAaAomv6db+PU8Mwx770wJ75FgQ3tYPy6VRBpsigIADAsmrNxxumcPkYMwql9WfHGgah0QAAQQOCBB2pXkqxD7cHjdySHxJMz7MShORzREwWMD6xnfOwmX2MqZej8f+cJv3A4MbSukIHELgUGFX4pXdrYVdvB9wSiFgvsChwq7EmcU++Jo65Qb7/wTJV3byYbf55wEzKFyB1q7qSLxyw5JJDrsO9f0BYXfhnkfM3EyBVD7fS1zm8KMNm6R64bL4koeEXWb2naoRQCB/BAi70uwlV3KkCcbmCOSZgPXhVWwZk87HHc2VHHnWX6aDQEyAK7s4FxAoXIFUwq7Yh1mxD7i5sqtwzxdmnl8CVtA1ZepMmTH99viShVzZlV89ZjYIHEog+QdZDrddLAD/ZPNnknxbE67s4hxDAAEEnBMg7MrAnntyZIDGLgjkiUBi2NW7R+eDfpK7tWUM8mTqTAOBghLgnl0F1W4mi0ALgUzCLu7Zw0mEgPkCrQVd1qySv/+Phd3W79bShbz+ze89M0Ag1aDLkkq8T9fer2q5ZxenDwIIIKCRAGFXBs1I/jCbn9rIAJFdEDBEwHq9W4/YEiTJ63fzfmBIIykTgTQFWvuwO/knu3k/SBOVzREwRKC117/1YffipcvlR1UX2LNIfj9o6/3BkKlTJgIFK9DWD6wl/p8/eZlSXv8Fe9ow8TwRaGvpwpeWvS8nndg1frVn4vtB8upPbY2VJ2RMAwEEENBWgLArw9ZY/wm+a/ose+8Rwwa1WJ83wyHZDQEENBQ43M1oY+XyfqBh4ygJgQwFrA+7rr1lWou9n3xkcjzwTny+b59e9v37OnaojG/P+0GG8OyGgAYCh3v9xz7MWrR0ebzSxPcG64ttvT9oMEVKQACBQwhYH14//sziFs+ecFyn+HKGye8B991xnVQNHxzfntc/pxYC5gokfv+eOIvYv/PJ3x8kfwYYC8DXrt9k7578/YG5MlSOAAIImCdA2GVez6gYAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEDggABhF6cCAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIICAsQKEXca2jsIRQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQIuzgHEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEjBUg7DK2dRSOAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCBA2MU5gAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAgggYKwAYZexraNwBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABwi7OAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAWMFCLuMbR2FI4AAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIEHZxDiCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCBgrQNhlbOsoHAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAgLCLcwABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQMBYAcIuY1tH4QgggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAoRdnAMIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAALGChB2Gds6CkcAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEECDs4hxAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBAwVoCwy9jWUTgCCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggABhF+cAAggggAACCCCAAAIIIJAjgX37a+XGyQ/L7RPHy8B+p+boqBwGAQQQQAABBBBAAAEEEMhvAcKu/O4vs0MAAQQQQAABBBBAwAiBFWs2yLW3TDuo1uuvHC63TRxvfz0WFI0fOUSqhg82Yl7JRRJ2Gdk2ikYAAQQQQAABBBBAAAHNBQi7NG8Q5SGAAAIIIIAAAgggUAgCVtg1ZepMmTH9dundo7M95Y1bd8jEOx6Um64ZZWy4RdhVCGcvc0QAAQQQQAABBBBAAAGnBQi7nO4Ax0cAAQQQQAABBBBAAAFpLexKvpIr+aqo2N+tMOyFJe/KoqXLbcnEq8Fao7WO9eCMOfZSglbAtnP3l/ZmTz4yOb604PzFb8jyVevknknXic9baj8f2++xabdKxw6VEtvmO6f1lqmP/tXepm+fXmI9/8Tsv8njzyy2vzZi2KD4OLGaf3zFxfLEsy/K2vWb7G3uu+O6FoFebLvY862NkTjvxOc5nRBAAAEEEEAAAQQQQACBQhMg7Cq0jjNfBBBAAAEEEEAAAQQ0FGgt7Er+2qHCri/27o9fERa7GmzqlBsOeU+s2JKJiQGRFVzNWbjMDqoSg6y2wq67ps+KB1X+QKPc/cAsO3SLhVexrw0acJodZsXmYLUgdqzkmltb6vChGXNk1569dmgWCAbt+34lzlvDllISAggggAACCCCAAAIIIJAzAcKunFFzIAQQQAABBBBAAAEEEDiUwKHu2RW7UsoKoA4VdllXaA3sd6o9dHK41Nrxkq/QsraxAqf/nDpTfjflBnsZxXSu7EoMxFrbL/FrsaAqsWbr+FaYZT2s+5NZ22/Zvit+r7Lk+o46stIOu5LH4OxCAAEEEEAAAQQQQAABBApVgLCrUDvPvBFAAAEEEEAAAQQQ0EigtSu7rPISr7iy/p4Y8rR2BZSpYVdiIPbYnxbEl0BMbNEJx3Wyr2Aj7NLoxKUUBBBAAAEEEEAAAQQQ0EKAsEuLNlAEAggggAACCCCAAAKFLXCosCsx0DrpxC4FE3ZZZ4N1lVdrj9ZCvsI+e5g9AggggAACCCCAAAIIFLoAYVehnwHMHwEEEEAAAQQQQAABDQQOFXYl3s8qn8Ou5GUMl69aZ9+fy+ctPag7hF0anLCUgAACCCCAAAIIIIAAAloJEHZp1Q6KQQABBBBAAAEEEECgMAUOFXZZIdD7azbIY9NutWFytYxhcj2xgMmqwarFuodYW/fnigVVbd2zK/lYsYBv+NCz41d3WcszWssb/njCJQc5FOYZw6wRQAABBBBAAAEEEEAAgW8ECLs4GxBAAAEEEEAAAQQQQMBxASvwufaWaQfVMWLYoPgVTslXNGXznl1WIVZIddf0WXZNffv0kh9fcbE88eyLSsKutes3xecauxdX7x6d41+LzS1xu+uvHG6HX1zZ5fjpSgEIIIAAAggggAACCCCgmQBhl2YNoRwEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAIHUBQi7UrdiSwQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAc0ECLs0awjlIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIpC5A2JW6FVsigAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAghoJkDYpVlDKAcBBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQCB1AcKu1K3YEgEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAQDMBwi7NGkI5CCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACqQsQdqVuxZYIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAKaCRB2adYQykEAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEhdgLArdSu2RAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQ0EyAsEuzhlAOAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIBA6gKEXalbsSUCCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggIBmAoRdmjWEchBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBFIXIOxK3YotEUAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEENBMg7NKsIZSDAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCQugBhV+pWbIkAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIKCZAGGXZg2hHAQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAgdQFCLtSt2JLBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABzQQIuzRrCOUggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAgikLkDYlboVWyKAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCGgmQNilWUMoBwEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAIHUBwq7UrdgSAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEBAMwHCLs0aQjkIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAKpCxB2pW7FlggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAApoJEHZp1hDKQQABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQSF2AsCt1K7ZEAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBDQTICwS7OGUA4CCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggEDqAoRdqVuxJQIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAgGYChF2aNYRyEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEUhcg7Erdii0RQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQ0EyDs0qwhlIMAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIJC6AGFX6lZsiQACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAgggoJkAYZdmDaEcBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQACB1AUIu1K3YksEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAHNBAi7NGsI5SCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCKQuQNiVuhVbIoAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIaCZA2KVZQygHAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAgdQHCrtSt2BIBBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQEAzAcIuzRpCOQgggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAqkL/H/DsbWYmwaT7AAAAABJRU5ErkJggg==",
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = px.line(data_frame=bio.system_snapshot(), y=[\"A\"], \n",
" title= \"Initial System State\",\n",
" color_discrete_sequence = ['red'],\n",
" labels={\"value\":\"concentration\", \"variable\":\"Chemical\", \"index\":\"Bin number\"})\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"id": "2ba3f95a-2d63-4a11-a150-c5f1ea4a7716",
"metadata": {},
"source": [
"#### Now do 4 rounds of single-step diffusion, to collect the system state at a total of 5 time points: t0 (the initial state), plus t1, t2, t3 and t4"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "28e46cc9-4ddf-4d9c-8ab2-20c08bb3106e",
"metadata": {},
"outputs": [],
"source": [
"# All the system state will get collected in this object\n",
"history = MovieArray()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "122c2273-3f40-4cd0-8ed6-77c03f0f6ed7",
"metadata": {},
"outputs": [],
"source": [
"# Store the initial state\n",
"arr = bio.lookup_species(species_index=0, copy=True)\n",
"history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "45c018bc-5f69-4da2-b8c2-4c412e10bc5e",
"metadata": {},
"outputs": [],
"source": [
"# Do the 4 rounds of single-step diffusion; accumulate all data in the history object\n",
"for _ in range(4):\n",
" bio.diffuse(time_step=delta_t, n_steps=1, delta_x=delta_x , algorithm=algorithm)\n",
"\n",
" arr = bio.lookup_species(species_index=0, copy=True)\n",
" history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "a77b0423-fd2b-4219-8dc0-4f2fa65813be",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(5, 300)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Now, let's examine the data collected at the 5 time points\n",
"all_history = history.get_array()\n",
"all_history.shape "
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "128ff22d-0ce2-45ed-8ec0-2a067a613990",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(5, 300)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compute time derivatives (for each bin)\n",
"df_dt_all_bins = np.apply_along_axis(np.gradient, 0, all_history, delta_t)\n",
"df_dt_all_bins.shape"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "f32bd8a9-b062-4db9-9e20-cd230464e6f8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Let's consider the state at the midpoint in time (t2)\n",
"f_at_t2 = all_history[2] # The middle of the 5 time snapshots\n",
"f_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "398588e7-2066-4171-b95c-089b83126034",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Computer the second spatial derivative (simple method)\n",
"gradient_x_at_t2 = np.gradient(f_at_t2, delta_x)\n",
"second_gradient_x_at_t2 = np.gradient(gradient_x_at_t2, delta_x)\n",
"second_gradient_x_at_t2.shape"
]
},
{
"cell_type": "markdown",
"id": "8849a5a0-4a36-4e9a-b0d6-34823d147bed",
"metadata": {},
"source": [
"#### Now, let's use the computed values to see how the left- and right-hand side of the diffusion equation compare\n",
"at all bins (except 2 at each edge), at the middle simulation time t2"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "8071dbaf-2d87-44e6-b0fc-b361d6094a35",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lhs = df_dt_all_bins[2] # t2 is the middle point of the 5\n",
"lhs.shape"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "3909d680-2ea2-43ca-9094-bc370ee61c43",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rhs = diffusion_rate*second_gradient_x_at_t2\n",
"rhs.shape"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "19514cd4-ceb5-450d-8c3e-9b43456f02c0",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.023163289760024783"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"num.compare_vectors(lhs, rhs, trim_edges=2) # Euclidean distance, ignoring 2 edge points at each end"
]
},
{
"cell_type": "markdown",
"id": "258e81c8-0720-4a70-a7e0-af06167490d8",
"metadata": {},
"source": [
"The above number is a measure of the discrepancy from the perfect match (zero distance) that an ideal solution would provide. \n",
"\n",
"Notice how vastly closer to zero it is than the counterpart value we got with the very coarse simulation in experiment \"validate_diffusion_1\""
]
},
{
"cell_type": "markdown",
"id": "5c069ce7-2e2b-4f13-bf2f-32509c865871",
"metadata": {},
"source": [
"# VARIATIONS on the BASELINE\n",
"We'll now tweak some parameters, and observe the effect on the estimate of the discrepancy from the exact diffusion equation\n",
"\n",
"The baseline distance was **0.023163289760024783**"
]
},
{
"cell_type": "markdown",
"id": "3002d2cd-5142-483b-86db-126c3f73413f",
"metadata": {},
"source": [
"## Variation 1 : reduce spatial resolution\n",
"Expectation: greater discrepancy"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "e3918a7f-10a6-4192-8e00-bfd1b7d6627e",
"metadata": {},
"outputs": [],
"source": [
"n_bins = 100 # Reducing the spatial resolution"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "1523c8e7-ac8a-42d9-bc0a-9b9c6ea00f56",
"metadata": {},
"outputs": [],
"source": [
"# Initialize the system\n",
"bio = BioSim1D(n_bins=n_bins, chem_data=chem_data)\n",
"\n",
"# Initialize the concentrations to 2 superposed sine waves\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=1, amplitude=10, bias=50)\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=2, amplitude=8)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "0b2fd9b1-ed93-440d-b67e-ab9a660bef7e",
"metadata": {},
"outputs": [],
"source": [
"history = MovieArray() # All the system state will get collected in this object\n",
"# Store the initial state\n",
"arr = bio.lookup_species(species_index=0, copy=True)\n",
"history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "020e44ca-3dae-47f5-bdb2-377d7efe7ee2",
"metadata": {},
"outputs": [],
"source": [
"# Do the 4 rounds of single-step diffusion; accumulate all data in the history object\n",
"for _ in range(4):\n",
" bio.diffuse(time_step=delta_t, n_steps=1, delta_x=delta_x , algorithm=algorithm)\n",
"\n",
" arr = bio.lookup_species(species_index=0, copy=True)\n",
" history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "0aa5578c-4454-48f4-be62-5fe58ea7c5e1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(5, 100)"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Now, let's examine the data collected at the 5 time points\n",
"all_history = history.get_array()\n",
"all_history.shape "
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "03a64e4c-6b8a-4732-afa8-199028190340",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(100,)"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compute time derivatives (for each bin) : simple method, using central differences\n",
"df_dt_all_bins = np.apply_along_axis(np.gradient, 0, all_history, delta_t)\n",
"\n",
"# Let's consider the state at the midpoint in time (t2)\n",
"f_at_t2 = all_history[2] # The middle of the 5 time snapshots\n",
"f_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "4998952c-d051-4d3e-83c2-f443d3e67e47",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(100,)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Computer the second spatial derivative (simple method, using central differences)\n",
"gradient_x_at_t2 = np.gradient(f_at_t2, delta_x)\n",
"second_gradient_x_at_t2 = np.gradient(gradient_x_at_t2, delta_x)\n",
"second_gradient_x_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "885a0f51-eaa9-43f8-ae5e-fe8a6bd7bcc9",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0705630110277808"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compare the left and right hand sides of the diffusion equation\n",
"lhs = df_dt_all_bins[2] # t2 is the middle point of the 5\n",
"rhs = diffusion_rate*second_gradient_x_at_t2\n",
"\n",
"num.compare_vectors(lhs, rhs, trim_edges=2) # Euclidean distance, ignoring 2 edge points at each end"
]
},
{
"cell_type": "markdown",
"id": "76a989a9-311c-4e58-9358-5a51b88f705e",
"metadata": {},
"source": [
"#### The discrepancy value has indeed gone up from the baseline 0.023163289760024783"
]
},
{
"cell_type": "markdown",
"id": "976f3b8d-bf0e-4a23-978f-0f2064d654f6",
"metadata": {},
"source": [
"## Variation 2 : increase spatial resolution\n",
"Expectation: smaller discrepancy"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "aab88c73-3afa-48f6-959e-1e9043350cc6",
"metadata": {},
"outputs": [],
"source": [
"n_bins = 900 # Increasing the spatial resolution (from the baseline of 300)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "e996032a-63bd-4d41-8cce-b97a158c82b0",
"metadata": {},
"outputs": [],
"source": [
"# Initialize the system\n",
"bio = BioSim1D(n_bins=n_bins, chem_data=chem_data)\n",
"\n",
"# Initialize the concentrations to 2 superposed sine waves\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=1, amplitude=10, bias=50)\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=2, amplitude=8)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "6e26a594-77d6-49d8-9e35-ad2e5d03bdff",
"metadata": {},
"outputs": [],
"source": [
"history = MovieArray() # All the system state will get collected in this object\n",
"# Store the initial state\n",
"arr = bio.lookup_species(species_index=0, copy=True)\n",
"history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "75f95eee-5e93-41a4-8d6a-6663b455faa5",
"metadata": {},
"outputs": [],
"source": [
"# Do the 4 rounds of single-step diffusion; accumulate all data in the history object\n",
"for _ in range(4):\n",
" bio.diffuse(time_step=delta_t, n_steps=1, delta_x=delta_x , algorithm=algorithm)\n",
"\n",
" arr = bio.lookup_species(species_index=0, copy=True)\n",
" history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "2e75b83d-ad36-441e-be9e-e4c7321eb1b4",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(5, 900)"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Now, let's examine the data collected at the 5 time points\n",
"all_history = history.get_array()\n",
"all_history.shape "
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "1c5accca-8be9-4410-8cf7-5122774e7306",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(900,)"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compute time derivatives (for each bin) : simple method, using central differences\n",
"df_dt_all_bins = np.apply_along_axis(np.gradient, 0, all_history, delta_t)\n",
"\n",
"# Let's consider the state at the midpoint in time (t2)\n",
"f_at_t2 = all_history[2] # The middle of the 5 time snapshots\n",
"f_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "d103fc3f-1578-4d7f-88de-87a7fc896c5a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(900,)"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Computer the second spatial derivative (simple method, using central differences)\n",
"gradient_x_at_t2 = np.gradient(f_at_t2, delta_x)\n",
"second_gradient_x_at_t2 = np.gradient(gradient_x_at_t2, delta_x)\n",
"second_gradient_x_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "fab3cfc4-407f-4db0-af21-dda78212933f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.007721522026370068"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compare the left and right hand sides of the diffusion equation\n",
"lhs = df_dt_all_bins[2] # t2 is the middle point of the 5\n",
"rhs = diffusion_rate*second_gradient_x_at_t2\n",
"\n",
"num.compare_vectors(lhs, rhs, trim_edges=2) # Euclidean distance, ignoring 2 edge points at each end"
]
},
{
"cell_type": "markdown",
"id": "5fbf609b-e96a-4281-bfa2-f2f0f43fae72",
"metadata": {},
"source": [
"#### The discrepancy value has indeed gone down from the baseline 0.023163289760024783"
]
},
{
"cell_type": "markdown",
"id": "27746f62-229e-4915-9b1d-4dac6190eade",
"metadata": {},
"source": [
"## Variation 3 : reduce time resolution\n",
"Expectation: greater discrepancy"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "a143cba9-2e3c-4b0a-a357-7eca88cc454f",
"metadata": {},
"outputs": [],
"source": [
"n_bins = 300 # Restoring the baseline number of bins\n",
"delta_t = 0.02 # Reducing the time resolution (from baselin 0.01)"
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "3fd3d708-dc6f-4387-b49a-090813c52d4f",
"metadata": {},
"outputs": [],
"source": [
"# Initialize the system\n",
"bio = BioSim1D(n_bins=n_bins, chem_data=chem_data)\n",
"\n",
"# Initialize the concentrations to 2 superposed sine waves\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=1, amplitude=10, bias=50)\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=2, amplitude=8)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "cf1d4c86-17be-4c53-b7d6-c14fc80a1604",
"metadata": {},
"outputs": [],
"source": [
"history = MovieArray() # All the system state will get collected in this object\n",
"# Store the initial state\n",
"arr = bio.lookup_species(species_index=0, copy=True)\n",
"history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "b31a59ae-8d90-4cd1-beb6-0a2c518502b2",
"metadata": {},
"outputs": [],
"source": [
"# Do the 4 rounds of single-step diffusion; accumulate all data in the history object\n",
"for _ in range(4):\n",
" bio.diffuse(time_step=delta_t, n_steps=1, delta_x=delta_x , algorithm=algorithm)\n",
"\n",
" arr = bio.lookup_species(species_index=0, copy=True)\n",
" history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "50f13dcb-759d-410d-b7a3-925e475a927c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(5, 300)"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Now, let's examine the data collected at the 5 time points\n",
"all_history = history.get_array()\n",
"all_history.shape "
]
},
{
"cell_type": "code",
"execution_count": 38,
"id": "e2f93e1b-4dd9-4163-96a7-30c6552fce91",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compute time derivatives (for each bin) : simple method, using central differences\n",
"df_dt_all_bins = np.apply_along_axis(np.gradient, 0, all_history, delta_t)\n",
"\n",
"# Let's consider the state at the midpoint in time (t2)\n",
"f_at_t2 = all_history[2] # The middle of the 5 time snapshots\n",
"f_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 39,
"id": "149b6566-22e0-4994-82d4-adc67c4fb08e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Computer the second spatial derivative (simple method, using central differences)\n",
"gradient_x_at_t2 = np.gradient(f_at_t2, delta_x)\n",
"second_gradient_x_at_t2 = np.gradient(gradient_x_at_t2, delta_x)\n",
"second_gradient_x_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 40,
"id": "9c62a79a-3188-4ffc-9b79-abd2a69e12fa",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.04453006107638132"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compare the left and right hand sides of the diffusion equation\n",
"lhs = df_dt_all_bins[2] # t2 is the middle point of the 5\n",
"rhs = diffusion_rate*second_gradient_x_at_t2\n",
"\n",
"num.compare_vectors(lhs, rhs, trim_edges=2) # Euclidean distance, ignoring 2 edge points at each end"
]
},
{
"cell_type": "markdown",
"id": "c79b51f0-da17-47cf-b4bb-e423951e121a",
"metadata": {},
"source": [
"#### The discrepancy value has indeed gone up from the baseline 0.023163289760024783"
]
},
{
"cell_type": "markdown",
"id": "4ff62063-0f8b-4d28-8b47-44bcf0396874",
"metadata": {},
"source": [
"## Variation 4 : increase time resolution\n",
"Expectation: smaller discrepancy"
]
},
{
"cell_type": "code",
"execution_count": 41,
"id": "967b3cb4-7d89-4114-86bb-8dfe76f6cb71",
"metadata": {},
"outputs": [],
"source": [
"delta_t = 0.005 # Increasing the time resolution (from baselin 0.01)"
]
},
{
"cell_type": "code",
"execution_count": 42,
"id": "4ae7056d-6c3d-4db4-93ca-c9688fef52f7",
"metadata": {},
"outputs": [],
"source": [
"# Initialize the system\n",
"bio = BioSim1D(n_bins=n_bins, chem_data=chem_data)\n",
"\n",
"# Initialize the concentrations to 2 superposed sine waves\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=1, amplitude=10, bias=50)\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=2, amplitude=8)"
]
},
{
"cell_type": "code",
"execution_count": 43,
"id": "81a3bd58-8ef4-43fb-9477-b27ab0160d2c",
"metadata": {},
"outputs": [],
"source": [
"history = MovieArray() # All the system state will get collected in this object\n",
"# Store the initial state\n",
"arr = bio.lookup_species(species_index=0, copy=True)\n",
"history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 44,
"id": "e32ba5c9-185c-4d84-8c6d-415a55cf98e7",
"metadata": {},
"outputs": [],
"source": [
"# Do the 4 rounds of single-step diffusion; accumulate all data in the history object\n",
"for _ in range(4):\n",
" bio.diffuse(time_step=delta_t, n_steps=1, delta_x=delta_x , algorithm=algorithm)\n",
"\n",
" arr = bio.lookup_species(species_index=0, copy=True)\n",
" history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 45,
"id": "b92026f5-82e8-4d4d-b376-a9e9e46d0c00",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(5, 300)"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Now, let's examine the data collected at the 5 time points\n",
"all_history = history.get_array()\n",
"all_history.shape "
]
},
{
"cell_type": "code",
"execution_count": 46,
"id": "4fdb6d6e-df25-4a1d-a886-0125f82e897c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compute time derivatives (for each bin) : simple method, using central differences\n",
"df_dt_all_bins = np.apply_along_axis(np.gradient, 0, all_history, delta_t)\n",
"\n",
"# Let's consider the state at the midpoint in time (t2)\n",
"f_at_t2 = all_history[2] # The middle of the 5 time snapshots\n",
"f_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 47,
"id": "8db28e0f-b63d-43d9-8913-87c62c87739e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Computer the second spatial derivative (simple method, using central differences)\n",
"gradient_x_at_t2 = np.gradient(f_at_t2, delta_x)\n",
"second_gradient_x_at_t2 = np.gradient(gradient_x_at_t2, delta_x)\n",
"second_gradient_x_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 48,
"id": "01e1ec43-d7b9-4b6b-a57e-2026c400c0ff",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.011808914990091101"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compare the left and right hand sides of the diffusion equation\n",
"lhs = df_dt_all_bins[2] # t2 is the middle point of the 5\n",
"rhs = diffusion_rate*second_gradient_x_at_t2\n",
"\n",
"num.compare_vectors(lhs, rhs, trim_edges=2) # Euclidean distance, ignoring 2 edge points at each end"
]
},
{
"cell_type": "markdown",
"id": "f03489cb-5a2c-490f-a48d-50b3cbe75edd",
"metadata": {},
"source": [
"#### The discrepancy value has indeed gone down from the baseline 0.023163289760024783"
]
},
{
"cell_type": "markdown",
"id": "2e42b48c-0c11-4a33-80d7-4c51bbf1fac3",
"metadata": {},
"source": [
"## Variation 5 : Use a better (higher-order) method to numerically estimate the SPATIAL derivatives\n",
"Expectation: presumably smaller discrepancy (unless dwarfed by errors from approximations in the simulation)"
]
},
{
"cell_type": "code",
"execution_count": 49,
"id": "69401e31-a5f6-4faf-b23f-024a73a6e805",
"metadata": {},
"outputs": [],
"source": [
"n_bins = 300 # Restoring the baseline number of bins\n",
"delta_t = 0.01 # Reducing the baseline time resolution"
]
},
{
"cell_type": "code",
"execution_count": 50,
"id": "f732c9a0-4b95-4bcf-821e-0221ba7170e9",
"metadata": {},
"outputs": [],
"source": [
"# Initialize the system\n",
"bio = BioSim1D(n_bins=n_bins, chem_data=chem_data)\n",
"\n",
"# Initialize the concentrations to 2 superposed sine waves\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=1, amplitude=10, bias=50)\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=2, amplitude=8)"
]
},
{
"cell_type": "code",
"execution_count": 51,
"id": "98b3b3bf-3ed7-4867-895b-3ae39bfb0424",
"metadata": {},
"outputs": [],
"source": [
"history = MovieArray() # All the system state will get collected in this object\n",
"# Store the initial state\n",
"arr = bio.lookup_species(species_index=0, copy=True)\n",
"history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 52,
"id": "38387573-8c4c-4b8f-9e28-f5d3196cc3c6",
"metadata": {},
"outputs": [],
"source": [
"# Do the 4 rounds of single-step diffusion; accumulate all data in the history object\n",
"for _ in range(4):\n",
" bio.diffuse(time_step=delta_t, n_steps=1, delta_x=delta_x , algorithm=algorithm)\n",
"\n",
" arr = bio.lookup_species(species_index=0, copy=True)\n",
" history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 53,
"id": "27729fc9-a692-43ee-9ac3-fef64ade8c0e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(5, 300)"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Now, let's examine the data collected at the 5 time points\n",
"all_history = history.get_array()\n",
"all_history.shape "
]
},
{
"cell_type": "code",
"execution_count": 54,
"id": "03150f7c-590e-4531-bea4-e18d8df7ffa3",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compute time derivatives (for each bin) : simple method, using central differences\n",
"df_dt_all_bins = np.apply_along_axis(np.gradient, 0, all_history, delta_t)\n",
"\n",
"# Let's consider the state at the midpoint in time (t2)\n",
"f_at_t2 = all_history[2] # The middle of the 5 time snapshots\n",
"f_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 55,
"id": "dc283a71-7228-486f-8cbc-09d1a2e51c68",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 55,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Computer the second spatial derivative (MORE SOPHISTICATED METHOD, using 5-point stencils)\n",
"gradient_x_at_t2 = num.gradient_order4_1d(arr=f_at_t2, dx=delta_x)\n",
"second_gradient_x_at_t2 = num.gradient_order4_1d(arr=gradient_x_at_t2, dx=delta_x)\n",
"second_gradient_x_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 56,
"id": "92b9fc7a-659f-4a62-a0a4-5fe8de345d0e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.006831207619477847"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compare the left and right hand sides of the diffusion equation\n",
"lhs = df_dt_all_bins[2] # t2 is the middle point of the 5\n",
"rhs = diffusion_rate*second_gradient_x_at_t2\n",
"\n",
"num.compare_vectors(lhs, rhs, trim_edges=2) # Euclidean distance, ignoring 2 edge points at each end"
]
},
{
"cell_type": "markdown",
"id": "abdddfd0-e02b-46fc-a3b4-2b9581740a0b",
"metadata": {},
"source": [
"#### Not surprisingly, the discrepancy value (in part caused by poor numeric estimation of spatial derivatives) has indeed gone down from the baseline 0.023163289760024783"
]
},
{
"cell_type": "markdown",
"id": "38642e67-a039-416d-9945-c0a020d26490",
"metadata": {},
"source": [
"## Variation 6 : Use a better (higher-order) method to numerically estimate the TIME derivatives\n",
"Expectation: presumably smaller discrepancy (unless dwarfed by errors from approximations in the simulation)\n",
"\n",
"(Note: we revert to the original method for the *spatial* derivatives)"
]
},
{
"cell_type": "code",
"execution_count": 57,
"id": "ba76ad1c-97e7-463c-ac2c-30a283db42dc",
"metadata": {},
"outputs": [],
"source": [
"# Initialize the system\n",
"bio = BioSim1D(n_bins=n_bins, chem_data=chem_data)\n",
"\n",
"# Initialize the concentrations to 2 superposed sine waves\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=1, amplitude=10, bias=50)\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=2, amplitude=8)"
]
},
{
"cell_type": "code",
"execution_count": 58,
"id": "9d0a38bf-6c59-4dee-87e8-6aa2e149fb47",
"metadata": {},
"outputs": [],
"source": [
"history = MovieArray() # All the system state will get collected in this object\n",
"# Store the initial state\n",
"arr = bio.lookup_species(species_index=0, copy=True)\n",
"history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 59,
"id": "47078710-eb53-44c5-83bc-406dbbd5173b",
"metadata": {},
"outputs": [],
"source": [
"# Do the 4 rounds of single-step diffusion; accumulate all data in the history object\n",
"for _ in range(4):\n",
" bio.diffuse(time_step=delta_t, n_steps=1, delta_x=delta_x , algorithm=algorithm)\n",
"\n",
" arr = bio.lookup_species(species_index=0, copy=True)\n",
" history.store(par=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 60,
"id": "f8d28717-210e-4adc-8d53-ceeb84b3bb26",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(5, 300)"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Now, let's examine the data collected at the 5 time points\n",
"all_history = history.get_array()\n",
"all_history.shape "
]
},
{
"cell_type": "code",
"execution_count": 61,
"id": "028d732d-f97b-4103-ae94-85a60cfce4dc",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compute time derivatives (for each bin) : MORE SOPHISTICATED METHOD, using 5-point stencils\n",
"df_dt_all_bins = np.apply_along_axis(num.gradient_order4_1d, 0, all_history, delta_t)\n",
"\n",
"# Let's consider the state at the midpoint in time (t2)\n",
"f_at_t2 = all_history[2] # The middle of the 5 time snapshots\n",
"f_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 62,
"id": "ca637810-beda-4abf-b232-473b17e1a76d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Computer the second spatial derivative (BACK TO SIMPLE METHOD, using central differences)\n",
"gradient_x_at_t2 = np.gradient(f_at_t2, delta_x)\n",
"second_gradient_x_at_t2 = np.gradient(gradient_x_at_t2, delta_x)\n",
"second_gradient_x_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 63,
"id": "845f75bd-25fe-4120-8975-e8dcd8e36afc",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.023351269997869635"
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compare the left and right hand sides of the diffusion equation\n",
"lhs = df_dt_all_bins[2] # t2 is the middle point of the 5\n",
"rhs = diffusion_rate*second_gradient_x_at_t2\n",
"\n",
"num.compare_vectors(lhs, rhs, trim_edges=2) # Euclidean distance, ignoring 2 edge points at each end"
]
},
{
"cell_type": "markdown",
"id": "fd1d0295-dfe2-4e7d-9457-cc2fea1a98f2",
"metadata": {},
"source": [
"#### Here, the discrepancy value has remained largely the same from the baseline 0.023163289760024783"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "30dffc16-7072-4a62-b606-642138990e5f",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
},
"toc-autonumbering": false,
"toc-showcode": true,
"toc-showmarkdowntxt": false
},
"nbformat": 4,
"nbformat_minor": 5
}