{
"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: Nov. 28, 2022"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "d142d450-f2af-42fa-9026-2b44919c5b40",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Added 'D:\\Docs\\- MY CODE\\BioSimulations\\life123-Win7' to sys.path\n"
]
}
],
"source": [
"# Extend the sys.path variable, to contain the project's root directory\n",
"import set_path\n",
"set_path.add_ancestor_dir_to_syspath(3) # The number of levels to go up \n",
" # to reach the project's home, from the folder containing this notebook"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "9897632a-191c-4cdc-9bc2-a1143c2f0d36",
"metadata": {},
"outputs": [],
"source": [
"from experiments.get_notebook_info import get_notebook_basename\n",
"\n",
"from src.life_1D.bio_sim_1d import BioSim1D\n",
"from src.modules.reactions.reaction_data import ReactionData as chem\n",
"from src.modules.movies.movies import MovieArray\n",
"from src.modules.numerical.numerical import Numerical as num\n",
"\n",
"import numpy as np\n",
"\n",
"import plotly.express as px\n",
"import plotly.graph_objects as go"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "0c291bcf-fd3a-443d-9df1-1b8ce9fefcae",
"metadata": {},
"outputs": [],
"source": [
"# We'll be considering just 1 chemical species, \"A\"\n",
"diffusion_rate = 10.\n",
"\n",
"chem_data = chem(diffusion_rates=[diffusion_rate], names=[\"A\"])"
]
},
{
"cell_type": "markdown",
"id": "82504138-da32-486a-930d-bc9419cb1dc2",
"metadata": {},
"source": [
"# BASELINE\n",
"This will be our initial system, whose adherence to the diffusion equation we'll test.\n",
"Afterwards, we'll tweak individual parameters - and observed their effect on the closeness of the approximation"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "bca9a060-f62d-433a-9924-1fb2314870f3",
"metadata": {},
"outputs": [],
"source": [
"# Parameters of the simulation run (the diffusion rate got set earlier, and will never vary)\n",
"delta_t = 0.01\n",
"n_bins = 300\n",
"delta_x = 2 # Note that the number of bins also define the fraction of the sine wave cycle in each bin\n",
"algorithm = None # \"Explicit, with 3+1 stencil\""
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "696122b0-c8c4-4d43-a4fc-319f95221602",
"metadata": {},
"outputs": [],
"source": [
"# Initialize the system\n",
"bio = BioSim1D(n_bins=n_bins, chem_data=chem_data)\n",
"\n",
"# Initialize the concentrations to 2 superposed sine waves\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=1, amplitude=10, bias=50)\n",
"bio.inject_sine_conc(species_name=\"A\", frequency=2, amplitude=8)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "338bd51c-6b12-4cf0-8526-797389fbc55c",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/vnd.plotly.v1+json": {
"config": {
"plotlyServerURL": "https://plot.ly"
},
"data": [
{
"hovertemplate": "Chemical=A
Bin number=%{x}
concentration=%{y}",
"legendgroup": "A",
"line": {
"color": "red",
"dash": "solid"
},
"marker": {
"symbol": "circle"
},
"mode": "lines",
"name": "A",
"orientation": "v",
"showlegend": true,
"type": "scatter",
"x": [
0,
1,
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
12,
13,
14,
15,
16,
17,
18,
19,
20,
21,
22,
23,
24,
25,
26,
27,
28,
29,
30,
31,
32,
33,
34,
35,
36,
37,
38,
39,
40,
41,
42,
43,
44,
45,
46,
47,
48,
49,
50,
51,
52,
53,
54,
55,
56,
57,
58,
59,
60,
61,
62,
63,
64,
65,
66,
67,
68,
69,
70,
71,
72,
73,
74,
75,
76,
77,
78,
79,
80,
81,
82,
83,
84,
85,
86,
87,
88,
89,
90,
91,
92,
93,
94,
95,
96,
97,
98,
99,
100,
101,
102,
103,
104,
105,
106,
107,
108,
109,
110,
111,
112,
113,
114,
115,
116,
117,
118,
119,
120,
121,
122,
123,
124,
125,
126,
127,
128,
129,
130,
131,
132,
133,
134,
135,
136,
137,
138,
139,
140,
141,
142,
143,
144,
145,
146,
147,
148,
149,
150,
151,
152,
153,
154,
155,
156,
157,
158,
159,
160,
161,
162,
163,
164,
165,
166,
167,
168,
169,
170,
171,
172,
173,
174,
175,
176,
177,
178,
179,
180,
181,
182,
183,
184,
185,
186,
187,
188,
189,
190,
191,
192,
193,
194,
195,
196,
197,
198,
199,
200,
201,
202,
203,
204,
205,
206,
207,
208,
209,
210,
211,
212,
213,
214,
215,
216,
217,
218,
219,
220,
221,
222,
223,
224,
225,
226,
227,
228,
229,
230,
231,
232,
233,
234,
235,
236,
237,
238,
239,
240,
241,
242,
243,
244,
245,
246,
247,
248,
249,
250,
251,
252,
253,
254,
255,
256,
257,
258,
259,
260,
261,
262,
263,
264,
265,
266,
267,
268,
269,
270,
271,
272,
273,
274,
275,
276,
277,
278,
279,
280,
281,
282,
283,
284,
285,
286,
287,
288,
289,
290,
291,
292,
293,
294,
295,
296,
297,
298,
299
],
"xaxis": "x",
"y": [
50,
50.544429428667165,
51.088179283950524,
51.630571063807565,
52.170928407051974,
52.70857815921861,
53.242851432961885,
53.77308466117989,
54.29862064106969,
54.81880956733467,
55.333010052784,
55.84059013458597,
56.34092826446227,
56.83341428113769,
57.317450363391465,
57.79245196208926,
58.25784870961231,
58.713085305138996,
59.157622374276286,
59.59093730158345,
60.01252503457715,
60.42189885785704,
60.81859113604284,
61.20215402426825,
61.57216014503327,
61.92820323027551,
62.26989872758102,
62.596884369518115,
62.90882070514144,
63.20539159277945,
63.48630465328596,
63.751291683005014,
64.00010902576838,
64.2325379033164,
64.44838470360496,
64.64748122653477,
64.82968488671305,
64.99487887293242,
65.1429722641268,
65.27390010164028,
65.38762341772012,
65.48412922022146,
65.5634304335874,
65.62556579624368,
65.67059971462322,
65.6986220741107,
65.70974800727217,
65.704117619809,
65.6818956747483,
65.64327123545478,
65.5884572681199,
65.51769020445475,
65.43122946538178,
65.3293569465877,
65.21237646686652,
65.08061318024517,
64.93441295294744,
64.77414170631202,
64.60018472683976,
64.41294594460129,
64.21284718129132,
64.00032736926742,
63.775841742961624,
63.53986100410003,
63.29287046221103,
63.03536915194446,
62.767868928764315,
62.49089354461371,
62.204977705185655,
61.910666110463616,
61.60851248022481,
61.29907856622386,
60.98293315279715,
60.66065104764711,
60.33281206458206,
60.00000000000001,
59.66280160491487,
59.32180555433006,
58.977601415768284,
58.63077861876623,
58.28192542714066,
57.93162791582594,
57.580468954074114,
57.229027196796395,
56.87787608580947,
56.527582862731656,
56.178707595252185,
55.83180221847259,
55.48740959299154,
55.146062581373926,
54.80828314461176,
54.47458146014833,
54.14545506299808,
53.82138801145301,
53.50285007882258,
53.19029597260686,
52.884164582453884,
52.58487825820056,
52.29284211924252,
52.00844339642252,
51.73205080756888,
51.46401396775535,
51.204662835292,
50.954307194392996,
50.71323617540191,
50.481717813388244,
50.2599986458607,
50.04830335027312,
49.846834421928385,
49.65577189281366,
49.47527309182776,
49.305472446787945,
49.146481328529276,
48.9983879373349,
48.86125723186071,
48.735130900642396,
48.62002737619715,
48.51594189165739,
48.42284657979783,
48.34069061424248,
48.269400392563504,
48.20887976090919,
48.159010279725436,
48.11965153006181,
48.09064145988159,
48.071796769724486,
48.062913337001035,
48.063766678128985,
48.07411244765571,
48.093686973444406,
48.12220782693885,
48.1593744274586,
48.20486867941727,
48.25835564129731,
48.31948422515935,
48.387887925409686,
48.46318557549796,
48.544982131167615,
48.63286947883482,
48.72642726762714,
48.82522376357119,
48.92881672437982,
49.03675429325235,
49.148575910068345,
49.263813238324204,
49.38199110613446,
49.502628459594334,
49.625239326778704,
49.74933379063347,
49.874418968999976,
50,
50.125581031000024,
50.25066620936653,
50.374760673221296,
50.497371540405666,
50.618008893865536,
50.736186761675796,
50.851424089931655,
50.96324570674764,
51.07118327562018,
51.17477623642881,
51.27357273237286,
51.36713052116517,
51.455017868832385,
51.53681442450203,
51.61211207459031,
51.68051577484064,
51.74164435870269,
51.79513132058273,
51.8406255725414,
51.87779217306115,
51.90631302655559,
51.92588755234429,
51.936233321871015,
51.937086662998965,
51.92820323027551,
51.90935854011841,
51.88034846993819,
51.840989720274564,
51.791120239090816,
51.7305996074365,
51.659309385757524,
51.57715342020218,
51.484058108342616,
51.37997262380286,
51.26486909935761,
51.13874276813929,
51.00161206266512,
50.85351867147074,
50.69452755321206,
50.52472690817225,
50.344228107186346,
50.15316557807162,
49.951696649726884,
49.74000135413931,
49.518282186611756,
49.2867638245981,
49.04569280560701,
48.795337164708016,
48.53598603224466,
48.26794919243113,
47.9915566035775,
47.70715788075749,
47.41512174179946,
47.11583541754612,
46.80970402739315,
46.497149921177424,
46.17861198854701,
45.85454493700193,
45.525418539851685,
45.19171685538825,
44.85393741862609,
44.51259040700847,
44.16819778152742,
43.82129240474783,
43.47241713726835,
43.122123914190546,
42.77097280320363,
42.41953104592591,
42.068372084174065,
41.71807457285935,
41.36922138123378,
41.02239858423173,
40.678194445669945,
40.337198395085146,
40,
39.66718793541796,
39.3393489523529,
39.017066847202855,
38.70092143377615,
38.3914875197752,
38.0893338895364,
37.79502229481435,
37.509106455386295,
37.23213107123569,
36.964630848055556,
36.70712953778898,
36.460138995899975,
36.22415825703838,
35.99967263073258,
35.787152818708684,
35.58705405539872,
35.39981527316025,
35.22585829368798,
35.065587047052574,
34.919386819754834,
34.7876235331335,
34.670643053412306,
34.568770534618224,
34.48230979554525,
34.41154273188011,
34.35672876454522,
34.31810432525169,
34.295882380191,
34.29025199272783,
34.3013779258893,
34.329400285376764,
34.374434203756316,
34.436569566412594,
34.515870779778524,
34.612376582279865,
34.726099898359706,
34.8570277358732,
35.00512112706758,
35.17031511328693,
35.35251877346523,
35.551615296395035,
35.76746209668359,
35.99989097423161,
36.24870831699498,
36.51369534671404,
36.794608407220544,
37.09117929485854,
37.40311563048187,
37.73010127241896,
38.07179676972448,
38.427839854966706,
38.79784597573175,
39.18140886395714,
39.57810114214295,
39.987474965422834,
40.409062698416534,
40.84237762572369,
41.286914694861004,
41.74215129038767,
42.20754803791071,
42.68254963660851,
43.16658571886228,
43.659071735537715,
44.159409865414,
44.666989947215995,
45.1811904326653,
45.701379358930296,
46.226915338820085,
46.75714856703811,
47.29142184078137,
47.82907159294801,
48.36942893619241,
48.911820716049476,
49.45557057133282
],
"yaxis": "y"
}
],
"layout": {
"autosize": true,
"legend": {
"title": {
"text": "Chemical"
},
"tracegroupgap": 0
},
"template": {
"data": {
"bar": [
{
"error_x": {
"color": "#2a3f5f"
},
"error_y": {
"color": "#2a3f5f"
},
"marker": {
"line": {
"color": "#E5ECF6",
"width": 0.5
},
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "bar"
}
],
"barpolar": [
{
"marker": {
"line": {
"color": "#E5ECF6",
"width": 0.5
},
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "barpolar"
}
],
"carpet": [
{
"aaxis": {
"endlinecolor": "#2a3f5f",
"gridcolor": "white",
"linecolor": "white",
"minorgridcolor": "white",
"startlinecolor": "#2a3f5f"
},
"baxis": {
"endlinecolor": "#2a3f5f",
"gridcolor": "white",
"linecolor": "white",
"minorgridcolor": "white",
"startlinecolor": "#2a3f5f"
},
"type": "carpet"
}
],
"choropleth": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "choropleth"
}
],
"contour": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"type": "contour"
}
],
"contourcarpet": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "contourcarpet"
}
],
"heatmap": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"type": "heatmap"
}
],
"heatmapgl": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"type": "heatmapgl"
}
],
"histogram": [
{
"marker": {
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "histogram"
}
],
"histogram2d": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"type": "histogram2d"
}
],
"histogram2dcontour": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"type": "histogram2dcontour"
}
],
"mesh3d": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "mesh3d"
}
],
"parcoords": [
{
"line": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "parcoords"
}
],
"pie": [
{
"automargin": true,
"type": "pie"
}
],
"scatter": [
{
"fillpattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
},
"type": "scatter"
}
],
"scatter3d": [
{
"line": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatter3d"
}
],
"scattercarpet": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattercarpet"
}
],
"scattergeo": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattergeo"
}
],
"scattergl": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattergl"
}
],
"scattermapbox": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattermapbox"
}
],
"scatterpolar": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterpolar"
}
],
"scatterpolargl": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterpolargl"
}
],
"scatterternary": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterternary"
}
],
"surface": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"type": "surface"
}
],
"table": [
{
"cells": {
"fill": {
"color": "#EBF0F8"
},
"line": {
"color": "white"
}
},
"header": {
"fill": {
"color": "#C8D4E3"
},
"line": {
"color": "white"
}
},
"type": "table"
}
]
},
"layout": {
"annotationdefaults": {
"arrowcolor": "#2a3f5f",
"arrowhead": 0,
"arrowwidth": 1
},
"autotypenumbers": "strict",
"coloraxis": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"colorscale": {
"diverging": [
[
0,
"#8e0152"
],
[
0.1,
"#c51b7d"
],
[
0.2,
"#de77ae"
],
[
0.3,
"#f1b6da"
],
[
0.4,
"#fde0ef"
],
[
0.5,
"#f7f7f7"
],
[
0.6,
"#e6f5d0"
],
[
0.7,
"#b8e186"
],
[
0.8,
"#7fbc41"
],
[
0.9,
"#4d9221"
],
[
1,
"#276419"
]
],
"sequential": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
],
"sequentialminus": [
[
0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1,
"#f0f921"
]
]
},
"colorway": [
"#636efa",
"#EF553B",
"#00cc96",
"#ab63fa",
"#FFA15A",
"#19d3f3",
"#FF6692",
"#B6E880",
"#FF97FF",
"#FECB52"
],
"font": {
"color": "#2a3f5f"
},
"geo": {
"bgcolor": "white",
"lakecolor": "white",
"landcolor": "#E5ECF6",
"showlakes": true,
"showland": true,
"subunitcolor": "white"
},
"hoverlabel": {
"align": "left"
},
"hovermode": "closest",
"mapbox": {
"style": "light"
},
"paper_bgcolor": "white",
"plot_bgcolor": "#E5ECF6",
"polar": {
"angularaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"bgcolor": "#E5ECF6",
"radialaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
}
},
"scene": {
"xaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
},
"yaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
},
"zaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
}
},
"shapedefaults": {
"line": {
"color": "#2a3f5f"
}
},
"ternary": {
"aaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"baxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"bgcolor": "#E5ECF6",
"caxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
}
},
"title": {
"x": 0.05
},
"xaxis": {
"automargin": true,
"gridcolor": "white",
"linecolor": "white",
"ticks": "",
"title": {
"standoff": 15
},
"zerolinecolor": "white",
"zerolinewidth": 2
},
"yaxis": {
"automargin": true,
"gridcolor": "white",
"linecolor": "white",
"ticks": "",
"title": {
"standoff": 15
},
"zerolinecolor": "white",
"zerolinewidth": 2
}
}
},
"title": {
"text": "Initial System State"
},
"xaxis": {
"anchor": "y",
"autorange": true,
"domain": [
0,
1
],
"range": [
0,
299
],
"title": {
"text": "Bin number"
},
"type": "linear"
},
"yaxis": {
"anchor": "x",
"autorange": true,
"domain": [
0,
1
],
"range": [
32.544724436364255,
67.45527556363574
],
"title": {
"text": "concentration"
},
"type": "linear"
}
}
},
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0IAAAFoCAYAAAB69ksNAAAgAElEQVR4Xu2dCZgcRdmAv+y92d0kJAL+oALh0HCJqIBA5FY8gAACQTkiiSER86ugEIgEiAQCCqgBEk7DJQQQwqEihEMJCKioBAhyKSj+IuTa7H1k/65aep2dzOx091RXdc+8/Tx5Atmqr6rfryfZd6v6q2F93iVcEIAABCAAAQhAAAIQgAAEyojAMESojLLNrUIAAhCAAAQgAAEIQAACmgAixIMAAQhAAAIQgAAEIAABCJQdAUSo7FLODUMAAhCAAAQgAAEIQAACiBDPAAQgAAEIQAACEIAABCBQdgQQobJLOTcMAQhAAAIQgAAEIAABCCBCPAMQgAAEIAABCEAAAhCAQNkRQITKLuXcMAQgAAEIQAACEIAABCCACPEMQAACEIAABCAAAQhAAAJlRwARKruUc8MQgAAEIAABCEAAAhCAACLEMwABCEAAAhCAAAQgAAEIlB0BRKjsUs4NQwACEIAABCAAAQhAAAKIEM8ABCAAAQhAAAIQgAAEIFB2BBChsks5NwwBCEAAAhCAAAQgAAEIIEI8AxCAAAQgAAEIQAACEIBA2RFAhMou5dwwBCAAAQhAAAIQgAAEIIAI8QxAAAIQgAAEIAABCEAAAmVHABEqu5RzwxCAAAQgAAEIQAACEIAAIsQzAAEIQAACEIAABCAAAQiUHQFEqOxSzg1DAAIQgAAEIAABCEAAAogQzwAEIAABCEAAAhCAAAQgUHYEEKGySzk3DAEIQAACEIAABCAAAQggQjwDEIAABCAAAQhAAAIQgEDZEUCEyi7l3DAEIAABCEAAAhCAAAQggAjxDEAAAhCAAAQgAAEIQAACZUcAESq7lHPDEIAABCAAAQhAAAIQgAAixDMAAQhAAAIQgAAEIAABCJQdAUSo7FLODUMAAhCAAAQgAAEIQAACiBDPAAQgAAEIQAACEIAABCBQdgQQobJLOTcMAQhAAAIQgAAEIAABCCBCPAMQgAAEIAABCEAAAhCAQNkRQITKLuXcMAQgAAEIQAACEIAABCCACPEMQAACEIAABCAAAQhAAAJlRwARKruUc8MQgAAEIAABCEAAAhCAACLEMwABCEAAAhCAAAQgAAEIlB0BRKjsUs4NQwACEIAABCAAAQhAAAKIEM8ABCAAAQhAAAIQgAAEIFB2BBChGFP+nTkL5FePPC0vPLYoxlEIDQEIQAACEIAABCAAAQiEJZAaEXrh5b/L0VPPHXR/799ktDx8+6Vh71l22HeSZPZ98De/l2+fc4WcccqxcsJRnx0Uz5eZy847RT6zzydDjRVWhPx5ZA8y6ZiD5bvTJ4YaO2jjr3zjfPn3f1ZF4hh0jGLb+Ryz49x+9bmyw3ZbDvyxiXs54OhT9bNxy+XfK3ba9IcABCAAAQhAAAIQSDCB1IiQz/AHC26TRYsfkChi4sdIogj53+xny5h/v5/bf3f54ezpxh8lE/JgfFIZAVWu1JW9qqaERQlc5nNg4l4QoTizSWwIQAACEIAABCCQHAJlKULZ+IdaESomVUFXhPzxh1r5UbHKTYR8ftkrP35Obrzj13r1xl+pQ4SKeVrpCwEIQAACEIAABMqLQMmIUKZ0+KsIfipzrSb425/ybUfzV2D8FZnMb8bVN9x/fv7VDZ6U7HHCilCurXnZgwy1IqbE4KIrbh20xS/XXP155ruPzHv1Y2bOIx/Pz3pbB9X4/pXN0P/zoFv9gvJTcQvdi7+ClHkf2Vsrc7VR7TPv1+fvx4m6PbO8/prhbiEAAQhAAAIQgEDyCJScCGV/45prq1P2nw21IpRPhNQ3/ZnvE/nfiGd+0xz0G/nM95+CFFZQorfLjtts8B5L9n2p/1dX5ntU2X2HWkXJ9X5UrlUaXyAyt+9lCmauPw8jfUG3BQ51L2qOaszM97xy8Rlqa5yf41xSHCRvyfv4MyMIQAACEIAABCBQvgRKToSyvyHNJTLFilCux8X/xj/znZWgIqTi5Vp5UX+e612oXPeUPb4vV7lWX1R/v/hCPnkYSg4Vv4/tuO3AVr188qCkK5fE5PvzXFxzrfTkW4UJuzXOZ54pNvnuxW+bKx/qfoKucpXvXzXcOQQgAAEIQAACEEgWAUTIy0fYFaGhxCVzpSOMCGU+Ftlb+9TXMuPmkhw11p+ef2WD1Z/svtmPXz55yCVbfl/VR11+ZbV88hD2z4f6aOQTxUzxLSRC2dva/PEy5SbfnIfKZbYYJusjzmwgAAEIQAACEIAABHIRQIQiiJAvKpnfQOeSqagilJ2oXJXTMmP7YpRru1kuqcre2pWrfHa+ktX+3DJXZcIKj4nKbL4YZa445ROhzK2Hud5/CiJC+d5B8nnk2qrIXzkQgAAEIAABCEAAAsklgAiFFKF8W6TiFKFc7+Vkys/b767WJcULvafizzFTYgqtCBWKqR5tFyKkxg36vlO+1a1cuRxqRSh7xS25H2tmBgEIQAACEIAABCBQiAAi5BEaakUl+5vofN9UFyNC2WWgs5OWT1b8LWpqRefg/XYbdOiqmo/683wHxPqCk2tLnRo/TEnxuERoqJLhQbcHqnvJtzKXS4Syt/35uRjqHaFCHzK+DgEIQAACEIAABCCQPAKI0Hs5yffNfLb45DrzJ7NCWpR3hDLff8lXgjvXtrfMcbP7+V/LLlaQvYoy1Df4uarGKVyqz3MrXi9YLCGsIOUSQFWmPF/BBdU+877z3UuusuKZzDO3xhV6N0rNJ/tcI9VHXX4BiuR9zJkRBCAAAQhAAAIQgEA2gdSIUOZ7Hv5NZG7xyvdT/yBV4/Q31C//XY6eeu4An6HOEco+e0i9H3LiUZ+Vb59zxaCiBmHeEcp1f2oyhc6pGaoCW64CA7mkIvt9oELnCKl5BXmvplgRUuPkK3CQr0pbvnvJZqE47DxurD73KLsSXPb7QEOdI+Q/MPkOfeWvHAhAAAIQgAAEIACBZBJIjQglE5/bWQ21euF2ZowOAQhAAAIQgAAEIACBZBNAhJKdnyFnl+9g1RTfElOHAAQgAAEIQAACEICAFQKIkBXM5gfxV4NyHfBpfjQiQgACEIAABCAAAQhAoLQIIEKllU/uBgIQgAAEIAABCEAAAhAIQAARCgCJJhCAAAQgAAEIQAACEIBAaRFAhEorn9wNBCAAAQhAAAIQgAAEIBCAACIUABJNIAABCEAAAhCAAAQgAIHSIoAIlVY+uRsIQAACEIAABCAAAQhAIAABRCgAJJpAAAIQgAAEIAABCEAAAqVFABEqrXxyNxCAAAQgAAEIQAACEIBAAAKIUABINIEABCAAAQhAAAIQgAAESosAIlRa+eRuIAABCEAAAhCAAAQgAIEABBChAJBoAgEIQAACEIAABCAAAQiUFgFEqLTyyd1AAAIQgAAEIAABCEAAAgEIIEIBINEEAhCAAAQgAAEIQAACECgtAohQaeWTu4EABCAAAQhAAAIQgAAEAhBAhAJAogkEIAABCEAAAhCAAAQgUFoEEKHSyid3AwEIQAACEIAABCAAAQgEIIAIBYBEEwhAAAIQgAAEIAABCECgtAggQqWVT+4GAhCAAAQgAAEIQAACEAhAABEKAIkmEIAABCAAAQhAAAIQgEBpEUCESiuf3A0EIAABCEAAAhCAAAQgEIAAIhQAEk0gAAEIQAACEIAABCAAgdIigAiVVj65GwhAAAIQgAAEIAABCEAgAAFEKAAkmkAAAhCAAAQgAAEIQAACpUUAESqtfHI3EIAABCAAAQhAAAIQgEAAAohQAEg0gQAEIAABCEAAAhCAAARKiwAiVFr55G4gAAEIQAACEIAABCAAgQAEEKEAkGgCAQhAAAIQgAAEIAABCJQWAUSotPLJ3UAAAhCAAAQgAAEIQAACAQggQgEg0QQCEIAABCAAAQhAAAIQKC0CiFBp5ZO7gQAEIAABCEAAAhCAAAQCEECEAkCiCQQgAAEIQAACEIAABCBQWgQQodLKJ3cDAQhAAAIQgAAEIAABCAQggAgFgEQTCEAAAhCAAAQgAAEIQKC0CCBCpZVP7gYCEIAABCAAAQhAAAIQCEAAEQoAiSYQgAAEIAABCEAAAhCAQGkRQIRKK5/cDQQgAAEIQAACEIAABCAQgAAiFAASTSAAAQhAAAIQgAAEIACB0iKACJVWPrkbCEAAAhCAAAQgAAEIQCAAAUQoACSaQAACEIAABCAAAQhAAAKlRQARKq18cjcQgAAEIAABCEAAAhCAQAACiFAASDSBAAQgAAEIQAACEIAABEqLACJUWvnkbiAAAQhAAAIQgAAEIACBAAQQoQCQaAIBCEAAAhCAAAQgAAEIlBYBRKi08sndQAACEIAABCAAAQhAAAIBCCBCASDRBAIQgAAEIAABCEAAAhAoLQKIUGnlk7uBAAQgAAEIQAACEIAABAIQQIQCQKIJBCAAAQhAAAIQgAAEIFBaBBCh0sondwMBCEAAAhCAAAQgAAEIBCCACAWAVKjJv1a2F2rC1xNEoLa6Qhrrq2Vlc2eCZsVUChGoGDZMNtmoVv69qqNQU76eMAKbblQn767tlN71fQmbGdMZisDophpp6+yVjq5eQKWIQFN9lZ7tuvaeFM26vKe62Zj68gbg8O4RIQPwESEDEC2GQIQswjY4FCJkEKblUIiQZeCGhkOEDIG0HAYRsgzcwHCIkAGIEUMgQhHBZXZDhAxAtBgCEbII2+BQiJBBmJZDIUKWgRsaDhEyBNJyGETIMnADwyFCBiBGDIEIRQSHCBkA5ygEIuQIfJHDIkJFAnTYHRFyCL+IoRGhIuA57IoIOYQfcWhEKCI4A90QIQMQWREyANFiCETIImyDQyFCBmFaDoUIWQZuaDhEyBBIy2EQIcvADQyHCBmAGDEEIhQRHCtCBsA5CoEIOQJf5LCIUJEAHXZHhBzCL2JoRKgIeA67IkIO4UccGhGKCM5AN0TIAERWhAxAtBgCEbII2+BQiJBBmJZDIUKWgRsaDhEyBNJyGETIMnADwyFCBiBGDIEIRQTHipABcI5CIEKOwBc5LCJUJECH3REhh/CLGBoRKgKew66IkEP4EYe2JUKHTZolY0aPkOsvPSPiTN11W77idZk4fY7ctmC27DRurLGJIEIGULIiZABiwBAVa9dK9Z//KDVP/04q//GmVP7rLZGeoc9KGNbbI32VVbJ+zBhZv9nm0rf9DlLzqd1k5WZbSV9tXcCRaeaaACLkOgPRx0eEorNz2RMRckk/+tiIUHR2rnqaEqGTTr1Inn52xaDbGD2qSR5fMl//mQsRWvLAMpk171qZO3OKTDh478iIEaHI6OLviAjFw3hYa4tUL3/OE59npeaPv9e/V/3tNXODVVVJz7bbSfdHdpDuXT723u+7amHiSh4BRCh5OQk6I0QoKKlktUOEkpWPoLNBhIKSSk47EyK0w76TJFN6/LtTcrTp+zaSC8+a6kSETFFGhEyRjCEOImQG6rDODql9+CGpv2+JVD//F6la8eIGgfsaGqXrE7tJ9447S89HxnnS8j7pa2gINIGK/7wtVa+/JrUvPu/Ff06GvfZqzn693qpR5wEHScdBn9O/s2oUCG/sjRCh2BHHNgAiFBvaWAMjQrHijS04IhQb2tgCFytCSnZeef2fAys/+Sbqrwipr/srR/nkKXNlKXM72vgJM2Tv3XaSZc8sl1Vr1umhph1/qHxw8030yo9/+X1yCUz2ypXqP2PyEZJrReuFxxbpkIhQbI9f8YERoeIY1vzhGWm47iqpu/duUTLkX0pAunfaWbo98en6qLdis+NHpWfc9sUN5vX23xFa9X8rpdqTLb3S9OrLouZR9crLolai/EutDrUffpS0HTdJyxeXOwKIkDv2xY6MCBVL0E1/RMgN92JHRYSKJWi/f7EipFaDDv3MnnrVZ6hLidCrf39Li4sSD3Upsdl27AcG3htSMrJyVbPcs2iu/vr86+6ShTfdK76QqPZKgHzR8b+evQVP9VUxsgUmW9rU1y+75g49vvrat7921MA7QGq++eKYyhLvCBkgiQiFh6iEZ/itN8vwRdfq1Rn/6t5lV+n4zOf0Ly0e3vY101ehYglKjGp/84jU333n4Ll582mbeJy0H3Us2+dMJyVAPEQoAKSENkGEEpqYAtNChNKZN0QofXkrRoR80QjyDk6ud4TOvOBqefHlN3JKi09Syc/Rh+yn5clfEfKlK9dKjYqpVozUu0mZX1fxVMGDIHP1Jez2+x7dIA7FEhL2jCNCwRNSu+w3UvvQr6XuoQek6uWXdEe1ytNxwGel8yDv116fDh4sYstCIpQZtvaxh6X2Ue/XY0ul+oXn9Zd6N32/dO53oPfrAOnc90CkKGIewnZDhMISS057RCg5uQgzE0QoDK3ktEWEkpOLoDNJigj5hQ1yzdtfRconQplyo1aJcgnMa2/8S2+f81eXco3jrzhlfk21Z2tc0KfJQTtEqAB0r6pb/X13y/Cbb9ArLf6l3vVpO3GytB9xlNX3cMKI0MBk1T388j6pv+1mqXvkoYFKdes32VTavBUitVJkYtueg8c3NUMiQqlJ1QYTRYTSmTtEKJ15Q4TSl7diREjdbZitcdnlszNXhHwRKiQq6h2h7BUhEyKk7mP3XccNbNPL3JaHCCX4uUaE8ien/o5bpeniCwaqva0fOVJvLWs9cYozcYgkQhm3qIouDPfua7gnRZkFHTo+f4isO/UMrwLdrgl+WtM7NUQovblDhNKZO0QonXlDhNKXt2JFqFCxBCU7+arG5doaN9TWtWJWhFRm8m2NyyVhiFBKnmVEaMNEqTLXI2eeqqvAqatnq62l5X9Pk/ajJ1pd/cn1CBUrQpkxVYEFtUqk3nfyCz2oSnNrz73Qmeil5GMTepqIUGhkiemACCUmFaEmggiFwpWYxohQYlIReCLFipC/KpRdAc6XC7+QQqF3hFQcv3Jb5qqQkqXdd91enwNUjAipd3vUHFataR6ocOcXS1BFErIlSa0QqYutcYEfJTcNEaH/clcy0Hj5j6Tx0ou1GKiqa83f+760HXtcLIUPomTcpAj541esXOlVvlsoDVddLurQV1XkQd1z85nniNo+x1U8AUSoeIauIiBCrsgXNy4iVBw/V70RIVfko49rQoQyJSZzJpmrO0FEKF+czKpxUbfG+UUO/Op1/jz9OSrhuvfBJwemr95L8ivWsTUu+vMVe09EqB9x/c8Xy/Cf3aTfA1Ln/bR9+Xjv1wleCeyPxp6DMAPEIUL++KqgQv2dt2kWlf96S6+EtX/pGP2rZ+ttw0yTtlkEEKH0PhKIUDpzhwilM2+IUPryZkqE0nfn7mdM+WwDOSh3Ear+0x9l+C03aAka1tUpnQd+xlsNOUHaD+uvUZ+0K04R8u+15snHZfidiz0pWizD2lql+2Mfl7YjPSE6aqI+BJYrPAFEKDyzpPRAhJKSiXDzQITC8UpKa0QoKZkIPg9EKDgr0y0RIQNEy1WEKtas8eTnRu/9mBt10QC1+qFWgdq9VSBVYjqplw0R8u+9TlWa81aH6u+5S/9R574H6NWhti9NTMxWwaTmKXteiFBaMrXhPBGhdOYOEUpn3hCh9OUNEXKXM0TIAPtyFCF16KgqhqCKBahLVUxbe/7F0vuhLQwQjTeETRHKFKIR55w1UD1PnZfUfN4FVJgLkWpEKASshDVFhBKWkIDTQYQCgkpYM0QoYQkJMB1EKACkmJogQgbAlpMIqQIII74/WxquXajP0lHiowRIiVBaLhcipNl4vBpuXuSVE58rqgS3KqjQOmWaV1Bhtn6nimtoAohQep8QRCiduUOE0pk3RCh9eUOE3OUMETLAvlxEqPr552Sj6Sf1n53jfRPfcsq3vHNzTk/dN/HOROi9Z01VlWu48sfS9OMf9svkZpvL2nmXpkomDXxsQodAhEIjS0wHRCgxqQg1EUQoFK7ENEaEEpOKwBNBhAKjMt4QETKAtBxEaPitN8nI73xTl8Tu3nFnWb3g+tSek+NahPxHTgnlqG9/fWB7Ydvxk/TqGqtDuT+UiJCBv6wchUCEHIEvclhEqEiAjrojQo7AFzEsIlQEvCK7IkJFAlTdS12Emi6cI02XzNOk9Dfr3upFX22dAXJuQiRFhPy7b1h0raj3h4a1tuithquvvkG6PrGbGzgJHhURSnByCkwNEUpn7hChdOYNEUpf3hAhdzlDhAywL1URqnzjb9I4/zJR36ivHzVKWr7xbWn1tsP1VVcboOYuRNJESJFQZy81/PQaqbv/Hln/vo2l9aSp+pf6b65+AohQep8ERCiduUOE0pk3RCh9eUOE3OUMETLAvhRFqOaPv5fGyy+TuvuW6INAlQSp1aBSuJIoQvobfa+AgpKhhuuvloqV70rHIRO0DHWO37cUsBd9D4hQ0QidBUCEnKEvamBEqCh8zjojQs7QRx4YEYqMruiOiFDRCEtva1zdg7/SK0E1v1smXXvsqSWo4+AvGCCVjBBJFSGfjjpzaLgnRLXLfiM922ynZajtq19L/UpcsdlHhIol6K4/IuSOfTEjI0LF0HPXFxFyxz7qyKUuQjvsO0m22XJzuWfR3KiIYuuHCBlAW0orQuqA1MbLfyRVL78kHV84VEtQ1yd3N0ApOSGSLkKKVNUrf9UrQ2qFSFWWU6txrSedLN07fTQ5IC3PBBGyDNzgcIiQQZgWQyFCFmEbHAoRMgjTUqhSFqH5190lSx//o6xa0yxXXvht2WncWEtUgw2DCAXjNGSrkhChvj69FU79qli5UlpPnCytM74tPVsm64E1kC5Jgwj59zn8xuu1EKnS5aqAglodaj/6yyYwpC4GIpS6lA1MGBFKZ+4QoXTmDRFKX95KWYQOmzRLDhz/cfnTC6/Ipu/bSC48a2qiEoQIGUhH2kWo4t139Fa4xit+pKvBtXgC1HLKN6WvaYQBOskLkSYRUvRqfv+0lqH6O27VpbVbvW1ySohUhblyuhCh9GYbEUpn7hChdOYNEUpf3oyJ0FNPiXR02Aewxx4idRtWE16+4nWZOH2O3LZgtrz2xr/kkoWL5fEl8+3Pb4gRESED6UizCFW9tEKvAg2/7Wbp3fwD/ZXhvjbdAJXkhkibCCmSw1rWDWyVq/zHm9J54Gc8IZoqHZ/9fHJBG54ZImQYqMVwiJBF2AaHQoQMwrQYChGyCNvQUMZE6GMfE/nznw3NKkSYP/1JZJddNujgb4vz3w1S7wopKUrS9jhEKESe8zVNqwjVPvFbafBWguqW/lq6P/oxLUHth3/JAJFkh0ijCPlEVa6GX3eV1D30gBZXtTrU5r07tH5Eaa7eZT5JiFCyP1dDzQ4RSmfuEKF05g0RSl/ejInQsceK/Pvf9gHceqvI+9+/wbj+trgZk4/QXzvp1IsStz0OEXovbcpS/Wva8YeKnzT1ZyqRr/79Lf3lXFUv0ihCNX94RkafcIwu2dyz7Ydl1c/ulJ6ttrb/4XEwYppFSOEa1tkhTRfM0VsZ1dW5z/6yesF1sn6TTR3QtDckImSPtemRECHTRO3EQ4TscDY9CiJkmmj88YyJUPxTDTyCvy0uu8PoUU2J2h5X9iLkJ2ruzCky4eC9N0iwsteVq5oHSv4pKRozeoRcf+kZA23TJkL1994to6ZP1t9Qq2+iV914m373pFyutIuQnyd1COtGXh6VzCoJUjKk8lmqFyKU3swiQunMHSKUzrwhQunLWymKUPa2OD8rauEh3/fcLjJX9iKUvWyXnYTxE2bIadOOGZCkJQ8s2+BlrzSJkJKgjaaeqEsyt06aImvnXerVaq5y8ew5G7NUREgBrPzXWzqfNU89qfO47rtnybrTZjpjG+fAiFCcdOONjQjFyzeu6IhQXGTjjYsIxcs3juilKELq++ejD9lv0A4rxU4tMKgrc0EhDqZBY5a9CCkzVct0q9asG2Dmv8iVWe3Cf7Er15+lRYTq775DryAoCVp7/sXSOu0bQZ+TkmpXSiKkE+Plc+S5Z0nDwsv1/3Z8/hC9OlRqq3yIUHo/hohQOnOHCKUzb4hQ+vJWiiKUliyUtQjl2hZ35gVXy70PPikvPLZIgopQKpJ9550i6iU675tmuewykW99KxXTZpIhCNx8s8h0r+JfS4vIRz4icscdIjvuGCIATSEAAQhAAAIQgED5EECE3qtvnlnKz9+/uPUWmw3UP0/zilDdL++T0Sd9pexXgvyPdcmtCGX8fVW14kVdBKPqb6/pFaE186+S9kMPL4m/0VgRSm8aWRFKZ+5YEUpn3lgRSl/eWBFyl7OyFiGFPddLW5l/lusdoVnzrtUrRv6V5K1x9T9fLCMuOE8q3/i7rPvOmbJu5tnunraEjFzKIqQQKxlqvOpyGX7zIn1ArtoC2eL9Wr/xJgnJQLRpIELRuCWhFyKUhCyEnwMiFJ5ZEnogQknIQrg5IELheJlsXfYipF7aeuX1fw6U8lNb45Y9s3zg/9NcNa7uviUy4sI5UvXyS/qMoHVnzZa+mlqTz08qY5W6CKmkqANYG713hho8IapYvVrajzxGy1D3xz6eypypSSNCqU2dIELpzB0ilM68IULpyxsi5C5nZS9CCn3mOUG56pun8RwhdeBm09xzpfr556T1a9Nl3ZnnlMWhm0E+SuUgQj6H4Ytv0UUUqpf/Rbp2/5S0nvyN1G6VQ4SCPN3JbIMIJTMvhWaFCBUilMyvI0LJzMtQs0KE3OXMiQip7WaZVdoybz9zy5k7LOFGTtrWOHW+TJO3EqQOTW07/qvSfJYnQSnfFhUuI0O3LicRUiRqnnxcrw6pd8V6N9u8f6ucJ0RSWWkSa+yxEKHYEcc2ACIUG9pYAyNCseKNLTgiFBva2AIjQrGhLRjYugjlOpC04CwT3iBJIqTOk1ESVPvEb6X96C9L85mzpfeDH0o4QbvTKzcRUnQr33xDvzfUcNUVGnbr1K9rGerdYku78IsYDREqAp7jroiQ4wREHB4RigjOcTdEyHECIgyPCEWAZqiLdRFK2omyJjgmReBUYeoAACAASURBVISq//RH752g86T2kaXSftgR3na42dKzzXYmbrGkYpSjCOkE9vUNvDdU+c9/SMfnvqhXhzr3+nQq8osIpSJNOSeJCKUzd4hQOvOGCKUvb4iQu5whQgbYJ0GEql98wVsJOk/qfnW/dHz2815hhHOke4edDNxd6YUoWxF6L5V199+jhajmqSeke8edtQy1TTwu8YlGhBKforwTRITSmTtEKJ15Q4TSlzdEyF3OrIuQ2hp34PiPy4zJR7i7a8Mjuxahqtde0dvh6pf8XDr3PUCvBHV9/JOG77J0wpW7CKlMVv/lT1qG6u+4VdaPGqWLKLRMnyF9jU2JTTQilNjUFJwYIlQQUSIbIEKJTEvBSSFCBRElrgEi5C4l1kVoyQPL5JKFiwfKU7u7dXMjuxQhtcVJSZCqDta1x15eYQRPgvYcb+7mSjASItSf1IqVK3V5bSVEw9pape24Sfq9oZ5x2ycy64hQItMSaFKIUCBMiWuECCUuJYEmhAgFwpSoRoiQu3RYFyH1jtBQF1Xjgj8MFe/8Rx+WOvymn+rzYZq9Etmd+x8YPECZtkSEBidePT9Khqr+ukI699lfb5XrOOjgxD0diFDiUhJ4QohQYFSJaogIJSodgSeDCAVGlZiGiJC7VFgXIXe3Gt/ILlaEKpqb9TtBDdcskO7td9DnBKmX37kKE0CENmRU++hST4bmS+3DD0nP2G20DLWeNLUwTIstECGLsA0PhQgZBmopHCJkCbThYRAhw0AthEOELEDOMwQiZIC9bREa1tkhTRfMkcYrfiQ9W22tCyO0H/4lA3dSHiEQodx5rnr1ZWlYMF8abrhO+mpq+88b8n6t32TTRDwYiFAi0hBpEohQJGzOOyFCzlMQaQKIUCRsTjshQu7wOxEh9Z7QrHnXDrrruTOnyISD93ZHooiRbYuQeieo6ZJ5+nBMJUFpqPhVBF7jXRGh/EiVZCsZUlvlKt59xxPso7QQJaH4BiJk/KNgLSAiZA210YEQIaM4rQVDhKyhNjYQImQMZehA1kVo/nV3ycKb7pXbFsyWncaN1RNevuJ1mTh9jkw7/tBUVpOzKUJKgJQIrR8zxpOgc6X1xMmhk17uHRChwk9A/c8XaxlSZ1N1fXJ3XVWufcKRhTvG2AIRihFuzKERoZgBxxQeEYoJbMxhEaGYAccQHhGKAWrAkNZFaPyEGXL0IfttIDxKkG6/79FUVpOzJULDb71JRs04Wad2zfyrpO3Y4wOmmWaZBBChYM9D5b/eko1O+orU/OEZkaoqWXv+xdI6ZVqwzjG0QoRigGopJCJkCbThYRAhw0AthUOELIE2OAwiZBBmyFDWRUhVjcu1Dc7fLkfVuNwZrH3itzLmSK8YQk9P/zek3nYlrmgEEKEQ3LznbeTMU6VhUf9WViXfa+ddIn0NjSGCmGmKCJnh6CIKIuSCevFjIkLFM3QRARFyQb24MRGh4vgV09u6CLEiFD5dVStelPd98QCpWLtWWr75HWk+e074IPQYIIAIhX8Y1GrkyJmnybDWFunecWdZdeNi6f3QFuEDFdEDESoCnuOuiJDjBEQcHhGKCM5xN0TIcQIiDI8IRYBmqIt1EeIdoXCZq/jP27LxgXuL2qbUfujhsvr6W8IFoPUGBBChaA9F9fPPyegTjpHKN9+Q9SNHyuqrb5TOAw6KFixCL0QoArSEdEGEEpKIkNNAhEICS0hzRCghiQgxDUQoBCzDTa2LkJo/VeOCZbHy3/8nI86eKfV336HPCGr+/jzp2bK/wARXdAKIUHR21S8slwaviIJaIVIypIootE6bIetHjIgeNGBPRCggqAQ2Q4QSmJQAU0KEAkBKYBNEKIFJKTAlRMhdzpyIkLvbjWfkuIoljJh9pjRe+WPp3mVXWetJUNen0llePB7q0aMiQtHZqZ5qi2bDVZd7QjRf1MG+bV8+QQtR9w47Fhe4QG9EKFa8sQZHhGLFG1twRCg2tLEGRoRixRtLcEQoFqyBgiJCgTAN3SgOEWr0znIZMXum9HqHWaqVoPYjjjYwU0IoAoiQmedg+M9u1EJU/cLz0jl+X13Ao+OznzcTPEcURCg2tLEHRoRiRxzLAIhQLFhjD4oIxY7Y+ACIkHGkgQNaEyFVLU6dE6TOEBrqomqcSP2Sn2sJUu8FNZ93gbSc8q3ACaVhYQKIUGFGQVvUPv6Y3ipX9+tfSs9WW2sZap3cX+Ld9IUImSZqLx4iZI+1yZEQIZM07cVChOyxNjUSImSKZPg41kQo/NTS08PkilDN07/TElTzx9/rbyrXzpnn7UWqSA+MFMwUETKbpKq/vaZlqOG6q6Svuvq994a+Ib3v/x+jAyFCRnFaDYYIWcVtbDBEyBhKq4EQIau4jQyGCBnBGCmIdRHKd44QB6qKrsY10pOguvvvkfbDjpBmT4J6N/9ApMTSKT8BRMj80zGsu/u994YuF1Xko33CkVqIuj65u7HBECFjKK0HQoSsIzcyICJkBKP1IIiQdeRFD4gIFY0wcoDEiFC5H6iqvpFUK0EN1yyQrt320BLU9YndIieWjoiQi2dAbetU7w3V/P5p/fwqGWo//EtGpoIIGcHoJAgi5AR70YMiQkUjdBIAEXKCvahBEaGi8BXVOTEidOYFV8uyZ5bL40vmF3VDLjqb2BrX+JNLZMScs6X3gx/S2+E6Dpng4lbKYkxWhOJNs9rWqbbKqbLvvZu+v/+9IU+I+mpqihoYESoKn9POiJBT/JEHR4Qio3PaERFyij/S4IhQJGxGOlkRoVznBuWa/dyZU2TCwekrEV2sCNXfcauM9M4LqlizWktQ69SvG0kuQXITQITifzLUQcCN6r0h79ewrk5pPWmqFqKesdtEHhwRiozOeUdEyHkKIk0AEYqEzXknRMh5CkJPABEKjcxYBysilDnbfO8IGbsjB4GKEaHaZb/Rh6ZWL/+LtMw4VZrPOd/BHZTXkIiQvXw3/PQaLUNVr70iHQcdLK3TZ0jnp/eLNAFEKBK2RHRChBKRhtCTQIRCI0tEB0QoEWkINQlEKBQuo42ti5DR2SckWFQRUt8cKgmqe/BX0v6lifrQ1PUbb5KQuyrdaSBCdnNbt/TXWoZqH3tYurffQW+Ta/vKiaEngQiFRpaYDohQYlIRaiKIUChciWmMCCUmFYEngggFRmW8ISJkAGkUERrW1qolqOGG66Rzr0/rQ1O7d97FwGwIUYgAIlSIkPmvV720Qhqvmi/Db1okfU0jpOW994bWjxoVeDBEKDCqxDVEhBKXkkATQoQCYUpcI0QocSkpOCFEqCCi2BpYF6HlK16XidPn5L2hcjlQtemHF0rTvO/rdyaUBHV89vOxJZnAgwkgQm6eiGGtLf3vDXlV5SpWrZK2icfp94a6d9w50IQQoUCYEtkIEUpkWgpOChEqiCiRDRChRKZlyEkhQu5yZl2Exk+YIXvvtpPsvuv2csnCxQNV4g6bNEsOHP9xmTH5CHc0Io4cdkVo+C036FLZw7q6tQS1TpoScWS6RSGACEWhZq5P/e0/00JU/dyfpXPvfbQMdRz8hYIDIEIFESW2ASKU2NQMOTFEKJ15Q4TSlzdEyF3OrIuQXyxh6y02k6+fedmACKnKcpli5A5J+JHDiFDtI0u9Q1PPELVVaN1pM2XdmbPDD0iPogggQkXhM9K55nfLtAzV/eJe6d1iq/6tcl+bPmRsRMgIeidBECEn2IseFBEqGqGTAIiQE+xFDYoIFYWvqM7OREiVyVZS5G+FK4cDVatWvOhJ0EypfXSptH35BH1oaph3JIrKNJ0HCCBCyXgYKv/xpvfeUH+JbamslBaviIJaHerdbPOcE0SEkpG3KLNAhKJQc98HEXKfgygzQISiUHPbBxFyx9+6CKktcNtvt4VceNZUyfzvUj9QtWLNGr0dbvjPbpTO/Q7U5wX1jNveXebLeGREKFnJb1wwX783VPnPf0j7oYd7MjRDunbbY4NJIkLJyluY2SBCYWglpy0ilJxchJkJIhSGVjLaIkLu8mBdhLJvVa0K+ddtC2bLTuPGuqMRceQgW+OaLpwjTZd48vORcZ4EXSSd+x8YcTS6FUsAESqWoPn+dfffo7fK1Tz1hHTt+gm9MtR+xNGDBkKEzHO3FRERskXa7DiIkFmetqIhQrZImxsHETLHMmwk5yIUdsJJbF9IhBoWXatLZffVVOvtcFHOUEnifad1TohQMjNX/Zc/aRmqv+NWfZ6Wfm/I+9VXW6cnjAglM29BZoUIBaGUvDaIUPJyEmRGiFAQSslqgwi5y4d1EfKLJah3hErlGkqEqp9/Tt73hQNFlQ5WhRFUgQQutwQQIbf8hxy9p0dGeKunjT/+oW6mtsqtnXeprN9kU0QowWkrNDVEqBChZH4dEUpmXgrNChEqRCh5X0eE3OUEETLAPp8IVaxcKRsftLdUvvmG/oZu9fW3GBiNEMUSQISKJRh//7pf3icbTZ+sf4CwfuRIaT7/YunwCoxsslGt/HtVR/wTYASjBBAhozitBUOErKE2OhAiZBSnlWCIkBXMOQexLkJpPi8oX5pyipD3k+0xR35Rap/4rT4w8t1fLJW+hkZ3mWbkAQKIUDoehsp/vSUjv/O/Uvfgr/SEu/bZX6oWXC7/3mTLdNwAsxwgkHQRUsJd9bfXper116Tyb69J1d9fl4q33/ZEvDVSFtUKpqp+2Ps/m0nPth+W7p121quaabsQobRlrH++iFD68oYIucuZdRFavuL1QecHubt1cyPnEqGRM0+VhmsXyvoxY+Sdh5ZJ74e2MDcgkYoigAgVhc96Z/XO0MhzzpKK/7wtUlUlLad8S9adejo/WLCeiegDJkGEKtau7Zcc9UsJzz/flKpXXtb/rZ+tmC+1sqmlyPvBWO9WY/v/+yPbJ/rfBkQo5ocipvCIUExgYwyLCMUIt0Bo6yKUWSUu19z8c4XcIQk/crYIDb/1Jhk142T9TdvKxUuk0/tJNldyCCBCyclF0Jmob2KbfjBX/3BBvNVW9dP25rPnSPtRxwYNQTuHBGyKkBKd6uXPeYdWv6hXeXz5UVuV811KUnq32lp61K+xW0vvllt5WzJHSZ/351EutR1arWhWrF7lydZf9XzyyZbaKdCz7XbStcuu0rX7p6T7E7vpeSThQoSSkIXwc0CEwjNz3QMRcpcB6yLk7lbjGzlThGr+8IyMOexgGdbZoV/ybp0yLb6BiRyJACIUCZvzTrpq3FuvSM/UaV6Z7Sf1fLr22FPW/OAnnMnlPDtDTyBOEVLSobYg1z6yVP+eTzjU6rwSDC0822wrPd4qvS8/6mtxX0rmtRT9+VmpetVbifIO2K5+6QXJJWhqnu1HHKXPnFPPuKsLEXJFvrhxEaHi+LnojQi5oN4/pnURylc1bv51d8nt9z0qjy+Z745GxJF9EVL/AG+y7x76H+K2Y4+XNfOvihiRbnESQITipBtf7Mzy2dnb5dQPHNZ9d5YurMCVPAImRUgJRe1vHpGa3z7q/f6o3uqWeanVwm5vdUVtQVMrLb78JPXZUCKkhEgJUs3Tv5MaJXPePfqXup8Or9hO+xcnWJciRCh5n6UgM0KEglBKVhtEyF0+EiNCSx5YJrPmXStp3RpX+c9/yIjZM6X+3rul45AJ3qGp86T3gx9yl1lGzksAEUrnw5F9jlDVX1fI8NtuEbUVteLdd/RBrO0Tj5O2iV+RvuEN6bzJEp11MSJU+e//k+pn/yA1f/qj/r36T3+QiubmAVJKdLq93Kv8d3/s4/p3tS05rZd6lmt+94TUeocL1zy5zNtW95eBW+neYSfp2mu8dH5qb+nac2/vHdT3xXqbiFCseGMLjgjFhja2wIhQbGgLBk6MCJ15wdWy7Jnl1leEfAHLJpUpZKrS3at/f0s32WbLzeWeRXMHNVcrQuqdIPUNmaoM9J/HnkplhaCCT0uJNECE0pnIfAeqqi1Go777v2yXS3Baw4qQWh1Rqz61v31M51VtNfYvdcBul/ceTden99XvX6p3a9IsPoXSpu6/9tGlUn/XHYNWv9R2vo7PH+KtFB0unQccVChMpK8jQpGwOe+ECDlPQegJIEKhkRnrYEWE8slG9l3MnTlFbB+0quZ2ycLFeQXspFMvkpWrmgfkR0nRmNEj5PpLzxiY/tqLLhNVJU79Y/zu/Uv1P9JcySWACCU3N0PNLJ8I+X3YLpfcvBYSIbW9TW1zq3nice9dnwcHbQ1Td6W2uSnp6fz0ft72sE+VbcVAJf31v7xX6u5bIuqwbv9S2+davv5NafvqFFGiaOpChEyRtBsHEbLL28RoiJAJitFiWBGhzKnle0co2vSL71VIhMZPmCGnTTtmQNA2aP/YYyIHeT+N8ypZrbnsSmk7flLxkyJCrAQQoVjxxha8kAipgbOry6n3QlpP/oa0TprCKm1smSkcOFuE1Hsxtb9VKz797/moggeZlzpuoHOvT0vn/gfq39N4Bk9hKsW1UFXplPw33HDdAD/FqXXyyV6RnulG3pdDhIrLkaveiJAr8tHHRYSisyu2p3URKnbCpvvnWq3yt8WpM48mTp8jty2YLTuNG6uH3uDPNt5Y5N13dXU4VSWOK/kEEKHk5yjXDIOIkN9P/eR8xIXnSd0v79N/pH5KrqpwtZ58il5d4LJLYNN6keYHH5XqpQ/pLW+ZqxlqJkpY9TY3b8Wn03sHRp2xwxWcgHrOmy6eO8B14AcA3g8BiikSgQgFz0GSWiJCScpGsLkgQsE4xdGq7EUoG2rmVrhAIuSV9O0dv4+03f+rkt6nHsfD5ypmZeUwqamqlPbOHldTYNwIBIZ5n7WGukppaQ+et4rn/iJ1Z82UykcfHhhRfV67TpkhPZ/7PJ/ZCHkI2qXiNe/g0od+LZWe/FQ9/huRlpb/dq2rk95P7i493mq6ykfvbrsHDUu7IQio57z2onlSqXirHwA0Nkr3pJOk69TvSN+m7w/Nrr62Urp7+qSnd33ovnRwR6CmukIP3tVN3txlIdzITcOrw3WgtTECTkRIbTdbtWZdzptwXTXOlx81j0AitGiRtOxzQKR/ZIxlkUChCCBCoXAlpnEUEfInX/HiC1Jz1QKpvvlGkY7+F+/Xb7GFdH99hnQff4J3cOaoxNxnWicybO0aqXzkEal67BGp9ASo4o3B293W7/xRTz6/IL2f+pT0ePIjngxxxUNAiVDt3PMHhEix7poyNbQQIULx5CfuqIhQ3ITNx0eEzDMNGtG6COUqNhB0sjbaZZfxzvWOUHaZ78wDVW3MkTGKI8DWuOL4ueodZmtcvjlWe0JU++tfSN2vfynq8GN19Wy9rXR89vPScfDnvZLE413dXurGrVizRqr/+Hup+eMzUvPs7/V/V6xePXAfPR8eJ10f/6R0e7+a9t1L3tlqnPSu70vdfaZ5wnUPeM+6V1ih3vs1rK1Vet//P/p4h3bvV5Bnna1x6cw+W+PSlze2xrnLmXURSlqxBCU6mYe4ZotakKpxiJC7BzjKyIhQFGru+5gQof9+l97jVd+6Txqu/PGAEKmvdXzmc/o9IvW+CteGBFRRg7pHHnrvMNNHBld386pmdu6+53/LWmdUzyxUNQ7W8RJQ72Q1XnqRfuZVYR9V4bTtqGN1pbmecdvnHRwRijcvcUVHhOIiaz6uKhzTcN1CabrofPPBiRiIQNmLUOYZQYrY7ruOG1QaW/1ZkHOEAtGmUSIIIEKJSEPoSRgVoYzR1Zk1jVddIfV339H/TaJ3qW8OWydP835yfrh3aOWY0HMtlQ7DWlukVpe0VkUOHpWqV/466NbUYaad+3gFDjyBVEUO+hoac946IpSMJ0KVKW+89GIZ7lWbGxCiY4+T5jPPyVmZDxFKRt7CzgIRCkvMUXvv35sxx0zQBWSkj9VyR1kQ6yKkpOLA8R+XGZOPcHXPxsdlRcg40lgDIkKx4o0teFwi5E+44j9vez+Zu0oabvqpqP/Wl1rl8L7R79jvIH1oZTlUMxs4zPSRpVL79JMDcqhwKNHR4rP/QdLh/VJlroNciFAQSvbaKCFquOLH0nDzIp1fldfm8y7QZeYzL0TIXk5MjoQImaQZX6wR55wljVf8SNQ5YJVv/TO+gYg8JAHrIlTo3J405gsRSlfWEKF05WtAVLyqcZtsVCv/XtVf7CCua1hnhy67XX/3nVL78EOi/t+/1Dkt/Qd77qt/V/+ApflS5y5VP/8XqVr+nNQ++bjUeOKjtmpkXuqAaH3P+x3Yf1i0J4dhL0QoLDE77ZUQjZh1utQ96FU99a6uPfb0zsO7YkD4ESE7eTA9CiJkmqj5eOoMsI2mT9ZHO7z7i6Wy8QF7mR+EiIEIWBch9Y7QUJfrqnGBqGU1QoSiUHPXBxFyx76YkeNeEco1NyUFdQ/+0jv48zG9fWFgpei9xmqFSK+QqPNvvN/zbQ0r5r5N9FUyp85Wql7xglS9tEKqX3pBqj35yb4fNZaSOy0+6jDT/T9T1Dk0/twRIRNZjC9G/b13y4jvnS7qkFYluutOPUNavnmabPS+EdLW2SsdXb3xDU5k4wQQIeNIjQZUfxdvfPC+orYer7nsSmk7fpJQLMEo4lDBrItQqNmlpDEilJJEvTdNRChd+fJn60KEskmpf8CUENU96hUMeOp3+h+ygcv7BrJrl12l21s16f7w9tI7dqyWip4PetvHIqyiRM2SFh7vXZ6ql7zfvZfkq155eYN3e/zY6v2n7o/sID077SxdH/2YdHnFDoJudwszP0QoDC03bdXqYNP3z5aGRdfqCSjJX3/lldKy+96IkJuURB4VEYqMLvaO6nP2vgP3FrUaqwRIiZC6EKHY0ecdABEywB4RMgDRYghEyCJsg0MlQYQG3Y73boV6h6bmvWICNV7RBb/YQvZtqy11SooGfn3wQ/3//T/9f6bOMVo/cqTupn8q/17RhormtV5ltjX94To6pfI//x4Irdv19krFqpV6O5sqZKCkJ3Mrn99YrVT1bLuddI/bQbp32El/k9vtyY+al40LEbJB2cwYNU89KaO+fcqAPHeeeJKsnj3XyMqgmRkSpRABRKgQIUdfV8URvnKk3nLd7f3Q7J0HHhv4IRki5Cgn3rBORCizCtvcmVNkwsF7i9oyl6timzs0wUdGhIKzSkJLRCgJWQg/h8SJUNYtqNUhtd1MFRuo/uuLWmgq33xT/z5o5Sj8rQfv4a08aeHxRKfHEx4lPur/XRd5QISCpzAJLZVMN/74EmnySm4rKVfCvHbuxdJ++FFJmB5zKEAAEUrmI9J04RxpumSe/jy9s3TZoHdMESF3ObMuQpnn9GQeVjr/urvk9vseHXSmjzss4UZGhMLxct0aEXKdgWjjJ12EhrortWJT+Y833pMj7/f/+5f+/yrvXB5/RceXJfWPZF9dnQ63fsRI7yfxo/pDe5KTWZxB/3dlpfdeUoOsH+1tcdvpo1p61Mu3SbsQoaRlJNh8xvzrdamcPFmqfv+07qDO2Vo779JYtk8GmxGtghBAhIJQsttGvYe30Ulf0X+Pr/z5/d5xA58eNAFEyG4+MkezLkJq5ee2BbNlp3FjJVOEVDW5WfOuFYoluHsYymVkRCidmU6zCKWTuLlZI0LmWNqMpKrGtb+7Wiqvv16G37JIF9roGbuNtB13orR9ZVJZn7FlMw9hx0KEwhKLt33d0l97xUjOkKpXX5bm2d+Xlv89bYMBEaF4czBUdOsipOTnygu/vYEIsSLk7iEot5ERoXRmHBFKZ97UrBGhdOYus3y22uKpKsupn2yrS73jsOr6W1gdSmBqEaHkJCWzOEL7UcfK6gXX5ZwcIuQuZ9ZF6MwLrpZlzyzXW+D8FaGtt9hMJk6fI4d+Zk+58Kyp7mhEHJmtcRHBOeqGCDkCX+SwiFCRAB12R4Qcwi9i6FznCKkztkbOPFVv89QHsZ49R1qnTCtiFLqaJoAImSYaMV5mcYQdd9bnBeU7YgERisjYQDfrIqTm7G+Dy5z/tOMPlRmTjzBwS/ZDIEL2mRczIiJUDD13fREhd+yLHRkRKpagm/75DlRVP+VWMqQOhVSXendozYLrqSznJk0bjIoIJSMRA8URvIqg73rFEXq22jrvxBAhdzlzIkLubjeekRGheLjGFRURiotsvHERoXj5xhkdEYqTbnyx84mQP2L93XdoIVLFQNT5U6tuXCzd3k++udwSQITc8lej1z34Kxn95SP7iyPc8nPpPOCgISeFCLnLmXUROunUi+TpZ1dsUBSB8tnuHoJyGxkRSmfGEaF05k3NGhFKZ+4KiZC6q5pnnpKGa66U+rvvlF7vfKzWr31dWqZ+3eohwumkG9+sEaH42AaJrM4JGnm2Vxzh5Zdk3Rnfk3XfPatgN0SoIKLYGlgXIfVe0NGH7LfBNjiKJcSWYwJnEUCE0vlIIELpzBsilN68BREhdXcV/3lbGq++Uhq8X8PaWqX1pKnS6slQzzbbpffmUzxzRMhd8lRlRSVBtY8u9aorTpK1518kfY1NBSeECBVEFFsD6yKkVn78Q1Qz74ry2bHlmMCIUEk8A4hQetPIilA6cxdUhPy7U4UURs2YKuodInWIr6oq1zNu+3TefIpnjQi5S57aDqe2xan3gd557Hd5iyNkzxARcpcz6yLEipC7ZDNyPwFWhNL5JCBC6cwbK0LpzVtYEVJ3WukdEjz6hGOk+vnn9DeBa+ddIm3HHp9eCCmcOSLkJmmNP/6hjPj+bP3cv/PAY6F+CIAIucmZGtW6CKktcAtvunfgUFU1ieUrXtfls9NaOY5iCe4e4CgjI0JRqLnvgwi5z0HUGbAiFJWc235RREh/Y9HZISNmnS4Ni67VN6BEaO0Pfyx9tXVub6hMRkeE7CdaF0fwfgCgrpWLl0jnPvuHmgQiFAqX0cbWRUjNPlf57Fzb5YzeaYzBEKEY4cYQGhGKAaqFkIiQBcgxDYEIxQQ25rBRRcifliqvPeo735RhrS26mtzqn3pbDllnvQAAIABJREFU5YYoIRzz7ZRNeETIbqrV6uf7vnCgfs6bz7tAWk75VugJIEKhkRnr4ESEjM0+IYEQoYQkIuA0EKGAoBLWDBFKWEJCTAcRCgErQU2LFSF1K1UrXpTRU0/Qv6stQ2t+dIW0H35Ugu6y9KaCCNnLqSoUsvGBe+sDhtXK55r5V0UaHBGKhM1IJ0TIAEZEyABEiyEQIYuwDQ6FCBmEaTkUImQZuKHhTIiQlqHXXpGGq66Qhuuvlr6aGl1iW1WV6938A4ZmSphMAoiQneeh6vVXZcS5s0QVCVFy33zO+dL7gQ9GGhwRioTNSCcnIqQKJqxasy7nDbzw2CIjN2YzCCJkk3bxYyFCxTN0EQERckHdzJiIkBmOtqOYEiE9756e/hLb3plDlf94Uzq+cKiWoc69Pm37tkp+PEQo/hRXrFkjI847S4bftEi/D6QkqHvnXSIPjAhFRld0R+sidNikWTJm9Ai5/tIzip58UgIgQknJRLB5IELBOCWtFSKUtIwEnw8iFJxVkloaFaH3bkz99FydN1S77DfSvcOOenVInbfCZY4AImSOZb5II84/Rxp/9AP97puSoM79DixqUESoKHxFdbYuQvnOESrqLhx3RoQcJyDk8IhQSGAJaY4IJSQREaaBCEWAloAucYiQuq3qF1/wZOgKGX7zIukb3qBXhlq8X+s32TQBd53+KSBC8eaw8cqfeKtBs6R30/drCWo/sr9aXDEXIlQMveL6IkLF8dO9ESEDEC2GQIQswjY4FCJkEKblUIiQZeCGhotLhPzp6XNXLpyjt82pLXKrr7kBGTKQO0TIAMQ8IWoffkjGfOVI/cyqwgimzshChOLLWaHI1kVIbY07cPzHZcbkIwrNLTVfR4RSkyo9UUQoXfnyZ4sIpTNvataIUDpzF7cIKSo1Tz0pG009UVfd6t1sc1l142Lp3mXXdAJLyKwRoXgSUfXKX70KceN1mex1p82UdWfONjYQImQMZehA1kVInSF0ycLF8viS+aEnm9QOiFBSM5N7XohQuvKFCKUzX5mzRoTSmUMbIqTIqBLEY758pFT/+Vl96OraeZdK2/GT0gktAbNGhMwnQZfJPnhfqXzzDV0hTq1emrwQIZM0w8WyLkLqHaGhLqrGhUsgrcMTQITCM0tCD1aEkpCFaHNAhKJxc93Llgjp+/S2Go2ceao0LLpW/68SoTU/+IlXe7vKNYbUjY8ImU2ZWgEac+QXpeYPz0jXHnvKyp/fr4Xd5IUImaQZLpZ1EQo3vXS0ZkUoHXnyZ4kIpStf/mwRoXTmTc0aEUpn7qyK0HuIGq5ZII0LL5fKN/4m7YceLq3TZkjXbnukE6CjWSNC5sAPW9eszwpquOE66Ry/b3+Z7Bi2biJC5nIWNhIiFJZYjvaIkAGIFkMgQhZhGxwKETII03IoRMgycEPDuRAhNfW6X90vDZ4M1T7xW+na9RPSOn2G3o7EFYwAIhSMU5BWTRecJ02XXqRLvTfPniudBxwUpFvoNohQaGTGOjgRIfWe0Kx5/cvf/jV35hSZcPDexm7MZiBEyCbt4sdChIpn6CICIuSCupkxESEzHG1HcSVC6j6rl/9FGhfMl/rbf6bLFLdO+4ZeHeqrrraNIXXjIUJmUqZkfMSc78n60WP6y2QfdayZwDmiIEKxoS0Y2LoIzb/uLll4071y24LZstO4sXqCy1e8LhOnz5Fpxx+aympyiFDB5yxRDRChRKUj8GQQocCoEtcQEUpcSgJNyKUIqQlWrFrlrQzN11vlhrW1SuuUaVqIerbs/96BKzcBRKj4J6P+jlu9s4K+5z2DK72VoPP1cxfnhQjFSXfo2NZFaPyEGXL0IfttIDxKkG6/79FUVpNDhNw9wFFGRoSiUHPfBxFyn4OoM0CEopJz28+1CPl3P/zG6/XqkCpf3HHwF/RWOXXuEBciFMczUPvIUk+CzpLqF56XdaeeIevOOieOYQbFRIRiR5x3AOsipKrG5doG52+Xo2qcu4ehXEZGhNKZaUQonXlTs0aE0pm7pIiQoqcOsmz0VodqH10q3TvurGWo7ZivpBNszLNmRSg64Oq//EkXR6h9/DFpPXGyNJ87V/qaRkQPGLAnIhQQVAzNrIsQK0IxZJGQoQggQqFwJaYxIpSYVISeCCIUGlkiOiRJhBSQqpdWSONV82X4TYu89zZGS+vJ35AWT4j6hjckgldSJoEIRctE5T/e9FaCZkn9kp9LxyETZK0nQb1bbBUtWMheiFBIYAabWxch3hEymD1CRSKACEXC5rwTIuQ8BZEngAhFRue0Y9JESMEY1rJOvzOk3h2qWLNG2k44SVrUe0PbfcQpqyQNjgiFz4Z6B02Xyb7+aunac7xeCVIVC21diJAt0huOY12E1BSoGucu4YzsbbGorpDG+mpZ2dwJjhQRQIRSlKysqSJC6cxdEkXIJ6m2ym009QSpWLtWerb9sKy68Tb9O5cIIhT+KWi6cI40XTJP+hoa5Z2lj1t/lhCh8Dkz1cOJCJmafFLiUCwhKZkINg9EKBinpLVChJKWkeDzQYSCs0pSyySLkOJU+eYbMvqEY6T6+ef0N7Crr7lBOj7zuSQhdDIXRCgcdlUme+T3Tvf2XlbJysVLpHOf/cMFMNAaETIAMWII6yJ00qkXydPPrpDsogiqiMLuu46T6y89I+KtuOuGCLljH2VkRCgKNfd9ECH3OYg6A0QoKjm3/ZIuQorOsNYWGfWdb4oqd6yudafNlHXfPUt/U1uuFyIUPPP1994tG53UX3RDibSrg3sRoeA5M93SughRLMF0CokXlgAiFJZYMtojQsnIQ5RZIEJRqLnvkwYR8inpn+qf6wlQT490HnCQrL76Rlk/cqR7iA5mgAgFg64laOqJ+plZe/7FsZ8VNNSsEKFgOYujlXURonx2HGkkZhgCiFAYWslpiwglJxdhZ4IIhSWWjPZpEiFFrPaJ38pGXztRKv7ztvRstbWs/uktutR2uV2IUOGMJ0mC1GwRocI5i6uFdRFiRSiuVBI3KAFEKCipZLVDhJKVjzCzQYTC0EpO27SJkCJX8/TvpOGqy0V9o9szdhv9U/7Wk6YmB6qFmSBCQ0NuuPpKabr4fF1oQ22jXHf6LAtZGXoIRMhdCqyLEOWz3SWbkfsJIELpfBIQoXTmTc0aEUpn7tIoQop05Vv/1OW1VZntvuoaffiqKrG9fuNN0pmIkLNGhPIDa7hmwWAJUu+TJeBChNwlwboIqVulfLa7hDMyIpTWZwARSmvmEKG0Zi6tIuTzVu8NNXqrQ+qgzPYjj9GHr3bvsmta0xF43ohQblRJlSA1W0Qo8ONtvKETETJ+F4YCnnnB1XLvg0/KbQtmy07jxg5EPWzSLHn172/p/99my83lnkVzB41I1ThDCbAUhhUhS6AND4MIGQZqMRwrQhZhGxwq7SKkUNTdf4+3MjRfap56Urr22NNbGZohHV88zCCl5IVChDbMScO1CwevBH3nzEQlDhFylw5E6D32apXqp7f9SgtPpgipct8rVzUPyI+SojGjRwwq840IuXuAo4yMCEWh5r4PIuQ+B1FngAhFJee2XymIkCJY/ednpXHBfKn/+WLp/eCHpOVk770hb6tcqV6I0ODMNlx31WAJ8kqsJ+1ChNxlBBF6j72qZqcEaOL0OYNESBV3OG3aMTLh4L11SyVMlyxcLI8vmT+QNUTI3QMcZWREKAo1930QIfc5iDoDRCgqObf9SkWEFMWKd9/RMtTg/RrW3aW3ybV6QtS7+QfcQo5hdETov1DTIEFqtohQDB+EgCERIQ+UWuX56sTPydZbbDZIhJaveH0DMcr1Z4hQwKctIc0QoYQkIuQ0EKGQwBLUHBFKUDJCTKWURMi/7YafXqNlqOr1V6X90MO9laEZ0rXbHiGoJL8pItSfo4brr5ami86XYc1rpUVVhzv1jMQmDxFyl5qyFyH1XtDb767WW92yJSeoCLlLHyNDAAIQgAAEIBCKwC9+IXLZZSIPPyyy++4i3/qWyMSJoULQOOEErrxS5NxzRdasETnnHJFZ7ktkJ5xY2U6vrEUoe5tbVBFiRShdnx9WhNKVL3+2rAilM29q1qwIpTN3pbgi5Gei+sUXdInt4T+7UXo3fb8usa22yvVVV6czWRmzLvcVIbXq568EqXOCWr59euJzyoqQuxSVvQjNmndtTvrTjj9UZkw+QnK9I6T6vPDYooF+iJC7BzjKyIhQFGru+yBC7nMQdQaIUFRybvuVsggpshXNzf3nDan3htY1S+uUafrdod4ttnILvsjRy1mEGhZd2y9B6rDU0z0J+tZ3i6RppzsiZIdzrlHKWoSygeTaCkfVOHcPZ1wjI0JxkY03LiIUL984oyNCcdKNL3api5BPbvgtN+gS21UrXpSOg7+gV4c69/p0fGBjjlyuIpRWCVKPAyIU84diiPCIUAacXCKkvsw5Qu4e0DhGRoTioBp/TEQofsZxjYAIxUU23rjlIkKKYtXfXpPRJxyjZWj9mDGy+uobpHOf/eMFHFP0chQhdXjuyHPPEunpkbXnX5y68uiIUEwfhgBhEaEAkAo1YWtcIULJ+joilKx8BJ0NIhSUVPLaIULJy0mQGZWTCCkew1pbZNSMk6X+3rs9M6oS9X7JugSeOVMod+UmQk0XzpGmS+ZpLGmUIDVvRKjQUx3f1xEhA2wRIQMQLYZAhCzCNjgUImQQpuVQiJBl4IaGKzcR8rE1XvEjGfH92Xp1QZXYXjP/KulraDRENf4w5SJCWly/802pv+NWLa5rfvATaTt+UvyAYxgBEYoBasCQiFBAUEM1Q4QMQLQYAhGyCNvgUIiQQZiWQyFCloEbGq5cRUjhq/3NI7LR1BOlYuVK6Rm3vay6/hbp2fbDhsjGG6YcRKjyzTf0Vsbq55/Tkrr6mhuk4zOfixdsjNERoRjhFgiNCBlgjwgZgGgxBCJkEbbBoRAhgzAth0KELAM3NFw5i5BCWPO7ZV4Rhcul7hf3Ss9HxkmLd/hq23GTDNGNL0ypi9DwxbfoQ3GVBHXud6DOS+cBB8UH1EJkRMgC5DxDIEIG2CNCBiBaDIEIWYRtcChEyCBMy6EQIcvADQ1X7iKkMA7r7JCR3var4bfepKm2fPM7uixzX22dIcrmw5SyCKl3gZp+cIHetqhWgNRKUJq2LebLNiJk/nMQNCIiFJTUEO0QIQMQLYZAhCzCNjgUImQQpuVQiJBl4IaGQ4T+C1KVZh4581T9DbjaKrd6wfXSvePOhkibDVOKIjSokIWHa93ps/SvUrkQIXeZRIQMsEeEDEC0GAIRsgjb4FCIkEGYlkMhQpaBGxoOERoMsvrPz8pGXztRl9r2q8q1fONbiVsdKjURyixtXgrvA+X6eCJChv7SihAGEYoALbsLImQAosUQiJBF2AaHQoQMwrQcChGyDNzQcIjQhiDVysSIc84StUKkrp6ttpa18y5N1DsqpSRCtQ8/5BWtOEEq1q7VxSpW3XhbaopWhPkYIkJhaJltiwgZ4IkIGYBoMQQiZBG2waEQIYMwLYdChCwDNzQcIpQfZP09d3nvDd0otUsf1KtDbccer3917baHIfrRw5SCCCnhbPQKIjQsnC8Va9Z4ZbG/Ki3TZ0jPdh+JDibBPREhd8lBhAywR4QMQLQYAhGyCNvgUIiQQZiWQyFCloEbGg4RKgDSe19InTnUdOnF+jBWtW1LFVJoOfkbWo5cXWkXIbX6M/LbX0/9wbZh8o8IhaFlti0iZIAnImQAosUQiJBF2AaHQoQMwrQcChGyDNzQcIhQMJDqTJuR3ztd6n55n+6giimowz279tgzWADDrdIsQjVPPSkbfX2yKKbrR46U1VffmKhth4ZTNRAOEYqLbOG4iFBhRgVbIEIFESWqASKUqHQEngwiFBhV4hoiQolLSaAJIUKBMA00qnvwV7qynPomXl1qq1zzuRfI+jFjwgUqsnUaRUitAjX9YK40XLtQV+ZTFflW/9Q7xNZ7B6scLkTIXZYRIQPsESEDEC2GQIQswjY4FCJkEKblUIiQZeCGhkOEwoNU5w41XXyB3jKnvqFXqxrrzv6+tKqDWC1tl0uVCHmMht96s4w4/2ypWLmyvxrfqWfoX7Z4hc+y+R6IkHmmQSMiQkFJDdEOETIA0WIIRMgibINDIUIGYVoOhQhZBm5oOEQoOsjaR5bqYgr1d9+pg3R8/hC9QtTxuS9GDxqwZ1pEqO6+JTL8tpul7te/1HfWftgR0n7MV/RBqeV2IULuMo4IGWCPCBmAaDEEImQRtsGhECGDMC2HQoQsAzc0HCJUPMj6O27V7w/5qx2tk6bo1Y71m2xafPA8EZIuQhX/eVtGeiXIFRt19X5oC2k+e460H35UbEySHhgRcpchRMgAe0TIAESLIRAhi7ANDoUIGYRpORQiZBm4oeEQITMg9fsv3z9bGm5epLfL9dXWSdtXp0jL1FO0BJi+kipCattgw8LL/1tlz+PQ8s3T9C/FpJwvRMhd9hEhA+wRIQMQLYZAhCzCNjgUImQQpuVQiJBl4IaGQ4QMgXwvTNXfXtOHsfrV5dQ7MO3elrlWT4hMVphLogipQhIjZp0uioG61OpP83kXSO9mm5uFnNJoiJC7xCFCBtgjQgYgWgyBCFmEbXAoRMggTMuhECHLwA0NhwgZApkVpvr556ThuoVSf/ttolZJ1NW9y67e+UOn9G8PK7KoQpJESN3rCG81rPbhh/rv06sGt3bepUbFL54s2Y2KCNnlnTkaImSAPSJkAKLFEIiQRdgGh0KEDMK0HAoRsgzc0HCIkCGQecKo94YaFl3jSdFVot6bUZd6d6h18snSOulrkctuOxchb/tf7W8elcafXCK1T/y2/768EuLN3/u+VzDiuKJFL96suImOCLnhrkZFhAywR4QMQLQYAhGyCNvgUIiQQZiWQyFCloEbGg4RMgSyUBhPHOrvu9sruf1jqf7zs7q1emem/eiJ3ntE39AHtIa5XImQmvvwO2+T+rvu+K/YeeXD2zypa/nf03Qpca7cBBAhd08GImSAPSJkAKLFEIiQRdgGh0KEDMK0HAoRsgzc0HCIkCGQIcLU/OEZabjqCi1GqrCCuro+sZu0f2miLsEd5J0amyJU+a+3vDLhN0n9nYul6pW/Dtxpz7YfltYTJ0vb8ZOkr6ExBIHybIoIucs7ImSAPSJkAKLFEIiQRdgGh0KEDMK0HAoRsgzc0HCIkCGQEcIowVBb5oZ7W+dU1Tn/UiLUtfue0rXXeC1I6p2b7CtuEVJb+lTBh3pv9cff+qbmoLa/qXeclLSpuXEFJ4AIBWdluiUiZIAoImQAosUQiJBF2AaHQoQMwrQcChGyDNzQcIiQIZBFhBnW2qKlo+5X93vv3TwySIq0fHjbzbr2+rQnR5/qFySv6EJTU38p6nXt/StKxVxq/KpXXpbqFS9I1Usr9BxUAYSByyvsoA5AbZt4nHQecFDZl8GOyhoRikqu+H6IUPEMBREyANFiCETIImyDQyFCBmFaDoUIWQZuaDhEyBBIg2HU9rPaJx6Xmqd/5/16UirffGNQdLUNbf2uH5fe3XaX9v/5gN5K19fQoNuo9456Mw5yXb/ppvrPVOW6irfflormtXp7W/Xzy6XKE5/ql17cIL4fR634dBwyQa8AqZUgruIIIELF8SumNyJUDL33+iJCBiBaDIEIWYRtcChEyCBMy6EQIcvADQ2HCBkCGWMYJUI1f3xGapQcee8XDVqtMTCulqexY6Xbe+en5yPbD2zJK/cDUA2gHRQCETJNNHg8RCg4q7wtESEDEC2GQIQswjY4FCJkEKblUIiQZeCGhkOEDIG0GEa9TzTy909IxeuvSs/Lr3oHmL4+UHRBrfxUvlemW01JrQLpc4y87W1q5UitFvV+cAvvvaOdRBU76Nl2O/07V/wEEKH4GecbAREywB4RMgDRYghEyCJsg0MhQgZhWg6FCFkGbmg4RMgQSMth4i6WYPl2ymI4RMhdmhEhA+wRIQMQLYZAhCzCNjgUImQQpuVQiJBl4IaGQ4QMgbQcBhGyDNzAcIiQAYgRQyBCEcFldkOEDEC0GAIRsgjb4FCIkEGYlkMhQpaBGxoOETIE0nIYRMgycAPDIUIGIEYMgQhFBIcIGQDnKAQi5Ah8kcMiQkUCdNgdEXIIv4ihEaEi4Dnsigg5hB9xaEQoIjgD3RAhAxBZETIA0WIIRMgibINDIUIGYVoOhQhZBm5oOETIEEjLYRAhy8ANDIcIGYAYMQQiFBEcK0IGwDkKgQg5Al/ksIhQkQAddkeEHMIvYmhEqAh4DrsiQg7hRxwaEYoIzkA3RMgARFaEDEC0GAIRsgjb4FCIkEGYlkMhQpaBGxoOETIE0nIYRMgycAPDIUIGIEYMgQhFBMeKkAFwjkIgQo7AFzksIlQkQIfdESGH8IsYGhEqAp7DroiQQ/gRh0aEIoIz0A0RMgCREBCAAAQgAAEIQAACEIBAugggQunKF7OFAAQgAAEIQAACEIAABAwQQIQMQCQEBCAAAQhAAAIQgAAEIJAuAohQuvLFbCEAAQhAAAIQgAAEIAABAwQQoYgQD5s0S179+1u69zZbbi73LJobMRLd4iKw5IFlMmvetRuEf+GxRQN/Rh7joh8+7vIVr8vE6XPktgWzZadxYwcFKJSnQl8PPxt6hCGQL3d8BsNQtNf2pFMvkqefXTFowMy/F9UXCn2mCn3d3t2Uz0iF8sbnrXyeBe7UHAFEKAJL9ZfRylXNA/Kj/kEYM3qEXH/pGRGi0SUuAuofhUsWLpbHl8zPOQR5jIt8+LjjJ8yQVWvW6Y7ZIlQoT4W+Hn429AhDYKjc8RkMQ9JeW5WzzL8Xz7zgaln2zPKBPyv0mSr0dXt3Ul4jFcobn7fyeh64WzMEEKEIHNVfRqdNO0YmHLy37l3oL58IQ9DFAIFCeSGPBiAbDJFvVaFQngp93eAUCZWHwFArQkP9MILcJeORys5fobwU+noy7qr0Z5GdN/7NK/2cc4fmCSBCIZnm+gd/qC09IcPT3CCBXNsE/O0f5NEgaEOhguYks50aOns7HZ9HQwkJESbM1jg+gyHAWmo6/7q75Pb7HtUrQoU+h3zmLCUlwDCZefN/KJu9HZzPWwCQNClrAohQyPQX+kci+92GkOFpHiOBzO0c5DFG0BFDB80JIhQRcIzdgsonn8EYkxAxtJ+7uTOn6F0OhT6HiFBE0Ia7ZectV3g+b4ahE64kCSBCIdNa6B8JRCgkUIvN/dypn5CRR4vgAw4VNCeIUECgFpsFFSE+gxaTEmAoPx/Tjj9UZkw+Qvco9DlEhAKAjblJrrzlGpLPW8yJIHxJEECEIqQx1/5otRydXXUnQmi6xEjA3yrn54k8xgg7Qugw7whlft7IYwTYhrsEFSE+g4bBFxHOz0WuKo2FPlOFvl7EtOhagMBQecvuyueNxwkChQkgQoUZbdCCijkRoDnokl1hJ7u6H3l0kJQhhsz3zXShPBX6erLusjRnM5TEZlYn4zOYjPwXeqm+0Geq0NeTcZelN4tCeePfvNLLOXcUPwFEKCJjzlCICM5it8wcqWF333XcBiXOyaPFhAwxVGYJZtVs9KimQeV9C+Wp0NeTcZelOYuhcsdnMHk596U118z894TU1wp9pgp9PXl3nu4ZBckbn7d055jZuyGACLnhzqgQgAAEIAABCEAAAhCAgEMCiJBD+AwNAQhAAAIQgAAEIAABCLghgAi54c6oEIAABCAAAQhAAAIQgIBDAoiQQ/gMDQEIQAACEIAABCAAAQi4IYAIueHOqBCAAAQgAAEIQAACEICAQwKIkEP4DA0BCEAAAhCAAAQgAAEIuCGACLnhzqgQgAAEIAABCEAAAhCAgEMCiJBD+AwNAQhAAAIQgAAEIAABCLghgAi54c6oEIAABCAAAQhAAAIQgIBDAoiQQ/gMDQEIQAACEIAABCAAAQi4IYAIueHOqBCAAAQgAAEIQAACEICAQwKIkEP4DA0BCEAAAhCAAAQgAAEIuCGACLnhzqgQgAAEIAABCEAAAhCAgEMCiJBD+AwNAQhAAAIQgAAEIAABCLghgAi54c6oEIAABCAAAQhAAAIQgIBDAoiQQ/gMDQEIQAACEIAABCAAAQi4IYAIueHOqBCAAAQgAAEIQAACEICAQwKIkEP4DA0BCEAAAhCAAAQgAAEIuCGACLnhzqgQgAAEIBCSwGGTZsmY0SPk+kvPCNmT5hCAAAQgAIENCSBCPBUQgAAEUk7gzAuulnsffHKDuzj0M3vKhWdN1X++5IFlMmvetTJ35hSZcPDeqbxjRCiVaWPSEIAABBJLABFKbGqYGAQgAIFgBJQILXtmuTy+ZP5Ah+UrXpeJ0+fItOMPlRmTjwgWKOGtEKGEJ4jpQQACEEgZAUQoZQljuhCAAASyCeQSIdVm/IQZsvduO+lVIV+MblswW3YaN1Z8qVDtnn52hQ45elTTIJnKRTpIP9Vm++22GFiNUnFOOvUiWbmqWe5ZNFeH9eemBG7VmnX6z5S0fXDzTfTKlX/581X/H2Rsfyz/ntT/F4qR+XWeLghAAAIQKB8CiFD55Jo7hQAESpRALhGaf91dsvCmewckIJcIvfr3twatGCk52XbsB4Z8B0fJSKF+QUVICZAvIf58M2VMxVGXL0+5xs5uky1cftwXHlukY+WKUaKPBbcFAQhAAAIFCCBCPCIQgAAEUk4g3ztCmVKRb0Uos/CAivPiy28MiEcuLLm2p2X3CypC/mqVGid7furPsgUv19j+u09KqNSltgNmr/AowTv6kP30FkG216X8YWf6EIAABAwSQIQMwiQUBCAAARcE8m2NU6sjaouYWg1JiwhlFnNQqzm33/fowHa9XBLj35fqp67MbXWZufDflUKEXDyhjAkBCEAgmQQQoWTmhVlBAAIQCEwgnwipADvsO0lvf9t3z10GrZYEWdnJNYEIzlX0AAAEGElEQVQg/YpZETIhQv42uKDzDwyahhCAAAQgUFIEEKGSSic3AwEIlCOBfCKUWTnOtghln/eTr1iCX947c2XHL+8dZEXI3xqXueo1VIlwVoTK8RPCPUMAAhDITQAR4smAAAQgkHIC+UTILwxge2tc9nx8Wdlmy803qBpXrAipFa/M85IytwP6aVXz2X3X7fX5SYhQyh92pg8BCEDAIAFEyCBMQkEAAhBwQSBpxRIUA1WgwC+LrQRIrRDlKp8dVoRUxbrMK1OC/D/3ZSizXWbVuOzVKhc5Y0wIQAACEHBPABFynwNmAAEIQAACEIAABCAAAQhYJoAIWQbOcBCAAAQgAAEIQAACEICAewKIkPscMAMIQAACEIAABCAAAQhAwDIBRMgycIaDAAQgAAEIQAACEIAABNwTQITc54AZQAACEIAABCAAAQhAAAKWCSBCloEzHAQgAAEIQAACEIAABCDgngAi5D4HzAACEIAABCAAAQhAAAIQsEwAEbIMnOEgAAEIQAACEIAABCAAAfcEECH3OWAGEIAABCAAAQhAAAIQgIBlAoiQZeAMBwEIQAACEIAABCAAAQi4J4AIuc8BM4AABCAAAQhAAAIQgAAELBNAhCwDZzgIQAACEIAABCAAAQhAwD0BRMh9DpgBBCAAAQhAAAIQgAAEIGCZACJkGTjDQQACEIAABCAAAQhAAALuCSBC7nPADCAAAQhAAAIQgAAEIAABywQQIcvAGQ4CEIAABCAAAQhAAAIQcE8AEXKfA2YAAQhAAAIQgAAEIAABCFgmgAhZBs5wEIAABCAAAQhAAAIQgIB7AoiQ+xwwAwhAAAIQgAAEIAABCEDAMgFEyDJwhoMABCAAAQhAAAIQgAAE3BNAhNzngBlAAAIQgAAEIAABCEAAApYJIEKWgTMcBCAAAQhAAAIQgAAEIOCeACLkPgfMAAIQgAAEIAABCEAAAhCwTAARsgyc4SAAAQhAAAIQgAAEIAAB9wQQIfc5YAYQgAAEIAABCEAAAhCAgGUCiJBl4AwHAQhAAAIQgAAEIAABCLgngAi5zwEzgAAEIAABCEAAAhCAAAQsE0CELANnOAhAAAIQgAAEIAABCEDAPQFEyH0OmAEEIAABCEAAAhCAAAQgYJkAImQZOMNBAAIQgAAEIAABCEAAAu4JIELuc8AMIAABCEAAAhCAAAQgAAHLBBAhy8AZDgIQgAAEIAABCEAAAhBwTwARcp8DZgABCEAAAhCAAAQgAAEIWCaACFkGznAQgAAEIAABCEAAAhCAgHsCiJD7HDADCEAAAhCAAAQgAAEIQMAyAUTIMnCGgwAEIAABCEAAAhCAAATcE0CE3OeAGUAAAhCAAAQgAAEIQAAClgkgQpaBMxwEIAABCEAAAhCAAAQg4J7A/wMsbDTzzs2LmgAAAABJRU5ErkJggg==",
"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(pars=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(pars=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(pars=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(pars=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(pars=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(pars=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(pars=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(pars=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(pars=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(pars=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(pars=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(pars=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(pars=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(pars=bio.system_time, data_snapshot=arr, caption=f\"State at time {bio.system_time}\")"
]
},
{
"cell_type": "code",
"execution_count": 60,
"id": "f8d28717-210e-4adc-8d53-ceeb84b3bb26",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(5, 300)"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Now, let's examine the data collected at the 5 time points\n",
"all_history = history.get_array()\n",
"all_history.shape "
]
},
{
"cell_type": "code",
"execution_count": 61,
"id": "028d732d-f97b-4103-ae94-85a60cfce4dc",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compute time derivatives (for each bin) : MORE SOPHISTICATED METHOD, using 5-point stencils\n",
"df_dt_all_bins = np.apply_along_axis(num.gradient_order4_1d, 0, all_history, delta_t)\n",
"\n",
"# Let's consider the state at the midpoint in time (t2)\n",
"f_at_t2 = all_history[2] # The middle of the 5 time snapshots\n",
"f_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 62,
"id": "ca637810-beda-4abf-b232-473b17e1a76d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(300,)"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Computer the second spatial derivative (BACK TO SIMPLE METHOD, using central differences)\n",
"gradient_x_at_t2 = np.gradient(f_at_t2, delta_x)\n",
"second_gradient_x_at_t2 = np.gradient(gradient_x_at_t2, delta_x)\n",
"second_gradient_x_at_t2.shape"
]
},
{
"cell_type": "code",
"execution_count": 63,
"id": "845f75bd-25fe-4120-8975-e8dcd8e36afc",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.023351269997869635"
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compare the left and right hand sides of the diffusion equation\n",
"lhs = df_dt_all_bins[2] # t2 is the middle point of the 5\n",
"rhs = diffusion_rate*second_gradient_x_at_t2\n",
"\n",
"num.compare_vectors(lhs, rhs, trim_edges=2) # Euclidean distance, ignoring 2 edge points at each end"
]
},
{
"cell_type": "markdown",
"id": "fd1d0295-dfe2-4e7d-9457-cc2fea1a98f2",
"metadata": {},
"source": [
"#### Here, the discrepancy value has remained largely the same from the baseline 0.023163289760024783"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "30dffc16-7072-4a62-b606-642138990e5f",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
},
"toc-autonumbering": false,
"toc-showcode": true,
"toc-showmarkdowntxt": false
},
"nbformat": 4,
"nbformat_minor": 5
}