{
"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: July 14, 2023"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "30dffc16-7072-4a62-b606-642138990e5f",
"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": "a1bc9cbe",
"metadata": {},
"outputs": [],
"source": [
"\n",
"from src.life_1D.bio_sim_1d import BioSim1D\n",
"from src.modules.chemicals.chem_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"
]
},
{
"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": "iVBORw0KGgoAAAANSUhEUgAABb8AAAFoCAYAAAB38YZnAAAAAXNSR0IArs4c6QAAIABJREFUeF7s3Ql4VNX9xvF3ss8kIYALyCK7gIq1tSht3e2KBZRWxR2xaNUuWot/qLXWukC1Lt20lIpoxSoqbpVqa5XihlKtdQFcQFDBpcqShJkks/2fc8PEJGSZhHtn7vK9z2MjmTvn/s7ndwP2ncO5oXQ6nRYHAggggAACCCCAAAIIIIAAAggggAACCCCAAAI+EggRfvuom0wFAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAwBIg/OZGQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEPCdAOG371rKhBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQIv7kHEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBHwnQPjtu5YyIQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAHCb+4BBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAd8JEH77rqVMCAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQIDwm3sAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAwHcChN++aykTQgABBBBAAAEEEEAAAQQQQAABBBBAAAEEECD85h5AAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQ8J0A4bfvWsqEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBAi/uQcQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEfCdA+O27ljIhBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAcJv7gEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAAB3wkQfvuupUwIAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAgPCbewABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEDAdwKE375rKRNCAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQIPzmHkAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBDwnQDht+9ayoQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEECL+5BxBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQR8J0D47buWMiEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABwm/uAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAHfCRB++66lTAgBBBBAAAEEEEAAAQQQQAABBBBAAAEEEECA8Jt7AAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQMB3AoTfvmspE0IAAQQQQAABBBBAAAEEEEAAAQQQQAABBBAg/OYeQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEPCdAOG371rKhBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQIv7kHEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBHwnQPjtu5YyIQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAHCb+4BBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAd8JEH77rqVMCAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQIDwm3sAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAwHcChN++aykTQgABBBBAAAEEEEAAAQQQQAABBBBAAAEEECD85h5AAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQ8J0A4bfvWsqEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBAi/uQcQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEfCdA+O27ln46oevmLtLzL63WTXMuUK+qSh/PlKkhgAACCCCAAAIIIIAAAggggAACCCCAAAItBQi/s7wjVry0WlPPn9Pi7MsvmqbJ4w/NcoTG0zZvrdE5M6/Xnv1212UzpilcVqI16zfq7Iuu1bmnT9phPBNgL3n8Oc29+kING9SvS9fqavi9eMkyXXL1/B2useCGmRq7/6guXTvbk7taY7bj2nVerK5Bl14zXw//c3mLIceMHtriQ4VMXw/cf5R+dPbx3bq8HWN068K8CQEEEEAAAQQQQAABBBBAAAEEEEAAAR8KEH53oamZIPSdjR91ezW1G8PvzLxefPXNHUJ2E07f/Jcl6k7Qnw2tm8PvzIcSn9t3RNMHFc0/wPh409YmLzuCazvGyMaccxBAAAEEEEAAAQQQQAABBBBAAAEEEAiCAOF3F7psR/jd1uU6WvndhfJ2ODXbYDmzqr29Fd6mvrfefk9fO/zAnSmnzfdmW6PtF85iwM5qe3Tp8xo+ZIC1It+O4NqOMbKYFqcggAACCCCAAAIIIIAAAggggAACCCAQCAHC7y60uaPwOxOUXjVrun4ye55eWbXWGvnoo8a1uWo4sz1GJvh+/8NPWlRy5onjre0zWgew7W3D0fo6ZrDOwtvMBc12Jzfe+kBWW6t0dG7r19qqtXmdmVXlzSe+R59dWtTReiuW9jyPn3C4NUzzbVsyq9VbXyebVeyZ2s2Yme1p2rtVMqF1pueZ8zK11tXXW1vdtH4902NzfmdjmO1xzNGZRxduZ05FAAEEEEAAAQQQQAABBBBAAAEEEEDA1wKE311ob2fht9kepPle0G2t5G3re53t+d38oZWmhmtuulMnT/5y0x7g7QW12YbfmZXfbQXorXnaW52cqaHv7r2t0L6tmtp6b3s1tmfd+vzmoXHzULt5SNz6+9kG/V3Z8qWjVdvmtdm/WahZPzi56cGj2d4bGf9sPbpwO3MqAggggAACCCCAAAIIIIAAAggggAACvhYg/O5CezsLv5uH1JlhTQi76KGlTXuE72z43V65JsCeNXtei1XT2YbfZsy2HnbZehV2e3My32+9dUp7gb4xXLb8paYtVNqrsa35mOu0HjfjaVZ+N3/4aFe/35ZrV1bZd2fLkmzujUxd2Xp04XbmVAQQQAABBBBAAAEEEEAAAQQQQAABBHwtQPjdhfZ2N/xuvtLYrvA7m7C6K+F3hqGtbVhah+BtBdvmWh98tKlpi5DmK7Lb20vcXLO9GluPl6mv9Qrzrobc7Z3f2W3QlnfzlfLZhN9tbfPS2d8UyNSVrUdn8+B1BBBAAAEEEEAAAQQQQAABBBBAAAEEgiJA+N2FTrsh/M6ErB9v2tpilffOrvxujyETdH9u3xEt9r5uHsZu/PBjnX3RtZo9a7rG7j+qaai29rFuHva2F363t+K6eY2Z/bJzFX639smE4ZktVToKvzOGu/auavobAGa8bFd+d8WjC7czpyKAAAIIIIAAAggggAACCCCAAAIIIOBrAcLvLrTXDeF3ew+cdCr8bi+gzgS6JvB+d+NHWv7Cyk4fDNnW3uJdXfndul35Cr9bh90dhd/tzTHb8DvTg+Yr67tw23IqAggggAACCCCAAAIIIIAAAggggAACgRQg/O5C250KvzvaiqN1cNqVPbKz3fbk0aXPa/iQAU0P0GxO0t7DNJt//52NH6n1ntsmHI/V1WvfkUNaCLeuqXUAnDm5vZA/V+G3md+9D/9L3zr6MIXLSna4S1pv/dJ6O5bMG9rzM6+3nnt7Y2TOzfZBnV24pTkVAQQQQAABBBBAAAEEEEAAAQQQQAAB3woQfnehtU6F3x0FpK3D4tYPljTlZ77Xem/ubMPvzBYezfewNuNm6nrx1TdbbLGSIctct/VWJub1tvYFb2t1dFvz6ezapl5zmAdcOrXyu/lWI5mtTTLz7mgrmPYeeto6uM6YZ7MNTFc8unA7cyoCCCCAAAIIIIAAAggggAACCCCAAAK+FiD8zrK9mZC2+enNQ9GOtrbo7IGXZszW+2Nn9rRua9zWtZhzDzloP82aPa9FSJ1t+N08QG/NkamjLabOHh7ZmVlmzOYPk2wd4Lf1oMm2HhLZeuV5V0Px9m6Dth5Sac5t6yGerffmbv5hQut5mHvHHIseWtpiH/COxjDnd+aR5e3MaQgggAACCCCAAAIIIIAAAggggAACCPhegPDb9y12boLZbk3iXAWMjAACCCCAAAIIIIAAAggggAACCCCAAAIItC1A+M2d0S2Bjh7w2K0BeRMCCCCAAAIIIIAAAggggAACCCCAAAIIIGCjAOG3jZhBGopV30HqNnNFAAEEEEAAAQQQQAABBBBAAAEEEEDAewKE397rGRUjgAACCCCAAAIIIIAAAggggAACCCCAAAIIdCJA+M0tggACCCCAAAIIIIAAAggggAACCCCAAAIIIOA7AcJv37WUCSGAAAIIIIAAAggggAACCCCAAAIIIIAAAggQfnMPIIAAAggggAACCCCAAAIIIIAAAggggAACCPhOgPDbdy1lQggggAACCCCAAAIIIIAAAggggAACCCCAAAKE39wDCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAr4TIPz2XUuZEAIIIIAAAggggAACCCCAAAIIIIAAAggggADhN/cAAggggAACCCCAAAIIIIAAAggggAACCCCAgO8ECL9911ImhAACCCCAAAIIIIAAAggggAACCCCAAAIIIED4zT2AAAIIIIAAAggggAACCCCAAAIIIIAAAggg4DsBwm/ftZQJIYAAAggggAACCCCAAAIIIIAAAggggAACCBB+cw8ggAACCCCAAAIIIIAAAggggAACCCCAAAII+E6A8Nt3LWVCCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAoTf3AMIIIAAAggggAACCCCAAAIIIIAAAggggAACvhMg/PZdS5kQAggggAACCCCAAAIIIIAAAggggAACCCCAAOE39wACCCCAAAIIIIAAAggggAACCCCAAAIIIICA7wQIv33XUiaEAAIIIIAAAggggAACCCCAAAIIIIAAAgggQPjNPYAAAggggAACCCCAAAIIIIAAAggggAACCCDgOwHCb9+1lAkhgAACCCCAAAIIIIAAAggggAACCCCAAAIIEH5zDyCAAAIIIIAAAggggAACCCCAAAIIIIAAAgj4ToDw23ctZUIIIIAAAggggAACCCCAAAIIIIAAAggggAAChN/cAwgggAACCCCAAAIIIIAAAggggAACCCCAAAK+EyD89l1LmRACCCCAAAIIIIAAAggggAACCCCAAAIIIIAA4Tf3AAIIIIAAAggggAACCCCAAAIIIIAAAggggIDvBAi/fddSJoQAAggggAACCCCAAAIIIIAAAggggAACCCBA+M09gAACCCCAAAIIIIAAAggggAACCCCAAAIIIOA7AcJv37WUCSGAAAIIIIAAAggggAACCCCAAAIIIIAAAggQfnMPIIAAAggggAACCCCAAAIIIIAAAggggAACCPhOgPDbdy1lQggggAACCCCAAAIIIIAAAggggAACCCCAAAKE39wDCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAr4TIPz2XUuZEAIIIIAAAggggAACCCCAAAIIIIAAAggggADhN/cAAggggAACCCCAAAIIIIAAAggggAACCCCAgO8ECL9911ImhAACCCCAAAIIIIAAAggggAACCCCAAAIIIED4zT2AAAIIIIAAAggggAACCCCAAAIIIIAAAggg4DsBwm/ftZQJIYAAAggggAACCCCAAAIIIIAAAggggAACCBB+cw8ggAACCCCAAAIIIIAAAggggAACCCCAAAII+E6A8Nt3LWVCCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAoTf3AMIIIAAAggggAACCCCAAAIIIIAAAggggAACvhMg/PZdS5kQAggggAACCCCAAAIIIIAAAggggAACCCCAAOE39wACCCCAAAIIIIAAAggggAACCCCAAAIIIICA7wQIv33XUiaEAAIIIIAAAggggAACCCCAAAIIIIAAAgggQPjNPYAAAggggAACCCCAAAIIIIAAAggggAACCCDgOwHCb9+1lAkhgAACCCCAAAIIIIAAAggggAACCCCAAAIIEH5zDyCAAAIIIIAAAggggAACCCCAAAIIIIAAAgj4ToDw23ctZUIIIIAAAggggAACCCCAAAIIIIAAAggggAAChN/cAwgggAACCCCAAAIIIIAAAggggAACCCCAAAK+EyD89l1LmRACCCCAAAIIIIAAAggggAACCCCAAAIIIIAA4bcN98DGT2I2jMIQCCCQEehZUaKGeFLR+iQoCCBgo8AuPUpVG4urPp6ycVSGQgCBPr3K9PHWeiVTaTAQQMBGgX67hMX/17IRlKEQkFQQknbvFdYHm8gxuCEQ2FkB8+cUh/sFCL9t6BH/QWYDIkMg0EyA8JvbAQFnBAi/nXFlVAQIv7kHEHBGgPDbGVdGDbYA4Xew+8/s7RUg/LbX06nRCL9tkCX8tgGRIRAg/OYeQMBxAcJvx4m5QEAFCL8D2nim7bgA4bfjxFwggAKE3wFsOlN2TIDw2zFaWwcm/LaBk/DbBkSGQIDwm3sAAccFCL8dJ+YCARUg/A5o45m24wKE344Tc4EAChB+B7DpTNkxAcJvx2htHZjw2wZOwm8bEBkCAcJv7gEEHBcg/HacmAsEVIDwO6CNZ9qOCxB+O07MBQIoQPgdwKYzZccECL8do7V1YMJvGzgJv21AZAgECL+5BxBwXIDw23FiLhBQAcLvgDaeaTsuQPjtODEXCKAA4XcAm86UHRMg/HaM1taBCb9t4CT8tgGRIRAg/OYeQMBxAcJvx4m5QEAFCL8D2nim7bgA4bfjxFwggAKE3wFsOlN2TIDw2zFaWwcm/LaBk/DbBkSGQIDwm3sAAccFCL8dJ+YCARUg/A5o45m24wKE344Tc4EAChB+B7DpTNkxAbeG3yteWq1r5y7STXMuUK+qSsfm7/TAm7fW6JyZ1+vCs4/X2P1HdftyhN/dpvv0jYTfNiAyhCsECjZvVsHmT9T4ddP2r58otMn8+yYVmK9btn/dfm6opqbN2lO9eysdKVc6ErG+psorlA6HG79X3vhPKhKRzPe3n5M5P7J7b8VLw4qWVShdWalUjx5Kl1e4wogiEPCyAOG3l7tH7W4WIPx2c3eozcsChN9e7h61u1WA8NutnaEuLwrkK/yO1TXo0mvm6+F/Lm/BdvlF0zR5/KHKZ/i9eMkyLXpoqS3BO+G3i34qCL9d1AxKaVMgtK1WxW+8rsJ31qnwvfdU8N67Ktrwrgr+95EKPv64MdDevNn1eonhI5QYNETJQYOV3P41MWiwEoOHWiE5BwIIdCxA+M0dgoAzAoTfzrgyKgKE39wDCNgvQPhtvykjBlcgH+H3mvUbdfZF12r8kQfpR2cf34RvguJZV83TjHOnaNPmalZ+N7stWfltw88o4bcNiAxhi0DhRx+q6I3XVfTGahW++YaK31hl/brw/Y1ZjW+tsu7Vu/Gf3r2V6rn9a6/eSrf6dapXL6V67SLztfURikUVipp/tqlgW1QmfA/FYo1fo9us1wqsf49K5rzaxu8XmH83QX19TNparfTWrSqo3qqCj//Xaf2pnj2VHDxUiYGDlBzcGI43BuWDlBg6vNP3cwICQRAg/A5Cl5ljPgQIv/OhzjWDIED4HYQuM8dcCxB+51qc6/lZINfhd2bFd9/de7cIvlsbZ1Z+m+1CZs2ep/c//MQ6ZcENM1tsH5IJ0jOvn3ni+KZxM6uuzzjh67rlrkf0yqq11hhmdfln9hluBfBtjWtWfi9/YaUumzFN4bIS6z2trzNm9FBrZbg5zLYmmbHNr9uqgW1PXPBTRPjtgiYErISidWtV9OYbjUH36pUqeutNFb2+UgXV1e1KWKumh41QcsBAJQfsqWT//kr2H9AUYKd22901ij0rStQQTypan2yqyWyvUvT2GhWtX6fCd9arcN1a62vRurdVtObNTmtP9utvrRCP7zNGyb1GKj56b8X32Y8V453KcYKfBAi//dRN5uImAcJvN3WDWvwkQPjtp24yF7cIEH67pRPU4QeBXIffmRB59qzpHe6BbcLvqefP0dFHjWsKoVtvR2LGunj2PF05a7qGDeqn1sF6Jvw2fcrsHZ4ZNxNem/3EW2+x0jr8bqvmR5c+r+FDBqh3z0rdcuffdM7px1hBeeaax0843Nq+hW1PXPRTQvjtomb4rBSzgrpkxfMqWbFcRStfU9Fbb6j4tVfanWW6LKzEXiOVGLGXEqP2VmKvUYoP30uJkd1/MEA+SNsKvzuro3DjBmtbl6J161T47noVWkH5ehWuf1uFH7zf/ocCQ4er4XMHKP75A9Uwdpzin/lsZ5fidQQ8K0D47dnWUbjLBQi/Xd4gyvOsAOG3Z1tH4S4WIPx2cXMozXMCuQ6/TdBsVnLPvfpCK7Bu72hrz+/WYfd1cxdp8MC+VsicOZq/z3yv9cMm2wqjW3+vdfhtrmOO5lu0dNRo8/51735gnU/47aIfCcJvFzXD46WYLT9Kn/yXSpc+rpLnn1Xxyy+1OSOzJUlixCjF9xpprWJOjByt+IiRSg7c0+MCjeV3J/zubOLWSvm311gfHlgfJLy+SsUrX93hbenSMjV87vOKH/B5xceOU/1BX1Bq1906G57XEfCEAOG3J9pEkR4UIPz2YNMo2RMChN+eaBNFekyA8NtjDaNcVwt4Nfzu12fXNh+YabBbb0nSfMuRrobfZjzzYM7jJhze7kr11luimPdkVqzX1dfvEMB354Zgz+/uqLV6D+G3DYgBHqLkhRUqXfpPlT7+mEqee2YHCbM9R8O4Lymx736KjxylxMi929xn20+EToTf7fkYc7O6vnjFcpX8+3kVfvjBDqeaDxUaDjhQ8QMPUsPnD7LCcQ4EvChA+O3FrlGzFwQIv73QJWr0ogDhtxe7Rs1uFyD8dnuHqM9LArkOv7uy7cm1cxc1bVdiTJuv/M6E3+MO2LvFyu/m9tkE3eb8jlZ+dxZ+m1Xel1w9v8Ve5M1XjhN+u+ingfDbRc3wQCnmoZRljzyskif+qdJ//XOHfbrj++6n+sOPVMOXDrW24TAPcgzakcvwu7Vt4XvvquTFFSp+/jlruxnz4URbR8PYgxq3Sfn8WDUc+AUl++4RtDYxXw8KEH57sGmU7AkBwm9PtIkiPShA+O3BplGy6wUIv13fIgr0kECuw+/OHniZ2Ut70+ZqdRR+my1TOtuOxI7w2+zj3dF12tp6hfDbpT8AhN8ubYxLygrV16nkmadU+tjfVbrsCRWveq1FZeZBlPUm6D7iy9bXVK9eLqk8f2XkM/xua9YlK55Tyb+fU7HZf/3fz8nsL976SO7Rz1oVHh97oBoO+qIaDhibP0CujEA7AoTf3BoIOCNA+O2MK6MiQPjNPYCA/QKE3/abMmJwBXIdfhvpzOrv8Uce1GIfbRMkL3n8OWs/8GzC78zDKy+/aFrT6m8TeGceQNnWqutsAvHWe35nrrPghplNW59kQvoHHnlKH3y0qemhnJm5fW7fEdb3WPntop8twm8XNcMlpZh9pUvNyu6lj1n7dzc/UrvsqvrDjlT9UV9V/SGHKdmvv0uqdk8Zbgu/W8uYB2iaQLz4+caV4WZ/9taHWbFff8RXVPeNo62vfKjhnvsryJUQfge5+8zdSQHCbyd1GTvIAoTfQe4+c3dKgPDbKVnGDaJAPsJv45xZAf7wP5c3sWf26u5VValsHnjZPEh//8NPmsbJhOHZBN3mTZ098NKckwnAMxfJ1FpWWtpi73Gz1/d+ew/TyyvXEH677QeK8NttHcl9PQXV1dYWJqV//5vK/v6ICj75uKmIdGWl6r94iBoOP1L1hxyh+KjRuS/QY1d0e/jdFmfJi/9uXB2+/FmVPrVUBZs2tTjN7Nted/QExb45ScmBgzzWEcr1iwDht186yTzcJkD47baOUI9fBAi//dJJ5uEmAcJvN3WDWrwukK/w2+tuua6fB17aIE74bQOiB4cwq3/Dd9+psr/9dYeVv2ZFd8NhR1kru9n+ouvN9WL43XqWZmW49WHI4/9Q8X//0+Jls6973dfGq+6bkxQf85muA/EOBLopQPjdTTjehkAnAoTf3CIIOCNA+O2MK6MGW4DwO9j9Z/b2ChB+2+vp1GiE3zbIEn7bgOiRIcz+3eEH71d40UJrW5PmR92Xv6a6Y76t2IRJSpdXeGRG7izTD+F3c9nCDz9Q2UP3q/Qfj6jsn39vgZ7sP0Cxbx6j+vETVP+lQ9zZEKryjQDht29ayURcJkD47bKGUI5vBAi/fdNKJuIiAcJvFzWDUjwvQPjtjRYSftvQJ8JvGxBdPkTJ8qcVvusORe5dpFB0W1O18f0/p+ipZyg2cTJ7OtvYQ7+F381pQrGotQ986aNLFF7yYIvtUcw+4XXfmKA6E4Qf+WWlS8tsVGUoBCTCb+4CBJwRIPx2xpVRESD85h5AwH4Bwm/7TRkxuAKE397oPeG3DX0i/LYB0YVDFG54T5E7b1f4L7eraN3apgqTe/RT9Lgpip16hhJDhrmwcu+X5Ofwu0V30mmV/Pt5lT76sMoeeVjFq1c1vZyOlKvuyC+rfvxE1X31GzLBOAcCOytA+L2zgrwfgbYFCL+5MxBwRoDw2xlXRg22AOF3sPvP7O0VIPy219Op0Qi/bZAl/LYB0SVDmFXd4QfvswLv0meelNJpq7J0WdhajRs98RTVH3akVFDgkor9WUZgwu9W7TMfuJjV4KVmH/lnn1YoHm88o6hI9V88WHUmCB8/Qcl+/f3ZeGbluADht+PEXCCgAoTfAW0803ZcgPDbcWIuEEABwu8ANp0pOyZA+O0Yra0DE37bwEn4bQNinoco/dcTity+QGVLHpLZ1ztzNBz4BUVPnarYxGPZxzuHPQpq+N2cOFRTo7Klj6nswftV9tgjMr/OHGa7HROCm+12EsNH5LAzXMrrAoTfXu8g9btVgPDbrZ2hLq8LEH57vYPU70YBwm83doWavCpA+O2NzhF+29Anwm8bEPMwRNGaNxW+c6Eid9wm80DCzJEcMFDRE05W9JSpSg7cMw+VcUnC7x3vgdKnn1Tp35eobMlfVfT2mk+D8L33Vd3k4xT91gncr/zodCpA+N0pEScg0C0Bwu9usfEmBDoVIPzulIgTEOiyAOFYijEpAAAgAElEQVR3l8l4AwLtChB+e+PmIPy2oU+E3zYg5miIgupqhRfdofC9d6lkxXNNV02HI4pNPk6xb09R/SGH5agaLtOeAOF3x/dG8aqVKvv731S2eJGKX3ul6eSGsQcpdtyJih5/ktIVFdxgCOwgQPjNTYGAMwKE3864MioChN/cAwjYL0D4bb8pIwZXgPDbG70n/LahT4TfNiA6PETps08pcvMfFb7/nhZXqj/4UMVOPE2xScda+3pzuEOA8Dv7PhS9+YbC99yp8OK7W6wIjx17nKInnab6I47KfjDO9L0A4bfvW8wE8yRA+J0neC7rewHCb9+3mAnmQYDwOw/oXNK3AoTf3mgt4bcNfSL8tgHRoSHCD9yrit9cp+L//qfpCsk9B1lbmphgMNl3D4euzLA7I0D43T29kuVPK/KX2xW+d5FCdTFrkNQuuyp60qmKnjZNiSHDujcw7/KNAOG3b1rJRFwmQPjtsoZQjm8ECL9900om4iIBwm8XNYNSPC9A+N2yhZu31uicmddrz36767IZ0xQuK3FFjwm/bWgD4bcNiDYPUf7nW1T+2+tVtPatppGjU05W9NQz1HDQF22+GsPZLUD4vXOioW21Cj+wWJHbb1XJ8882Ddbwuc8rdvxJ1j+pHj127iK825MChN+ebBtFe0CA8NsDTaJETwoQfnuybRTtcgHCb5c3iPI8JUD43bJdK15arbsfWqrq2qhmnDtFwwb1c0U/Cb9taAPhtw2INgxhVrqW3zZf5b+5ToUfvG+NmBg6XNFp0xU94RSlevWy4SoMkQsBwm/7lM0HQOHbFyhy1x0tHuwam/QtxU46VXVHfdW+izGS6wUIv13fIgr0qADht0cbR9muFyD8dn2LKNCDAoTfHmwaJbtWgPC7ZWuum7tIhxy0n5587mUNHthXk8cf6oreEX7b0AbCbxsQd2II8xDL8j/+XuV/vFEFmz6xRjKru2sv/D/VHfmVnRiZt+ZLgPDbGfmyxx5V+I4/K/zg4qYLJHfvo9gJJ1lbASWGjXDmwozqGgHCb9e0gkJ8JkD47bOGMh3XCBB+u6YVFOIjAcJvHzWTqeRdIG/h9yefSK+8kvv577KLNGZMm9c1W57M/s1CzfrByXrr7Q3WCnC3bH1C+G3DrUL4bQNiN4Yo+ORjlf/+BlXcMk+hmhprBLOKtXbGT9Tw+QO7MSJvcYsA4beznTAfEkUW3aHwwltVvGpl08Xin/msoieeqti3pyjVs6ezRTB6XgQIv/PCzkUDIED4HYAmM8W8CBB+54Wdi/pcgPDb5w1mejkVyFv4/dBD0sSJOZ2rdbEJE6QHH2zzumbLE7Pi+0dnH6/M3t8Xnn28xu4/Kvd1troi4bcNLSD8tgGxC0MUbtygil//ytrPOFRfJxUUKPbNY1Q7Y5bio/fpwkic6lYBwu/cdabkPy8ofMdtCi9epIKtW60Lp0tKVPfVbyh24qmq+/LXpMLC3BXElRwVIPx2lJfBAyxA+B3g5jN1RwUIvx3lZfCAChB+B7TxTNsRgbyF3888I118sSNz6nDQL3xBuuqqHU6J1TXo0mvm67gJhzeF3WYLFHOYMDzfB+H39g6YTyimnj/H+tWY0UN105wL1Kuq0vr14iXLdMnV861/P/qocTss2yf8zs1tXPTWm6q44WpF7rlLSiSULi5W7NsnqPZH/6fEkGG5KYKr5ESA8DsnzC0uEmqoV9lDDyiycIFKn/yXlE5br6d23U3R4060Hhab2Gtk7gvjirYKEH7byslgCDQJEH5zMyDgjADhtzOujBpsAcLvYPef2dsrkLfw295p7PRoa9Zv1NkXXav3P2zcijhztM5Xd/pC3RyA8FuSCb6vnbuoReCd8Wz9WlufXBB+d/Puy/Jtxa/8VxXX/VLhhx+UUimly8KKnnK6an9woZL9+mc5Cqd5SYDwO7/dKtzwniLmIZl33q7Cd99pKiY+5jPWQzKj3z6RB8jmt0Xdvjrhd7fpeCMCHQoQfnODIOCMAOG3M66MGmwBwu9g95/Z2ytA+N3oaRYNL39hZYvFwm2tBrdXP/vRAh9+m31oZl01TzPOnaJhg/rtIGfC7uZPKG0rKCf8zv6G68qZpc8+pYrrrlbpE49Zb0tXVqp22lnadu4Pldpl164MxbkeEyD8dknD0mmVLluqyMJbVfbwg43bDJmfxeJi1X/l69b+4HVf+bpUVOSSgimjMwHC786EeB2B7gkQfnfPjXch0JkA4XdnQryOQNcFCL+7bsY7EGhPgPBbyoTc4w7YW5PHH9qCyoTi6979IO9bnwQ+/G5raf6ZJ463GtNWA835F8+epytnTW8Kyz/Y1BgIcdgjUPa3v6r8+qtV/MIKa0Dz4L3oeedr23fOUbpHD3suwiiuFqgqL1ZDIqlYfcrVdQapuFB1tcJ3/6XxIZn//U/T1M3PZ+yk0xQ7ZaoSI0cHicSTc+1VWaJtdXE1xBu3teFAAIGdFAilpXRIu/Us1SfVDUql+NnaSVHejkALgb69y8T/13LvTZFWWiGF3FsglbUpEApJu/Us00ebyTG4RRDYWQHz5xSH+wUCH36bldx3P7S0aWl+5omkx084XN84ctwOG7a3FX6ntu+N6/52u7vC0F13SVdeqdCrrzYW2q+f0jNmKH3WWVI47O7iqc5WgZD5LzKTJ5j/4XCfwMqVCs2dq9Dtt0ubNzfVlz74YOmcc5Q+8UT31UxFlkBBKGRt587PFjcEAvYIJBJpFRWFtv9s8ZNljyqjIPCpgPlzi/+v5d47IpFMq6iQ8Nu9HWq7MtMx8/+3+NnyWueo140C5s8pDvcLEH63Cr9NyzJ71cz8/ima89vb1XzpflvhN9ue7NyNHrnjNlVcf42K3l5jDZQYNES1P7pI0ZNP37mBebdnBdj2xDutCz90n4pf/LfC99/btD94qvcuip54iqJnTFdi8FDvTCYAlbLtSQCazBTzIsC2J3lh56IBEGDbkwA0mSnmXIBtT3JOzgV9LMC2J95obuDDbxNmX3PjnZr9k+nqVVVpda35njTs+e3cjRy5c6Eqr75Che+sbwy99xqpmgtnKfat4527KCN7QoDw2xNt2qHIsn88osht82W2Lsoc9YcdoejU6YpNOMabk/JZ1YTfPmso03GNAOG3a1pBIT4TIPz2WUOZjisECL9d0QaK8IkA4bc3Ghn48Duzr3ff3Xtb+3xntj258OzjNXb/UWr9gEsThpvDnJs5WPndtZs9fP89qrz6ShW98XpT6F39sytU9/WjuzYQZ/tWgPDb260t3LhBkVtvVvlt81Xwv4+syST79FX0lKnWavBk3z28PUEPV0/47eHmUbqrBQi/Xd0eivOwAOG3h5tH6a4VIPx2bWsozIMChN/eaFrgw2/Tpkzg/cqqtVbXLr9oWosnlJqV4JdcPd967eijxjXtD0743bWbvOTfz6vHxTNUsv1Blokhw1RzyWWKTZzctYE42/cChN/+abHZFiWy4E8q/dcTTZOq+9p4Rc88W3VHfsU/E/XITAi/PdIoyvScAOG351pGwR4RIPz2SKMo01MChN+eahfFulyA8NvlDdpeXl7D79ahc3OyMaOH6qY5FzRtReJmTlZ+d9ydwnfXq8dll8is+DZHcvc+qv3Jpdp2ylQ3t5Xa8ihA+J1HfIcubfb0j9zyJ0XuuFUFW7Y0/l6w5yBFT/+Otp18mlK77ubQlRm2uQDhN/cDAs4IEH4748qoCBB+cw8gYL8A4bf9powYXAHCb2/0Pq/hd1tbiHiDrWWVhN9tdy1UW6vKX12l8j/epFBDvdLhiGrP+6Fqz/+x0mVhL7aamnMkQPidI+g8XMb8XhB+4D5FbpmnkueftSpIl5SobvxERadNV/0XD8lDVcG5JOF3cHrNTHMrQPidW2+uFhwBwu/g9JqZ5k6A8Dt31lzJ/wKE397ocd7Cb7Pqe9ZV8zTj3CkaNqifN7TaqZLwuxVMMqnyW29W5ZzLVbDpEykUUvT4E1X986uU2m13T/ea4nMjQPidG+d8X6Xo9dUqv/kPitz9F4VqaqxyzINvo1O/o+gJpyhVVZXvEn13fcJv37WUCblEgPDbJY2gDN8JEH77rqVMyAUChN8uaAIl+EaA8NsbrST8tqFPhN+fIpb94xH1uHRW08MsGw4Yq63X/lbxffezQZohgiJA+B2UTjfOMxTdpsi9i6zV4MUvv2R9z/ztkNix37YekNnwuc8HC8TB2RJ+O4jL0IEWIPwOdPuZvIMChN8O4jJ0YAUIvwPbeibugADhtwOoDgyZt/DbzMVsezJ4YN8WD5d0YI6OD0n4LSvsrvrxD1T6zJOWt/Uwy59fqdjREx335wL+EyD89l9Ps51RyYv/tkLw8H33KFQXs94WH/MZKwSPfvsEpSPl2Q7FeW0IEH5zWyDgjADhtzOujIoA4Tf3AAL2CxB+22/KiMEVIPz2Ru/zGn6vWb9RCxc/phnnTFG4rMQbYm1UGeTwu+CTj1V5+c9UfsdtUiqlVM+eqv3xLG0787tKFxd7tqcUnl8Bwu/8+rvh6gVbtyqyaKEVhJsP18yRrqxU9LgTrd9fEiNHuaFMz9VA+O25llGwRwQIvz3SKMr0nADht+daRsEeECD89kCTKNEzAoTf3mhV3sJvs+f3OTOv1yur1rYpNWb0UN005wL1qqp0vWQQw+9QfZ0qfv9rVfz6WoW21VpB97ap01X7fz+1AnAOBHZGgPB7Z/T8917zN0oit/xJZQ8/oFBDgzXBhrHjFJ12lmKTjlW6pNR/k3ZoRoTfDsEybOAFCL8DfwsA4JAA4bdDsAwbaAHC70C3n8nbLED4bTOoQ8PlLfx2aD55GTZQ4Xc6rcg9d6ny8ktUuHGD5V339aNVfcUvlRg8NC/+XNR/AoTf/uupHTMyD9At//MCRW69WYXvrLOGTPXureiUUxWdeqYSQ4fbcRlfj0H47ev2Mrk8ChB+5xGfS/tagPDb1+1lcnkSIPzOEzyX9aUA4bc32kr4bUOfghJ+l7ywQlUXfl/Fr75sqcVH7209zLLhwC/YoMgQCHwqQPjN3dChQDqtsices7ZEKfv736RkUgqFVH/o4YpOna7YN74pFRWB2IYA4Te3BQLOCBB+O+PKqAgQfnMPIGC/AOG3/aaMGFwBwm9v9D7v4feKl1Zr6vlzWmgtuGGmxu7vnf1c/R5+F767Xj0umanwXx+w+pTsu4dqfnqZoiecbAVOHAjYLUD4bbeof8cr/OB9RRb8SZHbF8j8u/V7VJ++ip58uqKnn6lk/wH+nXw3Zkb43Q003oJAFgKE31kgcQoC3RAg/O4GGm9BoBMBwm9uEQTsEyD8ts/SyZHyGn6b4PvauYta7O1tHoJ59kXX6tzTJ2ny+EOdnLttY/s1/A7V1Kjyl1eofP5ca5/ddDii2u+dr9ofXqh0Wdg2PwZCoLUA4Tf3RJcFEgmFH12iyIJ5Kl36uJROS4WFqvvy1xQ9Y7rqjvoqH9ZJIvzu8p3FGxDISoDwOysmTkKgywKE310m4w0IdCpA+N0pEScgkLUA4XfWVHk9MW/hd6yuQZdeM1/HTTh8h1XeJhS/+6GlumzGNIXLSvIKlM3FfRd+JxKquOWPqrj6ShVs3mwFRtETTlL1pVcqtdvu2ZBwDgI7JUD4vVN8gX9z0fq3rS1RInf8WWafcHMk9xyk6GnTtO2UqUrtultgjQi/A9t6Ju6wAOG3w8AMH1gBwu/Atp6JOyhA+O0gLkMHToDw2xstz1v4vXlrjWZdNU8zzp2iYYP6tdAyq7+vufFOzf7JdPWqqnS9pJ/C77K//VU9fn6xita8abmb/bzNvt5mf28OBHIlQPidK2l/X8f8jZWyh+5X+S3zVLL8aWuy6eJi1Y2faK0Grz/YG3+7yM4uEX7bqclYCHwqQPjN3YCAMwKE3864MmqwBQi/g91/Zm+vAOG3vZ5OjZa38JuV3061tHvjFq9epR7/d4FKn15mDZAYNEQ1l16u2MTJ3RuQdyGwEwKE3zuBx1vbFCh6fbW1N3j5X25TqLa28fe5IcMUnTZd0ZNOV6qqKhByhN+BaDOTzIMA4Xce0LlkIAQIvwPRZiaZYwHC7xyDczlfCxB+e6O9eQu/Dc/iJcu06KGl7Pmdx3vFbAlQeeXPVX7rzVYVqZ49VXvRxao967w8VsWlgy5A+B30O8C5+YfqYorcc5e1LUrxf//TdKHolJMVnTpdDZ8/0LmLu2Bkwm8XNIESfClA+O3LtjIpFwgQfrugCZTgOwHCb9+1lAnlUYDwO4/4Xbh0XsNvU6fZ33vq+XNalLzghpk77APehTnl/FSvbntS8YffquKaq1Swdatltu3Ms1Uz61IrAOdAIJ8ChN/51A/OtUv+84LCt96syN13KlRfZ008vs8YRad+x3rOQTpS7jsMwm/ftZQJuUSA8NsljaAM3wkQfvuupUzIBQKE3y5oAiX4RoDw2xutzHv47Q2mjqv0Wvhd/NKL6nXumSp643VrYta+3tf9TvFRo/3QDubgAwHCbx800UNTCNXUKHLXQkUWzJPZAsoc6YoKRb89RdEzz1Z89D4emk3HpRJ++6aVTMRlAoTfLmsI5fhGgPDbN61kIi4SIPx2UTMoxfMChN/eaCHhtw198kr4bfa57XHFpSr/003WrJP9+qv6yqsVm3CsDQoMgYB9AoTf9lkyUtcESp99ytobPHzvoqY3NhwwVtEzzpLZGsXrB+G31ztI/W4VIPx2a2eoy+sChN9e7yD1u1GA8NuNXaEmrwoQfnujc4TfNvTJC+F32T8eUdWF31fhxg1SQYG2TTtL1T/9hbW6kQMBtwkQfrutI8GrxzwPIbLwNkVuvVlF69ZaAKlevRSdcqqiZ3xHiaHDPYlC+O3JtlG0BwQIvz3QJEr0pADhtyfbRtEuFyD8dnmDKM9TAoTf3mhXzsPvzVtrdM7M63XGCV/XLXc9oldWNYYKrY8xo4e2eBCmmzndHH4XfPw/VV10gcIPLrYIE3uN1OYbb1Z8/8+5mZTaAi5A+B3wG8BN00+nVbr0cWtLlPCjS6REQgqFVH/woYqeMV1135igdHGxmyrusBbCb8+0ikI9JkD47bGGUa5nBAi/PdMqCvWQAOG3h5pFqa4XIPx2fYusAnMefmdYTAg+66p5mnHuFA0b1K+FlnkI5t0PLdVlM6YpXFbiekm3ht/lty9Q5aWzrAdapkvLVPvjmar93gWeCmpc33wKdESA8NsRVgbdSYHCD95X5M+3WP9Yf4vGbB+1ex9FTz7dekhmsv+AnbyC828n/HbemCsEU4DwO5h9Z9bOCxB+O2/MFYInQPgdvJ4zY+cECL+ds7VzZFeG32vWb9Q1N96p2T+Zrl5VlXbO15Gx3BZ+F61/Wz3PO0sly5+25tsw7kva8vs/KjFoiCPzZ1AE7BYg/LZblPFsFUgmVfb3vylyyzyVPfGYlE5b20nVHfVVRaedZX01v3bjQfjtxq5Qkx8ECL/90EXm4EYBwm83doWavC5A+O31DlK/mwQIv93UjfZrcWX4vXjJMi1/YSUrv7t4D4XicVX87npV/GqOQvV1SvXsqepfzFH0pNO6OBKnI5BfAcLv/Ppz9ewFCt99R5Fb/qjyhbep4JOPrTcmB+6p6GnTtO3UM5TadbfsB8vBmYTfOUDmEoEUIPwOZNuZdA4ECL9zgMwlAidA+B24ljNhBwUIvx3EtXHonIffZlX32Rddq/c//KTdaezRZxfNvfrCHbZDsXHetg7lhpXfxS+9qF7nnqmiN1635hY79jht/eV1SvXexda5MhgCuRAg/M6FMtewUyDU0KCyhx9QZP48lT77lDW02Qvc7Alu9gavP+QwOy/X7bEIv7tNxxsR6FCA8JsbBAFnBAi/nXFl1GALEH4Hu//M3l4Bwm97PZ0aLefhd2YiHe357dRknRo3n+F3KLpNPX7xM5XPnyulUkr2668tv/6D6o84yqnpMi4CjgsQfjtOzAUcFCh6601Fbp6ryKKF1jMXzJEYNkLRqWcqOuVUpXr1cvDqHQ9N+J03ei7scwHCb583mOnlTYDwO2/0XNjHAoTfPm4uU8u5AOF3zsm7dcG8hd/dqtalb8pX+F32j0dUdeH3Gx+8Vlio2unnquYnP1M6Uu5SKcpCIDsBwu/snDjL3QKhupjC992jyPw/quQ/L1jFmocPx46ZrOjU6WoYe1DOJ0D4nXNyLhgQAcLvgDSaaeZcgPA75+RcMAAChN8BaDJTzJkA4XfOqHfqQoTfO8XX+OZch98FH/9PVRddoPCDi63rx0fvoy033az4vvvZMBuGQCD/AoTf+e8BFdgrUPzqy42rwe9dJPM3dqzfu/fe19oSJXrciUpXVNh7wXZGI/zOCTMXCaAA4XcAm86UcyJA+J0TZi4SMAHC74A1nOk6KkD47SivbYPnNfzuaP/vMaOH6qY5F6hXVaVtk3VqoFyG3+W3L1DlpbOsv0qfLgurZuYlqj3n+9bKbw4E/CJA+O2XTjKP1gKh2lpF7lqoyII/qXjVa9bL6fIKRb99gqLf+a71YaaTB+G3k7qMHWQBwu8gd5+5OylA+O2kLmMHVYDwO6idZ95OCBB+O6Fq/5h5C79jdQ269Jr5GnfA3vrMPsO1cPFjmnHOFIXLSnTd3EU65KD9NHb/UfbP2IERcxF+F61/Wz3PO0sly5+2ZlB/8KHa8ts/KjlwTwdmxJAI5FeA8Du//lw9NwIlzz+ryC3zFH7gPoUa6q2LNhwwVtEzzlLs2G9ZW6TYfRB+2y3KeAg0ChB+cycg4IwA4bczrowabAHC72D3n9nbK0D4ba+nU6PlLfxu/sBLM7lrbrxTs38y3VrpveKl1br7oaW6bMY0Kwx3++Fk+B2Kx1Xxu+tV8as5CtXXKdV7F1VfcbWix5/odhbqQ6DbAoTf3abjjR4UKNi8WZE7brNWgxe9vcaaQapnT0VPOEXRM89SYuhw22ZF+G0bJQMh0EKA8JsbAgFnBAi/nXFl1GALEH4Hu//M3l4Bwm97PZ0azRXhd++elZr9m4Wa9YOTrfDbbIfSPAx3avJ2jetU+F380ovqde6ZKnrjdatUE3hXX3WtFYpwIOBnAcJvP3eXubUrkE6rdNlSRRbMU/hvf5USCetU8zd9zN7gdeMnKl1cvFOAhN87xcebEWhXgPCbmwMBZwQIv51xZdRgCxB+B7v/zN5eAcJvez2dGi1v4XfzbU8mjz/U2upk8MC+Mv++eMkyLX9hZWBXfpuHofX4xc9U/qebrL4nBwzUll//QfWHHeHUfcC4CLhKgPDbVe2gmDwIFH74gSJ/vkWR2+arcOOGxj8L+g9oDMGP/Iri++3fraoIv7vFxpsQ6FSA8LtTIk5AoFsChN/dYuNNCHQoQPjNDYKAfQKE3/ZZOjlS3sLv1pMy26CcM/N6vbJqrfbos4vmXn2hhg3q5+TcbRvbzpXfZY89qqoLzlPh+xut+mp+8CPVXnSx9XBLDgSCIkD4HZROM89sBMoeedjaEsX8+ZA5EsNHKPatE1R73g+VjpRnM4x1DuF31lSciECXBAi/u8TFyQhkLUD4nTUVJyKQtQDhd9ZUnIhApwKE350SueIE14TfrtDoZhF2hN8F//tIVTMvVPiBe60q4vt/Tlt+fZPi+4zpZlW8DQHvChB+e7d3VO6cgPlQNLz4boUXL1Lxf//TdKH4vvupbuKxik2arMSwER0WQPjtXH8YOdgChN/B7j+zd06A8Ns5W0YOrgDhd3B7z8ztFyD8tt/UiRHzFn43f+ClV1Z4t9eAnQ2/y/98i3r8bKZCNTVKhyOqueQy1Z51nhP9ZkwEPCFA+O2JNlFkHgXMgzHD992jsvvvVfHKVz8NwkeNVt2xx6nu6EmKjxq9Q4WE33lsGpf2tQDht6/by+TyKED4nUd8Lu1bAcJv37aWieVBgPA7D+jduCThdzfQWr+lu+F30fq31fO8s1Sy/GlryPojjrL29k72629DVQyBgHcFCL+92zsqz71A0bq1Ct+/WGUPLlbxyy81FZDYa6RiEyerbtJkxUfvY32f8Dv3/eGKwRAg/A5Gn5ll7gUIv3NvzhX9L0D47f8eM8PcCRB+5856Z66Ut/DbFG0ecnnIQftp7P6jdmYOeX9vd8LvyhuuUeUVl1q1p3bdTVtn/0qxY4/L+1woAAE3CBB+u6EL1OBFgcJ31ily/2KV/vV+lbz470+D8GEjFJt4rMpOPEHVI/dRfTzlxelRMwKuFSD87l5rQnUxFVRXK1RdrYKaaoVqa1WwZbNCNdUqqN5q/Tq0eXPja5nvmXOtf7Zaf2vSjGEO87cn08VFUlGx0sXFUnGx0kXma9H2r9u/X1ikdEmxVNTs+9Z55vVm798+TuZ7ipQrVVmpdI8qpXr0UKpXL6X69FVi0JDuTZ53ZSVA+J0VEych0CUBwu8ucXEyAh0KEH574wbJa/i9Zv1GLVz8mGacM0XhshJviLVRZVfC75IXVqjqh99V8epV1kjRk05T9S/mKNWzp2fnT+EI2C1A+G23KOMFUaDw3XcUfmCxyh5+QCUrnmsiSA4Zquikyao/epIaPntAEGmYMwK2CwQx/A5tq/00uDZBdCymgs2fNH3PCqu3bFGoZmvj92prVLB1a2PQXb1VBZs+sb0P+Row1bu3kn32UHKPfkrt0e/TrwMGKrl7n8Zf77pbvsrz9HUJvz3dPop3qQDht0sbQ1meFCD89kbb8hZ+mz2/z5l5vV5ZtbZNqTGjh+qmOReoV1Wl6yWzCb/NypUeV1yq8j/dZM0nMXS4tt7we9V/8RDXz48CEci1AOF3rsW5nt8FCjduUPjBxapY8oAKnnnm0yB8z8GKTTpWdd88Rg0HjPU7A/NDwDEBv4bf5vkChevWqWjtWypcv05Fb74us21f0Ruv22aZ6r2LUj2qlDarqauqlK6otFZWpyurrMUh6coeja9Zq613kQoLbbl2qL5OijNrAQQAACAASURBVMcViselZEKheEJKbP+1+X4iLsXN9+ON399Wuz24NyvVt6pg6xYVbnhP5mHE2R5JE4b3G6Bkn75K9t1DqX79G8Pyfv0av9d/gNJl4WyHC8R5hN+BaDOTzLEA4XeOwbmcrwUIv73R3ryF397gya7KzsLvssceVdUF5zX9x3HNBRep5uKfZzc4ZyEQQAHC7wA2nSnnRMDs+R19+x0VLr5XZQ/c1/TMCXPx5MA9ra1R6r45SQ1jx+WkHi6CgF8EvBp+mwC4aO1amS2Tit58Q4Xr1jaG22vWWN/r7DCrmZsH1ykTVFeaoLpK6V69lK6oaPx3872e23/dLMz2S9Bb+MH7Kvjgfeu/9c0/BdbXDSr88AOZDx/Nr83q92wOE/hbAbkJyvsPUHLQEMX3GaPEqNHWCvKgHYTfQes4882FAOF3LpS5RlAECL+90em8hd9m5fesq+ZpxrlTNGxQy/+QW/HSat390FJdNmOaJ7ZDaS/8LvjfR6qaeaHCD9xr3Q0NYw/S1utvVHzUaG/cHVSJQJ4ECL/zBM9lfS/Q+oGXBR//T+G/PqCy++9R6VPLmuZvHrxsPSxzwiQ1HPRF37swQQR2VsDN4bfZF7v4rUywvU6Fa95sXM1tVnV/8H6HUzerkxNDhik5eIiSw4YrYb6aMHb4XkpXuv9vZ+5sX+18v9mb3FopbgXijQF5wcYNKvxgowo/+KAxLH/v3Q4vaVa/x/fZT/F991Nyr5GK77Ov4nuPsT5k8OtB+O3XzjKvfAoQfudTn2v7TYDw2xsddWX4bfYCv+bGOzX7J9M9u+1JZOGt6vGzmdbehuY/SKt/doW2TTvLG3cFVSKQZwHC7zw3gMv7VqB1+N18ogWbNyv80H2NQfiypU0vWXvZ7t5H9V/+muqP+prqDznMtz5MDIHuCuQ7/LYCVbNq++21Knx7+1fr12tkfrY7OpJ7DlZi2DDrwY3JIcOUGD5C1veGDlW6tKy7JLxvJwTMXunFr76s4tdeVtHK11T0xmoVv/qKQrFom6OaDyzjo/ZWfN8xSo7ax1poE99v/52owD1vJfx2Ty+oxD8ChN/+6SUzyb8A4Xf+e5BNBa4MvxcvWablL6zM2crv6+Yu0s1/WdLC6/KLpmny+EOt75l6Lrl6vvXvRx81boe6mq/8NnsiVp1/nkqfedI6v+7rR2vrr35j7evHgQAC2QkQfmfnxFkIdFWgo/C7+VgmeClb8qDK7r9XZY//o8VlzAe6206dZq0ErT/4MCVGjupqGZyPgO8EchF+W1uTZFZsr1urwrVrVGT24X7rTZlVxR0dib1GWuF2YsRIJQeZYHu4koMHW6u6ObwjYG1J8/pqFb/2iopee7UxFF/5arsTSIzYS/G9Rikx5jONW6eMHGX13ksH4beXukWtXhEg/PZKp6jTCwKE317okpTz8Nus6j77omv1/oftP+F9jz67aO7VF+6wHYpTpCb8NsePzj5+h0uYLViunbuo6eGbbZ2bCb8rr52jytm/sMYwq+Sqr75BsW9OcqpsxkXAtwKE375tLRPLs0C24XfzMs2WCWVL/6mSf/5d4UceltkqpfmR2mVX1R9yuBoOO0L1Bx9KmJbnHnP5/AjYFX6bfbfNau2i7eF24bq3VWRC7rVvdTgxswVJYvBQ6+cvMXSYUoOHWF+t1dz9B+QHhavmTMCsEjeheNGq16xgvHj1ShW++06712/47AGKH3CgGg48SA1fPMTVi3QIv3N2G3GhAAkQfgeo2UzVcQHCb8eJbblAzsPvTNUd7flty8y6MEhH4bd5bfDAvk2rwFuH4eYyH/99map++F0Vr15lXXXb6Weq5tIrrafScyCAQNcFCL+7bsY7EMhGoDvhd+txi1/5r0qXPq6SJ5eqdPkzCkW3tTjFPJCt/vAj1XDokdYWKfzNp2w6wzleF8g2/DY/L2altlnBa/bdNvtvW6t5166x9oPu6DAfNCWGjVBi0CAlh42w9t5ODDGB91CZ1zgQaC4Q2lbbGISvWqnC7avFi1/9r7UlY+vDfEBinu/QcOjhqv/SIa76EJPwm/saAfsFCL/tN2XE4AoQfnuj93kLv93E03rbk8yWJ7G6Bl16zXyNO2DvpvDbrFy/ePY8XTlreuPK9PPPl379a2s65q8Wbvn1TWo48Atumh61IOA5AcJvz7WMgj0iYEf43XqqJSueU+lT/1LJsidU8vxzCtXXtTjF/BV7E4JbocrBhxHSeeReocyuCTQPvws++djae9vaf9uE3FbAvd5avd36b060FUJagfbgoY0Bt9maxITcw0coHSnvWlGcjUAbAmZ/eGs/8eXPqGT5Myp99qkdw/CBe6r+0CPUcPBhef8Qk/Cb2xgB+wUIv+03ZcTgChB+e6P3hN+t+pTZlmX2rOnad9RQK/w+bsLhGrt/456mO4TfoZD1/Yaf/Vz1My/2RtepEgGXC5SVFCqZTCmeTLu8UspLp9MKbf99EA33C0RKC9UQTymRcu5nq/DJZSpc+oSK/vWECp95egeU1KjRShxxpJKHHq7k4YcrXdXT/XBUiEArgYL33lPB22sUeusthdatU8nba5Res0ahN99QqLa2Q6/U8OFKDRmm1LBhSptV23vtZW1TkmL/fO6zPAkULn9WRcuWquCJJ1T03LNSXcsPMVMj9lLi0MOUOvRwJY44Quldd8tZpZXhItXEEjm7HhfqmkA6LfGfgV0zc8PZJsEoDxeplp8tN7SDGjwuYP6c4nC/QF7Db7P1yTkzr9crq9buIDVm9NCmfbZzzZjZ6uQbR47rfOX3FVdo28TJSg0fkesyuR4CvhWwwu9UWvFEyrdz9MvETITa+BEghxcEwqVFakgklczVB0uxmIqefUYFSx9XkQnFn39uB6bkZ/ZX8tDDlPj28Up+fqwXGKkxIAIFa95Swbp1Ml9Da9eo4M03VPD22yp4fXXHAmVlSu41UulBg6z/PkybYHv4iMaAe/DggOgxTS8LFD61/UPMpUtV+Gw7H2IefoSSh5l/DnP0Q8zKSLFqonEvc/q69nQoZB4i5us5+nFy5gOL8rJi1cb42fJjf5lTbgXMn1Mc7hfIa/jd0V7b+aRrvs93Nnt+Zx54mc+auTYCfhJg2xM/dZO5uEnAiW1PujI/swdt6bNPq+RfT6j06WUqfvmlHd4e33tfNYz7ohoOOVyxCcd0ZXjORaBLAmaLnqK1a1X4zjqZB00WrlvbuP/2mjXW9zo6Uj17ymzpY+27PXSYKvYZqS1991SDecBkn75dqoOTEXCzQCgWVclzy1Xy5BMqffJfKnnx3zv+vj3mM9a2VtaDj8d9SemKCtumxLYntlEyEAJNAmx7ws2AgH0CbHtin6WTI+Ut/HbLAy9NHUv+uVwnT/6K5dx6W5PWD7hsK7An/HbyFmXsIAoQfgex68w5FwL5Dr9bz9E8eM3aL9wEKs8+ZT2cra2j4bMHKDFylBKj91V8nzFKjBrNgzRzccP44BqhmhoVv/m6tfd20frGB0yaB02a/bfN3scdHSbEtgLuwYOt/bcTgwY3ht3D99rhoebZPvDSB6RMIeACZluf0uXbP8R86l8yD0FufTR87vONz3o45AjrAcg7cxB+74we70WgbQHCb+4MBHZeIPKXPyt8z50qXfr4zg/GCI4LBD78zjzU8uF/Lm/CXnDDzKY9vs03Fy9Zpkuunm+9fvRR43TZjGkKl5U0nU/47fh9ygUCJkD4HbCGM92cCbgt/G49cStUef5ZFT+/XCXm68svqWDLljZ9Ur16KT5qH8X33U/JUaMVH723FYyny+1bcZizxnChbguE6mIqfOcdFW7coKJ316vg3XcaHzRpVnGvW6uCzZs7HDu552Alhm5/wOTQYY0PmjQPnBwyVOmycNZ1EX5nTcWJPhMwP2Pmb/KUPLVMpU8+oaJW2wKZnyMThNcf9VXVfW28kgP37JIA4XeXuDgZgawECL+zYuIkBNoUKHzvXfU8/5xPQ2+2fvLEnZK38NvotN5SxBNibRRJ+O3VzlG3WwUIv93aGeryuoDbw++2fAs/+lBFq1ep+LWXVbTyNRW9bv79VZktK9o6kgMGKj5qb8X3HaOkCcdHjlJ8zGe83rpA1m8++Ch8f4MKPnhfRRveU8H7G1Ww4T0r6A5tfyCf2YLBBOAdHYm9RiqxfcV2ctBgJYaZ7UpM6D3cNlfCb9soGcjjAgUf/68xDF+2VGWP/0OF777TYkbmA8u6iccqNmmyEsM6f2YS4bfHbwjKd6UA4bcr20JRHhAo/9Mf1OPyn8ls5ZjabXdtufa36n3q8R6onBLzGn6bLUYWLn5MM86Z0mIltdfaQvjttY5Rr9sFCL/d3iHq86qAF8Pv9qyL3l5jrTA0W6VYofgbq1W86rV2W2P2Ek+MGKm68ROU6ruHEgMHyqz65ciPgPlQo2DjBivILjSh9vsbVGiC7Q/eV+GGxu+ZvYazOcxK0sTAQUr266fkwEFKmWB7iFnFPUTJ/gOyGWKnzyH83mlCBvCpgPl9ufSxv6vs0SUqWd7y4Znx0fuobtJk1X3zGMVHjW5TgPDbpzcG08qrAOF3Xvm5uAcFzHZ5Pb//XZU894xVffSk01R9+S+VqqoSe357o6F5C7/NXtvnzLxer6xa26bUmNFDddOcC9SrqtL1koTfrm8RBXpMgPDbYw2jXM8I+Cn8bg/d7D9bvHqVCl9fqeJXX7ECcROqtnWkwxGldt1VoVhMyT36Nf7Tdw+l+g+wvib79beC8mSfPZTq3dszfc53oYXvrt8eYr/XGGybQPvddxu/t/E9mb8ums1htkswPbBCbdOb/gOU6tuvsS/m1+Y1lzxckvA7m45yDgKyHiobfuA+lT24WMUvvdhEkhixl2ITJ1thuPmwMnMQfnPXIGC/AOG3/aaM6F+Byl9fq8rLL7EmaP4m4dbf/kH1XzykxZ9T/p29f2aWt/DbP4QS4befuslc3CBA+O2GLlCDHwWCEH631TfzVxPNCvHila9ZDz40Dzws/GCjCt9Z3+lDDzPjpUvLWgbi/fsr1afvp6G5CWMDspI8FN2msscetUJs80BJc5QuW2o9RLLgfx9l9aNjVspY4fUe24PsAQMbV+Rv/5rs209mX3evHITfXukUdbpJwHxQlgnCzRZGmcNshxKbeKzqJhyj3Y74Av9fy01NoxZfCBB++6KNTMJhAbOApuf3zlLxf/9jXanmBz9Szc+u2OGqrPx2uBE2DU/4bQMk4bcNiAyBQDMBwm9uBwScEQhq+N2ZpglxrVXJ72+0QvGQtUJ5w/aVyhsbt+CIbutsGOt1s0LcBLd1X/26VPzpw7ETI0cptVufrMZQQUjp4hKlS0qsMdIlxc3+vcR6TeZ7pWVtjlewaZO1F2HBtm0K1dZYK9sLqrcqZH69rVbmwaLW15qaludFoyqoqW46r8C8t7Y2u5qbnZXcvY/M3uvWhwNmlXa/AY0fElhh9x7Wa+3V3uWLueQNhN8uaQRleFbA/A2d8IOLVfbgfSpZ8dyn8xg2TDUTjlX90ZPU8NkDPDs/CkfATQKE327qBrW4UaDyqstUed0vrdLi+4zRlt/OVXy//dsslfDbjR3csaa8ht+xugZdes18PfzP5dqjzy6ae/WF6tdnV+t74w7YW5PHH+oJRcJvT7SJIj0kQPjtoWZRqqcECL+73y4TFDeG458G49Y+1WZbj48+bNyv+v2N3b+Ay9+Z6r2L0hUVSpWXK92jZ2PIb7aHGbCnkgMGKNl/YGPAPWCgy2fiTHmE3864MmowBczvpeEH7lXZg/er5PlnmxCSew5SbMIxqps4WQ0HjA0mDrNGwAYBwm8bEBnClwIlL6xQz+9NV9Gbb1jzq7nkctX88MIO50r47Y1bIa/h93VzF2nwwL76xpHjdM1Nd+rkyV/WsEH9tOKl1br7oaW6bMY0TzwIk/DbGzc7VXpHgPDbO72iUm8JEH4736/MCvKCTZ8oVL1VBWZ1dbVZTV0r1dcp1BBXKN4gNTRYX82vZX1tkOKNrzX/XuO5O77HGq+dw2wpki6vULq8XKmKSuur9WsrvDZfK5WurFQ6Um59z7yWMl/N9yORT99TUWk9yIejcwHC786NOAOB7gj0a9ii6lvvUOmD96v02ac+DcIHDFRswrGqm3iMGsaO687QvAeBwAoQfge29Uy8HYFQXUyVV/xcFX/4rXVGw7gvactvblJi6PBOzQi/OyVyxQl5C7/NAy9nXTVPM86dYq32bh5+r1m/UdfceKdm/2Q6D7x0xW1CEQjkVoDwO7feXC04AoTfwek1M82tAOF3br25WnAEmj/w0nyoGH7AbI2yWKXPPCUlkxaE2VIp9s1jVDfpWDUc+AUpFAoOEDNFoBsChN/dQOMtvhUoWf60ep1zpgrffcdaHFL986u07bRpWf9ZQvjtjVvDleE3K7+9cfNQJQJOCRB+OyXLuEEXIPwO+h3A/J0SIPx2SpZxgy7QPPxublGwebPC99+jsofusx64mzlSu+2u2KRvKXbMZGvlHgcCCOwoQPjNXYGArOfa9PjpRSq/fYHFUfe18dp63e+U7NO3SzyE313iytvJeQu/zYwXL1mm5S+s1KwfnKzfzr/P2vakd89KnTPzeh0/4XD2/M7bbcGFEcivAOF3fv25un8FCL/921tmll8Bwu/8+nN1/wq0F363CMK3bFHZkgcVXrxIpUsfb3rJPJfACsInH6f4/p/zLxIzQ6CLAoTfXQTjdN8JlP3jEVVdcJ71LJ/Urrtp6y+vs/686M5B+N0dtdy/J6/ht5muWeU99fw5LWa+4IaZGrv/qNxrdPOK7PndTTjehkA7AoTf3BoIOCNA+O2MK6MiQPjNPYCAMwLZhN8tgvDqaoUfXKyye+5U6VPLml5K7DVS0VOmKnrCyUrtsqszxTIqAh4RIPz2SKMo03aBgk2bVPXjH1h/TpgjesJJqr7yV0r17NntaxF+d5sup2/Me/id09k6dDHCb4dgGTawAoTfgW09E3dYgPDbYWCGD6wA4XdgW8/EHRboavjdvByzoi98390KL/qLil/5b9NLsQnHKHby6ar78tccrp7hEXCnAOG3O/tCVc4KmL8dVDXzRzIBeHLgntpyw02qP+yInb4o4fdOE+ZkgLyG39fNXaQPPtqky2ZMU7isxJpwrK5Bl14zX+MO2JttT3JyC3ARBNwnQPjtvp5QkT8ECL/90Udm4T4Bwm/39YSK/CGwM+F3c4Hi//5HkT/fosi9dylUU2O9lOy7h6JTTlHs1KlKDBriDzBmgUAWAoTfWSBxim8Eita+Za32tp4PUVCg2rPOU81Pf650WdiWORJ+28Lo+CB5C78zIfdxEw7fYYsTHnjpeN+5AAKuFiD8dnV7KM7DAoTfHm4epbtagPDb1e2hOA8L2BV+ZwhCdTFF7rlL4YW3qmTFc00y9V84WLFTpio26VjbAhEPs1O6zwUIv33eYKZnCZgHWlbO/oUq5v7O+nVixF7acuPNavjsAbYKEX7byunYYHkLvzdvrdGsq+ZpxrlTNGxQvxYTXLN+o6658U7N/sl09aqqdGzydg3Mtid2STIOAo0ChN/cCQg4I0D47YwroyJA+M09gIAzAnaH382rLHrrTUVuvVmRv9ymgi1brJfS5RXWHrCxE05WwwFjnZkUoyKQZwHC7zw3gMs7LlB++wJV/uKn1hYn5qiZ9TPVXDjTkesSfjvCavugeQu/Wfltey8ZEAHfCBB++6aVTMRlAoTfLmsI5fhGgPDbN61kIi4TcDL8zkw1FI+r7OEHFLltvkqf/JeUTlsvxUfvba0Gjx5/slK9erlMhnIQ6L4A4Xf37XinuwWKXl+tnt8/SyUv/tsqtGHcl7Tl9390dGsrwm933xNNf9an09v/dM9DvWZ7k1mz52nu1Rc2rf42q77PvuhanXv6JPb8zkNPuCQCbhAg/HZDF6jBjwKE337sKnNygwDhtxu6QA1+FMhF+N3crfCddYr8eYHMqsGC/33U9FJs4mTFTj5NdUd91Y/MzClgAoTfAWt4AKYbikVVOftyVdz4a2u25gPL6stmK3rSaY7PnvDbcWJbLpC3ld+Z6jNh9/sfftI0oQU3zNxhH3BbZuvQIGx74hAswwZWgPA7sK1n4g4LEH47DMzwgRUg/A5s65m4wwK5Dr+bppNMquwfj1gPySx77FEpmbReSvbrbz0kM3rqGUoO3NPh2TM8As4IEH4748qo+REo+9tfVTXrQhW+964UCmnbqWeo5tIrlaqqyklBhN85Yd7pi+Q9/N7pGbhgAMJvFzSBEnwlQPjtq3YyGRcJEH67qBmU4isBwm9ftZPJuEggb+F3M4PCjz5UZOGt1opwszLcOkIh1X/pEMVOOUOxiccoXVLqIjVKQaBjAcJv7hA/CBRueE89zz9XpU88Zk0nvve+2vK7Pyq+3/45nR7hd065u30xwu9u0336RsJvGxAZAoFmAoTf3A4IOCNA+O2MK6MiQPjNPYCAMwJuCL+bZpZOq3TZUmtv8LK/PaRQQ4P1klldGJt8vKKnTFX8M591BoJREbBRgPDbRkyGyr1AIqGKm36ryquvlNnuxPweXPOTn2vbGdOlgoKc10P4nXPybl0wr+H35q01Omfm9Xpl1dodih8zeqhumnOBelVVdmtiuXwT4XcutblWEAQIv4PQZeaYDwHC73yoc80gCBB+B6HLzDEfAq4Kv5sBFGzerMidf7a2RSl64/WmV8zKw9ip2x+SmaO/cp+PvnBNbwsQfnu7f0GuvmTFcvX8/tkqeutNiyF6wkmq/sUcpXbZNW8shN95o+/ShfMafl83d5FV7I/OPr5LRbvtZMJvt3WEerwuQPjt9Q5Sv1sFCL/d2hnq8roA4bfXO0j9bhVwa/jd3Kvk+Wet1eDhB+6zViGaw2yDUveNb1qrwesPP9LaJoUDAbcIEH67pRPUka1AwaZN6nHpLEXuvF1Kp5UYOUpbbrhRDWPHZTuEY+cRfjtGa+vAeQu/zarvWVfN04xzp2jYoH62TirXgxF+51qc6/ldgPDb7x1mfvkSIPzOlzzX9bsA4bffO8z88iXghfA7YxOqrVXk7r9Yq8GLX36piSw5YOCnD8nsPyBflFwXgSYBwm9uBs8IpNOK3HGbelx2sUwAni6vUM3/Xazas78nFRa6YhqE365oQ6dFEH53StT5CYTfnRtxBgJdESD87ooW5yKQvQDhd/ZWnIlAVwQIv7uixbkIZC/gpfC7+ayK//ufxodk3nWHQttqm16qO/Irin73ezJfORDIlwDhd77kuW5XBIpeX62eF5wn87drzBE79jhVX/FLJfv07cowjp9L+O04sS0XyFv4bao3254MHthXk8cfastk8jUI4Xe+5LmuXwUIv/3aWeaVbwHC73x3gOv7VYDw26+dZV75FvBq+J1xC9XFFL7/XkVuu6UpwDGvxffbX7U/+JFiEyfn5QFt+e4r18+vAOF3fv25escCZvuoyjmXq2Lu76VEQonBQ7X1+t+r/pDDXElH+O3KtuxQVF7D7zXrN2rh4sc045wpCpeVeEOsjSoJvz3bOgp3qQDht0sbQ1meFyD89nwLmYBLBQi/XdoYyvK8gNfD7+YNKHrzDUVuvVmRO25VQXW19VJyz8HaNu0sRU8+XalevTzfLybgDQHCb2/0KYhVhh+6Tz0uvkiFGzdY06+5+OequeAiV1MQfru6PU3F5S38Nnt+nzPzer2yam2bUmNGD9VNcy5Qr6pK10sSfru+RRToMQHCb481jHI9I0D47ZlWUajHBAi/PdYwyvWMgJ/C7wx6qL5O4cV3KzL/jyr5zwvWt80DMmMTj1H0jOlqOOiLnukPhXpTgPDbm33zc9WFG95Tz/PPVekTj1nTrD/iy9YDLZMeeE4C4bc37sy8hd/e4MmuSsLv7Jw4C4FsBQi/s5XiPAS6JkD43TUvzkYgWwHC72ylOA+Brgn4MfxuLlD86suK3DxXkXsXKRTdZr0UH723olOnK3r8SUpXun8hWNc6ytluECD8dkMXqMESSCRUeeNvVHHNVTLbnST3HKSts69V3dfGewaI8NsbrSL8tqFPhN82IDIEAs0ECL+5HRBwRoDw2xlXRkWA8Jt7AAFnBPwefmfUQrW1ity1UJEF81S8aqX17XSkXNHJxyn6ne8qvu9+zgAzaiAFCL8D2XbXTbpkxXL1/P7ZKnrrTaVLSlR73vmq/fFMpUvLXFdrRwURfnujXXkPv1e8tFpTz5/TQmvBDTM1dv9R3hCURPjtmVZRqEcECL890ijK9JwA4bfnWkbBHhEg/PZIoyjTcwJBCb+bN6bk+WcVuWWewg/cp1BDvfVSw2cPsLZEiU0+TumysOf6SMHuEiD8dlc/glZNwaZN6nHpTEXuXCil09aDLM0DLc2DLb14EH57o2t5Db9N8H3t3EUt9vY2D8E8+6Jrde7pkzR5/KGeUCT89kSbKNJDAoTfHmoWpXpKgPDbU+2iWA8JEH57qFmU6imBIIbfmQYVbN6syB23KbLgTyp6e4317VRVlaLHn6zomWcrMXyEp3pJse4RIPx2Ty8CVUk6rfKFt6rysotlfn9L7tFP1VderdjEyZ5mIPz2RvvyFn7H6hp06TXzddyEw3dY5W1C8bsfWqrLZkxTuKzE9ZKE365vEQV6TIDw22MNo1zPCBB+e6ZVFOoxAcJvjzWMcj0jEOTwu6lJ6bRKly21tkQJ/+2v1h655qj/wsGKnvEd1X3zGGvLAA4EshUg/M5WivPsEih6fbV6fv8slbz4b6moSLVnnaeaWZcoHY7YdYm8jUP4nTf6Ll04b+H35q01mnXVPM04d4qGDer3/+3dC3RU1b3H8f9kknklAQFbwYqAaHlYLZaFglpLtdYKIhQtWh+3CEWqtV0+Fl6weqltFcSr9dZWpSmvqlWjgEpBbasXFQWlWq8PfPDwgfKwKkhek2Qyc9c+YcbJSE/WWwAAIABJREFUMEkmw9kz55z9zVquanJm7/3//PdpyC+HPW0WrZ7+vvmO+2XONdOkR3fnv8kH4XeX9hwXI9CpAOF3p0RcgEBeAoTfebHxIgQ6FSD87pSICxDIS4Dwuy2bf+cOidy9yPrH/9GH1hfjvQ6U+vMulLqLpknLof3zcuZFZgkQfpvV72JW69++Tcrn/14qfn+btYym446X3bf+XmKD3HPMcWd+hN+dCTnj60ULv3ny2xkbgFUg4EQBwm8ndoU1eUGA8NsLXaQGJwoQfjuxK6zJCwKE3+13MfT4SutIlNA/nkhd1Pjt71gheHTMOC+0nxo0CRB+a4Jl2JSAOte78obZUr5kgfW5+IFfkj2/vFHqzz3fc0qE3+5oadHCb8WzbNUzUr1iNWd+u2OvsEoECiZA+F0waiYyTIDw27CGU27BBAi/C0bNRIYJEH533nD/1g8ksmSBlN+9SEo+/cR6QctBva03yKy/YLK09O7T+SBcYZQA4bdR7S54sRV33S6VN/1GfDU11tx1U6fLnmt/JYlK55/qkA8W4Xc+aoV/TVHDb1WuOt978uVz21S++LaZ+5wDXnia3Gfk2JPcrbgSgVwECL9zUeIaBLouQPjddTNegUAuAoTfuShxDQJdFyD87ppZ+OGHJLKwSoLPP5t6YcPYM6V+8o9FPRXOBwJKgPCbfaBDIPS3x6Tif26RwAvPW8M3njRaPp9zq6eOOMnmRvitYzfZP2bRw2/7Syr8iITfhTdnRm8LEH57u79UVzwBwu/i2TOztwUIv73dX6orngDhd372pZs2WkeiRP6yREr27LEGifUb0Po0+HkXSrxnr/wG5lWeECD89kQbHVNE4KX1Unn9talfuqn3Hvj8NzcZc/wS4bdjtmKHCylq+H3r/GrZ8fFncv2MKRIOtb5DdfIs8JHDh8rEMSe5QpHw2xVtYpEuEiD8dlGzWKqrBAi/XdUuFusiAcJvFzWLpbpKgPB7/9rla4xKePlSiSz6o6iASn0kAkGJjhtvnQ3eNPKE/ZuAV7tSgPDblW1z3KJL33lbuv3mvyS0aoW1NvXmu7VX/afUTZ4miUBrvmfCB+G3O7pctPCbN7x0xwZhlQgUQ4DwuxjqzGmCAOG3CV2mxmIIEH4XQ505TRAg/Lavy2VvviGRP90lkYceEF9drTVw8+AhUj95mtSfc75nz+O1T9A7IxF+e6eXxajEv2O7VN54vUQeuFekpUUS5RVSe+nPpfanl0uioqIYSyrqnITfReXPefKihd+7Pq+RWTdWyYxLz5WB/Q5us+DN72+Tm++4X+ZcM016dHf+ofg8+Z3zfuNCBHISIPzOiYmLEOiyAOF3l8l4AQI5CRB+58TERQh0WYDwu8tknb7AV1srkQfvk8iiKinb8Lp1fSIckYaJP5C6aZdI89eO7nQMLnC3AOG3u/tXrNWX7N4tFb+9Scr/NF/U3ypRT3fX/8dUqZkxy3rq29QPwm93dL5o4TdPfrtjg7BKBIohQPhdDHXmNEGA8NuELlNjMQQIv4uhzpwmCBB+6+1yYP0LEllcJeGHl1lhlvpoHvYN60iUhrMmSSIU1rsARi+KAOF3UdhdO6kv2iAVd/1Byn/3363vIVBSIvVnTZKaa2ZLS99+rq3LroUTftslqXecooXfqqz1r7wls+ZUyfx5V6We/lZPfU+/+ha59EfjOfNbb+8ZHQHHChB+O7Y1LMzlAoTfLm8gy3esAOG3Y1vDwlwuQPhdmAaW7NolkfvvlsjiBVK6eaM1abxbN2mYdL4VhMcGDS7MQpilIAKE3wVhdv8ksZiU37tEKubdIP6dO6x6ot85TWpm3yDNQ4a6vz6bKiD8tglS8zBFDb9Vbcmwe/vOT1OlLr5tpowY5p5vsBx7onmXMrxxAoTfxrWcggskQPhdIGimMU6A8Nu4llNwgQQIvwsEnTZN8NmnrSNRQo+tEF9zs/UV9caYKgSPjptg1BvZFV6/MDMSfhfG2bWzJBISXrFcKn/zSyndsqn1/wOGj5A9v54rTceOcm1ZuhZO+K1L1t5xix5+21vO/o2WPIpFjXL9jCkSDrW+Q+2yVc/IdfMWWv8+9pSRbb6mPkf4vX/uvBqBTAHCb/YEAnoECL/1uDIqAoTf7AEE9AgQfutxzWXUkk/+LeV3L5LInxeKf+sH1kvUub71510o9ZN/LLF+A3IZhmscKED47cCmOGRJweeflW7XzJCy11+1VqT+1seea38l0dPPcMgKnbcMwm/n9STbigi/96okg++VT65rE3Cro1lumV8td869wnrzzVvnV1uvuHL6pJQn4bc7NjurdI8A4bd7esVK3SVA+O2ufrFa9wgQfrunV6zUXQKE3w7oVzwuoSf/1vo0+D+eEInHRXw+aRx9stRddLFETxsj4vc7YKEsIVcBwu9cpcy5ToXd3WbPkuDT/2sV3XJIX6mZeZ3UTzrPOuObj/YFCL/dsTsIv/f2SYXa/fv2tv5r3UsbUk93Jz8/ccxJ1tcyw3D1OcJvd2x2VukeAcJv9/SKlbpLgPDbXf1ite4RIPx2T69YqbsECL+d1S//Rx9KZPGfJHLvEvF/vLM1JOtzsNRfMFnqfzRVWnr3cdaCWU1WAcJvNkZSoPS9LVL569kSfnSZSCIh8Z49pfaKq6Vu6k844ijHbUL4nSNUkS8j/BZp8zS3OuIkGX6r3sy+eaGMHD409eab6ozyX8ypkhtmTUu9SSfhd5F3MdN7ToDw23MtpSCHCBB+O6QRLMNzAoTfnmspBTlEgPDbIY3Isgx1JnBkwXwJrnkm9dXomHHWkSjRk0917sJZmRB+swmUQPcrL5PyP7ce75sIR6T2kp9J7c+vkkRFBUBdECD87gJWES81PvxWYfd7W3ekjjHJFn7/YNzo1BtwZgu/44lEEVvI1AjkJuDL7TKuQqBLAs0tCSnz81fhuoRWxIt9PhHrOxbftorYBab2kkBzS9z6/0Dr3uK+8lJrqcUhAtxbDmlEO8tobmmRsi2bRe68U3xLlojs2tV65YABIhdfLImpU0UOPNDZRRi7OvXToTnfuMyptJMNvWWL+G66SXxVVakLE5ddJolrrxX58peNvRv2p/AS9Y2KD8cLGB9+q2NNFty3ap9GqTe2nPmzC2Tu7fd0+uT39k8bHN9oFoiAm77h8+S3e/ar9a2e7/euaVjPyqDUNTRLYyzumjWzUAScLKACb/Uzz5cPCMmnexqlJe6m77ZOlmVtCLQK9OkZlu2f8bOWY/dDom18Grn/XoksWSCB9etSS244a5L1NHjjqBMdW4ZpCzPxyW/Tf1zxb/tIKmdfI+HlD35xb06cJLXX/JfE+h9m2i1ga719eoVtHY/B9AgYH35nsqY/+R0OBawjUdRZ4Jz5rWcDMioC2QQIv9kXCOgR4NgTPa6MigDHnrAHENAjwLEnelx1j1r25hsSWfhHiVTfJ766Wmu62KDBUnfRNKk/5wJJVFbqXgLjdyBgYvht6oYIPb5Syu+6ve3xRKefITXX/FKahww1lcXWujn2xFZObYMRfmfQZobfmW9wqcJw9XHl9EmpV3Lmt7b9ycCGChB+G9p4ytYuQPitnZgJDBUg/Da08ZStXYDwWzux1gl89XUSWVptnQ1e9vqr1lyJYEgaJp4t9VOmS9Mxw7XOz+DZBQi/vb0zfA31ErnvHim/6/dSumVTqtjG0SfLnmt/Jc3DvuFtgAJXR/hdYPA8pyP87iT8Vl9Wgfh181rfCEAdh3L9jCmingpPfhB+57n7eBkC7QgQfrM1ENAjQPitx5VRESD8Zg8goEeA8FuPazFGDby0XiKLqyS8fKn4oq1H2TR//RipV0+Dn32OJEIcHVCovhB+F0q6sPP4d2yXyB//YL2JZcnu3dbkiUBQ1NFDtZddYf3tCz7sFyD8tt9Ux4iE3zaoEn7bgMgQCKQJEH6zHRDQI0D4rceVUREg/GYPIKBHgPBbj2sxRy35/HOJPHCPRBZVSenGd6ylxLt1k4ZJ50ndRRcT0BWgOYTfBUAu4BRlr7wsFXf8TsKPLhOJxVrvqQO/ZB0zVPfjn0i8F286q7MdhN86de0bm/DbBkvCbxsQGQIBwm/2AALaBQi/tRMzgaEChN+GNp6ytQsQfmsnLuoEweeetULw0MpHxNfcbK2l6bjjpW7KxRIdN0ESgS/+tnVRF+qxyQm/PdDQeFzCKx+V8jtvl8CLa1MFNQ85Uuou+Zk0nH0u90+B2kz4XSDo/ZyG8Hs/AdXLCb9tQGQIBAi/2QMIaBcg/NZOzASGChB+G9p4ytYuQPitndgRE5R88m8pv2eJRP68QPwfvG+tKd6zl9Sfd6F1LEqs3wBHrNMriyD8dm8n1RvIlt+9SMr/eEfqXhGfT6KnfFfqLvm5NH7r2+4tzqUrJ/x2R+MIv23oE+G3DYgMgQDhN3sAAe0ChN/aiZnAUAHCb0MbT9naBQi/tRM7boLQU39vfRr8sb9aa0tUVkrjiJHSOPZMaThzosR79HDcmt22IMJvt3VMxL/1A+sp7/K/LBFfbW3rvREMSf0PL7Ce9I4NPMJ9RXlkxYTf7mgk4bcNfSL8tgGRIRAg/GYPIKBdgPBbOzETGCpA+G1o4ylbuwDht3Zix05gvXnfoiqJ3H+P+D/6MLVO6wnXn14uzcO+YZ0VzkfXBQi/u25WrFeoI03K7/pD63neez9aeveR+qnTrTPy4wccUKylMe9eAcJvd2wFwm8b+kT4bQMiQyBA+M0eQEC7AOG3dmImMFSA8NvQxlO2dgHCb+3Erpgg8MLzEl72oISXPygln32WWnPj6JMlOna8NIybYL3BHx+5CRB+5+ZUzKsi1fdJpOoOCfzrpdQymr92tNRddoXUn31OMZfG3BkChN/u2BKE3zb0ifDbBkSGQIDwmz2AgHYBwm/txExgqADht6GNp2ztAoTf2oldN0Fw9VMSemSphFcsl5Ldu1Prbzp2lETHjZeGMyZIS99DXVdXIRdM+F1I7dznUmffRx66XyJLFkjpxndSL2w4Y7zUT/+pNI46MffBuLJgAoTfBaPer4kIv/eLr/XFhN82IDIEAoTf7AEEtAsQfmsnZgJDBQi/DW08ZWsXIPzWTuzqCYLPPi2hFcsl/PBDbZ4Ibz7q6xI98/vSMOZMiQ0a7OoadSye8FuHav5jBtc8I5HFf7L2cfIjESmXugsvkrqf/FRa+vbLf3BeqV2A8Fs7sS0TEH7bwEj4bQMiQyBA+M0eQEC7AOG3dmImMFSA8NvQxlO2dgHCb+3EnpkguHaNhB5ZboXh/p07UnXFDj9C1JOzjWPHS9Mxwz1T7/4UQvi9P3r2vLZ0yyYJP/AXUcG3OtYn+RH76iCpmzJd6s8533qzVz6cL0D47fweqRUSftvQJ8JvGxAZAgHCb/YAAtoFCL+1EzOBoQKE34Y2nrK1CxB+ayf23gSJhKg3CVRBePivD4t/20epGtVxKOppcHU8ijomRUpKvFd/DhURfueApOES9Sau4aXVEl76gJS9+kpqhkQwJA1nTpD6yT+WpuOO1zAzQ+oUIPzWqWvf2ITfNlgSftuAyBAIEH6zBxDQLkD4rZ2YCQwVIPw2tPGUrV2A8Fs7secnCKx/QUKP7g3Ct36Qqle9QWbDmHESHTdBmk78liTKyjxvkSyQ8LtwrS7ZtUvCjyyV0NJqCa57TiSRsCZPBALSOPoUiY4/y9qHPOVduJ7YPRPht92iesYj/LbBlfDbBkSGQIDwmz2AgHYBwm/txExgqADht6GNp2ztAoTf2omNmiDw8j8l9OgyCT/6sPg/eO+LILxbN4l+93RpPGOCRE85VRLhiKddCL/1ttdXXyfhVX+V0EP3S/Dpp8TX3PxF4P2tkyU64WwCb70tKOjohN8F5c57MsLvvOm+eCHhtw2IDIEA4Td7AAHtAoTf2omZwFABwm9DG0/Z2gUIv7UTGztB4F8vSVC9WeaKR6T03c1tHKLfGyv15/2HRMeM86QP4beetoYeXynh6vsk9MQq8TVGU5NEv3OaRM86Rxq+N5YnvPXQF3VUwu+i8uc8OeF3zlTtX0j4bQMiQyBA+M0eQEC7AOG3dmImMFSA8NvQxlO2dgHCb+3ETCAiZa/9X+vRKI8ul9LNG9uYqLPBG797ukRPO12ahxzpCS/Cb/vaGHz6fyW8rFrCjy4TX01N28BbPeE99kwCb/u4HTkS4bcj27LPogi/begT4bcNiAyBAOE3ewAB7QKE39qJmcBQAcJvQxtP2doFCL+1EzNBhkDZhtcl9MgyCT32V1H/nv7RNGKktPQ5WJpO+KY0HX+ia8Nwwu/8t33Jnj0SfPpJCT6xSsreeN36xUnyw3rCm8A7f1yXvpLw2x2NI/y2oU+E3zYgMgQChN/sAQS0CxB+aydmAkMFCL8NbTxlaxcg/NZOzAQdCPi3fSShvz8uwX88IcHVT4mvob7N1fEDDpDGE06S5hO+KY0jT5Dmo4e5wpPwu2ttCj7/rARefEFKX39Vwg8/1ObF0ZNPlejESTzh3TVST11N+O2OdhJ+29Anwm8bEBkCAcJv9gAC2gUIv7UTM4GhAoTfhjaesrULEH5rJ2aCLggE1q+T4HNrpGztGgm+8Lz4amvbvDpRUSGNx39Tmk48SZpGnShNxwzvwuiFu5Twu2Pr4DOrJaAC7+eeleDaNftcrH7hET37HGkY931RvwDhw2wBwm939J/w24Y+EX7bgMgQCBB+swcQ0C5A+K2dmAkMFSD8NrTxlK1dgPBbOzET7IeAetPMwNo1Elj7nATWPSclu3a1DcNDYWkadbw0nThaGtX/HjtqP2az76WE319Y+qINElj7vNXH4JpnJPDi2n2gW/oeKo3qmJvjjpfoaWOk5aDe9jWDkVwvQPjtjhYSftvQJ8JvGxAZAgHCb/YAAtoFCL+1EzOBoQKE34Y2nrK1CxB+aydmAhsFyt7c0BqGP/es9b/+j3fuG4aPOFaa1NPh6p8Rx0oiELRxBbkNZXL4rY6uCa59TsrWPC3B59dI4JWXRWKxNnCx/odZZ7o3nXCSNH7r29LSu09usFxlpADhtzvaTvhtQ58Iv21AZAgECL/ZAwhoFyD81k7MBIYKEH4b2njK1i5A+K2dmAk0CpS+/ZZ1bEZAPVGswvCdO/aZrWH8WRL76iCJH9RbmoceKc1HHiWJ8gqNqxIxKfwue+M1Kdm9W8peeVlCDz8k6mn9zI/YoMGp42rUkSbxA7+k1Z/BvSVA+O2OfhJ+29Anwm8bEBkCAcJv9gAC2gUIv7UTM4GhAoTfhjaesrULEH5rJ2aCAgr4t75vvXGiOk9aheKl77y9z+zq3PCmo48RKfVL7ND+Ej/scGkeNFha+g2Q5sFDbFmt18Lvks8+ldLNm6x//B9ttZ7kDj75Nwm8/M+sXs1DjpSmE77Z+mT38SdKvNeBtrgyiJkChN/u6Dvhtw19Ivy2AZEhECD8Zg8goF2A8Fs7MRMYKkD4bWjjKVu7AOG3dmImKKJAySf/lsD6F6Rsw+tS+sbrUvbWG1kD8fQlthzSV2KHHyGxAQOlZeDhEjts7z+HH5FzJW4Mv9Wbi5ZtfFv8WzZL6bubpVT9+/vvSdnbb4qvpqbd2tX53LEhQ6V56FHSfNwo6wnveI8eOVtxIQKdCRB+dybkjK8TftvQB8JvGxAZAgHCb/YAAtoFCL+1EzOBoQKE34Y2nrK1CxB+aydmAgcKlG7aKKXvbRF1bIpfBb3q3zdtFP+HWztcrTqrOnbYQIkNPVIaTz3dujZeXi4qME8/ysPp4Xf4kaVS+u4W8W96R0q3qNrfEfV0d0cfzV87WmL9+lu/GIgPGCjN6hcEQ4+SeLduDuwwS/KSAOG3O7pJ+G1Dnwi/bUBkCAQIv9kDCGgXIPzWTswEhgoQfhvaeMrWLkD4rZ2YCVwmUPbWm+J//10pe/stKdmyyQqJSzdvFP+O7R2Hw0cPk3hFpXWNT0TKAqXSGIpY4XCiW3eJdz9AEhWVkujW7YvPqa9Vdrf+O96tu6gjWTI/fPV14quvF/W/JXX14qurFV80Kr7amtbP7f2ar65OpKlpn9eX7NolwdX/EP+HH4ov2tBxuH/EV8UK+NUxMAMGSmygegL+MGn5yiEu6yLL9ZIA4bc7ukn4bUOfCL9tQGQIBAi/2QMIaBcg/NZOzASGChB+G9p4ytYuQPitnZgJPCKggmN1hrj1lPi770rJB+9ZT4r7t30opR9u7fBoEKcQqDf6bDnkEIl9pa/EBg2Rlv4DJHb4VyU2YIC0HNrfKctkHQi0ESD8dseGIPy2oU+E3zYgMgQChN/sAQS0CxB+aydmAkMFCL8NbTxlaxcg/NZOzAQGCagzxkv2fC7+mj3S09ckuz/YaT2hrT6nzs327d7d+u979khJTevnUv+9Z4/4Gur30UoEQ5Ioj4gKrhORiMQj5ZIoL5dEOCIJ9e8R9bVy6/gV2XtN6+fLJa6+VtlN4r16SctX+kq8e3eDukGpXhEg/HZHJwm/begT4bcNiAyBAOE3ewAB7QKE39qJmcBQAcJvQxtP2doFCL+1EzOBgQJOP/PbwJZQsosFCL/d0TzCbxv6RPhtAyJDIED4zR5AQLsA4bd2YiYwVIDw29DGU7Z2AcJv7cRMYKAA4beBTadkbQKE39pobR2Y8NsGTsJvGxAZAgHCb/YAAtoFCL+1EzOBoQKE34Y2nrK1CxB+aydmAgMFCL8NbDolaxMg/NZGa+vAhN82cBJ+24DIEAgQfrMHENAuQPitnZgJDBUg/Da08ZStXYDwWzsxExgoQPhtYNMpWZsA4bc2WlsHJvy2gZPw2wZEhkCA8Js9gIB2AcJv7cRMYKgA4behjads7QKE39qJmcBAAcJvA5tOydoECL+10do6MOG3DZyE3zYgMgQChN/sAQS0CxB+aydmAkMFCL8NbTxlaxcg/NZOzAQGChB+G9h0StYmQPitjdbWgQm/beAk/LYBkSEQIPxmDyCgXYDwWzsxExgqQPhtaOMpW7sA4bd2YiYwUIDw28CmU7I2AcJvbbS2Dkz4bQMn4bcNiAyBAOE3ewAB7QKE39qJmcBQAcJvQxtP2doFCL+1EzOBgQKE3wY2nZK1CRB+a6O1dWDCb1s5GQwBBBBAAAEEEEAAAQQQQAABBBBAAAEEEEDACQKE307oAmtAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQsFWA8NtWTgZDAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQcIIA4bcTusAaEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBGwVIPzOk3PZqmfkunkLrVePPWWkXD9jioRDgTxH42UImCmw+f1tMv3qW2T7zk9TAEcNOUzunHuF9OheKQ3RJpl980JZ+eQ66+u/vnqKTBxzkplYVI1ADgLqnrr5jvtlzjXTrHso+dHZvbT+lbdk8uVzrcvT78EcpuQSBIwQaO/eunV+tSy4b1Ubg/TvVdxbRmwPisxDIPPeyfwz3q7Pa+SSmb+V197cYo2++LaZMmLY4NRM/CyWBzovMUKgo3sr875SIH0O6iXz510lA/sdbPlwbxmxTSgSAeMECL/zaLn6QeaW+dWpgE59g1EfV06flMdovAQBcwVUmPCLOVVyw6xpqT9wpWuk31vJP6xdNX1Smx9+zNWjcgS+EEj/YSZbeN3RvZR5H6ofeta9tIFf6rLBEBCRrtxbmWDcW2whBLILqF/I3rnkYbno3NOtX9QmH4aYM2ua9We85C9sRw4faj30kHkv8bMYOwuB/O6tzn6e4t5iZyGAgFcFCL/z6KwKEfr37Z16AjXzm0QeQ/ISBIwU6Cj8Vn84m3Vjlcy49NxUMM4vmozcJhTdBYFsT6d2di+psPu9rTtSv8Dt7JdSXVgOlyLgGYGOnvxWRWZ7AIJ7yzPtpxDNAtnC7vS/xZT5dX4W09wQhveMQOa901n4zb3lmdZTCAIIZAgQfndxS2R+A1EvJyjoIiKXI7BXIPPYk/QnVrPdVzyRytZBoGOBbAFdZ/eSevouPbzr7AcjeoCAiQK5HnuSfnRD5i9subdM3DnUnItA5r2R7cGi5P10yY8mWEfiJZ8K52exXIS5xlSBzHsr89iT9CNPyDlM3SXUjYAZAoTfXexz8pvCD8aNTh29QPjdRUQuR6AdAfWDzY6PP7OOW9i285N9zi4m/GbrIJBf+J15Dnj6vaTC7/S/zURAxy5DYF+B9sLv9Cszj27IfIKOe4udhUB2gcxfFKnw+8EVq9scv5UZfvOzGLsJgc4FOvtbs+rPg9UrVlvHuYaCQesXS9xbnbtyBQIIuE+A8LuLPeM3ol0E43IEuiCQHi58trtmn/PACb+7gMmlRgrw5LeRbafoAgjkEn6rZaQH3jz5XYDGMIXrBdIffAiHAlY9PPnt+rZSgAMEst1bmctKPxrv4IMO5G9VOKBvLAEBBPQIEH7n4cpZWHmg8RIEchBIDxfU5Zz5nQMalyCQJsCZ32wHBPQI5BN+c+a3nl4wqncE2gvnMu83zvz2Ts+ppDACuQTfaiWZ7wtDzlGY/jALAggUXoDwOw9z3gU5DzRegkAWgSdWvyiHDzik3Te0TH9qjr8uzhZCoHOBXN6UL/Neyjy6i79h0bkzV5gn0N4vllY9uU7On3iqBZJ5L3FvmbdPqDh3gY6OY8j2Bpi/mFMlN8yaZv2ZkZ/FcnfmSvMEOrq31L2jPkYMG2z9b+af+bi3zNsvVIyAKQKE33l2Wn2juG7eQuvVY08Z2eZMujyH5GUIGCeg/oA1+fK5qboz76XkDz8rn1xnXZP+RmLGYVEwAh0IZL6Bkbp06g/HyJXTJ1mv6uxeSr8X0994FnQETBfo6N7KvK+U1eLbZqZCBfXf3Fum7yDqzyaQ7b7K/Jkq85rMe4ufxdhbCOwr0Nm9pd5TafrVt8j2nZ9aL872Zz7uLXYWAgh4UYDw24tdpSYEEEAAAQRGDXraAAAK6UlEQVQQQAABBBBAAAEEEEAAAQQQQMBwAcJvwzcA5SOAAAIIIIAAAggggAACCCCAAAIIIIAAAl4UIPz2YlepCQEEEEAAAQQQQAABBBBAAAEEEEAAAQQQMFyA8NvwDUD5CCCAAAIIIIAAAggggAACCCCAAAIIIICAFwUIv73YVWpCAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQMFyD8NnwDUD4CCCCAAAIIIIAAAggggAACCCCAAAIIIOBFAcJvL3aVmhBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAcMFCL8N3wCUjwACCCCAAAIIIIAAAggggAACCCCAAAIIeFGA8NuLXaUmBBBAAAEEEEAAAQQQQAABBBBAAAEEEEDAcAHCb8M3AOUjgAACCCCAAAIIIIAAAggggAACCCCAAAJeFCD89mJXqQkBBBBAAAEEEEAAAQQQQAABBBBAAAEEEDBcgPDb8A1A+QgggAACCCCAAAIIIIAAAggggAACCCCAgBcFCL+92FVqQgABBBBAAAEEEEAAAQQQQAABBBBAAAEEDBcg/DZ8A1A+AggggAACCCCAAAIIIIAAAggggAACCCDgRQHCby92lZoQQAABBBBAAAEEEEAAAQQQQAABBBBAAAHDBQi/Dd8AlI8AAggggAACCCCAAAIIIIAAAggggAACCHhRgPDbi12lJgQQQAABBBBAAAEEEEAAAQQQQAABBBBAwHABwm/DNwDlI4AAAggggAACCBRfYNmqZ2TdSxvk+hlTJBwKFH9BrAABBBBAAAEEEEAAAQ8IEH57oImUgAACCCCAAAIIeFGgIdoks29eKCufXNemvD4H9ZL5866Sgf0Otj6vguPqFavlzrlXSI/ula6kIPx2ZdtYNAIIIIAAAggggIDDBQi/Hd4glocAAggggAACCJgqkAy/e3+5p1w5fVKK4db51fLiK2+5OuzO7Cnht6m7nLoRQAABBBBAAAEEdAoQfuvUZWwEEEAAAQQQQACBvAXaC78zn/TODI6T/33GqaPkkpm/tebPfFo8c1HJuUYOHyrvbd0hC+5bZV0y9pSRqaNIdn1eY4131fRJMmLYYOvr6a+bOOYkSV5z0Tnfk0UPPC6vvbnFuu7XV0+Rrx95uEy/+hbZvvNT63OLb5uZGie55qOHDpQ5t99rff2oIYftE/Cr666btzC1/GxjpNed/vW8G8ELEUAAAQQQQAABBBBwqQDht0sbx7IRQAABBBBAAAGvC2QLv7N9Llv4rQLiqT8ck3piXD0tvuPjz9o9Uzv9iJVkYJwMsieNGy3pwXYu4bfqTfIYlvWvvCWTL5/bJsxWn7tlfnXqmmSorUJyNZf6yFxzZp2b399mhelzZk2zQvTkGOl1e32PUB8CCCCAAAIIIIAAAh0JEH6zPxBAAAEEEEAAAQQcKdDemd9qsekhcXtPfqe/eWRm2JxZcOYT3MmvqwBafahjV7ry5Hd6QJ7tdZmfy3bsiQq3fzGnSm6YNU16HlAps26skhmXnps66zwZkCfXx9EpjtzGLAoBBBBAAAEEEECgiAKE30XEZ2oEEEAAAQQQQACB9gXaO/Yk84lsr4bf6QF5zx7d2hyZkq6WfNKb8Ju7CQEEEEAAAQQQQACBtgKE3+wIBBBAAAEEEEAAAUcKtBd+q8WmB72PPbVO1r20IXWkSbYQ2I1PfmeG38mnwAf2Ozhrvwi/HbmNWRQCCCCAAAIIIIBAEQUIv4uIz9QIIIAAAggggAAC7Qt0FH6nn4ft1fA789iTzDfbzJQj/OZuQgABBBBAAAEEEECgrQDhNzsCAQQQQAABBBBAwJEC7YXfyTd6vPRH4603hyzUsSftvdmmenPN5BnkuZzvrbA7O/M721wq8F/11Asyf95VqXO/1RPtW7d9nNXBkU1lUQgggAACCCCAAAIIFFCA8LuA2EyFAAIIIIAAAgggkLtAR294ufi2mTJi2GBrsEKF3+mh9WtvbrHmnvWz8+XVDZtl5PChVgC9P+G3CtHTP5Jnead/TtWafl2fg3qlwnCe/M59b3ElAggggAACCCCAgBkChN9m9JkqEUAAAQQQQAABBBBAAAEEEEAAAQQQQAABowQIv41qN8UigAACCCCAAAIIIIAAAggggAACCCCAAAJmCBB+m9FnqkQAAQQQQAABBBBAAAEEEEAAAQQQQAABBIwSIPw2qt0UiwACCCCAAAIIIIAAAggggAACCCCAAAIImCFA+G1Gn6kSAQQQQAABBBBAAAEEEEAAAQQQQAABBBAwSoDw26h2UywCCCCAAAIIIIAAAggggAACCCCAAAIIIGCGAOG3GX2mSgQQQAABBBBAAAEEEEAAAQQQQAABBBBAwCgBwm+j2k2xCCCAAAIIIIAAAggggAACCCCAAAIIIICAGQKE32b0mSoRQAABBBBAAAEEEEAAAQQQQAABBBBAAAGjBAi/jWo3xSKAAAIIIIAAAggggAACCCCAAAIIIIAAAmYIEH6b0WeqRAABBBBAAAEEEEAAAQQQQAABBBBAAAEEjBIg/Daq3RSLAAIIIIAAAggggAACCCCAAAIIIIAAAgiYIUD4bUafqRIBBBBAAAEEEEAAAQQQQAABBBBAAAEEEDBKgPDbqHZTLAIIIIAAAggggAACCCCAAAIIIIAAAgggYIYA4bcZfaZKBBBAAAEEEEAAAQQQQAABBBBAAAEEEEDAKAHCb6PaTbEIIIAAAggggAACCCCAAAIIIIAAAggggIAZAoTfZvSZKhFAAAEEEEAAAQQQQAABBBBAAAEEEEAAAaMECL+NajfFIoAAAggggAACCCCAAAIIIIAAAggggAACZggQfpvRZ6pEAAEEEEAAAQQQQAABBBBAAAEEEEAAAQSMEiD8NqrdFIsAAggggAACCCCAAAIIIIAAAggggAACCJghQPhtRp+pEgEEEEAAAQQQQAABBBBAAAEEEEAAAQQQMEqA8NuodlMsAggggAACCCCAAAIIIIAAAggggAACCCBghgDhtxl9pkoEEEAAAQQQQAABBBBAAAEEEEAAAQQQQMAoAcJvo9pNsQgggAACCCCAAAIIIIAAAggggAACCCCAgBkChN9m9JkqEUAAAQQQQAABBBBAAAEEEEAAAQQQQAABowQIv41qN8UigAACCCCAAAIIIIAAAggggAACCCCAAAJmCBB+m9FnqkQAAQQQQAABBBBAAAEEEEAAAQQQQAABBIwSIPw2qt0UiwACCCCAAAIIIIAAAggggAACCCCAAAIImCFA+G1Gn6kSAQQQQAABBBBAAAEEEEAAAQQQQAABBBAwSoDw26h2UywCCCCAAAIIIIAAAggggAACCCCAAAIIIGCGAOG3GX2mSgQQQAABBBBAAAEEEEAAAQQQQAABBBBAwCgBwm+j2k2xCCCAAAIIIIAAAggggAACCCCAAAIIIICAGQKE32b0mSoRQAABBBBAAAEEEEAAAQQQQAABBBBAAAGjBAi/jWo3xSKAAAIIIIAAAggggAACCCCAAAIIIIAAAmYIEH6b0WeqRAABBBBAAAEEEEAAAQQQQAABBBBAAAEEjBIg/Daq3RSLAAIIIIAAAggggAACCCCAAAIIIIAAAgiYIUD4bUafqRIBBBBAAAEEEEAAAQQQQAABBBBAAAEEEDBKgPDbqHZTLAIIIIAAAggggAACCCCAAAIIIIAAAgggYIYA4bcZfaZKBBBAAAEEEEAAAQQQQAABBBBAAAEEEEDAKIH/B9Oe7InklJrwAAAAAElFTkSuQmCC",
"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": "c75683c7",
"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
}