diff --git a/examples/plots/02.00.GS.gif b/examples/plots/02.00.GS.gif deleted file mode 100644 index ead7415..0000000 Binary files a/examples/plots/02.00.GS.gif and /dev/null differ diff --git a/examples/plots/02.01-trained_GS.gif b/examples/plots/02.01-trained_GS.gif deleted file mode 100644 index 12e0c7f..0000000 Binary files a/examples/plots/02.01-trained_GS.gif and /dev/null differ diff --git a/examples/plots/02.02-trained_GS.gif b/examples/plots/02.02-trained_GS.gif deleted file mode 100644 index 9000540..0000000 Binary files a/examples/plots/02.02-trained_GS.gif and /dev/null differ diff --git a/examples/plots/02.03_gridsize.gif b/examples/plots/02.03_gridsize.gif deleted file mode 100644 index 6fa892c..0000000 Binary files a/examples/plots/02.03_gridsize.gif and /dev/null differ diff --git a/examples/plots/02.04-DNS.gif b/examples/plots/02.04-DNS.gif deleted file mode 100644 index 25c5515..0000000 Binary files a/examples/plots/02.04-DNS.gif and /dev/null differ diff --git a/examples/plots/02.04-LES.gif b/examples/plots/02.04-LES.gif deleted file mode 100644 index bcad4bd..0000000 Binary files a/examples/plots/02.04-LES.gif and /dev/null differ diff --git a/examples/plots/02.04-NNclosure.gif b/examples/plots/02.04-NNclosure.gif deleted file mode 100644 index 5e47eee..0000000 Binary files a/examples/plots/02.04-NNclosure.gif and /dev/null differ diff --git a/examples/plots/03.01_Burgers.gif b/examples/plots/03.01_Burgers.gif deleted file mode 100644 index 61822e8..0000000 Binary files a/examples/plots/03.01_Burgers.gif and /dev/null differ diff --git a/examples/plots/03.02_burgers.gif b/examples/plots/03.02_burgers.gif deleted file mode 100644 index 499044d..0000000 Binary files a/examples/plots/03.02_burgers.gif and /dev/null differ diff --git a/examples/plots/04.01_2DBurgers.gif b/examples/plots/04.01_2DBurgers.gif deleted file mode 100644 index c998b2d..0000000 Binary files a/examples/plots/04.01_2DBurgers.gif and /dev/null differ diff --git a/examples/plots/derivatives_CPUbenchmark_1024_1024_10.png b/examples/plots/derivatives_CPUbenchmark_1024_1024_10.png deleted file mode 100644 index ecd75c4..0000000 Binary files a/examples/plots/derivatives_CPUbenchmark_1024_1024_10.png and /dev/null differ diff --git a/examples/plots/derivatives_CPUbenchmark_1024_1024_200.png b/examples/plots/derivatives_CPUbenchmark_1024_1024_200.png deleted file mode 100644 index 6c6b0f5..0000000 Binary files a/examples/plots/derivatives_CPUbenchmark_1024_1024_200.png and /dev/null differ diff --git a/examples/plots/derivatives_CPUbenchmark_64_64_10.png b/examples/plots/derivatives_CPUbenchmark_64_64_10.png deleted file mode 100644 index 4f84cec..0000000 Binary files a/examples/plots/derivatives_CPUbenchmark_64_64_10.png and /dev/null differ diff --git a/examples/plots/derivatives_CPUbenchmark_64_64_1000.png b/examples/plots/derivatives_CPUbenchmark_64_64_1000.png deleted file mode 100644 index 586496c..0000000 Binary files a/examples/plots/derivatives_CPUbenchmark_64_64_1000.png and /dev/null differ diff --git a/simulations/NavierStokes_2D/NS_solvers.ipynb b/simulations/NavierStokes_2D/NS_solvers.ipynb deleted file mode 100644 index 234a860..0000000 --- a/simulations/NavierStokes_2D/NS_solvers.ipynb +++ /dev/null @@ -1,1029 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Comparison of solvers for 2D Navier-Stokes equations\n", - "In this notebook we compare the *implementation and performance* of different solvers for the 2D Navier-Stokes equations.\n", - "We consider the following methods:\n", - "* `IncompressibleNavierStokes.jl`\n", - "* SciML, scpecifically `DifferentialEquations.jl`\n", - "* `CoupledNODE.jl`: an in-house wrapper for SciML neural closures." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Setup and initial condition" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "using GLMakie\n", - "import IncompressibleNavierStokes as INS\n", - "T = Float32\n", - "ArrayType = Array\n", - "Re = T(1_000)\n", - "n = 256\n", - "lims = T(0), T(1)\n", - "x, y = LinRange(lims..., n + 1), LinRange(lims..., n + 1)\n", - "setup = INS.Setup(x, y; Re, ArrayType);\n", - "ustart = INS.random_field(setup, T(0));\n", - "psolver = INS.psolver_spectral(setup);\n", - "dt = T(1e-3)\n", - "trange = (T(0), T(10))\n", - "savevery = 20\n", - "saveat = savevery * dt;" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## IncompressibleNavierStokes.jl solver\n", - "We know that mathematically it is the best one.\n", - "First, call to get the plots.\n", - "Look at the generated animation of the evolution of the vorticity field in time in plots/vorticity_INS.mkv." - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Iteration 100\tt = 0.1\tΔt = 0.001\tumax = 2.94716\n", - "Iteration 200\tt = 0.2\tΔt = 0.001\tumax = 2.60922\n", - "Iteration 300\tt = 0.3\tΔt = 0.000999987\tumax = 1.93674\n", - "Iteration 400\tt = 0.399998\tΔt = 0.000999987\tumax = 1.81722\n", - "Iteration 500\tt = 0.499997\tΔt = 0.000999987\tumax = 1.66464\n", - "Iteration 600\tt = 0.599996\tΔt = 0.000999987\tumax = 1.59957\n", - "Iteration 700\tt = 0.699995\tΔt = 0.000999987\tumax = 1.3922\n", - "Iteration 800\tt = 0.799993\tΔt = 0.000999987\tumax = 1.32724\n", - "Iteration 900\tt = 0.899992\tΔt = 0.000999987\tumax = 1.15229\n", - "Iteration 1000\tt = 0.999991\tΔt = 0.000999987\tumax = 1.25634\n", - "Iteration 1100\tt = 1.1\tΔt = 0.00100005\tumax = 1.24546\n", - "Iteration 1200\tt = 1.2\tΔt = 0.00100005\tumax = 1.10215\n", - "Iteration 1300\tt = 1.3\tΔt = 0.00100005\tumax = 1.12894\n", - "Iteration 1400\tt = 1.40001\tΔt = 0.00100005\tumax = 1.16041\n", - "Iteration 1500\tt = 1.50001\tΔt = 0.00100005\tumax = 1.15735\n", - "Iteration 1600\tt = 1.60002\tΔt = 0.00100005\tumax = 1.11329\n", - "Iteration 1700\tt = 1.70002\tΔt = 0.00100005\tumax = 1.04566\n", - "Iteration 1800\tt = 1.80003\tΔt = 0.00100005\tumax = 0.972415\n", - "Iteration 1900\tt = 1.90003\tΔt = 0.00100005\tumax = 0.98953\n", - "Iteration 2000\tt = 2.00004\tΔt = 0.00100005\tumax = 0.971748\n", - "Iteration 2100\tt = 2.10003\tΔt = 0.000999928\tumax = 0.928539\n", - "Iteration 2200\tt = 2.20002\tΔt = 0.000999928\tumax = 0.901102\n", - "Iteration 2300\tt = 2.30002\tΔt = 0.000999928\tumax = 0.877316\n", - "Iteration 2400\tt = 2.40001\tΔt = 0.000999928\tumax = 0.849998\n", - "Iteration 2500\tt = 2.5\tΔt = 0.000999928\tumax = 0.831636\n", - "Iteration 2600\tt = 2.59999\tΔt = 0.000999928\tumax = 0.820711\n", - "Iteration 2700\tt = 2.69999\tΔt = 0.000999928\tumax = 0.800575\n", - "Iteration 2800\tt = 2.79998\tΔt = 0.000999928\tumax = 0.769892\n", - "Iteration 2900\tt = 2.89997\tΔt = 0.000999928\tumax = 0.757282\n", - "Iteration 3000\tt = 2.99996\tΔt = 0.000999928\tumax = 0.762677\n", - "Iteration 3100\tt = 3.09996\tΔt = 0.000999928\tumax = 0.771079\n", - "Iteration 3200\tt = 3.19995\tΔt = 0.000999928\tumax = 0.770211\n", - "Iteration 3300\tt = 3.29994\tΔt = 0.000999928\tumax = 0.751899\n", - "Iteration 3400\tt = 3.39994\tΔt = 0.000999928\tumax = 0.714522\n", - "Iteration 3500\tt = 3.49993\tΔt = 0.000999928\tumax = 0.66334\n", - "Iteration 3600\tt = 3.59992\tΔt = 0.000999928\tumax = 0.664664\n", - "Iteration 3700\tt = 3.69991\tΔt = 0.000999928\tumax = 0.666573\n", - "Iteration 3800\tt = 3.79991\tΔt = 0.000999928\tumax = 0.66811\n", - "Iteration 3900\tt = 3.8999\tΔt = 0.000999928\tumax = 0.668276\n", - "Iteration 4000\tt = 3.99989\tΔt = 0.000999928\tumax = 0.666109\n", - "Iteration 4100\tt = 4.09989\tΔt = 0.000999928\tumax = 0.66184\n", - "Iteration 4200\tt = 4.19988\tΔt = 0.000999928\tumax = 0.656757\n", - "Iteration 4300\tt = 4.29987\tΔt = 0.000999928\tumax = 0.650323\n", - "Iteration 4400\tt = 4.39986\tΔt = 0.000999928\tumax = 0.640125\n", - "Iteration 4500\tt = 4.49986\tΔt = 0.000999928\tumax = 0.625802\n", - "Iteration 4600\tt = 4.59985\tΔt = 0.000999928\tumax = 0.609068\n", - "Iteration 4700\tt = 4.69984\tΔt = 0.000999928\tumax = 0.592071\n", - "Iteration 4800\tt = 4.79983\tΔt = 0.000999928\tumax = 0.585452\n", - "Iteration 4900\tt = 4.89983\tΔt = 0.000999928\tumax = 0.580601\n", - "Iteration 5000\tt = 4.99982\tΔt = 0.000999928\tumax = 0.572566\n", - "Iteration 5100\tt = 5.09981\tΔt = 0.000999928\tumax = 0.563368\n", - "Iteration 5200\tt = 5.19981\tΔt = 0.000999928\tumax = 0.554177\n", - "Iteration 5300\tt = 5.2998\tΔt = 0.000999928\tumax = 0.545931\n", - "Iteration 5400\tt = 5.39979\tΔt = 0.000999928\tumax = 0.538927\n", - "Iteration 5500\tt = 5.49978\tΔt = 0.000999928\tumax = 0.534137\n", - "Iteration 5600\tt = 5.59978\tΔt = 0.000999928\tumax = 0.533798\n", - "Iteration 5700\tt = 5.69977\tΔt = 0.000999928\tumax = 0.532659\n", - "Iteration 5800\tt = 5.79976\tΔt = 0.000999928\tumax = 0.529582\n", - "Iteration 5900\tt = 5.89976\tΔt = 0.000999928\tumax = 0.524232\n", - "Iteration 6000\tt = 5.99975\tΔt = 0.000999928\tumax = 0.516467\n", - "Iteration 6100\tt = 6.09974\tΔt = 0.000999928\tumax = 0.506342\n", - "Iteration 6200\tt = 6.19973\tΔt = 0.000999928\tumax = 0.494262\n", - "Iteration 6300\tt = 6.29973\tΔt = 0.000999928\tumax = 0.480726\n", - "Iteration 6400\tt = 6.39972\tΔt = 0.000999928\tumax = 0.466241\n", - "Iteration 6500\tt = 6.49971\tΔt = 0.000999928\tumax = 0.4514\n", - "Iteration 6600\tt = 6.5997\tΔt = 0.000999928\tumax = 0.436702\n", - "Iteration 6700\tt = 6.6997\tΔt = 0.000999928\tumax = 0.422441\n", - "Iteration 6800\tt = 6.79969\tΔt = 0.000999928\tumax = 0.408869\n", - "Iteration 6900\tt = 6.89968\tΔt = 0.000999928\tumax = 0.402614\n", - "Iteration 7000\tt = 6.99968\tΔt = 0.000999928\tumax = 0.403815\n", - "Iteration 7100\tt = 7.09967\tΔt = 0.000999928\tumax = 0.404441\n", - "Iteration 7200\tt = 7.19966\tΔt = 0.000999928\tumax = 0.404621\n", - "Iteration 7300\tt = 7.29965\tΔt = 0.000999928\tumax = 0.404509\n", - "Iteration 7400\tt = 7.39965\tΔt = 0.000999928\tumax = 0.404261\n", - "Iteration 7500\tt = 7.49964\tΔt = 0.000999928\tumax = 0.40398\n", - "Iteration 7600\tt = 7.59963\tΔt = 0.000999928\tumax = 0.403784\n", - "Iteration 7700\tt = 7.69962\tΔt = 0.000999928\tumax = 0.4037\n", - "Iteration 7800\tt = 7.79962\tΔt = 0.000999928\tumax = 0.403781\n", - "Iteration 7900\tt = 7.89961\tΔt = 0.000999928\tumax = 0.404021\n", - "Iteration 8000\tt = 7.9996\tΔt = 0.000999928\tumax = 0.404338\n", - "Iteration 8100\tt = 8.09964\tΔt = 0.0010004\tumax = 0.404655\n", - "Iteration 8200\tt = 8.19968\tΔt = 0.0010004\tumax = 0.404873\n", - "Iteration 8300\tt = 8.29972\tΔt = 0.0010004\tumax = 0.404899\n", - "Iteration 8400\tt = 8.39976\tΔt = 0.0010004\tumax = 0.404711\n", - "Iteration 8500\tt = 8.4998\tΔt = 0.0010004\tumax = 0.404326\n", - "Iteration 8600\tt = 8.59984\tΔt = 0.0010004\tumax = 0.403743\n", - "Iteration 8700\tt = 8.69989\tΔt = 0.0010004\tumax = 0.403037\n", - "Iteration 8800\tt = 8.79993\tΔt = 0.0010004\tumax = 0.402284\n", - "Iteration 8900\tt = 8.89997\tΔt = 0.0010004\tumax = 0.401533\n", - "Iteration 9000\tt = 9.00001\tΔt = 0.0010004\tumax = 0.400853\n", - "Iteration 9100\tt = 9.10005\tΔt = 0.0010004\tumax = 0.400321\n", - "Iteration 9200\tt = 9.20009\tΔt = 0.0010004\tumax = 0.399901\n", - "Iteration 9300\tt = 9.30013\tΔt = 0.0010004\tumax = 0.399562\n", - "Iteration 9400\tt = 9.40017\tΔt = 0.0010004\tumax = 0.399248\n", - "Iteration 9500\tt = 9.50021\tΔt = 0.0010004\tumax = 0.398874\n", - "Iteration 9600\tt = 9.60025\tΔt = 0.0010004\tumax = 0.398393\n", - "Iteration 9700\tt = 9.70029\tΔt = 0.0010004\tumax = 0.397708\n", - "Iteration 9800\tt = 9.80033\tΔt = 0.0010004\tumax = 0.396734\n", - "Iteration 9900\tt = 9.90037\tΔt = 0.0010004\tumax = 0.395421\n", - "Iteration 10000\tt = 10.0004\tΔt = 0.0010004\tumax = 0.393757\n" - ] - } - ], - "source": [ - "(state, outputs), time_ins2, allocation2, gc2, memory_counters2 = @timed INS.solve_unsteady(; setup, ustart, tlims = trange, Δt = dt,\n", - " processors = (\n", - " ehist = INS.realtimeplotter(;\n", - " setup, nupdate = 10, displayfig = false,\n", - " plot = INS.energy_history_plot),\n", - " anim = INS.animator(;\n", - " setup, path = \"./simulations/NavierStokes_2D/plots/vorticity_INS.mkv\",\n", - " nupdate = savevery),\n", - " #espec = realtimeplotter(; setup, plot = energy_spectrum_plot, nupdate = 10),\n", - " field = INS.fieldsaver(; setup, nupdate = savevery),\n", - " log = INS.timelogger(; nupdate = 100)\n", - " )\n", - ");" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Second, call to time it (without spending time in animations ~ 5% extra)" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "_, time_ins, allocation, gc, memory_counters = @timed INS.solve_unsteady(;\n", - " setup, ustart, tlims = trange, Δt = dt);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## SciML\n", - "### Projected force for SciML" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "import DifferentialEquations: ODEProblem, solve, RK4\n", - "# create some cache variables so that we do not have to allocate them at every time step\n", - "F = similar(stack(ustart));\n", - "cache_F = (F[:, :, 1], F[:, :, 2]);\n", - "cache_div = INS.divergence(ustart, setup);\n", - "cache_p = INS.pressure(ustart, nothing, 0.0f0, setup; psolver);\n", - "cache_out = similar(F);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "* Right-hand-side out-of-place" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "create_rhs_op (generic function with 1 method)" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "function create_rhs_op(setup, psolver, cache_F, cache_div, cache_p, cache_out)\n", - " function right_hand_side(u, p = nothing, t = nothing)\n", - " u = eachslice(u; dims = 3)\n", - " INS.apply_bc_u!(u, t, setup)\n", - " INS.momentum!(cache_F, u, nothing, t, setup)\n", - " INS.apply_bc_u!(cache_F, t, setup; dudt = true)\n", - " INS.project!(cache_F, setup; psolver, div = cache_div, p = cache_p)\n", - " INS.apply_bc_u!(cache_F, t, setup; dudt = true)\n", - " return stack(cache_F)\n", - " end\n", - "end" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Test the forces" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "_Note:_ Requires `stack(u)` to create one array" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "F_op = create_rhs_op(setup, psolver, cache_F, cache_div, cache_p, cache_out);\n", - "F_op(stack(ustart), nothing, 0.0f0);\n", - "prob_op = ODEProblem(F_op, stack(ustart), trange);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "* Right-hand-side in-place" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [], - "source": [ - "function create_rhs_ip(setup, psolver, cache_F, cache_div, cache_p, cache_out)\n", - " function right_hand_side(du, u, p = nothing, t = nothing)\n", - " u = eachslice(u; dims = 3)\n", - " INS.apply_bc_u!(u, t, setup)\n", - " INS.momentum!(cache_F, u, nothing, t, setup)\n", - " INS.apply_bc_u!(cache_F, t, setup; dudt = true)\n", - " INS.project!(cache_F, setup; psolver, div = cache_div, p = cache_p)\n", - " INS.apply_bc_u!(cache_F, t, setup; dudt = true)\n", - " du[:, :, 1] = cache_F[1]\n", - " du[:, :, 2] = cache_F[2]\n", - " nothing\n", - " end\n", - "end\n", - "\n", - "temp = similar(stack(ustart))\n", - "F_ip = create_rhs_ip(setup, psolver, cache_F, cache_div, cache_p, cache_out);\n", - "F_ip(temp, stack(ustart), nothing, 0.0f0)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Solve the ODE using `ODEProblem`. We use `RK4` (Runge-Kutta 4th order) because this same method is used in `IncompressibleNavierStokes.jl`." - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [], - "source": [ - "prob = ODEProblem{true}(F_ip, stack(ustart), trange);\n", - "sol_ode, time_ode, allocation_ode, gc_ode, memory_counters_ode = @timed solve(\n", - " prob, RK4(); dt = dt, saveat = saveat);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Test the difference between in-place and out-of-place definitions" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@time sol = solve(prob, RK4(); dt = dt, saveat = saveat);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In place: `46.664913 seconds (8.02 M allocations: 500.621 MiB, 0.07% gc time)`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@time sol = solve(prob_op, RK4(); dt = dt, saveat = saveat);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Out of place: `53.804647 seconds (8.36 M allocations: 84.888 GiB, 2.96% gc time)`" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## CNODE" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "WARNING: Method definition (::Type{CoupledNODE.Grid})(Int64, DataType, Union{Float32, Float64}, Union{Float32, Float64}, Union{Float32, Float64}, Int64, Int64, Int64, Int64, Union{Nothing, Array{Float32, 1}, Array{Float64, 1}}, Union{Nothing, Array{Float32, 1}, Array{Float64, 1}}, Union{Nothing, Array{Float32, 1}, Array{Float64, 1}}, Int64, Any) in module CoupledNODE at /Users/luisaorozco/Documents/Projects/DEEPDIP/CNODE_dev/src/grid.jl:9 overwritten on the same line (check for duplicate calls to `include`).\n", - "ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.\n" - ] - } - ], - "source": [ - "import DiffEqFlux: NeuralODE\n", - "import CoupledNODE: create_f_CNODE\n", - "f_dns = create_f_CNODE((F_op,); is_closed = false);\n", - "import Random, Lux;\n", - "Random.seed!(123);\n", - "rng = Random.default_rng();\n", - "θ_dns, st_dns = Lux.setup(rng, f_dns);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Define the problem and solve it" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [], - "source": [ - "dns = NeuralODE(f_dns, trange, RK4(), adaptive = false, dt = dt, saveat = saveat);\n", - "sol_node, time_node, allocation_node, gc_node, memory_counters_node = @timed dns(\n", - " stack(ustart), θ_dns, st_dns)[1];" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Comparison\n", - "Bar plots comparing: time, memory allocation, and number of garbage collections (GC)." - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAMgCAIAAABwAouTAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdZ1gU1/s38LPLFnoV6VIUBBULRayIYEVFxajYe69Ro0ZjiTFqoomKJpbYNVGDYEeMKKKxIjYEUUA6LE3pddl5Xsz1m/8+i65GFpdlvp8XXjtnzpxzz47LvWfmzA6HoigCAADAVlxlBwAAAKBMSIQA/+fEiRPTpk2Li4tTdiAscvLkyWnTpsXExCg7EGAvJEJo+qZMmSL4BK9evbpz586hQ4cyMjKUHTKL3Lt379ChQ2lpacoOBNiLp+wAABqcg4ND9+7dmcV37949e/ZMR0fH1dVVupqmpmbr1q29vLwMDQ2/eIzs5eDg4OXl1axZM2UHAuzFwWQZYJvw8PC+ffu6u7s/fPhQ2bEAgPJhRAjwf9LS0vLy8hwcHHR0dOiS58+fi8ViFxeXioqK0NDQtLQ0c3PzIUOGaGpq0hUePHjw6NEjQoiPj4+jo2PdNsVi8e3bt2NjY2tqauzt7fv06aOurv6J8SQmJt65cyc3N7d58+aOjo6dO3fmcDjSFVJTU2/evJmTk9O8eXNPT087OzvptSUlJa9fvzYyMrKxsUlLS7t27VppaamLi0vPnj3pCpWVlaGhoSkpKaampoMHD9bV1WW2lUgkT5480dDQaNOmTUFBwZUrV3JycqytrX19fZl9ZxQXF9+9ezctLa2kpMTCwsLb27t58+bSFSoqKuLi4vT09Fq1apWXl3f16tXs7Oy+fft27NgxPT09NzfX3t5euvfMzMwHDx6kp6fzeDwTExMPDw8rKyuZTl++fHnnzp13796ZmZn5+PiYmZlJr83Pz09NTTU3NzczM0tISLhx40Z5eXnbtm379OnD5eKSEPz/KACWuXbtGiHE3d297qrZs2cTQq5evcqUmJubCwSCR48eWVpaMp8aS0vLV69eFRcXDxo0iClUU1PbvXu3TIORkZEtW7aU/sRZWlpGRkZ+NMj8/Pxhw4bJfFo7duzIVKipqVmwYIGamhqzlsvlzpw5s6qqiqlz8+ZNQsjkyZO3bt0qXXPEiBHV1dUPHjwwNzdnCs3MzF6+fMlsW1paSghp27ZtUFCQlpaWdPwPHjyQDnX58uUCgUA6TqFQuGXLFuk6L168IIQMGjRo//79QqGQrrZ582aKohYsWEAIuXTpElN5/fr1PJ7sd/QDBw4wFUpKSkaOHCm9ViAQrFu3TiKRMHUOHjxICPn+++9XrFgh/e2hW7duRUVFH33/gVWQCIF1/msiVFNTs7S0nDx58vXr1+/evTt27FhCSM+ePUePHt2uXbugoKDHjx8HBgaqq6sLBII3b94w20ZFRamrq2tpaf3www8PHz589uzZtm3bNDU1tbW1ExIS5ERYXl7eoUMHQoi3t3doaGhSUtK9e/cCAwMHDhzI1Fm4cCEhxNHR8cKFC0lJSaGhoe3atSOETJkyhalDJ8IWLVro6Ojs2LEjKirqwoULrVq1IoSsXbvW2Nh46tSp9E6NGTOGEOLp6clsSydCQ0NDTU3N1atXv3r16tWrV6tWreJwOEZGRiKRiKk5derUBQsWXLp06cWLF8+fP//jjz/oLw1BQUFMHToRWlhYaGhofPvtt1evXo2IiLh//z5VJxGGh4cTQhwcHM6dO/fmzZuEhITr169//fXXJ0+epCtIJJLBgwfTKe3GjRuJiYmnTp2ie9y4cSPTI50IbW1tjY2N9+zZExUVdenSpfbt2xNCFi9eLOfNBxZCIgTW+a+JkBAyffp0pkQsFltbWxNCzM3NS0pKmPKlS5cSQnbs2MGUuLu7czicy5cvS3dx/PhxQsikSZPkRLh582ZCyIABA8Ri8XsrJCUlcblcDQ2N9PR0pjA3N5c+o/vs2TO6hE6EhJDw8HCm2q1bt+jCWbNmMYU1NTUtWrQghDAZjk6EhJD58+dLdz137lxCyNKlS+XEHxMTw+Vye/TowZTQiZAQUnfQLJMIV61aRQgJDg7+UON0prSwsCgtLWUKnz9/Tr8h+fn5dAmdCPl8flxcnPT7xuPxLC0t5QQPLIRz5QAft2TJEua1mppajx49CCHTp0/X1tZmynv16kUISU5Ophfj4+OjoqJcXFx8fX2lmxo3bpyent6VK1fkdPfXX38RQjZs2CB9PlPauXPnJBLJhAkTpE/YGhsbT5s2jRBy9uxZ6codO3b08fFhFrt160afePz666+ZQh6PR0+sZeJn0AleZjEkJERO/O3atbOzs4uKiqL+/7l4hoaGM2bMkLMhIcTAwIAQQp99fW8FuuuFCxdKn7B1dnYeNGgQfR1XuvLgwYOdnJyYRTs7u1atWmVlZVVVVckPA1gFk2UAPoLL5cpc5zM2NiaEODg41C3MycmhFx8/fkwIqaqqWrlypUyDAoEgNze3vLy87qwTQkhtbW1sbCyXy6XP470Xfct/p06dZMrpG0JiY2OlC2XiVFNTMzQ0LCgokJlZQ8efm5srXainp2djYyNdYmdnp6urm5yczMRfVVW1e/fu4ODg1NRUkUgkkUiYyqWlpcy0I0JIq1atZK4m1jVy5Mjvv//+559/Pn/+vK+vb+/evX18fKTfKDn7fvHiRZkfQ2jdurVMNRMTk/j4+Nzc3Lqzb4C1kAgBPoLH48n8+aanHWpoaNQtZMYxhYWFhJDExMT9+/fXbdPAwKCiouK9ibC0tFQikTRr1oyZVPLeOoQQExMTmXK6pKSkRLqwbi9cLpfP5/P5/LrxS6cx8r/sKKN58+bFxcUlJSWampoSicTX1/fGjRu2trbDhg1r3rw5PSd2165dmZmZYrFYesNPuVnQ2tr64cOHa9euvXLlyvbt27dv366hoTFjxoxNmzbRQ0B632VmpX5o32WOEalzmAAIEiFAA6FHQqNGjTp69Oh/2lBbW1tNTe3t27eVlZUfutGCbpwZfTJEIhEhRPo+hHqSGSDScnJyOBwOHUNoaOiNGzd8fHzCwsKkp3ru3r277oYy9358iJOTU1BQUGVl5f3798PDww8ePBgYGFhaWkpf9vti+w7sgWuEAA2CPnd37969/zr4UFNTa9++vUQiefr06Yfq0BNEo6OjZcqjoqKYtQpRXFycmJgoXZKQkFBSUmJnZ0cPNOkgv/rqK+ksKBKJMjMz69m1urq6l5fXxo0bHz9+zOfz6ekz5AvuO7AHEiFAg2jXrp27u3tCQgI9jpHBzMl8rwkTJhBCvvvuu5qamvdWGD58uJqa2okTJ9LT05lCkUh0+PBhDofz1Vdf1S/2/88vv/wivbht2zZCCNMFfe5U5pdC161b99nnHuu+M0ZGRgKBoLq6ml6ku969e7d0zadPn165ckVLS0tmahLAp8CpUYCGcuDAge7du8+aNevRo0e+vr62trZ5eXkJCQmnT5+2sLCg76N4r7lz5546der69eteXl7Lli1zcHDIz8+Pjo6+fPny9evXCSEtWrRYsmTJ1q1bvb29N2/e7OTkRN/kV1paOnv2bOl5kvVkZGR0/PhxDQ0N+n6Po0eP7t+/v3nz5suWLaMreHp6crnc7du3m5qa9unTp6Sk5MCBAydPnjQxMal79vJTzJs3Lz09fcyYMY6Ojqampunp6du3by8rK5swYQJ9ZtXT09Pf3z8kJMTb23v9+vXW1taPHj1auXIlRVHr16/X19dX1L4DeyARAjSU9u3b37lzZ+bMmfv27du3bx9TbmxsXPdXY6QJhcJ//vlnzpw5p06d8vf3Z8qlfzp8y5Ytampqv/zyC/MbKzweb/HixVu3blXgLpiamm7evHn8+PHbt2+nS2xtbUNCQphpL05OTrt37168ePGiRYvokmbNmgUHB3/33XeflwhbtGhx6tSpiIgIpoTD4YwcOfK3335jSk6cODF37tyjR48yv+yjqan5008/MekZ4D/Bj24D61RUVGRkZKirq9edQJ+Xl1dYWGhhYcHMtExNTZVIJLa2ttLVCgoKioqKTExMpG9lq6qqyszM1NbWrjuhMS4u7vHjx6WlpcbGxi1atHBxcfnQDYIy0tPT79y5U1RUZGho2KZNm7Zt28pUyM3NjYyMLCgoMDAw8PT0lPm9zcrKyqysLB0dHZnJn5+yU2VlZdra2m3btn3x4kVRUVF4eHh+fn6LFi28vb3rTmfNzs6+f/8+fU9C7969NTQ0MjMzq6qqbGxs6Fma1dXVGRkZmpqapqamMtvWfc/Ly8sfPXqUlpZWXl5ubm7evn17+mZ/GampqXfu3CkuLjYxMenVq5fMM0NKSkry8vIMDAzoGxOlQ62oqGjRokXdX3ED1kIiBID3kE6Eyo4FoGFhsgwAALAaEiEAALAazpIDwHsIBIItW7bgwfHABrhGCAAArIZTowAAwGpIhAAAwGpIhAAAwGpIhAAAwGpIhAAAwGpIhAAAwGpIhAAAwGpIhAAAwGpIhAAAwGpIhAAAwGpIhAAAwGpIhAAAwGpIhAAAwGrsTYTV1dXKDgEAAJSPpc8jjI2NnTdv3s2bN5UdSONSUVGhoaGh7CigoeD4yhcWFrbiUChXoMJvkUQi4XJVdXhDUVRrXcnpPdu+fNeNNBGmpqYGBwc/f/5cR0fnq6++6tWrF7Pq5MmToaGhJiYmCxYssLa2pgtFItGOHTsyMjJ8fHwmT57M4XDkt19eXl5eXt6AO6CaxGKxskOABoTjK9/Tp0/jyjTE/TYqOxC2KslL2tiRKCMRNtLvDtu2bYuPj+/Vq1eLFi2GDBny999/0+V79uxZvXq1r68vRVE9evQoLS0lhNTU1PTq1evdu3d+fn7btm3bvHmzUmMHAABV0khHhDt37mQG+GVlZSdOnBg1ahRFUb/88suOHTv8/PzGjBlz//79kydPzpgx4/z581wud+/evRwOx8LCwt/ff9myZQKBQLm7AAAAKqGRjgilT3Pn5uYaGxvTL5KSkry8vOjyXr163bt3jxBy7969Xr160adDu3btWlRUlJSUpISgAQBABTXSESHj/v37J06ciIqKIoSIRCKBQKCrq0uvMjY2fvLkCV1uY2NDF3K5XENDw+zsbCcnJznNvnv3LjEx0d/fnykZP358v379GmgvVEVZWdlHL6+C6sLxlQ8zyZWPougLXgqkrq7O430k0zXqRBgbGzt8+PCjR4+2bt2aEKKuri4Wi2tra9XU1AghlZWVmpqadHlNTQ2zVVVVFV0uh46OjpGRUUBAAL3I4XA8PDw+ulWTV1tbizehCcPxlY/P5ys7BOAo/L/op0yjbbyJ8NWrV/379//111+HDx9Ol5ibmxNCMjMzW7RoQQhJT0+3sLAghFhaWsbHx9N1iouLCwsL6XI5eDyegYHBqFGjGnAHVBCXy1XdudfwUTi+8mG4rHycT8pbCtdIPxWJiYl9+vRZt27dmDFjmEIdHZ2+ffseP36cEFJUVHThwoURI0YQQkaMGHHt2rXc3FxCyJ9//unm5mZlZaWsyAEAQLU00hHhqlWr8vLytmzZsmXLFkKIk5PTpUuXCCGbNm3y9fX9999/X79+3bt3b/r+wvbt248bN87Nzc3Z2TkqKio4OFjJ0QMAgOpopIlw165ddAqkMfdCuLi4JCQkREdHGxkZOTs7S9efO3dudna2i4uLvr7+lw4XAABUViNNhCYmJh9apaOjw9xBIc3JyUn+TFEAAIC6Guk1QgAAgC8DiRAAAFitkZ4aVUWXLl1avGo9RSk7jnqQUBIuR4W/G7WwMI0Iu6TsKABAxSARKkxiYmK6VqvqPl8rOxC2qi7P/G2osoMAANWDRKhIHD0TYu2i7CjYqrJE2REAgEpS4fNgAAAA9YdECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArIZECAAArMZTdgDvR1FUYmJidHR0cXHxzJkzmfLQ0NCMjAz6ta6ubkBAAP26pKTkwIEDIpGod+/eAwYMUELEAACgmhrpiPDff//t1q3bzp07Z8+eLV0eGBj4999/R0dHR0dHx8TE0IUSicTb2/vWrVvm5uazZs3as2ePMkIGAACV1EhHhN26dcvLy4uLi2vXrp3MqmnTpo0ZM0a65MqVKwUFBffu3ePxeG3atJk2bdrMmTPV1NS+YLwAAKCqGmkilJPGIiIisrOzHR0dBwwYwOVyCSGRkZHe3t48Ho8Q0rt3b5FIlJSU5ODg8OXCBQAAldVIE+GHtGzZksPhZGRk7Nmzx9bW9sqVK2pqatnZ2VZWVnQFHo9nZGSUlZUlPxEWFxenpKRMnz6dKQkICOjRo0d9YhOLxYSi6tMC1BdFKisrlR1E41VZWcnn85UdReMlFouVHQLrUZTCP8J8Pv+jJwhVLBH+9ttv9Ivvv//e0dHx3LlzI0aM4PF4tbW1TJ2ampqPftqFQqGWlpabmxtTYmFhUc+/EVwul3A49WkB6otD8IdeDj6fj/dHDvoMEygVR+H/RT/lsKpYImTo6Oh07NjxzZs3hBBzc/O0tDS6vLy8vLCw0NzcXP7mQqHQ2NhYZiZOPeFT1Bjg2rAcampqeH/kwEdY+TjK+Qir0oGvra1lzl3k5eVFRUW1adOGEDJ48OB//vmnpKSEEHLu3DknJydbW1tlBgoAAKqjkY4IS0pK/P39y8rKKIrq27evoaHh6dOnCwoKOnbs2KNHD4FA8M8///Tv39/X15cQ0rVr1169evXo0cPd3f38+fOHDx9WdvgAAKAyGmki1NDQ2LJlC7NInzU2NjYODw+PiYmRSCQrVqxwdnZmKpw+ffrmzZsZGRmrV6/GcBAAAD5dgyTC0tLS/Px8HR0dIyOjz2uBx+O5urrKFHI4nDZt2tCnQ+uu6t279+f1BQAAbKawa4RisfjcuXNjxowxNzfX0dGxtbVt1qyZtrZ2v379du7c+e7dO0V1BAAAoEAKGBFKJJIjR46sX78+IyOjbdu2vr6+1tbWBgYGpaWlIpHo8ePH33777bfffjt79uw1a9YYGBjUv0cAAABFUUAi/Oabb/7888/58+ePHz/exsamboWSkpKQkJDffvutc+fOCQkJ9e8RAABAURSQCEeNGvXDDz9oamp+qIKOjs6kSZMmTpx48+bN+ncHAACgQApIhB4eHp9SDfNZAACgEWrYG+rv3r175syZnJycBu0FAADgsykyEVZWVpqZmZ0/f55eXLduXffu3UeOHNm2bVvm2YEAAACNiiIT4aNHj/Ly8vr160cIKSws3Lx58/Tp02NiYtq2bbt+/XoFdgQAAKAoikyEIpHIxMREQ0ODEHLt2jWxWLxhw4Z27dotXLjw/v37CuwIAABAURSZCNXV1SsrKymKIoRcvnzZycnJzMyMEKKjo4Mb6gEAoHFSZCJs165dUVHRoUOHYmNjz549S/8iNiEkJSXFxMREgR0BAAAoiiJ/a9TGxmbhwoX0Y9+NjY0XLVpEl1+4cEH6EbgAAACNh4J/dPvXX38dN25cRkZGz549DQ0NCSEURQ0fPrxTp06K7QgAAEAhFP/0CVdXV+kHR3A4nGnTpim8FwAAAIVQwDXCxMREiUTyKTVfv35d/+4AAAAUSAGJ8Ndff+3YsePx48fLysreW0Eikdy8eXPEiBHe3t717w4AAECBFHBqdOvWrVu3bp0zZ87cuXO9vLw8PDxatGihq6tbUVGRk5Pz+PHjyMjI9PR0f3//yMjI+ncHAACgQApIhFpaWuvXr1+8ePHRo0dPnz69YcOGmpoaZq2Dg4O/v/+sWbMcHR3r3xcAAIBiKWyyjL6+/qJFixYtWlRWVpaampqXl6evr29ubm5sbKyoLgAAABRO8bNGtbS02rRpo/BmAQAAGkLDPoYJAACgkUMiBAAAVkMiBAAAVkMiBAAAVkMiBAAAVlP8rFFCSH5+fkJCgra2trOzc0O0DwAAoCgKHhHm5OQMGjTI2Ni4W7dumzZtIoRQFNWhQ4fAwEDFdgQAAKAQikyEEolk6NChL168OHz48MyZM+lCDoczaNCgkJAQBXYEAACgKIpMhFFRUQ8ePLh8+fLkyZOtrKyY8vbt2+O5EwAA0DgpMhGmpKQYGxu3a9dOplxbW/vdu3cK7AgAAEBRFDlZxtDQ8O3btyUlJTo6OtLlT58+tbCw+K+tFRcXP3nyJDc3d+TIkdLlsbGx165da968ub+/v7q6Ol1YU1Nz7ty5jIwMLy+vTp061WcvAACAVRQ5IuzWrZu+vv6SJUvEYjGHw6ELnz179uuvvw4dOvQ/NXXv3j1jY+Np06aNHj1aujw0NNTT0zMjI+Pw4cO9e/cWi8WEEIqiBg0atHPnTpFINGDAgJMnTypqjwAAoMlT5IhQS0trz549Y8eOvX79upaWVmVlZd++fSMjI+3s7L777rv/1FSnTp0KCwuTk5NlTrSuX7/+p59+mj59ulgsdnZ2vnjx4vDhwyMjI2NjYxMTEzU0NLp06bJixYqAgAAmEwMAAMih4NsnRo4ceefOHTc3t3fv3uXn56enpy9evPj+/fsGBgb/qR11dXUNDQ2Zwrdv30ZFRdGDSx6PN3DgwLCwMELI1atX+/btS9f39fVNTk5+8+aNgnYIAACaOMXfUN+5c+e///5b4c0SQrKysng8XrNmzehFc3Nz+pH3mZmZlpaWdKFQKDQ0NMzMzGzZsqWcpsrLy7Ozs3/88UempG/fvh06dKhPePR5WlAmilRVVSk7iMarqqpKIBAoO4rGCx9h5WuAjzCfz+dyPzLka5BflmkgMmc7KYqiSzgcDkVRdcvlq62tLSwslN5KcZECAIDKUHAirKysPHz48N27d0UikXS5s7Pzr7/+Ws/GTU1NxWJxQUEBPSgUiUSmpqaEEDMzM6a76urqd+/emZmZyW9KU1PT0tJy69at9QxJGo+nSt8qmiYOEQqFyg6i8aqursb7Iwc+wsqnpI+wIq8R1tTU9O3bd+7cuY8fP1ZgswwjIyMXF5dLly4RQmpra8PCwvr160cI6du3b3h4OD2gvnr1qrW1tZ2dXUMEAAAATY8ivwE9ePDg33//PXfu3H+9WaKu0tLSqVOnFhcXUxQ1atQofX39/fv3E0LWrVs3derU169fP336VFNTk+7I29vb3t5+wIAB3bt3P3DgwLZt2z56RhgAAICmyESYl5enpaVV/yxICBEIBPR99NOmTSOEMDfO+/n53bhx49q1a6NHjx45ciSfzyeEcDicsLCwoKCg7Ozsixcvuru71z8AAABgCUUmQg8Pj+rq6uTkZFtb23o2xSTCutq3b9++ffu69ceNG1fPTgEAgIUUeQrR3Nx8x44dAQEBDx48qK6uVmDLAAAADUTB19J8fHwqKiq6dOkiFAo5Ujw9PRXbEQAAgEIo8tRoSUmJl5dXTU3NjBkzLC0tpWfBMje8AwAANCqKTIQREREikSg+Pr5169YKbBYAAKDhKPLUKJfL1dLScnBwUGCbAAAADUqRidDT01NdXf3ff/9VYJsAAAANSpGnRgUCwbp16wICApYuXerm5qalpcWs0tbWxvlSAABohBSZCNPS0hYuXEgIWbp0qcyqnj173rp1S4F9AQAAKIQiE6GlpeW1a9feu0pfX1+BHQEAACiKIhOhpqZmnz59FNggAABAQ8OPUwMAAKspYES4bdu2q1ev7tixQ19ff/Lkye+to5DnEQIAACicYk6NUhRFP+H9Q895x/PfAQCgcVJAIpw/f76zs7OlpaW+vn54eHj9GwQAAPhiFHCNMDc3d8CAAS9evKh/UwAAAF8YJssAAACrIRECAACrKew+wtLS0sLCwg92w+Npa2srqi8AAABFUVgiHDhwoJy1PXr0uH37tqL6AgAAUBSFJcLRo0ebm5t/aK2dnZ2iOgIAAFAghSXC+fPn9+jRQ1GtAQAAfBmYLAMAAKyGRAgAAKymgEQoFAq7du2qq6tb/6YAAAC+MAVcIzQxMbl792792wEAAPjyFPk8QoCmzczKRpSRquwo2IvD5R4/dmzcuHHKDgSaGiRCgE/17m0B2f2WCPHTEMohOLM0Pz9f2VFAE4TJMgAAwGpIhAAAwGpIhAAAwGoqdo0wISGhuLiYfi0UCtu1a0e/pijq9u3bIpGoa9euVlZWygsQAABUjIolwgULFmRmZpqamhJCzMzMjh07RpePGzfu6dOnbm5uc+fOPX78uPxfAAcAAGCoWCIkhKxatWrMmDHSJQ8fPrx27VpCQoK+vv6xY8dWrlyJRAgAAJ9I9a4RpqSkREZGZmVlMSUXLlzo16+fvr4+IWTEiBGxsbEpKSlKiw8AAFSKio0I+Xz+xYsXr169Gh0dPWfOnJ9//pkQkpmZyVwX1NLSMjAwyMzMtLGxkdNOdXV1Xl7e3r17mRIvLy97e/v6xCaRSOqzOShEbW1twzVONVzT8GkkEknDHWJ8hJWPUvxHmMvlcjgc+XVULBEGBwcLBAJCSGJioqur64ABA7y9vcViMZf7f0NbHo9XU1Mjv53KysqysrJHjx4xJS1btpSfOz9KIpEQCn8qlYoiHz30oMIoqra2tuEOMRJhI0Ap/Pjy+Xw1NTX5dVQsEdJZkBDSqlWr7t27P3782Nvb29TUNC8vjy4Xi8Vv3741MzOT346urq6Njc2BAwcUGBuPxyMf+94BDYtD1NXVG7J5UCoOh8/nN9wh5vFU7O9hE8ThNOhH+ENU7xohraqqKj4+3tzcnBDi6ekZERFBD6hv3bplbGzcsmVLZQcIAACqQZW+AeXn548fP97Ly0soFAYFBenp6fn7+xNCfH19165dO2bMGC8vr23bti1fvhzf7AAA4BOp0ohQT09v3LhxhYWFWVlZs2bNevDgAT2IVlNTu3nzpoeHR0JCQmBg4MKFC5UdKQAAqAxVGjnx+fwJEya8d5Went7SpUu/cKVndOkAACAASURBVDwAANAEqNKIEAAAQOGQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWQCAEAgNWaTiJMTEyMjIwsLi5WdiAAAKBKmkgiXLp0aa9evTZu3Ghvb3/37l1lhwMAACqDp+wAFODFixcHDx58+fKlmZnZrl27li1bhlwIAACfqCmMCM+cOdOvXz8zMzNCyMSJEx88eJCZmansoAAAQDU0hRFhenq6jY0N/VpPT09fXz8tLc3CwkLOJtXV1SKR6KeffmJKbGxsDA0N6xPGq1evavOyycvr9WkEPl91pUQiuXbtWsP1IJFISHwkEag3XBcgR21e6qtX4oY7xAkJCZL8UnyElaaiuFZcq/Dja2dn17JlS/l1OBRFKbbXL2/cuHG2trYbN26kF83NzU+cOOHt7S1nk1u3bvn6+goEAqbEysrKyMioPmEUFhYmZebWpwXlklBUTU2NUOo9UTkaAn4be9uGaz824U1ltbjh2m9oVdXVfD6fy+EoO5DP19Kiub6+fgM1Xl1d/U7YnBBVfX+qa6pLSkqMDOv1d0y5hNWF+nwFpyR/f//58+fLr9MURoSmpqYFBQX0a4lE8vbtW/o0qRyenp6lpaUNHxoAADR2TeEaoYeHx+3bt+nXDx8+1NHRsbOzU25IAACgKppCIhw2bFhlZeX8+fNDQkJmzpy5YMECoVCo7KAAAEA1NIVrhISQzMzM7du3Z2Vl9e7de9q0aVxuU0jwAADwBTSRRAgAAPB5MHICAABWQyIEAABWU1u/fr2yY4Av6tmzZ4QQbW3tgoKCu3fvampqamtr06v+/fdfbW1tDQ0NQkhRUVF4ePj9+/fz8/NNTEwEqnx/YROWk5MTHh7+6NGjoqIiU1NTHu/9N0QFBwe/fv3a0dGxuLj49u3bNTU1zZo1o1fV1NREREQUFxebmpoSQu7cuWNkZITDrVxVVVWRkZG3b9/Oz883NjamZ//du3cvOzub+amQ9PT05ORk+qgRQmpqaiIjIyMjI/Py8iwtLZn/CbGxsc+fP09NTa2urtbW1pb+H5KamhoTE5P2P5mZmVZWVl92RxsNClima9eu+/fvpygqLCyMEDJs2DBmlY2NTXh4OEVRUVFRxsbGQ4cOnT17to+Pj5eXl9LChQ+7ePGioaHhqFGjZs6c2b179wkTJnyo5t69ew8cOEBRVFRUFJfL7dChA7MqJCSEy+UOHjyYXuTz+Y8fP27oyEGOqKgoa2trFxeXKVOm9OvXz8LC4tatWxRFtW3bVk1NLS4ujq4WGBg4aNAg+nVCQoKjo2OnTp2mT5/u4eFhY2Pz/PlzetWECRNsbGz69Onj5uamr6/fq1evFy9e0Ku+++675s2be/3PwIEDv/i+NhZN4YZ6+GxmZmZ37969f/9+ly5dpMt//PHH6dOnb9q0iV4sKytTRnTwEWvXrt20adOsWbPoRenDVF5enpiYqKGhYW9vTwiZMGECs0pfX5/D4Tx9+rRjx46EkMOHD3t5eX3RuOHDiouLhwwZMmfOnLVr19IlSUlJRUVF9Ou2bduuX7/+9OnT0ptIJBJ/f/+ePXvu3buXnjC/YsWKoUOHxsfH0yP7gICAzZs3E0LKy8u/+eYbHx+fly9fGhgYEEJ69+596tSpL7mDjROuEbKatrb2ypUrly9fLlNeXl5eVVXFLGppaX3ZuOCTlJeXV1dXM4vMYTp48KCdnd3y5csnTpy4ZMkSQsi6deuY3yAkhEyaNOno0aOEkNzc3MePHw8cOPDLBg4fdObMGYFAsGrVKqakZcuWLi4u9Otvvvnm+vXrDx48kN7kwYMHr1+/3rJlC3Pb2Lp16woLC0NDQ2Ua19TU3LFjh1AoRPKTgUTIdvPmzcvIyLh69ap04bJly/744w9nZ+f58+dfuHBBIpEoKzyQY/Xq1d98842Hh8eSJUtu3LhBF8bHx3/99dcRERFhYWH37t1jhvXSJkyYcPr06erq6uPHjwcEBHzoyiJ8ec+fP2/fvv2Hjoient4333zDDBZpr169srW1lX5mgKamppOTU0xMTN0W+Hx+p06dYmNj6cXbt2/7/c+KFSsUtx8qBomQ7QQCwZo1a1asWCGd7fr27fvmzZtly5ZVVlZOmTJlyJAhFO43bXwmTJiQmJg4ZcqUvLw8Pz+/OXPmEEKuX7/eu3dvJycnuo66+nuelWFkZOTh4REaGnrs2LGJEyd+0aBBrtraWj6fL6fCokWL4uLimO89hJCqqqq6m/D5fLH4/T8Qz+fza2tr6detWrWa8z/Dhw+vX+wqDIkQyMSJE2tqas6cOSNd2KxZs0mTJh04cCA6OjosLOzJkyfKCg/ksLS0nD179vHjx69cubJ37978/PzS0tJPOZU9ZcqUlStX8vn89u3bf4E44RPZ29vHx8fLqaCurr569eqVK1dKb5KcnFxZWcmU1NbWvnr1ysHB4b0txMXF0VeOCSFmZmYD/0dmogCrIBECUVNT++GHH9asWcN8hSwvL2fWGhgYcLlcjio/u6epkj5M9EPEOBxOx44d7969K32J970GDhxob2+/bNmyhg0R/qORI0empqZKT4cpLS198+aNdJ1p06YVFhaeO3eOXuzWrZuuru6OHTuYCgcPHqyqqhoyZIhM4xRF7dmzJyUlZfTo0Q22ByoJ1waAEEL8/f23bt36+vVretHPz69Zs2bu7u5cLvfPP//s3r07xg2NkJubm4eHR8eOHauqqg4ePDhu3DgjI6P+/ft36NBhwIAB48ePLy4u5vF4CxYsqLstn8+/ePHie5vdsmULc5fh2LFju3fv3oD7AP8/MzOzY8eOTZo06fLlyy4uLhkZGSEhIYGBgdJP1OHz+Rs2bBgzZsygQYMIIerq6idPnvTz84uJifHw8Hj27Nnff/99+vRpXV1duv7NmzdXrlyZlZX19OnTd+/ehYSEMDcjRkdHz5s3j2n5l19+ee+59CYPvzXKOufPn2/durWjo2NGRsadO3eY74axsbF37twZNGiQhYVFdnb21atXk5KSuFxuu3bthg0bJv+6BShFcnLytWvXUlJShEJh586dBwwYQA/cJRJJcHDws2fPtLW1hw8f3rp16wcPHnC5XHd397y8vNDQ0EmTJkm3ExMTk5mZOWDAAELIgQMHmAtIhBBPT0/mciN8MRkZGcHBwVlZWaampr6+vq1btyaEBAUFde7c2dramhAikUiOHDlibm5OHzVCSE5OzqlTp7Kyspo3bz5q1Cjm1viIiIiEhARCiJaWVsuWLd3c3JiZOI8ePYqOjpbud+rUqez8pCMRAgAAq+EaIQAAsBoSIQAAsBoSIQAAsBoSIQAAsBoSIQAAsBoSIQAAsBoSIQAAsBoSIQAAsBp+Yg2g0UlOTv77778JIZMmTTI1NZVe9fvvv5eUlHTq1Klfv36f17i9vX1AQMAPP/yggEAJGTNmTHFx8eXLlxXSGoBSYEQI0Oi8fv165cqVK1euPHTokHR5XFzcvHnzVq5cyfzgsnxhYWEcDkfmySGFhYUVFRWKCrW0tLSkpERRrQEoBRIhQCPVqVOnI0eOSP8I4uHDh9u2bSsQCJQYFUDTg1OjAI3UlClTFi1adO/evW7duhFCxGLxn3/+uWTJktWrV8vUvHTp0tmzZ9++fWtpaTl9+vQOHToQQuLi4oKCggghR44c+eeffwghEydONDMzozdJSEj47bffsrKyWrduPXfuXKacEJKfn79v377o6GiBQODp6Tl16lTpJxK8efNm586d6enprVq1WrhwYQO/BwBfAkaEAI2Uvb19ly5djh49Si9euXIlLy9v3Lhx0nUoipo6daqfn19OTo6VldX9+/fd3d0vXLhACHn79i392IGXL19GR0dHR0czzy+MiYnp0aNHZmamUCjcvn177969a2pq6FUpKSkdO3b89ddf9fX1CSFff/21p6dnWVkZvfbly5dubm7BwcGmpqZJSUmdO3fOysr6Im8GQEOiAKCRCQsLI4RcuXJl3759urq6ZWVlFEX5+/sPGTKEoiiBQDBnzhy65p9//kkI+fvvv+nF2trakSNHmpqaVlVVURR15coVQsjjx4+lG2/WrJlAIIiOjqYXQ0NDCSHBwcH04rBhw7S0tF69ekUv0rNg1q1bRy/279/fxMQkOzubXtyzZw8hpGfPng31RgB8ERgRAjReAQEBYrH43LlzBQUFly9flnmOICHk6NGjHTp0GDlyJL3I5XIXLlwoEolevXolp1lfX18XFxf6df/+/fl8fnx8PCFEIpGEhYVNmDDBwcGBqdm1a1f6FGtNTU14ePi0adOYiazTp083MTFR3O4CKAeuEQI0Xrq6usOGDTt69Gh+fr6WltbgwYNlKsTHx5eUlLi5udGLNTU17969I4QkJSU5Ozt/qFkbGxvmNZfL1dPTKywsJIRkZmZWVla2bdtWunLbtm2PHz9OCElJSamtrW3Tpg2zisfjMSkTQHUhEQI0apMnTx4wYEBSUtLYsWOFQqHMWoqi7O3tly1bJlPu6uoqp001NbX3lldWVtJtynTB5XKZraSfX08IkUgk9FoA1YVECNCo+fj4WFpaJiUlTZkype5aJyen169f+/v7vze38Xg8QohEIvnEvlq0aCEUCp89eyZd+OzZM3t7e3otn89/8eIFs6q6uvrVq1dOTk6fvjsAjRC+ygE0alwu9+HDh0lJScxVPWmzZs1KSUnZsGGDdLa7e/cu/cLc3JwQIv96oTShUOjv7//XX3/FxMTQJUFBQY8ePRo7diwhhMfj+fr6Hjp0KCUlhV4bGBiYn5//mTsG0GhgRAjQ2MmZkOLv7//tt9/+8MMPQUFBLi4uZWVlsbGxCQkJ9OlNJycnT0/PyZMnb9q0SV1d/dixY9JX+N7rl19+iYqK6tatW58+fcrLy69fv967d+/FixfTa7dv396tWzc3N7d+/fqJRKI3b954eHgocE8BlIIjcz0AAJQuIyMjNDR04MCBVlZWddceOHDA0dGxR48eTMmTJ09CQkIyMjL09PTs7e0HDx5sbW1Nr6qtrb1x40Z6erpYLB4+fLixsfGxY8ccHBy6dOnCbH7s2LHWrVszKa2srOzo0aPMDfWjRo2SPu+ak5Ozd+/e1NTUVq1azZgxIzo6urq62s/Pr0HeCIAvAokQAABYDdcIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1XjKDgDgP0tOTn706BEhxNTUtGfPnnUrPHjwIC0tjRDSoUMHBweHLx0fNKSuXbu+ePEiJydHU1NT2bFAE4FECKonIiJi2rRphBA9Pb3s7GwNDQ3ptRKJZOTIkenp6YSQbdu2LV26VDlRQsMoKysrLS2lKErZgUDTgVOjoKo0NTWLioouXLggU37t2rX09HQMF5qq77777rfffhMKhcoOBJoOJEJQVYMHD9bQ0Dh69KhM+eHDh7lcrp+fn1KigoY2atSouXPn8ng4mwUKg/9MoKr09fX9/PyCgoIyMjIsLS3pQnqM2L9/fysrq/du9eTJk8OHDz9//ry6utrOzm78+PEDBgyQrrBly5aioqLNmzdHRET88ccfKSkpZmZmM2fO7N+/PyEkIyMjMDAwKiqKoqi+ffsuWbJE5sQsRVFnzpwJDg5OTU1VV1d3c3ObPXt2y5YtpeusWrVKR0fn22+/vXTp0vHjx9PS0nx9fZ2cnB49ehQQENCxY0eZmPfu3ZuSkjJ//nxmN98rNjb2yJEjT548KSsrMzU19fDwGD16tK2tLVMhLS1tz549Dx8+LCsrs7KyGjZs2JgxY7jc//s2HB4eHh4ePnLkSAsLix07dty/f19NTc3Hx2fx4sWampq1tbVHjhw5f/58bm6uo6Pj8uXL27Rpw2wbHx9/5MiRrl27+vj4BAYG3rx5s6yszNnZecmSJTKXaRMSEi5fvhwVFZWeni6RSFq1ahUQECBzFO7du3f+/HlfX19nZ+fdu3ffvn27uLj4yJEjjo6OO3fuzM7O/uGHH/h8Pl05Kytr3759z58/z87O1tXVtbCw8PT0HD16tPRZgRcvXuzfv//Zs2disdjW1nbMmDGDBg2S7vH06dNPnjyZPXs2IWTHjh3R0dFqamo9evRYsmSJoaGhnLcdmgIKQNUcPHiQEDJz5szQ0FBCyJYtW5hVv//+OyHk1KlT33zzDSFk27Zt0huuX7+ew+FwOBxra2snJyd6VDF37lzpOtbW1lwud9euXRwOx9DQkM49HA7n2LFjT58+NTY2FggErVq1ov8KDxw4UCKRMNtWVlYOHjyYECIQCDp27GhjY0MIEQqFZ86cke6Cx+NZWVktX76cEKKhoWFhYTFmzBj6HO/o0aNldlYkEvH5fEtLS7FYLOc92bx5M53SrKysOnXq1KxZM0LI+PHjmQpXrlyhEwNdQV1dnRDi4+NDX2+jrVu3jhCydu1aU1NT6d308fGpqqqiM0eLFi309fUJIbq6uvHx8cy2Fy9eJIRMnDixXbt2GhoaHh4ednZ29A5ev35dOlR6cyMjow4dOtjY2NBhL1q0SLrOrl27CCFLly6l30NTU1Ntbe179+5RFOXs7EwIYcKOi4szMjIihJiZmXl4eLRp04bezZiYGKa1P/74Q01NjRDSqlWrDh060K/Hjx8v/ZZOmDCBEPLrr7/q6+tL77uTk1NJSYmcdx6aACRCUD1MIqytrbW0tHRwcGBWde7cWU9Pr7y8vG4iPH78OCGkXbt2z549o0vS09M9PDwIIcePH2eqWVtb03+mL1y4QJecO3eOw+GYmJg4OTktWbKksrKSoqiMjAx7e3tCyD///MNsu3LlSkJI586dMzMz6ZK//vqLz+dramomJiYy1Xg8Ho/HMzIyCgsLo/8WFxYW1tbW2tjYCASCnJwc6Z3dtGkTIeT777+X84acPHmSEGJiYnLjxg2m8NmzZ0FBQfTrzMxMPT09NTW1gwcP0iU5OTmenp4y3wPoRMjn82fPnl1eXk5RVHZ2dtu2bQkhXl5eHTp0eP36NUVRtbW18+fPJ4SMGTOG2ZZOhGpqat27d8/NzaUL6dPUJiYmhYWFTM0NGzbQ7dCeP3/eunVrQoh08HQi5PP5w4YNy8rKoiiqurqaDkkmEY4bN44Qsnv3buYbSVVVVUhIiEgkohejo6N5PJ6mpuaVK1fokqSkJEdHR5n/HnQi5PP5a9asqaqqoiiqoKCge/fuhJCtW7fKefOhCUAiBNXDJEKKolasWEEIuX//PkVRsbGxhJA5c+ZQFCWTCJk0Iz2IoSgqMTGRx+O5uroyJXQi3L59u3Q1Omd069ZNevxHjz7XrVtHL5aUlGhra3O53NjYWOltZ82aRQhZuHAhU0KPRH/77TeZ/frxxx8JIT///DNTQp85VFNTS01N/dC7QdchhISHh3+ozpo1awgh48aNky5MTk7m8XhCoTAvL48uoRNhp06damtrmWpHjhyhM9yLFy+YwsLCQoFA0KJFC6aEToQcDicuLk66lxEjRhBCAgMDPxQbRVHXrl1jDiiNToRGRkZFRUUylWUSYdeuXQkhcgZtY8eOJYSsWbNGuvDmzZv0QJMZFNKJ0M/PT7paZGQkIWT48OFygocmAJNlQLVNnTqVw+HQU2boP9mTJ0+uWy0mJiYlJcXV1ZUefDBatmzp4ODw5MmT8vJy6fLhw4dLL9IXw4YNG8bhcJhCJycnQkhqaiq9+PDhw9LS0m7duklfOSOEzJgxgxBy/fp1mZDGjx8vUzJjxgyhULhv3z7qf/cGXLt2LTExcfDgwS1atPjQO/D69evExEQ7OzsfH58P1aF7nz59unShjY1N3759q6qqbt26JV3u5+cnfeGQ3k17e3t6aEjT09OzsLDIzMwUi8XS23p4eND1GVOmTCF1dl8ikbx69SoiIiIoKCgoKIh+D+Pj42XC9vPz09XV/dBO0ehx+bJly7Kyst5b4caNG6TOvvfq1at169YikejFixfS5XTaZnTo0IEQkpKSIj8GUHVIhKDaHBwcPDw8/vrrr5KSkhMnTrRu3bpz5851qyUmJhJCYmJiWtaRnJwskUjy8/OZyhwOR2auDf3nWGauCl1YXFxML9J/zenBmTT6LzWTL2mGhoZ1/8QbGxsPGzYsKSkpIiKCLtm/fz8hZObMmXLegaSkJPp9kFPnQ7HRJfSPDzBkki4dZ91MrKurW1tbW1ZWJl1IXxesWyK9+0FBQTY2No6Ojt7e3qNGjRo1ahSdpUpKSmS2pS8Qyvfdd9+Zm5vv27fP0tLS2dl58eLF9CVJem11dbVIJBIKhXUnT9EzmOTvu56eHpE6xNBUYdYoqLzJkyfPnj177ty52dnZixcvfm+dyspKQoilpSV9krMu6RmGHA5HekjEoCdZfEhtbS0hhJ6EIk0oFHI4nJqaGulCLS2t9zYyZ86c06dP//HHH97e3jk5ORcuXLCysqInrH4IPSaTmbz6ibHRW8nE9hn7zqh7ex/dBTNwjIiIGD16tJGR0datW3v27GlgYMDj8YqKilxcXKg698h/ys2g9vb2sbGxx48fv3jx4t27d3fu3Llz587evXufP39eR0eH7pc+BO8NTGbfP3E3oYlBIgSVFxAQsGTJkhMnTvB4vLrnG2kmJib0v/v27WugMJo3b04IoX/RRlpaWhpFUXQAH9WrV6927dqFhITk5uYeOnSopqZm5syZ8v86m5mZEULevHkjp46xsbFIJEpLS6Nnk0rHRv735ihERkaGTIlMF/RsnQMHDgwdOpSp8+TJk/p0qq+vv2DBggULFlRXV9+8eXPlypURERGbN2/etGmTpqamlpZWcXFxUVERPbyTCYw+asByODUKKk9PT2/06NEGBgZDhw41Nzd/bx0PDw91dfWoqKiCgoIGCsPd3Z3L5d6+fbuwsFC6/Pz584SQLl26fGI7s2bNqq6uPnLkyKFDh3g8Hn2NTY6OHTvq6enFxsbKnOWTRvdOz2dhlJWVhYeH/6fYPurevXvv3r2TLqFvcaFn5xJCRCIRIUT6ciMh5MqVKwrpXSAQ9OvXb+/evYSQqKgouvC9+56amvr06VN1dfVOnToppGtQaUiE0BQcOnTo7du3Z86c+VAFHR2dOXPmlJeXT5w4sbS0VHpVXl7eP//8U/8YzMzMhgwZUlxcTE9QpAvT0tJ++eUXQgg9d/RTTJgwQVtb+4cffkhMTBwyZIiFhYX8+nw+f/78+WKxeOLEidKX2SQSCX35kBAyY8YMDocTGBjIlBBCNmzYkJ+f37t3bwX+Lnlpaen333/PLCYmJu7du5fH4zEzmOgrc3R2pL18+fLnn3/+7B7PnTtXUVEhXRIXF0ekhnr0ZKUNGzYw34Fqa2uXL19eU1Mzbty4D52jBlbBqVFgi02bNj1//jw0NLRVq1be3t4tW7YsKChITEyMjIzs379/v3796t/Frl27Hj16tHv37ujo6L59++bl5Z0+ffrt27fz58/39vb+xEb09PTGjh1LT5P5xPS5du3a+/fvX79+3d7eftiwYSYmJhkZGTdv3uzWrRt996S7u/u33367adMmV1fX0aNHm5qaRkRE3L59u3nz5vT4SVG8vb3/+OOPp0+f+vj45OTk/Pnnn0VFRRs3bqRnDBFCJk6cePDgwSVLlty5c6d9+/Zv3rw5derUoEGDgoKCPq/HBQsWlJaW+vj42NnZaWtrx8XFnT17VkNDY8mSJXSF0aNHnzt37tSpU+3atRs1apSOjs7ly5efPn1qb2//008/KWa3QdUp894NgM9y/vx5V1fXjRs3yqmzY8cOV1dX6TvlKYqqqakJDAxs3749PXVCTU3Nyspq2rRpkZGRTJ3Bgwe7u7vLtLZz505XV1fpe+cpinr58qWrq+uyZcukCzMyMiZNmsTMCHVyctq3b5/03YcURXXu3NnX11dO8PRdbra2ttL388lXVVX1008/MflGKBT27Nnz3Llz0nWOHj1K34RHCNHS0ho7dmxycrJ0hX379rm6ujK/JEBLTk52dXVdsGCBTI8BAQGurq7MDXz0ucf58+c/ePCAOddqZWX1+++/y2x49uzZtm3b0lNyzMzMNmzYUFRU5OrqKn17/unTp+sePtrIkSM7dOhQUVFBL27bts3Ly0tHR4fuUV9ff8iQIY8fP5beRCwWb9u2jfmhO0NDwzlz5jB3T9LWrFnj6uoqsyFFUa6urkOHDq0bBjQlHApPMwH2qaqqKi0t1dfXb6BZghRFvX37VkND4/MegrFx48Y1a9b8+OOPq1at+q/blpeXV1ZWyvl5zIqKivLycgMDg/fODv1sly5dGjJkyPz58+l74SsrKysqKgwMDD5Uv6qqqqqq6qO3CX66srIyiUTCZMT3Ki0tra6uxm+HggycGgU2EgqFDfocHw6HQ/8A5mcoKSn57bfftLS0Pv2yojRNTU352VdDQ0P+vRYKoa6uXvduDWkKPwSfcrVPW1tbgT1Ck4FECNBYnDlzJiEh4cKFCyKRaOXKlZ+dSgHgP0EiBGgsjh49eunSJS6XO27cOOm5lwDQoHCNEKCxEIlEFRUV5ubmqvj49fLycpFIpKurK3PPPkDjh0QIAACshhvqAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCA1ZAIAQCAJ7TAlAAAIABJREFU1ZAIAQCA1ZAIAQCA1ZAIAQCA1ViaCCsqKm7evKnsKAAAQPlY+himR48ezZ079+HDh8oOpHEpLS3V1tZWdhTQUHB85duxY8eKP85z7dyVHcjnoyiKw+EoO4rPRNVU6seHiVJef/muWfqEenam/4/C29K04fjKV1lZKbHpXD10o7IDYauSvLKoU0rpmaWnRgEAAGhIhAAAwGpIhAAAwGpIhAAAwGpIhAAAwGpIhAAAwGpIhAAAwGoqmQgzMzNXr149ZMiQ8PBwpvDixYsDBgwYPnx4YmKiEmMDAADVopI31IvFYmdn58TExOzsbLokPT19+fLld+/eTUpKGjly5JMnT5QbIQAAqAqVHBFaW1sHBASYmJgwJefPnx8xYoSBgYGbm5tAIHj9Wgk/0gMAAKpIJRNhXVlZWWZmZvRrCwuLjIwM5cYDAACqookkQg0Njerqavp1ZWWlpqamcuMBAABV8fnXCDMyMg4dOpSVldWyZcvp06cbGBjUrRMcHBweHm5sbDxv3jzpM5kfUlNT8/Tp0+jo6JKSkm+++UZ6VUhIyLVr15o1azZv3jxTU1OZDR0dHcPCwujXr1+/btWq1efuFgAAsMtnjgjfvn3r7u6en5/v4+MTFRXVq1ev2tpamTq///770qVL3d3ds7Oze/ToUVVV9dFm79y5M3bs2JCQkNWrV0uX79u37+uvv3Z3d8/Nze3evXtOTs5PP/0UHR19+fLlXbt2EUL8/PweP368c+fOefPm9enTp1mzZp+3XwAAwDafOSJ8+PAhj8cLDAwkhAwZMkRLSys5OVl6HCaRSLZu3fr777/7+vpOmTKlU6dOwcHBY8eOJYSsX7/e09PT29ubrnn58uWYmJiVK1cSQry8vBISEmJiYlxdXWWa2r1795AhQ6ZOnerq6nrx4kU7O7vFixcTQoRCIf3vjRs3Ll682L9//8GDB380folEUlxcLH33RYcOHYyNjT/v3QAAANX1mYnQycmppKTk9evXDg4Od+/eNTY2Njc3l66QmZmZkpJCZzsOh+Pt7f3vv//SibB///7+/v4nTpzw8fEJCwubPn36hQsX5PQlEomSkpKYxOnt7R0VFbVv3z6ZagYGBhMnTvzE+AsLC7Ozszdt2sSUzJ4929fX9xM3b6rKyspU96me8FE4vvIx8wxAaSiqtLRUsU2qq6vzeB/JdJ+ZCK2trffv3+/q6qqvr19RUXH+/HmZ+SkikUhbW1tdXZ1eNDY2fvDgAf26a9euwcHB/v7+c+fO3bNnz+XLl11cXOT0JRKJNDQ0tLS0mKbi4+M/L2yGoaFh69atb9y4Uc92mhiKovAE8yYMx1c+gUBASJmyo2A3Dkcp/0U/MxHGxcXNnz//xIkT7u7uly5dGj169LNnz4yMjJgKQqGwpqaGWayurtbQ0GAWu3XrtmjRolWrVgUGBsrPgkxTFEXRX2ZlmgIAAKiPz5wsc/bs2S5dugwdOtTc3HzmzJkGBgbXrl2TrmBubl5VVZWXl0cvZmRkSJ87DQ0N3bVr1759+zZt2vTRYZmZmVltbW1OTs57mwIAAKiPz0yEzZs3T0xMFIvFhJCioqLMzEz67oj4+PioqChCSLNmzXr06PHXX38RQoqLiy9fvjxs2DB627CwsGnTpp0/f37mzJkhISFjx46VnrRSl6GhoaenJ91USUnJpUuXmKYAAADqiTt16tTz58/TKe3TjRs3Tk9Pz9XVdfLkyS4uLn369PHy8iKEHD16dOPGjXSdzZs3b9y48auvvurcubOXl1fPnj3p8levXl2+fNnd3Z0Q0rVr17NnzzI/k/327Vs3N7eAgACxWOzm5sbMXtm0adPmzZtHjBjRuXPn7t279+rVSxH7DgAAQDg2NjYpKSkWFharVq2aOXPmR2fXMCiKevTokUgksrW1bdeuHV2Ym5tbVVVlZWVFL+bn59+9e9fc3NzV1fVTpquJxeJnz54xiwKBwNnZmX5dUFBw584dMzMzNze3+s98i4qKmjdv3sOHD+vZThNTUlKio6Oj7CigoeD4yrdly5Y1t9+Jh29UdiBsVZKnvbFjSUHul++Zl5ycnJSU9Oeff27atOns2bMyl/rk4HA49KhOWvPmzaUXmzVr5ufn9x+i4fGk7yCUZmRk9J+aAgAA+BRcQkjLli3Xrl2bmJg4e/ZsZccDAADwRf3fZBl1dfURI0YoMRQAAIAvj3v//v3KykplhwEAAKAc3K5du7q4uIhEImVHAgAAoATcq1evisXin376SdmRAAAAKAG3X79+CxYsuHXrlrIjAQAAUAIuIaRZs2b5+fnKjgQAAEAJuISQ1NRUPIoPAADYiUsIuXLlSqdOnZQdCQAAgBLwKioqunXrNnToUGVHAgAAoAQ8DQ2NzZs3KzsMAAAA5eBWVFTs378/PT1d2ZEAAAAoAVdDQ2Pv3r0HDhxQdiQAAABKwCWE+Pn53bx5U9mRAAAAKAGXEGJjY5ORkaHsSAAAAJSASwgpLCzkcrkfrQoAAND0cCmKCg4O7tixo7IjAQAAUALeqFGj7t69i2uEAADATtzLly/v3bu3Z8+eyo4EAABACXhv3rwxNTVVdhgAAADKwTU1Nc3IyEhKSlJ2JAAAAErAI4Ts3LnzzZs3wcHByg4GAADgS+MRQqqqqgQCgbIj+Q8KCwsvXbr08uXLr776inluxtOnT/fv36+np7d8+XIDAwPlRggAAKqCe+PGjX/++cfa2lrZkfwHWVlZCQkJUVFRL168oEvevn07atSoKVOmtGnT5quvvlJueAAAoEK4Pj4+RUVF06ZNU3Yk/0GbNm2+//57R0dHpiQkJGTIkCHu7u4TJkzIz89PTU1VYngAAKBCuPfv309KSrK3t1d2JPWSkpJia2tLv7a1tU1JSVFqOAAAoDJ4Hh4eyo5BAdTU1CQSCf26trZWTU1NufEAAICq4NVn4/T09NDQ0Kqqqs6dO3fp0qVuhZiYmOvXr5uamg4fPlwoFH5Km8XFxY8fPy4oKBgxYoR0+YsXL8LDw01MTPz9/es2ZWdn9/jxY/r1mzdvmNEhAACAfJ//W9sXL17s0KHDvXv33rx58+OPP9atcO7cud69e2dmZu7fv79Pnz7MiE2Of//919jY+P+xd+cBNeX//8Df93a7ddv3lbQSGaJCtlChspYIUWNfsw0yGGPGvo4lYxjGHpPsVEQYe4WKmDZK0l66bbdu9/7+ON+5v/spUrm3U+7z8dc57/M+7/frOOnVOed93mfatGnjx48XL7906ZKTk9P79+///PPPQYMGlZSUhISEJCcnP3ny5NKlS4SQUaNGhYeHX79+fffu3WZmZsbGxk0+LgAAkCms0tJSFRWVxu7G5XL9/f2Dg4OHDBnyuTpr1qzZuXPnpEmTqqqqOnbsGBYW5uHhQQg5fPhw7969RUNd4uLiXrx4MXHiREKInZ1dcXFxSkqKnZ2deFM///zz9u3b/f39q6urbWxsrl27lp6ePmDAAEIINS5GXV39ypUrf/31l4aGxunTpxtyCFVVVWlpaaJVExMTFuurro8BAKA1YhkZGY0fP37mzJndu3dv+G5RUVG6urq2trahoaH6+vp9+vRhMBjiFT58+BAfHz9ixAhCCJvNHjp0aEREBJUI1dTUXFxcIiIibGxs4uLi3Nzc9u3bR+3F4XDq9pWbm/vs2TOqKXl5+aFDh/7zzz9BQUG1qllZWW3YsKGB8efl5b1+/drZ2VlU8uOPP/r4+DT0+L9RpaWldIcAUoTzWz8ej0d3CDJPKORyuZJtUlFRUV5evv46LBsbmwMHDhw4cKBnz57Lly8fPXp0Q5pOTU2tqakZMmRInz597ty50759+9DQUPFcmJWVpaysrK6uTq0aGho+f/6cWh4zZgyfzx86dOiOHTsWLFjw+++/jxw5sp6+srKyFBQUtLS0RE09efKkIUHWQ1dXt0uXLl/fzrdHVVWV7hBAinB+66GgoEBIOd1RyDYGg5YfUdajR48cHBw8PT2vXbvm6em5aNGiHTt2fHE3Pp//5s2blJQUU1PT0tJSMzOzmzdvuri4iCowGAyhUCi+i3ia9PHxSUlJGTt27M6dO+vPgrV2/FwJAABA0zCjo6OfPHkSGBh49+7dnTt3BgUF8fn8L+5mZGSkp6dnampKCFFRUbGxsUlKShKvYGBgUF5eLrrIzcnJEf/GRVxc3O+//75o0aLt27e/fv26/r4MDAx4PF5xcTG1mp2djc9lAACApDDt7e1FK4MHD66qqmrILVpq6GZubi4hpKqqKikpiXpjoaCg4MOHD4QQIyOjzp07X7lyhRBSXV0dERExePBgat/4+Hh3d/d9+/bt2LFj586dLi4uL1++rKcvfX39rl27Uk3x+XzxpgAAAL7S/4yTNDc3j4mJUVNT++JuhoaGs2bNGjJkiJeX161btzp16kQNH922bVtiYuLFixcJIWvWrJk9e/aLFy+io6O1tbXd3d2pfU+cOPH7779Tg1/GjBlTU1Pz999/r127lhDy8ePH6dOnf/z4saamZuzYsdra2r///jvV1PTp01+9ehUbG6uurj5s2DBJ/zsAAICMqv0kr1GuX78eFxdnaWk5YsQIajKXly9flpSUODo6UhWePXt269YtPT09b29vRUXFLzbI4/Go9wIpHA5HlPOeP39+8+ZNXV1db2/vTw4ubZTo6Oi5c+disEwtXC4Xgym+YTi/9du0adPqf4r4o9fRHYis4uaprLPlFuQ2f8+sO3fuODk5NW3nwYMH17pLaWNjI77arVs30WeSGkJBQcHb2/uTm2xtbW1tbZsQJAAAQD2YAwYMGD16dEVFBd2RAAAA0IC5a9euq1evUo/oAAAAZA0zICBg0aJF58+fpzsSAAAAGjAJIR07dvyaITMAAACtF1MoFF68eHHo0KF0RwIAAEAD5vDhw9PS0gICAuiOBAAAgAbMq1evxsfHL1++nO5IAAAAaMD6+++/CSH4ki0AAMgm1udeYAcAAJAFTLoDAAAAoNP/T4QPHz5cuHAhjaEAAAA0P+arV6/27t3bq1cvJycn0QflAQAAZASrU6dOSkpKPj4+x44da9++Pd3xAAAANCvW2bNnnZ2dNTQ06I4EAACABiwvLy+6YwAAAKANRo0CAIBMQyIEAACZhkQIAAAyDYkQAABkGhIhAADINCRCAACQaUiEAAAg01h0B9AU1dXVcXFxCQkJffv2tbKyogqzs7NDQkLU1dXHjx8vLy9Pb4QAANBaMDdt2pSTk0N3GI3z/PnzrVu3/v77748ePaJKKioqXF1d1dXV//333ylTptAbHgAAtCLMNWvWmJiYjB079tatW3QH01AODg5nzpzp1auXqOT8+fN9+vSZPHny+vXrnzx50upSOwAA0IXZuXPnqqqqS5cuOTs79+jRIzExke6QmiIpKaljx47Ucvv27VNSUuiNBwAAWgumtrb2uXPnuFzutWvXKioqPD096Q6pKaqrq1ms/3veyWKxqqur6Y0HAABaC9b169epJTc3Nw0Njd69e5eWlqqoqDRw/8rKyqysLAMDAyUlpbpbCwsLnzx5oq+v361bt4bHlJeXV1paamZmJl5YVFT05MkTXV3d7t27192lbdu2b9++pZbfvXvXtm3bhncHAACy7H9enzAxMfH29mYwGA3ff+HChVZWVpGRkXU3PXr0yNraet++fd7e3n5+fg1p7dmzZ23btjU2Nu7QoYN4+ZMnTzp06LB3795x48b5+vryeLzY2Njc3Ny3b9/GxcURQkaNGnX+/Pl3797duHFDTk7OwsKi4YcAAACy7H8SobGx8d9//62srNzAnW/fvp2cnGxubv7JrYGBgYGBgZcuXYqJiYmIiHjw4AFVHhUVlZeXJ6qWlZV17949atnExCQsLOzJkye1mlqxYsXSpUsvX74cExNz69atGzduHDhwQF1dPTMz8/jx44QQAwODP/74Y+HChadOnQoNDW1g/AAAAE1/j7C8vDwgIOD8+fNubm51txYUFNy5c+fMmTOEEA0NDQ8Pj3PnzvXu3ZsQ8vjx4wULFty8eVNXVzcrK2vQoEFz5szp27cvIURbW1tbWzshIUG8qaKioqioKCrhqaurDxs27ObNm3/88UetHgcMGDBgwIAGBp+Xl/fs2TNNTU1RycaNG319fRt++N+ksrKyRt0PgNYF57d+VVVVdIcg84TC0tJSyTapqKgoGkHyOU1PhD/++OPkyZM/dxMyMzNTQUFBX1+fWjUxMXn16hW1HBgYyGQynZycgoODJ0yY4O/vHxAQUE9H79+/Z7FYhoaGoqao26FfQ1dXt0uXLqI7ukwmU11d/Svb/AYIhcKGPx6GVgfnt35sNpuQMrqjkG0MBi0/ok1MhI8ePfrnn39E77PXxePx2Gy2aFVBQaGyslK0umzZsry8vG7duq1ZsyYwMLD+vng8nry8vOgv2VpNNZmcnJz4FSEAAMgmFnVVpKen16VLl4bvtmbNmjZt2qxevZoQkp+ff+zYMXl5efF7pAYGBlwut7KyUlFRkRCSl5cnuqQjhGRlZV2+fHnIkCGhoaFz5szR1dWtpy8DA4Py8vLy8nJqYGqtpgAAAL4G09XV1dXVde3atY3azc/Pr3fv3pqampqamnJycioqKhwOR7yCsbGxiYnJ7du3qdXbt29TDwgJITk5Oa6urv7+/mFhYb6+vk5OTtnZ2fX0ZWhoaGpqGhUVVbcpAACAr8SKiYkhhGhoaDRqtwkTJoiWDx065OnpSQ1UWbVqVWJi4rlz5+Tk5JYsWTJv3ryff/754cOHRUVF3t7eVP2VK1dOmzZt0aJFhJBly5bxeLyNGzfu2rWLEFJeXr5nz57s7GyBQLB582ZVVdU5c+Ywmcwffvhh/vz5RUVFjx8/zs3NHTt2rIQOHwAAZB3Lzs7uK5v44YcfOnXqRC07Ozt37dqVWg4ICNDT07t+/bqBgcH9+/epe6SEkP3794uP4Vm9ejWfzxdv0MDAYP369eIlc+fO1dXVDQ8P19fXf/DgwSdf3gcAAGgCRnFxsQwOmIyOjp47d27dFxZlHJfLVVVVpTsKkBac3/pt2rRp9T9F/NHr6A5EVnHzVNbZcgtym79nprGx8e7du5u/YwAAgJaA6eTktGDBgrNnz9IdCQAAAA2YV65c8fDwCAoKojsSAAAAGjAZDMaoUaNSU1PpjgQAAIAGTEJIcnKylpYW3ZEAAADQgHns2LHffvtt5MiRdEcCAABAAyY1R8zy5cvpjgQAAIAGzNOnT0dGRuIVdQAAkE2scePG0R0DAAAAbZhfrgIAAPDtYmppaWlpaU2ePJnuSAAAAGjAmjFjBiGkc+fOdEcCAABAA9amTZvojgEAAIA2eEYIAAAyDYkQAABkGis2NpYQoqGhYWFhQXcwAAAAzY1lb29PCPH09AwNDaU7GAAAgObGunHjBiFET0+P7kgAAABowHJxcaE7BgAAANpgsAwAAMg0JEIAAJBprAMHDhBCzMzMXF1d6Q4GAACgubFmzpxJCPH09GxdiZDL5aakpJiYmGhra1MllZWV//zzj5qaWs+ePemNDQAAWhFWYWEhIUReXp7uSBrhwYMHc+fO5XK5a9asmTRpEiGEz+e7uLg4OTklJSVZWVlt2LCB7hgBAKB1YGlqatIdQ6P17t372bNnAQEBopJr1661a9du/fr1NTU1VlZWy5Yt09DQoDFCAABoLZg8Hi8yMpK6Lmy94uPjqZkB5OTkrK2tk5KS6I4IAABaByabzZ4zZ87BgwfpjuSrlJWVKSgoUMscDofL5dIbDwAAtBZMBoPh5eVFzS/TWLm5ubm5ufVUEAqFGRkZJSUlTQ2voU0ZGBiIIsnJyTE0NPz6HgEAQBYwCSGWlpYZGRmN2i0yMtLMzMzGxqZr1642Njbx8fF166SkpHTq1MnZ2bldu3YbN25sSLP//vuvq6urlpYWm80WL09LS7OxsRk0aFC7du3WrVtXd0d3d/dz585VVlYmJycXFBRYW1s36nAAAEBmMQkhlZWVfD6/Ubupq6ufO3cuLy8vKytr8ODB1Gfua/nhhx88PDySk5OfPn26efPmxMREqjw5OZnH44mqlZWVpaWlUcscDmfGjBnHjh2r1dTSpUuHDBmSkpLy/Pnz7du337p1y97e/sGDB7t27aLe+rCyslqwYIGzs/OMGTNOnDjBZGKiAAAAaBAWISQ8PLxjx46N2s3BwYFaYDAYHh4ep06dqlWhpKTk6tWrycnJhBAzMzN3d/dTp05RF3NBQUFJSUnnzp1TVFQsKysbNmyYo6Mj9cKDiYmJiYlJQkKCeFOlpaWXLl16/fo1IaRdu3bDhg27ceNGTExMrR6nTp06derUBsafn5+fmJjYvXt3apXJZC5evHjEiBGN+kf49pSVlTEYDLqjAGnB+a1fVVUV3SHIPKGwtLRUsk0qKiqyWKz667CWLl165cqVCxcuNLmbEydODB06tFbhu3fvGAxGu3btqFUrKyvRZd/27dunTJkycuTI4OBgb29vCwuLT97tFMnMzBQKhWZmZtSqpaXlv//+2+RoKdra2qampuJDhCwtLVVUVL6y2dZOKBTiH+EbhvNbPzabTUgZ3VHINgaDlh9R1vbt23/88ceRI0c2bf8//vjj9u3bT548qVXO5XI5HI7oz08lJaWPHz9Sy3JycocPH548ebKpqem4ceMOHDhQ/1+pXC5XUVFRdLdTWVlZ1FSTMRgMJSUlOzu7r2wHAABaO1ZsbGy3bt2atvPx48fXrVsXFRVV93OGenp6XC6Xz+dT16RFRUX6+vqirTweLysrS19f//379zweT1FRsZ5e9PT0ysrKqqqqqBE0tZoCAAD4GswmZ8GzZ88uX748IiLC0tKy7tY2bdpoaWlFR0dTq48fP7a1taWWy8vLhw8fbmlp+erVKz09vVGjRlVWVtbTkZGRka6uruiiU7wpAACAr9TE0ZU3btyYMGHC9OnTX758GRISEhISIhAICCG7du2iZj5js9kzZsxYsmRJTEzMnj174uPjJ06cSO0bEBDQvn37AwcOsFisQ4cO6enpLV++nNrE4/FCQkKuX78uEAhCQkKuXLlCCJGXl585c+YPP/wQExMTFBT09OlTX19fCRw6AAAAISwLCwtCyJAhQ/bt29fw3crKykaNGvXq1atXr15RJZ6enoQQQ0NDKiMSQn7++eeNGzcGBAQYGhpGRkaKJv/csmWLpqYm9VxQTk7ur7/+Er0mX1lZGRISQrUWEhKira09bNgwQshPP/0kLy+/YMECfX39GzduaGlpSeTgAQAAGC4uLoSQvn37rlmzhu5gmk90dPTcuXPrjvGRcVwuV1VVle4oQFpwfuu3adOm1f8U8UfXN4gdpIibp7LOlltQ32xlUsJq2uRqAAAA3wbMwAIAADKN2dhZRgEAAL4lzIiICLpjAAAAoA2zR48edMcAAABAG2bXrl3pjgEAAIA2GCwDAAAyjRUYGEgI6dy5M6ZrAQAAGcTavHkzIcTT0xOJEAAAZBBLKBTSHQMAAABt8IwQAABkGhIhAADINCRCAACQaUiEAAAg05AIAQBApiERAgCATEMiBAAAmYZECAAAMg2JEAAAZBoSIQAAyDQkQgAAkGlIhAAAINOQCAEAQKYhEQIAgExDIgQAAJmGRAj/p6KiIi4uju4oQIri4uIqKirojgKgxflGEmFRUdHGjRvnzZsXGhpKdyytVUxMzIoVK+iOAqRoxYoVMTExdEcB0OJ8C4mwpqZmwIABCQkJXbt2Xb58+W+//UZ3RK2SUCgUCoV0RwFShFMM8EksugOQgCtXrpSXl584cYLJZFpaWvr6+s6bN4/F+hYODVqUnJycsrIyuqNousrKyqysrLS0NLoDabq2bdvKy8vTHQV8a76FbHHv3r2BAwcymUxCSP/+/fPz81NTUzt06NDMYdy/f3/t1l2t9w/uouKi1DfvXUeOpTuQpmtnbPDnvt3Sa9/KuhNRVGMwW+t9lLLysqkLl7PkWuv/el5JwW9bN82aNYvuQOBb01r/S4j78OGDiYkJtSwnJ6elpZWVlVV/IszMzExISOjYsaOoxNzc3NDQ8GvCSEpKepD0gbCVvqYRGgmJUMhQiXr2L92BNJlQ8eFDUlUuvQ5q5Dk19mMJs7VekQirq6pZrBpGa03kwrTHwcHB0nvM+fr1a0apmsLltVJqX9oEQkENny8vz6Y7kKbiVzLYStOmTZNsq05OTpMmTaq/zreQCNlsdk1NjWi1urpaQUGh/l369u07ePBgNTU1UUmnTp20tbW/JowePXr4fs3+dKuurs7KymrXrh3dgbRcPXr0oDuEr5Kenm5kZNSaby2aSLX1Hj16TJZqB1LG4/Fyc3Pbtm1LdyBfYeSPEm+yc+fOX6zzLSRCIyOjN2/eUMulpaXFxcXGxsb176Krq3vx4kXphwYAAC1da71JIm7kyJHXr18vKioihJw9e/a7777DZQ0AADTQt3BF6ODg4O7u7ujo2L179+vXr586dYruiAAAoNVgfDPvFT169OjDhw89e/Y0MjKiOxYAAGg1vp1ECAAA0ATfwjNCAACAJkMiBAAAmfYtDJaBRvntt9969uzp6OiYmJi4a9eusWPHOjs7U5t+/PHHyZMnW1tbE0LOnz9/4cKFwsLCNm3aeHh4DBs2jNao4RNqamqOHTsWGRlZWlpqamo6ZsyYfv36fbLmTz/9JC8vv3r16vT09A0bNvTs2XPKlCnUptzc3NWrV3/33Xfz5s0jhCxevHjJkiVffAEJpOru3bunTp16//69kZGRm5vbiBEjmEzm2rVrKysrN27cSNWJiop6+fIlddYIIdHR0YcPH87IyDA2Np48eXLfvn2p8r/++uvRo0dycnJGRkbt2rUbNGiQ6OReu3YtPDxc1KmCgsLWrVub8ShbEFwRypy///77xYsXhJB37979+eefc+fO5fPgPmQWAAAgAElEQVT51Kbg4OD3798TQvbv379w4cLBgwcvWrSoa9eud+/epTNi+IxVq1bt3LnT09Nz/vz5ZmZm9+/f/1xNNze3IUOGEELy8vKOHTu2evXqqqoqatPJkyeDg4MjIiKo1b179+bm5jZD8PA527ZtGz58eLt27ebMmePg4LBjx47Tp08TQkJCQrZs2SJ6AfrFixeiNBYcHDxw4MA2bdosWLCgQ4cO7u7u+/fvpzZFRUVlZGTY2tqy2ezw8HBLS8ugoCBq08OHDx88eNDxP80/LWXLgStCmWZubq6pqfnXX39Nnz5dvPzs2bOLFy+eOHEiIWTQoEE0RQdfcPbs2a1bt44aNYoQ4uLiIioXCATnzp17/vy5kpLS6NGjqakERcPilJSUevToce3aNWrHo0ePjhkzJi8vj44jgNqSkpICAwPv3LnTp08fqmTKlCnUS9KEkHHjxgUGBnp4eIh/VODjx4+zZs3av3+/r68vIWTw4MEdO3b09PQcOXIkNW2kra3tjBkzqMq3b992cXFxcnKi5luxtLScPXt2cx5gy4QrQpnGYDA2bdr0888/l5f/zxSdRkZGFy5cePXqFV2BQUMYGRkFBwe/fftWvFAoFI4ePXrv3r1mZmaKioqRkZGEkHPnzolPpeTn53f06FFCSGxsLIvF6tKlS/MGDp91/vx5W1tbURYkhDCZTNHsjz4+PhwO5+TJk+K7UHc+J0yYICpxd3c3NDS8evVq3fYHDBjQvXv38+fPU6sCgaDyP6KbBDIIiVDWOTk5derUad++feKFmzdvVlNT69Kli56enq+vb3JyMl3hQT0OHDiQk5NjYWHRpk2bWbNmffjwgRBy/fr1uLi4iIiIqVOnLl68eP78+XV3HDZsWExMTE5OzpEjR/z9/Zs7bvi8jIwMc3Pzz21lMBjr1q1bs2YNj8cTFaakpJiamjL/96Mo5ubm6enpn2zE0tIyMzOTWj5//rzOf9zd3SVxBK0SEiGQjRs3bty4UXT7hRBiaGh48eLFgoKC48ePl5SU9OvXj8vl0hghfFKHDh1u376dk5MTFBT04sULV1dXgUDw/Pnz3r171z/vPIvF8vb2Pnr06Llz53x8fJotYPgiJSWlkpKSeiq4u7ubm5v/8ccfohINDY26u5SUlCgpffpLOMXFxRwOh1r28vIq/Q9180A2IRECsbe3Hzhw4M6dO2uVq6mpDRky5OzZs0VFRXFxcbTEBl+ko6MzcuTIkydPvnz58v379yoqKg35erC/v//PP//s6Oioo6PTDEFCA/Xo0ePp06f1n8FNmzatW7dO9Ldp586d09PTMzIyRBXy8/MTExM/+bGUysrK2NjY1v4dFYlDIgRCCFm3bt3u3bs/fvxIrUZGRoqeGsbFxfH5/Nb9bZdvVFhYWHV1NbUcHR3N4XD09PScnZ1v3779+vVrqryysvKT+3bp0mX//v2//PJLM8UKDTNy5EgdHZ158+aJTtytW7eioqLE6/To0aNXr16icaFdu3YdMGDArFmzSktLCSE8Hm/evHk2NjYDBw6s1fjr16+9vLz09PS8vLykfyitCRIhEEKItbX1mDFjRHdHz5w506ZNm549e/bu3Xvw4MHbtm3DBz1aoN27dxsaGvbp08fe3n727NmHDh1SUFCwtrbesWOHk5PT0KFDe/fuvWrVqs/tPnny5E6dOtUtd3Z21vvP8ePHpXkEUBubzY6IiMjIyDAyMurTp4+pqemSJUvqXrVv2LCBetOJcubMGRaLZWJi0rdv3zZt2uTn51+8eFH01HD//v0WFhZKSkqDBw+2tLS8ffu26M75hQsX9MRQqVQGYa5RmcPlchUUFNhsdnV1dXl5ubq6OlVeVVVVVlamqqpKjcz++PFjeno6k8k0MzNTVlamNWT4rIKCgnfv3rHZbAsLC/HnguXl5cnJycrKypaWloSQiooKBoOhqKhYU1PD5XI1NDTEG+HxeNXV1SoqKoSQ4uJi8d8JysrKbHar/eJ5a/bhw4esrCx9ff02bdpQJSUlJRwOR/Rd5Y8fP8rJyVFnTXwXPT098fs3ZWVl1HBQFRWVWt9krqysrKioEC/R1NSU0uG0cEiEAAAg03BrFAAAZBoSIQAAyDQkQgAAkGlIhAAAINOQCAEAQKYhEQIAgExDIgQAAJmGRAggW+7cuSPBz4k8f/48JiZGUq0B0AKJEKDFefjwob29vb29/ZMnT8TL+Xz+wIED7e3tN23a1JB2srOzQ0JCiouLxQvHjBkj/u2Cr7R69erFixdLqjUAWiARArQ4JSUlsbGxCQkJf/75p3h5RETEvXv3nj17Jv6pgXo8f/587Nixb968kU6YAN8IJEKAFsrDw+PMmTOiz4AQQo4cOTJ06FBqMti6+Hx+Y7sQCASf21RTU1PP/Is1NTWN7QugxUIiBGihfHx8+Hz+hQsXqNXCwsLLly/7+fnVqpadne3v76+mpiYvL29gYLB27VoqS4WFhQUEBBBCJk6cSN1oTUxMFO21d+9eY2NjeXn5jh07nj9/XrzBS5cu2draslgseXn5gQMH1noEePDgQWNjYxaLZWZmduLECWkcOEAz+/SflgBAOzU1tZEjRx49enTChAmEkJMnTyorKw8fPly8TnFxcb9+/YRC4Z49e8zNze/fv7927dqSkpLt27fb2Nh4e3tv2LDB39/fzMyMEGJgYEDtFRoaqq+vv3PnTgUFhY0bN/r6+qalpenr6xNCLl++PHr0aBcXl6tXr1ZUVKxZs2bAgAExMTHW1taEkODg4BkzZvj4+MyYMSM7OzswMLC6upr6wAVAKyYEgBYmPDycEBIWFhYREcFkMjMyMoRCYffu3efPny8UCtls9uzZs6maP/30E4fDSUtLE+27efNmNptdWFgoFArDwsIIIU+fPhVvXEdHx9jYuLS0lFpNSUkhhBw8eJBa7dq1q6WlZVVVFbX64cMHDofj6+tLrbZv375Hjx4CgYBajY2NJYRQmRig9cKtUYCWy8XFxdjY+Pjx4y9evHj69Gnd+6LXrl3r0qWLUChM+4+xsXFVVdWrV6/qb1b0jUkLCwsVFZX09HRCCI/HS0hImDBhgujDdQYGBq6urrdu3SKElJSUJCUl+fj4MBgMamv37t07dOgg2UMGaH64NQrQcjGZzEmTJh05ciQvL8/GxsbOzq5WhaysrKysLAsLi1rl4p8vr6vWF88VFRV5PB4hJDMzUyAQGBsbi29t27bttWvXCCHv3r0jhBgZGdXaSu0L0HrhihCgRfPz80tJSdm/f//3339fd6uKisrQoUPr3urx9vZuQl9UgiwsLBQvLCwspD6DTm0tKioS31pQUNCEjgBaFCRCgBatffv248ePt7Gx8fX1rbu1T58+Dx48yM3N/eS+1P3Phl+xqaurGxsbR0ZGikp4PN6dO3ccHBwIIXp6etra2lFRUaKt2dnZ4iNRAVopJEKAlu7kyZMxMTHUqM5ali1bVlNT4+XlFR8fLxAIqqqqnj9/LprqxcLCQk5O7vjx44mJiWlpaQ3JiPPmzbt58+aWLVsqKiqKi4tnzpyZlZW1cOFCQgiDwZg1a1ZoaOjhw4f5fH5+fv6UKVOEn3/XEKC1QCIEaMWsra3Dw8Nzc3O7du2qrKzM4XC6d+9+//59aquRkdHmzZtPnTplY2NjYWERHx//xQaXLl0aEBCwYsUKdXV1bW3tkJCQ3bt3u7u7U1t/+uknT0/PqVOnqqqq6unp6enpubi4SPHwAJoFA3/QAbQ01dXVpaWlKioqotGb4oqKihQUFJSUlEQlAoHgxYsXmZmZ6urqlpaWda8duVwun89XVVVlsVjFxcUKCgocDke0tW5Jdnb2s2fPFBQU7O3t1dTUarWWmJiYnp5uaWlpZWVFvYahqqoqgcMGoAkSIQAAyDTcGgUAAJmGRAgAADINiRAAAGQaEiEAAMg0JEIAAJBpSIQAACDTkAgBAECmIRECAIBMQyIEAACZhkQIAAAyDYkQAABkGhIhAADINCRCAACQaUiEAAAg05AIAQBApiERAgCATEMiBAAAmYZECAAAMg2JEAAAZBoSIQBIUX5+fnh4+KtXr+gOBOCzkAgB/j8ej3f48GFvb+/OnTsbGRlZWFi4uLj89NNPL1++/GT96urqI0eOjBkzxsbGxtjYuEOHDp6enn/++SeXy23myFus6OhoNze3oKAgugMB+CwkQoD/8/DhQysrq6lTp549e7akpMTIyEhJSenx48e//vrrd999N3ny5Fr14+LiOnbs+P3334eGhpaVlRkaGjKZzMuXL0+fPt3Kyio2NpaWowCAxkIiBCCEkMePHw8aNOjdu3dTpkxJTU3NyMiIiYlJSEgoKCgICwvr169fVFSUeP3Xr1/3798/NTV17NixKSkpb9++jYmJefXqVW5u7q5duwQCQXJyMl3H0qK4uLjk5eVt3ryZ7kAAPoshFArpjgGAZtXV1dbW1mlpaStXrly3bl3dCkKhMCQkZOzYsaLVHj16xMTETJs27cCBAwwGo1b99+/f5+Xl2drafrFrHo+Xn5+vqKiora39yQpVVVW5ubkKCgq6urpfbCovL09PT4/NZouXZ2dny8vLf659Cp/Pz8nJUVVVVVNT+1ydioqKnJwcbW1tVVXV+iOhWuNwOFpaWvVUKywsrKio0NHRUVBQ+FyPeXl5qqqqmpqa9fdYUVGRm5urr6+vqKhYf02ATxACyLxTp04RQkxNTaurqxtSPzIykhCiqanJ5XKb3GlsbOywYcNEv7hVVVV9fHxyc3NFFdLT0318fDgcDlWhbdu2W7durampEW/ExsamR48eRUVFfn5+VDpRUFBYtGgRn88XCoV79+41NDSkdnd0dExOThbf18fHx9zc/M2bN0uXLqXyH4PBGDRoUGJioni1xMTEqVOnGhsbi35pdO7c+dSpU7UOp2/fvpaWlqWlpfPmzVNRUSGEODk5CYXC27dvm5ub//zzz6Kaubm5fn5+4tnU0tJyx44d4q3Fx8e7ubnJy8tTFaytrY8ePSpeIScnx9zc3MvL68OHD15eXiwWixDCZrMnTZr0NScFZBNLumkWoDUICwsjhIwbN476ffpFV69eJYSMHj2a+o3fBFeuXPH29ubxeG5ubo6Ojjwe79WrVxcvXlyxYgV15ffu3btevXp9+PChf//+gwcPLikpOXbs2NKlS+Pj448dOyZq582bNxwOZ/jw4VlZWdOnT2cymadPn965c6eSkhKHw1m3bp2Xl1fbtm2joqIePnw4bNiwhIQEUXbJyspKS0ubOXNmTEzMzJkz9fX1o6Kirl696uTk9PjxYzMzM6raP//8ExwcPHToUEtLS1VV1ZSUlDNnzkyYMKG0tHT69OmiSDIyMjIzM8eNG/fo0SN3d3djY2PqQrm8vDwtLS0vL4+qJhAI3N3dY2Ji+vfv7+zsrKys/O7du0ePHt28eXPRokVUndjY2IEDB3K53BEjRjg6OmZlZR09etTPz+/Nmzdr1qyh6tTU1KSlpXE4nIEDBxJCZs+ezWAwQkNDjx8/zuPxzpw507TzAjKK7kwMQL8uXboQQs6cOdPA+oMGDSKE7Nu3r2nd5efnq6urM5nM0NBQ8fKCgoKSkhJqedSoUYSQOXPmiLZmZ2e3bduWEHLx4kVRoZKSEiGkX79+osugpKQkFoulpKSkrq7+/PlzqrC6urp3796EkMuXL4v27d+/PyFETU0tJSVFVLhy5UpCyOjRo0Ul6enpxcXF4nEmJiaqqqrq6elVVVWJCk1MTAghlpaW4he1QqHw2rVrhJC5c+dSq8+ePSOEDBkyRCAQ1Dp2akEgEFCnY9u2beI9qqioMJnMuLg4qiQrK4v6DTZmzJjKykqqMCcnR0NDg8lkZmZmCgEaDINlAEhhYSEhpNZTtMrKygP/6927d9SmgoICQoiOjk7Tujt69OjHjx/9/Pw8PT3Fy7W0tKgbhrm5uZcuXVJRUfnll19EW/X19QMDAwkhBw8erNXgjh07RNemVlZWvXv3Li8vnzFjRteuXalCFovl6+tLCElMTKy174wZMywsLESrgYGBmpqaly5dEl3DmZiYqKuri+/SsWPHUaNG5ebmPn/+vFZra9asqf9ZZmVlJSFER0en1oNV0dPE6Ojo+Ph4MzOzBQsWiPc4c+ZMgUBw+PBh8b1YLNbevXtFjxj19PS8vLwEAsHr16/riQGgFtwaBSBMJpMQIvzfgWNcLnfmzJniJdeuXaOuyT5Zv+Hu379PCHF1df1chZiYGIFA4OjoWCs3Dx8+fO7cuU+ePBEvZLPZ3bp1Ey9p06YNIcTBwaFuYWZmZq2+hgwZIr6qoqLSt2/fy5cvx8bGDh06lCr8+PHjuXPnXr9+nZWVxePxCCEvXrygWqvVy4ABAz572IQQQjp37qyjo3Pq1Ckmk+nt7d2/f/9aWfbx48dUVLVuUw8bNmz79u21jt3U1FRfX79WCSFE9CcLQEMgEQIQHR2djIyM7Oxs8UJ1dfUbN25Qy8uWLaPu6VGoi54PHz40rbucnBxCiOghXF35+fmEECrpijMyMmKxWPn5+UKhUHRFpampKScnJ16NukKqdcFKDSWtqqqq1SaVIMVR/ebm5lKrDx8+HDlyZF5enrKyspGRETWspri4mBBSUVEhviODwRCNzfkcFRWVs2fP+vv7Hz9+/Pjx43Jycr179540adLUqVOpPy+oY68bFXXrVXSdSql79UkdJpWtARoIt0YBiL29PSGk7pWWy39qvQZA1X/06FHTuqN+WZeVlX2uAjWepVaaIYTweLyamhp5efm6L2w0Wd1eqMCobCoQCCZOnFhUVPTnn38WFBQkJSXFxMTExMSMGzeublMMBqNWSv4kJyenlJSUu3fvrlq1qmfPng8ePJgxY4aoQerYqTuo4srLy0VRAUgWEiEAGTlyJCHk9OnT9SQncdRIlosXL4oumxrF0tKSEJKUlPS5CtT9vZSUlFrl1PsP1FZJSU1N/WRJu3btCCGvX79+8+ZNv379pk6dKp6E6gm+IeTk5Pr16/frr7/ev3//5cuXRkZGZ8+epeaxo46u7nQEVAkVFYBkIRECEDc3N3t7+7y8vNmzZ9fU1HyxvoODw9ChQysqKqZMmVL3ZiMh5OnTp/fu3fvc7sOHDyeEHDp0iM/nf7JC9+7dtbW1o6Oja41GOXDgAKn34WIT/PXXX+Kr//7774MHD7S1tbt3704IEQgEhJBab6mnpKRQb1JKRIcOHZycnMh/D/YGDhzIYrEuXbpE3UAWkcaxA1CQCAEIg8EIDg7W0dE5fvz4wIEDw8PDRTcM8/LyDh06VHew5eHDh9u2bXv16tV+/fqFhYVRt/IEAkFsbOy8efMcHR3rDksR8fDw6Nu3b2xs7IQJE0QPJsvKykQDUxUUFBYvXkwI8ff3T0tLI4QIhcLjx48fPHhQSUlp4cKFEjz2yMjI3bt3UwkvOzv7+++/FwgECxcupO7ftm/fXklJKTIyMiIigqr/5s2b8ePHN7m7iIiIlStXil+GPnv27MaNG/Ly8tQYV2NjYz8/v4qKikmTJlFPBPl8/saNG8PDw42MjL7//vuvOViAT6P15Q2AFiQtLa1v376i/xoaGhqim4HKysrLly8vLS0Vr//+/XtnZ2dRfU1NTdFARwcHh1rTuNSSnZ3ds2dPQoicnJyVlVW7du2op2ui9+Sqq6snTJhACGEymdbW1tTYSA6Hc/78efF2lJSU9PX1azVOZYtbt26JF4aHhxNCZs6cKSqh3iPcsWOHvLy8vr5+ly5dqOQ3fPhwHo8nqrZz507qoMzMzDp16sRkMrt3706923Dy5ElRNRMTEyaTWfdIa71HePLkSao1NTW17777jhoxJCcnFxQUJNqlpKSEek2ezWZ36tSJml9NR0fn0aNHojrUe4SOjo61uqMmNd2/f/9n/+kB6sBcowD/4969e2FhYa9evSopKVFQUDAxMXF0dBw5cmStUf4iDx8+vHLlSmJiIpfL5XA4nTt3pi74vthRTU1NaGhoWFjYu3fvlJSUTE1Nhw0b5uzsLD7eJDw8PDQ09M2bN4qKit26dZs+fTo1eFJkx44dbDZ73rx54oXXrl1LSEgYP368eOU3b978/ffftra2ovclnJyc7t69m56enp+ff/DgwdTUVFVV1REjRkyaNIkawCly/fp1Kgw1NbVBgwZNmTIlOjr6wYMHw4cP79SpE1Xn999/53K5y5Ytq3WYaWlpp0+ftre3Hzx4MCGkoqLi7t27d+7cSU1NLSgo0NLSsrS0nDhxoo2NjfheAoEgJCTk6tWrWVlZysrKjo6O06dPF3+ZpLS0NCgoyMjIaNKkSeI7Pn78+Pbt20OGDGnIRK8AFCRCABklSoS1kiuArMEzQgAAkGlIhAAAINMwswyAjGrfvn1ZWVmtjxcCyCA8IwQAAJmGW6MAACDTkAgBAECmIRECAIBMQyIEAACZhkQIAAAyDYkQAABkGhIhAADINCRCAACQaUiEAAAg05AIAQBApiERAgCATGs1k24XFhYWFxebmZkxGIy6W9++fSsQCKhlFRUVPT295o0OAABaq1Yw6XZiYqKHh0dmZiafz6+urmaxPpG8tbS0VFRU5OXlCSHDhg3btWtXs4cJAACtUitIhIWFhenp6YqKip06daonEUZHR1tYWDR/eAAA0Kq1glujWlpaWlpaKSkp9Vd7+/Ytj8ezsLBQUFD4YptZWVn79+//5ZdfJBTjN6K6upq6qoZvEs5v/d6+fXv2cbKgpV8a1KempkZOTo7uKJrOQovtNdip+fttBYmwIdhs9qJFi3g8XmFh4cGDB0eNGlV//fT09ODg4M6dO4tKevXq1aZNGymH2dKVl5erqqrSHQVIC85v/YKDg1dffF7T53u6A/lKNXQH0FQVJSqn543Oz5Zsq0zml8eEfiOJMCkpSU1NjRBy4sSJyZMnp6ena2pq1lOfy+UWFBScOXNGVKKoqKijoyP1QFs2Ho+HK4ZvGM5v/fh8PkPblHR0pjsQWcXNI0RYWVkp2VbZbPYnH6iJ+0YSIZUFCSG+vr4LFix48eJFv3796qmvqalpaWkZGhraLNG1GjU1NUpKSnRHAdKC81s//JVAPwaDlh/Rb+09wpycnJKSEl1dXboDAQCA1qEVXBFWVVXt3LmzoKCAELJ161ZlZeWAgABCyMiRI+3s7H766ac7d+6EhIQ4ODhUVFQEBQV5eHhYW1vTHTUAALQOrSARUrS1tTdt2iReMm7cOCMjI0JI+/btDQwMoqKiOBzODz/8MHHiRJpiBACA1qcVJEI2m718+fK65RMmTKAWDA0NV61a1bxBAQDAN+Jbe0YIAADQKEiEAAAg05AIAQBApiERAgCATEMiBAAAmYZECAAAMg2JEAAAZBoSIQAAyDQkQgAAkGlIhAAAINOQCAEAQKYhEQIAgExDIgQAAJmGRAgAADINiRAAAGQaEiEAAMg0JEIAAJBpSIQAACDTkAgBAECmIRECAIBMQyIEAACZJvlE6O3tPWLEiOrq6lrlNTU1Xl5eEyZMkHiPAAAATSbhRBgdHX327Nlly5bJy8vX2iQnJxcYGHj69OnY2FjJdgoAANBkEk6EYWFhlpaWffv2/eRWBwcHGxubiIgIyXYKAADQZBJOhOnp6R06dKingrW1dXJysmQ7BQAAaDIJJ0IWi1VZWVlPhYqKCoFAINlOAQAAmkzCidDKyiomJqaiouKTW3k83uPHj9u3by/ZTgEAAJpMwolw9OjRZWVlK1eu/OTWtWvXFhYWenp6SrZTAACAJpNwIrSwsPjhhx927tw5YsSIGzduFBQUEEIKCwtv3rzp6em5cePGgICAjh07SrZTAACAJmNJvMX169fLyclt2bLl8uXL/9MTi7V06dINGzZIvEcAAIAmk3wiZDKZ69atmzFjxrlz5+Li4srKypSUlLp27Tp69GhTU1OJdwcAAPA1JJ8IKSYmJgsXLpRS4wAAAJIirUQowuVy+Xy+aFVZWZnNZku7UwAAgAaS/FyjK1euHD9+fFVVFbU6bNgwLTG+vr4S7xEAAKDJJHxFmJ6evmnTpkOHDolf9tnb2zs7OxNCCgoKDh06lJSUhFcJAQCghZBwIrxy5YqqqurEiRPFC/v27btp0yZCiEAguH79+pUrVxYvXizZfgEAAJpGwrdGExMTu3fvXvfTE//XGZNpY2OTmJgo2U4BAACaTMJXhGVlZSoqKuIlBw4cUFJSEq2qqKiUlpZKtlMAAIAmk/AVoZaW1rt378RLOnTo0LZtW9Hqu3fvtLW1JdspAABAk0k4Efbs2TM+Pv7169ef3JqWlhYTE9OzZ0/JdgoAANBkEk6EI0aMMDAwGD9+fFZWVq1N2dnZPj4+WlpaXl5eku0UAACgyST8jJDD4Zw8edLNza19+/ajR4+2t7fncDiVlZUxMTEXLlyoqqq6ePGisrKyZDsFAABoMsnPLDNgwICHDx8uXbr05MmTJ06coAoZDMaAAQO2bNlib28v8R4BAACaTCpTrNna2t64cSM/P//FixdFRUUaGhqdO3fW1dWVRl8AAABfQ4pzjero6AwYMEB67QMAAHw9yc81SpeioqLXr1+LT/ANAADwRa0gEV6+fNnOzo7NZvft2/dzdTZt2mRubu7t7W1lZZWQkNCc4QEAQKvWChJhmzZttmzZsnHjxs9VSE1NXb9+fUxMTEJCwtSpU/EdRAAAaLhWkAi7devm7Oysrq7+uQqnT592cXGxsLAghMyaNev27dvZ2dnNGCAAALRiEk6E3bt3v3PnDiHk8ePHFy5ckGzjn5Oenm5paUkt6+joqKmpZWRk1L+LUCgsLy+PFfPx40fpRwoAAC2OhEeNZmZmcrlcQsjt27ejo6NHjRol2fY/qbS0VE9PT7SqpKRExVCPgoKCt2/fTp8+XVSyaNGi0aNHSyvEVqKsrIzBYNAdBUgLzs7X4J8AABvVSURBVG/9RJ8TB9oIhRL/KoOioiKL9YVMJ+FE2K5du3Pnztna2lZUVFRVVRUVFdWtw+FwFBUVJdipvr6+qCOhUFhUVKSvr1//Ljo6Op06dXry5IkEw/gGCIXCWh8PgW8Jzm/92Gw2IWV0RyHbGAxafkQlnAiXL18+YcKEv/76i1rV0tKqW2fbtm1LliyRYKe2trZ79uyhluPi4thsNvW8EAAA4IsknAjHjBljb28fExMTHBz89u3bmTNn1q3Tq1evRrWZnp4eERFx79697OzsAwcOmJmZubq6EkK6dOnyxx9/ODo6jh07dsWKFevWrXN2dl6+fPm0adM4HI5kjgcAAL51kp9ZxtTU1NTUtKys7OXLlzNmzPj6Bj9+/BgbG8vhcJydnWNjYwUCAZUI+/btq6GhQQjhcDi3bt369ddfb9686erqunz58q/vFAAAZIS0pljz8/OTVFPUlV/d8n379omWra2tT548KakeAQBAdkjxPcLq6uqdO3c6Ozubmprq6uo6Ojru2rWrurpaej0CAAA0lrQSIZ/Pd3d3X7x48YcPH3r27Dlo0KCPHz8uXLjQw8OjpqZGSp0CAAA0lrRujQYHB9+6devo0aOTJk0Svbp0/Phxf3//4OBgX19fKfULAADQKNK6Irx165abm9vkyZPFX+CdNGnS0KFDo6KipNQpAABAY0krEVZVVVFDOmvR1NSsrKyUUqcAAACNJa1E+N1334WFhaWlpYkXvnnzJiwsrGvXrlLqFAAAoLGk9Yxw+vTpe/bssbOzmzJlCpX54uLijhw5oqioOG3aNCl1CgAA0FjSSoTa2tp37tyZP3/+zp07hUIhIYTJZA4ZMmT37t2fnHcNAACAFtJKhIQQS0vLsLCw4uLi1NRUQoiFhcUnnxoCtBahoaEFBQV0R9F0lZWVkp3vvvkNHTrUxMSE7ijgWyPFREjR0NCws7OTdi8AzWDiZH+hvRdDTp7uQJpIKCSt+itMwpSHWysrAwIC6A4EvjVST4QA35KqsTuJAr5kRA+Fs0uo5ywAkiXFKdYAAABaPiRCAACQaUiEAAAg06SVCPfs2bNq1SopNQ4AACAp0hosEx0dzePxpNQ4AACApEjrirB3795xcXH44hIAALRw0kqEU6dONTU19ff3rzXdKAAAQIsirUS4b9++R48enThxwsLCgsPhaIkJCgqSUqcAAACNJa1nhLa2trNmzfrkpi5dukipUwAAgMaSViJ0cnJycnKSUuMAAACS0hzvEX748KG0tLQZOgIAAGgsKSbC1NTU4cOHq6qqGhkZXblyhRDy7Nkze3v7nJwc6XUKAADQKNJKhLm5uf369fv3339/+uknU1NTqrBr1655eXmXLl2SUqcAAACNJa1nhAcOHJCTk4uOjlZXVz99+jRVyGQy7ezsEhISpNQpAABAY0nrivDly5eurq7q6uqEEIbYN9B0dXXz8vKk1CkAAEBjSSsRKisrFxcX1y1PTk7W1taWUqcAAACNJa1EOHDgwGvXrsXHxxOxK8Lr16/fvn3b1dVVSp0CAAA0lrQS4bhx4xwcHBwdHadPn56dnX358uXx48e7u7sPGjRo+PDhUuoUAACgsaSVCFksVnh4+Jw5c86fP5+ZmXnq1Knw8PCAgIBLly4xmfgIIgAAtBTSGjVKCFFWVt66deuWLVtycnIEAoG+vr6cnJz0ugMAAGgCKSZCCoPBoEbHIAsCAEALJMW7lIWFhQsWLDAzM1NSUlJSUjI3N1+4cGFhYaH0egQAAGgsaV0RFhUV9erV6+3bt8OHD/fz8yOEJCQk7Nu3Lyws7PHjxxoaGlLqFwAAoFGklQj37t2bk5Pz+PHjbt26iQqfP3/ev3//vXv3rlq1Skr9AgAANIq0bo0mJiYOHz5cPAsSQmxtbUeMGJGYmCilTgEAABpLWomwTZs2PB6vbjmPxzM2NpZSpwAAAI0lrUQ4Y8aMqKioa9euiReGhYVFRUVNnz5dSp0CAAA0loSfES5atKi8vJxa7tixo4eHh52dnY2NDSHk5cuXsbGxffv2TU9Pb9++vWT7BQAAaBoJJ8KLFy+Kz7WtqamZlpaWlpYmWn358mVSUhKmGwUAgBZCwolQlPMAAABaBUz7CQAAMk26U6yVlZW9ffs2OztbKBSKCq2srNq1ayfVfgEAABpIWolQKBQuX7587969FRUVtTZt27ZtyZIlUuoXAACgUaSVCA8dOrR169Y5c+a4ubkpKiqKb7KyspJSpwAAAI0lrUR47969wYMHBwUFSapBoVBYWVnJ4XAk1SAAAACR3mAZTU1NVVVVSbW2f/9+bW1tQ0PDQYMG5ebm1q1gbm6u9Z9p06ZJql8AAPjmSSsRTps27Z9//klJSfn6ppKSkpYtW3b37t2CggITE5OlS5fWrVNcXBwZGZmampqamrpr166v7xQAAGSEtG6N2tjYBAcH9+nTZ+TIkaampuJf5R00aJCDg0PDmzp+/Li7u3vnzp0JIcuWLbOzs9u/f3/de6Tq6uqampoSCR4AAGSHFL9HGBgYmJube/DgwVqbtm3b1qhEmJqaam1tTS1bW1tXVVW9f//e0tKyVrU+ffrU1NT07Nlz+/btHTp0+GKzNTU1RUVFolVVVVUWS7ovkwAAQAskrV/9u3btSkhIOHnypJub21deqBUXF6uoqFDLTCZTSUlJPIFRQkND7ezsKioq1q5d6+7u/vLly1pDVWvJy8uLj483NzenVhkMxubNm318fL4mzm9AaWkp3SG0aMIvVwFpEgp5PB6Xy5VS85/8YA40K6FQ4udXUVFRXl6+/jrSSoRJSUkjRoyYMGHC1zelo6NTUlJCLVdXV5eVlenp6dWqM3DgQEKImpranj17NDQ04uPje/ToUU+burq63bp1e/LkydeH942R4BCnbw+D7gBkHYOhoKAgvR9RBQUFQsql1Dg0CINBy68gaQ2Wad++vSh7faVOnTo9ffqUWn727Jm6urqRkdHnKvP5fIFAIP5IEgAAoB7SSoRz5sz5999/z549+/VN+fn53b17NzQ0NDMzc9WqVd9//z11nbty5cqjR48SQl68eBESEpKWlvby5Ut/f38zM7Pvvvvu6/sFAABZIK1bow8fPjQxMfH29u7YsWP79u3ZbLZo0+TJk4cNG9bwpgwNDUNDQ1evXl1QUDBkyJD169fXqsBgMA4fPvzjjz8qKir27NkzPDxcvDsAAIB6SCsRFhcXl5aW2tnZEUIyMzNrbWpsa66urnU/YSjKiDY2NmFhYU2NFAAAZJq0EqGfn5+fn5+UGgcAAJAUfI8QAABkmrSuCN++ffu5+dXwPUIAAGg5pJUIz5w5ExgY+MlN+B4hAAC0HFJ8Ruji4iJarampefPmzZYtW6ytrX19faXUKQAAQGNJKxEaGBgYGBiIl/To0cPNza1Tp07Pnz8fMmSIlPoFAABolGYdLKOmpta/f3+JvGUPAAAgEc09ajQjI4PBwJSNAADQUjTfqNHCwsILFy7cv39/wYIFUuoUAACgsZp11Kiqquqvv/7q7e0tpU4BAAAaq5lGjRJCtLS02rRp88XvQgEAADSn5hs1CgAA0AJhijUAAJBpEr4iPH369LZt2+qvs3DhQrxTDwAALYSEE6GioqKmpubntj579qygoCAnJ0eynQIAADSZhBPhqFGjRo0aVbf86dOnP/74Y0FBQceOHfv37y/ZTgEAAJpM6s8I3759O3PmTAcHh5cvX/7xxx/x8fEODg7S7hQAAKCBpDVqlBDy/v37X3755fDhw+rq6hs2bFiwYIGioqL0ugMAAGgCqSTCwsLCLVu27N69m8lkLlmyZMWKFerq6tLoCAAA4CtJOBGWl5fv2bNn06ZN5eXl/v7+a9euxduEAADQkkk4Ef7555+BgYFqamo//PBDu3btLl26VLdOr169unTpItl+AQAAmkYqt0ZLSko2bNjwua3btm1DIgQAgBZCwonQx8enT58+9ddp06aNZDsFAABoMgknQj09PT09Pcm2CQAAID2YaxQAAGQaEiEAAMg0JEIAAJBpSIQAACDTkAgBAECmIRECAIBMQyIEAACZhkQIAAAyDYkQAABkGhIhAADINCRCAACQaUiEAAAg05AIAQBApiERAgCATEMiBAAAmYZECAAAMg2JEAAAZBoSIQAAyDQkQgAAkGlIhAAAINOQCAEAQKYhEcL/yczM3Lt3L91RAAA0NxbdATTI+vXrDx48SAiZOXPmihUr6laIjIxctmzZ+/fvBw0atG/fPk1NzWaPkQiFwuLi4ubvV1Li4uIuXrw4a9YsugNpOmVlZTabTXcUANDKtIJEGBIS8ueff0ZGRgqFQldX1w4dOnh6eopXKCws9PLyOnTokLOz8+zZsxcuXHj06NHmj/PQoUOz5s6TV1Bq/q4lQiAU8Kv5Ru0s6A6kiYREqMxRLMj5QHcgANDKtIJEePDgwfnz51taWhJC5s+ff/DgwVqJ8NSpU926dRszZgwhZN26dZ07d96zZ4+amlozx1laWsoaOLNyzLZm7leyKukOoOkquSTQlO4gAKD1aQXPCF+9emVra0st29ravnr1qlaF169fiypYWlrKy8unpaU1a4gAANBqtYIrwoKCAtHlnbq6el5eXt0K1PWiqE5+fn79bWZkZERHRzOZ///vAA6H85WPl/h8fo1JD1bQ6K9phEZCoVAgEMjJydEdSJMJhar6Un08LFQ3YP3pSwhDel1IVU1NDZPJZDBabfxVFatW/fXzzz9LqX0+n09Me+G/MH2EVWw1if8Xnjp16rZtX7hR1woSoaamJpfLpZZLSkp0dHTqVigtLRWtfvz4UVtbu/42vby8srOzy8vLRSUaGhoSihcAAFqKhjwmawWJ0MrKKjExceDAgYSQxMRE8Ys/UYWwsDBqOSMjg8fjmZqafrFZfX19SUcKAACtTyt4Rujv7x8UFJSfn5+fnx8UFOTv70+Vz5kzh3peOH78+IcPH965c4fP569fv37EiBG0vD4BAACtUStIhJMnTx46dKiVlZWVlZW7u/vEiROp8nv37lHv7RkYGBw5cmTSpEkaGhppaWl79uyhNV4AAGhNGEKhkO4YAAAAaNMKrggBAACkB4kQAABkGhIh1KeiooLuEOALBAJBZeUXZgRKTk5OSUlpnnhAUsTf72og8RfJoOGQCGXO8OHDT548SQi5deuWlpbWL7/8Itpka2t79+5dQsjHjx/Hjx+voaFhYmKira29ZMkS2sKFz8vMzHRzc9PS0jI2NjY0NKznreEDBw4cPnyYEPL8+XMtLS0fHx/RpmfPnmlpaY0fP55aNTExefHihbQjh3qUlZXNnz9fW1tbX19fW1t74sSJGRkZhJDevXtbWlqKUt3BgwdF57G6unrFihU6Ojpt2rTR1NRcsGCB6G+jWbNmaWlp6ejoWFhYDB48OCgoqLq6mtq0ZcsWAzEWFq11nuGvh0QocwoKCqi/NKurq/l8/o4dO3JycqhNHz9+pP6TrFy5sqSk5N27d3l5ecnJyS4uLnRGDJ8xe/ZsY2PjnJycgoKC2NhYOzu7z9Vcu3btqlWrCDX/UU3NP//8k5ubS206cuSIlpaW6Ndrdna26BclND+hUDhq1Kj4+PjHjx9zudy0tLSePXvGxsYSQkpKSgoKCn777TeqZmVlpeiszZo169q1aw8fPiwuLo6Li4uOjp4wYQK1qby83M/PLzk5+ebNm3PmzPn99989PDwEAgEhhMvl9u7d+8V/njx5QscRtwhIhDLNwMDAy8trw4YNtcpfvHjh6uqqqqpKCNHS0nJzc6MjOviCFy9eeHh4KCgoEEKMjIyoSScIIampqePGjbO2tu7WrduJEycIIcePHw8ODqa2slissWPHnjp1ihBSVVV19uxZ8QtEoNf169cfPXr0999/UzOHqKurBwQEjB79f7O+LVu2bNu2baI/YiipqalHjhw5duyYlZUVIcTExOTEiRNXrlyJiYmhKigqKmpqapqamo4aNerOnTuPHz++cuUKtYnNZuv854sTcn3DkAhl3S+//HL06NFa05S7uLhs2LDhl19+uX37No/Hoys2qJ+Li8vixYu3bt36/9q705gmujUO4KeUVpFS1JbNSlHLoiJqkFarIiguNbhgjBilgKIRkRARJC6IJgTibjAxQIgBVIx7BIkBAr4KKpISEIl+EBSFyo4ssqgYnPth4tyGqtfXjXLn//s0fZ6Z0zMtzcPpTM959OjRwMAAHXz37t2CBQvmzJmj0WiysrImTZpECHnx4sXLly+ZA4OCguhvSrOzs5VKpf60hTBUCgsL5XL5t+a9mjZt2ooVK44fP64bfPz4sUQimTFjBhOZNGnS1KlTCwsL9VsQiURyufz+/fv0w9ra2gtf5OXl/b7zGGZQCNlOIpEEBQXpXikkhERHR586dUqj0axdu9bCwiIlJWWougffkZiYuHv37ry8vGXLlo0bNy4rK4sQkpWVZWdnt2vXLqFQKJVK586dq3+gi4sLl8utqKhIS0tjpmoCQ9DZ2WlhYfGdHeLi4lJSUrRaLRNpbGy0tLQctJuFhUVHR8dXW7Cysurq6qK3m5ub//mC/gKWnYbBXKPwp+3fv9/BweHJkydMhMPhqNVqtVo9MDCQnp4eHBysUqmkUukQdhL08Xi80NDQ0NDQ/v7+w4cP+/v7t7S01NXV6c/Hq2/z5s1HjhypqKhQqVRnzpz5C72FH2FlZVVZWfmdHSZMmODv7x8bGzt9+nQ6YmdnV19fP2i3+vr6bw0r6+rqPD096W2FQpGWlvarnR7+MCIEIhaLw8PDDx06pJ/icrmbN2/mcrn0fWtgmPh8/tatW7u7u1tbW21tbX/klxJ+fn6ZmZkbNmwwNsZ/wwZEpVJpNJrq6mrd4KDblw4ePHj9+nVmZVa5XP727ds7d+4wO5SXlz9//nzZsmX67VdXV2s0GpVK9Qf6PozhMwCEEBIREWFvb898YRIdHe3q6urm5sblcpOTk83NzV1cXIa2h6AvNDTU29t7xowZ/f398fHxjo6O48ePX7Nmzb59+xISEoKCguhbf5VKpf6xIpGoubnZxMREP1VWVsb8JTg5OdnY2PzZ0wAds2fPDggI8Pb2PnnypKura319fVpa2vz585k5lgkhlpaWO3bsOHnyJH07t42NTUxMjL+/f2JiokKhqKysDAkJ2blzp6OjI71/U1NTWVlZQ0NDRUVFQkLC9u3b582bR6daWlru3bvHtOzu7j6clzP8eSiErCOXyyUSCSFEJBIxF5DMzMyOHTuWkZExduxYQoizs3Nqamp0dLSRkZGLi0t+fr65uflQdhq+ZsqUKadOnaqrq+Pz+QqFIjc3l8PhmJmZFRUV7du3LykpSSAQREREKJVKBwcHeuQnFAo9PDzow5n3VCqVzpw5k9728vK6fPky8xQREREohH9ZSkpKSkrK0aNHGxoarK2tvb29V65cSQiZM2cOc/kwKiqqvLyceddiYmLs7OxOnDjR0NBgZWW1d+/ebdu20SlnZ+eCgoK9e/eamprKZLLMzEx3d3c6JZPJHj16FBcXxzz1rVu3Ro0a9fdO1WBg0m0AAGA1XCMEAABWQyEEAABWQyEEAABWQyEEAABWQyEEAABWQyEEAABWQyEEAABWQyEEAABWQyEEMDjv37+vqampqanRXyP3zZs3NTU1bW1tP914ZGTk9evXf62D/3XmzJlBqwIBDDsohAAGp6ioSCaTyWSya9eu6cbb29vt7e1lMtnBgwd/pJ2KigpfX99Xr17pBs+fP19SUvK7upqXl5ednf27WgMYEiiEAAbK0tIyPT1dN3Lp0iVTU9MfXy+iqanp2rVrnZ2dv79zAP9HMOk2gIHy8/M7ffq0Vqu1tbWlI+np6Rs3btRfJ7mvr6+oqKijo0MikcybN49eQKCnp6epqYkQUl9fT8+vLZFIRowYQR9CUVRxcXFDQ4OTkxOzsh2jqqqqvLycz+crlcpBk27TB2q1WplMJpfL/8B5A/x1FAAYmNzcXELI7du3bW1t4+Pj6eDTp08JIaWlpXw+PyQkhNk5MzNTLBZzuVxra2tCyPTp02traymKysjIGPRh12g0FEWJxeLg4GBPT08ej0cvwxQUFMS09vHjR39/f0KImZnZyJEjeTxebGwsk+3t7V26dCkhRCwW8/l8Hx8flUrl7u7+l14XgD8DX40CGCgjIyO1Wp2enk5RFCEkLS3N2dnZzc1Nd58nT574+vouX768vb29sbGxurr6w4cParWaEOLn55eTk0MIKS8vpz/tzAAuNTV14cKFPT09PT09UVFRqamppaWldCouLi4jIyMxMbGzs7Orqys4OJheBpbOxsTE3L1798aNG62trZ2dncbGxgUFBX/tBQH4U4a0DAPAV9AjwpycnKqqKg6H8/Dhw0+fPllbW584cYKiKN0RYWBgoI2NzYcPH5hj6aL16tUriqIGFUKaWCxWKBTMw+7ubg6Hk5CQQD+0srJSqVRMtr+/XyqVenp6UhT1+fNnc3NztVrNZFtbW01MTDAihOEO1wgBDJeDg4NSqTx37lxbW1tbW5vuMuW0kpISoVCYkJDARFpaWgghVVVVEyZM+Fazs2bNYrYFAsGYMWMaGxsJIR0dHc3NzQsXLmSyPB7P3d2dvi+0qampq6trwYIFTFYsFk+dOvUXzxFgyKEQAhi0TZs2RUVFabValUpFXwXU1d3dbWxsXFZWphtct26dQCD4TpuDViE3MjL6/PkzIaS9vZ0QMnr0aN3s6NGj+/r6CCEdHR2EEKFQOCjb39//b08KwKCgEAIYtPXr14eHh+fk5Hz1V/A2NjbGxsZXr1796rEcDudfPZetrS2Xy62trdUNvn79mr5xVCqVEkK0Wq1utra2dtBtpQDDDm6WATBoQqHw7NmzR44cWblypX529erVpaWlDx8+1A3SwztCiEgkIl9Gcj+Cz+e7urpevHjx48ePdKSuri4/P3/JkiWEEIFAMG3atIyMjIGBATpbXFz84sWLnzotAAOCESGAoduwYcO3Urt27bp58+by5cvDw8Pd3Nx6e3ufPXt24cIFelQ3ZcoUCwuLyMhIHx+fkSNHBgQE/M/R2+HDh1UqlZeXV2hoaF9fX3x8vKmp6YEDB+hsfHy8j4/P6tWrt2zZ0tzcHBsbSw8TAYY1jAgBDI5IJFq8eDE9ntO3aNGiyZMn09sCgaCwsDAsLOzq1at+fn67d+9+8OBBWFgYnTU1Nc3Pz1coFCUlJQUFBe/evSOEeHh4ODg46Dbo4eFhb29Pb3t5eeXn5/N4vODg4D179sycObO4uHjixIl0dtWqVVeuXHn9+nVgYGBycnJSUlJAQMCgX3QADDsciqKGug8AAABDBiNCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgNRRCAABgtf8AgYuaZdKb9vgAAAAASUVORK5CYII=", - "image/svg+xml": [ - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" - ], - "text/html": [ - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "import Plots\n", - "p1 = Plots.bar([\"INS\", \"SciML\", \"CNODE\"], [time_ins, time_ode, time_node],\n", - " xlabel = \"Method\", ylabel = \"Time (s)\", title = \"Time comparison\", legend = false);\n", - "#Memory allocation\n", - "p2 = Plots.bar([\"INS\", \"SciML\", \"CNODE\"],\n", - " [memory_counters.allocd, memory_counters_ode.allocd, memory_counters_node.allocd],\n", - " xlabel = \"Method\", ylabel = \"Memory (bytes)\",\n", - " title = \"Memory comparison\", legend = false);\n", - "#Garbage collections\n", - "p3 = Plots.bar([\"INS\", \"SciML\", \"CNODE\"], [gc, gc_ode, gc_node], xlabel = \"Method\",\n", - " ylabel = \"Number of GC\", title = \"GC comparison\", legend = false);\n", - "\n", - "Plots.plot(p1, p2, p3, layout = (3, 1), size = (600, 800))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Plots: final state of $u$" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using Plots\n", - "p1 = Plots.heatmap(title = \"\\$u\\$ SciML ODE\", sol_ode.u[end][:, :, 1], ticks = false);\n", - "p2 = Plots.heatmap(title = \"\\$u\\$ SciML CNODE\", sol_node.u[end][:, :, 1], ticks = false);\n", - "p3 = Plots.heatmap(title = \"\\$u\\$ INS\", state.u[1], ticks = false);\n", - "p4 = Plots.heatmap(title = \"\\$u_{INS}-u_{ODE}\\$\",\n", - " state.u[1] - sol_ode.u[end][:, :, 1], ticks = false);\n", - "p5 = Plots.heatmap(title = \"\\$u_{INS}-u_{CNODE}\\$\",\n", - " state.u[1] - sol_node.u[end][:, :, 1], ticks = false);\n", - "p6 = Plots.heatmap(title = \"\\$u_{CNODE}-u_{ODE}\\$\",\n", - " sol_node.u[end][:, :, 1] - sol_ode.u[end][:, :, 1], ticks = false);\n", - "Plots.plot(p1, p2, p3, p4, p5, p6, layout = (2, 3), size = (900, 600), ticks = false)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Plots: vorticity" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "vins = INS.vorticity((state.u[1], state.u[2]), setup)\n", - "vode = INS.vorticity((sol_ode.u[end][:, :, 1], sol_ode.u[end][:, :, 2]), setup)\n", - "vnode = INS.vorticity((sol_node.u[end][:, :, 1], sol_node.u[end][:, :, 2]), setup)\n", - "vor_lims = (-5, 5)\n", - "diff_lims = (-0.4, 0.4)\n", - "\n", - "p1 = Plots.heatmap(title = \"\\$\\\\omega\\$ SciML ODE\", vode,\n", - " color = :viridis, ticks = false, clims = vor_lims);\n", - "p2 = Plots.heatmap(title = \"\\$\\\\omega\\$ in SciML CNODE\", vnode,\n", - " color = :viridis, ticks = false, clims = vor_lims);\n", - "p3 = Plots.heatmap(title = \"vorticity \\$(\\\\omega)\\$ in INS\", vins,\n", - " color = :viridis, ticks = false, clims = vor_lims);\n", - "p4 = Plots.heatmap(title = \"\\$\\\\omega_{INS}-\\\\omega_{ODE}\\$\",\n", - " vins - vode, clim = diff_lims, ticks = false);\n", - "p5 = Plots.heatmap(title = \"\\$\\\\omega_{INS}-\\\\omega_{CNODE}\\$\",\n", - " vins - vnode, clim = diff_lims, ticks = false);\n", - "p6 = Plots.heatmap(title = \"\\$\\\\omega_{CNODE}-\\\\omega_{ODE}\\$\",\n", - " vode - vnode, ticks = false, clim = (0, 0.2));\n", - "Plots.plot(p1, p2, p3, p4, p5, p6, layout = (2, 3), size = (900, 600))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Divergence:" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "INS" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#div_INS = fill!(similar(setup.grid.x[1], setup.grid.N), 0)\n", - "#INS.divergence!(div_INS, state.u, setup)\n", - "div_INS = INS.divergence(state.u, setup)\n", - "maximum(abs.(div_INS))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "SciML" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "u_last_ode = (sol_ode.u[end][:, :, 1], sol_ode.u[end][:, :, 2]);\n", - "div_ode = INS.divergence(u_last_ode, setup)\n", - "max_div_ode = maximum(abs.(div_ode))\n", - "\n", - "using Printf\n", - "anim = Animation()\n", - "for (idx, (t, u)) in enumerate(zip(sol_ode.t, sol_ode.u))\n", - " ∇_u = INS.divergence((u[:, :, 1], u[:, :, 2]), setup)\n", - " title = @sprintf(\"\\$\\\\nabla \\\\cdot u\\$ SciML, t = %.3f s\", t)\n", - " fig = Plots.heatmap(∇_u'; xlabel = \"x\", ylabel = \"y\", title, aspect_ratio = :equal,\n", - " ticks = false, size = (600, 600), clims = (-max_div_ode, max_div_ode))\n", - " frame(anim, fig)\n", - "end\n", - "gif(anim, \"simulations/NavierStokes_2D/plots/divergence_SciML.gif\", fps = 15)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "CNODE" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "u_last_node = (sol_node.u[end][:, :, 1], sol_node.u[end][:, :, 2]);\n", - "div_node = INS.divergence(u_last_node, setup)\n", - "maximum(abs.(div_node))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**Conclusion:** While IncompressibleNavierStokes.jl guarantees a $\\nabla \\cdot u =0$ the other methods do not." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Animations\n", - "#### Animate SciML solution using `Makie`\n", - "! we can plot either the vorticity or the velocity field, however notice that the vorticity is 'flashing' (unstable)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "let\n", - " (; Iu) = setup.grid\n", - " i = 1\n", - " #obs = Observable(sol.u[1][Iu[i], i])\n", - " obs = Observable(INS.vorticity((sol_ode.u[1][:, :, 1], sol_ode.u[1][:, :, 2]), setup))\n", - " fig = GLMakie.heatmap(obs, colorrange = vor_lims)\n", - " fig |> display\n", - " for u in sol_ode.u\n", - " #obs[] = u[Iu[i], i]\n", - " obs[] = INS.vorticity((u[:, :, 1], u[:, :, 2]), setup)\n", - " #fig |> display\n", - " sleep(0.05)\n", - " end\n", - "end" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Animate SciML solution using `Plots.jl`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "function animation_plots(; variable = \"vorticity\")\n", - " anim = Animation()\n", - " for (idx, (t, u)) in enumerate(zip(sol_ode.t, sol_ode.u))\n", - " if variable == \"vorticity\"\n", - " ω = INS.vorticity((u[:, :, 1], u[:, :, 2]), setup)\n", - " title = @sprintf(\"Vorticity SciML, t = %.3f s\", t)\n", - " fig = Plots.heatmap(ω'; xlabel = \"x\", ylabel = \"y\", title, clims = vor_lims,\n", - " color = :viridis, aspect_ratio = :equal, ticks = false, size = (600, 600))\n", - " else\n", - " title = @sprintf(\"\\$u\\$ SciML, t = %.3f s\", t)\n", - " fig = Plots.heatmap(u[:, :, 1]; xlabel = \"x\", ylabel = \"y\", title,\n", - " aspect_ratio = :equal, ticks = false, size = (600, 600))\n", - " end\n", - " frame(anim, fig)\n", - " end\n", - " if variable == \"vorticity\"\n", - " gif(anim, \"simulations/NavierStokes_2D/plots/vorticity_SciML.gif\", fps = 15)\n", - " else\n", - " gif(anim, \"simulations/NavierStokes_2D/plots/velocity_SciML.gif\", fps = 15)\n", - " end\n", - "end" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Vorticity animation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "animation_plots(; variable = \"vorticity\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Velocity animation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "animation_plots(; variable = \"velocity\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Animate the difference in solution using `Plots.jl`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "anim = Animation()\n", - "for idx in 1:Int(ceil(trange[end] / saveat))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "the ODESolution saves the initial state too" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - " error_u = abs.(outputs.field[idx].u[1] - sol_ode.u[idx + 1][:, :, 1])\n", - " title = @sprintf(\"\\$|u_{INS}-u_{ODE}|\\$, t = %.3f s\", sol_ode.t[idx])\n", - " fig = Plots.heatmap(error_u; xlabel = \"x\", ylabel = \"y\", title,\n", - " aspect_ratio = :equal, ticks = false, size = (600, 600))\n", - " frame(anim, fig)\n", - "end\n", - "gif(anim, \"simulations/NavierStokes_2D/plots/u_INS-SciML.gif\", fps = 15)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Different initial conditions" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "t_range_bench = (T(0), T(1))\n", - "function bench_INS(ustart)\n", - " INS.solve_unsteady(; setup, ustart, tlims = t_range_bench, Δt = dt)\n", - "end\n", - "\n", - "function bench_ode(u0)\n", - " prob = ODEProblem{true}(F_ip, stack(u0), t_range_bench)\n", - " solve(prob, RK4(); dt = dt, saveat = saveat)\n", - "end\n", - "\n", - "using Statistics\n", - "total_diff, avg_diff, rel_error = [], [], [];\n", - "samples = 1:100;\n", - "for _ in samples\n", - " ustart = INS.random_field(setup, T(0))\n", - " state, _ = bench_INS(ustart)\n", - " sol_ode = bench_ode(ustart)\n", - " push!(total_diff, sum(abs.(state.u[1] - sol_ode.u[end][:, :, 1])))\n", - " push!(avg_diff, mean(abs.(state.u[1] - sol_ode.u[end][:, :, 1])))\n", - " #push!(rel_error, mean(abs.((state.u[1] - sol_ode.u[end][:, :, 1])/state.u[1]))) # zero division\n", - "end\n", - "\n", - "Plots.histogram(total_diff, bins = 20, label = \"Total error\",\n", - " xlabel = \"\\$\\\\sum(|u_{INS}-u_{ODE}|)\\$\")\n", - "Plots.histogram(avg_diff, bins = 20, label = \"Mean absolute error\",\n", - " xlabel = \"\\$\\\\langle |u_{INS}-u_{ODE}|\\\\rangle \\$\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "---\n", - "\n", - "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 1.10.3", - "language": "julia", - "name": "julia-1.10" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "1.10.3" - } - }, - "nbformat": 4, - "nbformat_minor": 3 -} diff --git a/simulations/NavierStokes_2D/plots/hist_mean_abs_error.svg b/simulations/NavierStokes_2D/plots/hist_mean_abs_error.svg deleted file mode 100644 index b354296..0000000 --- a/simulations/NavierStokes_2D/plots/hist_mean_abs_error.svg +++ /dev/null @@ -1,88 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ No newline at end of file diff --git a/simulations/NavierStokes_2D/plots/hist_total_error.svg b/simulations/NavierStokes_2D/plots/hist_total_error.svg deleted file mode 100644 index 795fd5d..0000000 --- a/simulations/NavierStokes_2D/plots/hist_total_error.svg +++ /dev/null @@ -1,101 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ No newline at end of file diff --git a/simulations/NavierStokes_2D/plots/velocity_SciML.gif b/simulations/NavierStokes_2D/plots/velocity_SciML.gif deleted file mode 100644 index 6f25c09..0000000 Binary files a/simulations/NavierStokes_2D/plots/velocity_SciML.gif and /dev/null differ diff --git a/simulations/NavierStokes_2D/plots/vorticity_INS.mkv b/simulations/NavierStokes_2D/plots/vorticity_INS.mkv deleted file mode 100644 index 78c96ce..0000000 Binary files a/simulations/NavierStokes_2D/plots/vorticity_INS.mkv and /dev/null differ diff --git a/simulations/NavierStokes_2D/plots/vorticity_SciML.gif b/simulations/NavierStokes_2D/plots/vorticity_SciML.gif deleted file mode 100644 index 0fa4a72..0000000 Binary files a/simulations/NavierStokes_2D/plots/vorticity_SciML.gif and /dev/null differ