From 62f05bcd04abafb3b4510dc9b28da8a53f429d69 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Fri, 3 Jul 2020 21:13:34 +0000 Subject: [PATCH] Rebuild model_inference/02-monte_carlo_parameter_estim.jmd --- .../02-monte_carlo_parameter_estim.html | 985 ++++++++++++++++++ .../02-monte_carlo_parameter_estim.md | 346 ++++++ .../02-monte_carlo_parameter_estim.ipynb | 204 ++++ .../02-monte_carlo_parameter_estim.pdf | Bin 0 -> 53826 bytes .../02-monte_carlo_parameter_estim.jl | 74 ++ 5 files changed, 1609 insertions(+) create mode 100644 html/model_inference/02-monte_carlo_parameter_estim.html create mode 100644 markdown/model_inference/02-monte_carlo_parameter_estim.md create mode 100644 notebook/model_inference/02-monte_carlo_parameter_estim.ipynb create mode 100644 pdf/model_inference/02-monte_carlo_parameter_estim.pdf create mode 100644 script/model_inference/02-monte_carlo_parameter_estim.jl diff --git a/html/model_inference/02-monte_carlo_parameter_estim.html b/html/model_inference/02-monte_carlo_parameter_estim.html new file mode 100644 index 00000000..9212b7e4 --- /dev/null +++ b/html/model_inference/02-monte_carlo_parameter_estim.html @@ -0,0 +1,985 @@ + + + + + + Monte Carlo Parameter Estimation From Data + + + + + + + + + + + + + + + +
+
+
+
+

Monte Carlo Parameter Estimation From Data

+
Chris Rackauckas
+ +
+ +

First you want to create a problem which solves multiple problems at the same time. This is the Monte Carlo Problem. When the parameter estimation tools say it will take any DEProblem, it really means ANY DEProblem!

+

So, let's get a Monte Carlo problem setup that solves with 10 different initial conditions.

+ + +
+using DifferentialEquations, DiffEqParamEstim, Plots, Optim
+
+ + +
+ERROR: Failed to precompile DiffEqParamEstim [1130ab10-4a5a-5621-a13d-e4788d82bd4c] to /builds/JuliaGPU/DiffEqTutorials.jl/.julia/compiled/v1.4/DiffEqParamEstim/nWq0E_XzcJh.ji.
+
+ + + +
+# Monte Carlo Problem Set Up for solving set of ODEs with different initial conditions
+
+# Set up Lotka-Volterra system
+function pf_func(du,u,p,t)
+  du[1] = p[1] * u[1] - p[2] * u[1]*u[2]
+  du[2] = -3 * u[2] + u[1]*u[2]
+end
+p = [1.5,1.0]
+prob = ODEProblem(pf_func,[1.0,1.0],(0.0,10.0),p)
+
+ + +
+ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
+timespan: (0.0, 10.0)
+u0: [1.0, 1.0]
+
+ + +

Now for a MonteCarloProblem we have to take this problem and tell it what to do N times via the prob_func. So let's generate N=10 different initial conditions, and tell it to run the same problem but with these 10 different initial conditions each time:

+ + +
+# Setting up to solve the problem N times (for the N different initial conditions)
+N = 10;
+initial_conditions = [[1.0,1.0], [1.0,1.5], [1.5,1.0], [1.5,1.5], [0.5,1.0], [1.0,0.5], [0.5,0.5], [2.0,1.0], [1.0,2.0], [2.0,2.0]]
+function prob_func(prob,i,repeat)
+  ODEProblem(prob.f,initial_conditions[i],prob.tspan,prob.p)
+end
+monte_prob = MonteCarloProblem(prob,prob_func=prob_func)
+
+ + +
+EnsembleProblem with problem ODEProblem
+
+ + +

We can check this does what we want by solving it:

+ + +
+# Check above does what we want
+sim = solve(monte_prob,Tsit5(),num_monte=N)
+plot(sim)
+
+ + +
+ERROR: UndefVarError: plot not defined
+
+ + +

nummonte=N means "run N times", and each time it runs the problem returned by the probfunc, which is always the same problem but with the ith initial condition.

+

Now let's generate a dataset from that. Let's get data points at every t=0.1 using saveat, and then convert the solution into an array.

+ + +
+# Generate a dataset from these runs
+data_times = 0.0:0.1:10.0
+sim = solve(monte_prob,Tsit5(),num_monte=N,saveat=data_times)
+data = Array(sim)
+
+ + +
+2×101×10 Array{Float64,3}:
+[:, :, 1] =
+ 1.0  1.06108   1.14403   1.24917   1.37764   …  0.956979  0.983561  1.0337
+6
+ 1.0  0.821084  0.679053  0.566893  0.478813     1.35559   1.10629   0.9063
+7
+
+[:, :, 2] =
+ 1.0  1.01413  1.05394  1.11711   …  1.05324  1.01309  1.00811  1.03162
+ 1.5  1.22868  1.00919  0.833191     2.08023  1.70818  1.39972  1.14802
+
+[:, :, 3] =
+ 1.5  1.58801   1.70188   1.84193   2.00901   …  2.0153    2.21084   2.4358
+9
+ 1.0  0.864317  0.754624  0.667265  0.599149     0.600943  0.549793  0.5136
+8
+
+[:, :, 4] =
+ 1.5  1.51612  1.5621   1.63555   1.73531   …  1.83822   1.98545   2.15958
+ 1.5  1.29176  1.11592  0.969809  0.850159     0.771088  0.691421  0.630025
+
+[:, :, 5] =
+ 0.5  0.531705  0.576474  0.634384  0.706139  …  9.05366   9.4006   8.8391
+ 1.0  0.77995   0.610654  0.480565  0.380645     0.809383  1.51708  2.82619
+
+[:, :, 6] =
+ 1.0  1.11027   1.24238   1.39866   1.58195   …  0.753107  0.748814  0.7682
+84
+ 0.5  0.411557  0.342883  0.289812  0.249142     1.73879   1.38829   1.1093
+2
+
+[:, :, 7] =
+ 0.5  0.555757  0.623692  0.705084  0.80158   …  8.11213   9.10669   9.9216
+9
+ 0.5  0.390449  0.30679   0.24286   0.193966     0.261294  0.455928  0.8787
+92
+
+[:, :, 8] =
+ 2.0  2.11239   2.24921   2.41003   2.59433   …  3.22292   3.47356   3.7301
+1
+ 1.0  0.909749  0.838025  0.783532  0.745339     0.739406  0.765524  0.8130
+04
+
+[:, :, 9] =
+ 1.0  0.969326  0.971358  1.00017  …  1.25065  1.1012   1.01733  0.979304
+ 2.0  1.63445   1.33389   1.09031     3.02672  2.52063  2.07503  1.69808
+
+[:, :, 10] =
+ 2.0  1.92148  1.88215  1.87711  1.90264  …  2.15079  2.27937   2.43105
+ 2.0  1.80195  1.61405  1.4426   1.2907      0.95722  0.884825  0.829478
+
+ + +

Here, data[i,j,k] is the same as sim[i,j,k] which is the same as sim[k]i,j. So data[i,j,k] is the jth timepoint of the ith variable in the kth trajectory.

+

Now let's build a loss function. A loss function is some loss(sol) that spits out a scalar for how far from optimal we are. In the documentation I show that we normally do loss = L2Loss(t,data), but we can bootstrap off of this. Instead lets build an array of N loss functions, each one with the correct piece of data.

+ + +
+# Building a loss function
+losses = [L2Loss(data_times,data[:,:,i]) for i in 1:N]
+
+ + +
+ERROR: UndefVarError: L2Loss not defined
+
+ + +

So losses[i] is a function which computes the loss of a solution against the data of the ith trajectory. So to build our true loss function, we sum the losses:

+ + +
+loss(sim) = sum(losses[i](sim[i]) for i in 1:N)
+
+ + +
+loss (generic function with 1 method)
+
+ + +

As a double check, make sure that loss(sim) outputs zero (since we generated the data from sim). Now we generate data with other parameters:

+ + +
+prob = ODEProblem(pf_func,[1.0,1.0],(0.0,10.0),[1.2,0.8])
+function prob_func(prob,i,repeat)
+  ODEProblem(prob.f,initial_conditions[i],prob.tspan,prob.p)
+end
+monte_prob = MonteCarloProblem(prob,prob_func=prob_func)
+sim = solve(monte_prob,Tsit5(),num_monte=N,saveat=data_times)
+loss(sim)
+
+ + +
+ERROR: UndefVarError: losses not defined
+
+ + +

and get a non-zero loss. So we now have our problem, our data, and our loss function... we have what we need.

+

Put this into buildlossobjective.

+ + +
+obj = build_loss_objective(monte_prob,Tsit5(),loss,num_monte=N,
+                           saveat=data_times)
+
+ + +
+ERROR: UndefVarError: build_loss_objective not defined
+
+ + +

Notice that I added the kwargs for solve into this. They get passed to an internal solve command, so then the loss is computed on N trajectories at data_times.

+

Thus we take this objective function over to any optimization package. I like to do quick things in Optim.jl. Here, since the Lotka-Volterra equation requires positive parameters, I use Fminbox to make sure the parameters stay positive. I start the optimization with [1.3,0.9], and Optim spits out that the true parameters are:

+ + +
+lower = zeros(2)
+upper = fill(2.0,2)
+result = optimize(obj, lower, upper, [1.3,0.9], Fminbox(BFGS()))
+
+ + +
+ERROR: UndefVarError: BFGS not defined
+
+ + + +
+result
+
+ + +
+ERROR: UndefVarError: result not defined
+
+ + +

Optim finds one but not the other parameter.

+

I would run a test on synthetic data for your problem before using it on real data. Maybe play around with different optimization packages, or add regularization. You may also want to decrease the tolerance of the ODE solvers via

+ + +
+obj = build_loss_objective(monte_prob,Tsit5(),loss,num_monte=N,
+                           abstol=1e-8,reltol=1e-8,
+                           saveat=data_times)
+
+ + +
+ERROR: UndefVarError: build_loss_objective not defined
+
+ + + +
+result = optimize(obj, lower, upper, [1.3,0.9], Fminbox(BFGS()))
+
+ + +
+ERROR: UndefVarError: BFGS not defined
+
+ + + +
+result
+
+ + +
+ERROR: UndefVarError: result not defined
+
+ + +

if you suspect error is the problem. However, if you're having problems it's most likely not the ODE solver tolerance and mostly because parameter inference is a very hard optimization problem.

+ + +

Appendix

+

This tutorial is part of the DiffEqTutorials.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqTutorials.jl

+
+

To locally run this tutorial, do the following commands:

+
using DiffEqTutorials
+DiffEqTutorials.weave_file("model_inference","02-monte_carlo_parameter_estim.jmd")
+
+

Computer Information:

+
+
Julia Version 1.4.2
+Commit 44fa15b150* (2020-05-23 18:35 UTC)
+Platform Info:
+  OS: Linux (x86_64-pc-linux-gnu)
+  CPU: Intel(R) Core(TM) i7-9700K CPU @ 3.60GHz
+  WORD_SIZE: 64
+  LIBM: libopenlibm
+  LLVM: libLLVM-8.0.1 (ORCJIT, skylake)
+Environment:
+  JULIA_DEPOT_PATH = /builds/JuliaGPU/DiffEqTutorials.jl/.julia
+  JULIA_CUDA_MEMORY_LIMIT = 2147483648
+  JULIA_PROJECT = @.
+  JULIA_NUM_THREADS = 4
+
+
+

Package Information:

+
+
Status `/builds/JuliaGPU/DiffEqTutorials.jl/tutorials/model_inference/Project.toml`
+[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.5.0
+[593b3428-ca2f-500c-ae53-031589ec8ddd] CmdStan 6.0.6
+[ebbdde9d-f333-5424-9be2-dbf1e9acfb5e] DiffEqBayes 2.16.0
+[1130ab10-4a5a-5621-a13d-e4788d82bd4c] DiffEqParamEstim 1.15.0
+[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.15.0
+[31c24e10-a181-5473-b8eb-7969acd0382f] Distributions 0.23.4
+[bbc10e6e-7c05-544b-b16e-64fede858acb] DynamicHMC 2.1.0
+[429524aa-4258-5aef-a3af-852621145aeb] Optim 0.22.0
+[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.41.0
+[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 1.5.2
+[731186ca-8d62-57ce-b412-fbd966d074cd] RecursiveArrayTools 2.5.0
+[f3b207a7-027a-5e70-b257-86293d7955fd] StatsPlots 0.14.6
+[84d833dd-6860-57f9-a1a7-6da5db126cff] TransformVariables 0.3.9
+
+ + +
+
+

+ Published from 02-monte_carlo_parameter_estim.jmd + using Weave.jl v0.10.2 on 2020-07-03. +

+
+
+
+
+ + + diff --git a/markdown/model_inference/02-monte_carlo_parameter_estim.md b/markdown/model_inference/02-monte_carlo_parameter_estim.md new file mode 100644 index 00000000..f5dba308 --- /dev/null +++ b/markdown/model_inference/02-monte_carlo_parameter_estim.md @@ -0,0 +1,346 @@ +--- +author: "Chris Rackauckas" +title: "Monte Carlo Parameter Estimation From Data" +--- + + +First you want to create a problem which solves multiple problems at the same time. This is the Monte Carlo Problem. When the parameter estimation tools say it will take any DEProblem, it really means ANY DEProblem! + +So, let's get a Monte Carlo problem setup that solves with 10 different initial conditions. + +````julia +using DifferentialEquations, DiffEqParamEstim, Plots, Optim +```` + + +```` +Error: Failed to precompile DiffEqParamEstim [1130ab10-4a5a-5621-a13d-e4788 +d82bd4c] to /builds/JuliaGPU/DiffEqTutorials.jl/.julia/compiled/v1.4/DiffEq +ParamEstim/nWq0E_XzcJh.ji. +```` + + + +````julia + +# Monte Carlo Problem Set Up for solving set of ODEs with different initial conditions + +# Set up Lotka-Volterra system +function pf_func(du,u,p,t) + du[1] = p[1] * u[1] - p[2] * u[1]*u[2] + du[2] = -3 * u[2] + u[1]*u[2] +end +p = [1.5,1.0] +prob = ODEProblem(pf_func,[1.0,1.0],(0.0,10.0),p) +```` + + +```` +ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true +timespan: (0.0, 10.0) +u0: [1.0, 1.0] +```` + + + + + +Now for a MonteCarloProblem we have to take this problem and tell it what to do N times via the prob_func. So let's generate N=10 different initial conditions, and tell it to run the same problem but with these 10 different initial conditions each time: + +````julia +# Setting up to solve the problem N times (for the N different initial conditions) +N = 10; +initial_conditions = [[1.0,1.0], [1.0,1.5], [1.5,1.0], [1.5,1.5], [0.5,1.0], [1.0,0.5], [0.5,0.5], [2.0,1.0], [1.0,2.0], [2.0,2.0]] +function prob_func(prob,i,repeat) + ODEProblem(prob.f,initial_conditions[i],prob.tspan,prob.p) +end +monte_prob = MonteCarloProblem(prob,prob_func=prob_func) +```` + + +```` +EnsembleProblem with problem ODEProblem +```` + + + + + +We can check this does what we want by solving it: + +````julia +# Check above does what we want +sim = solve(monte_prob,Tsit5(),num_monte=N) +plot(sim) +```` + + +```` +Error: UndefVarError: plot not defined +```` + + + + + +num_monte=N means "run N times", and each time it runs the problem returned by the prob_func, which is always the same problem but with the ith initial condition. + +Now let's generate a dataset from that. Let's get data points at every t=0.1 using saveat, and then convert the solution into an array. + +````julia +# Generate a dataset from these runs +data_times = 0.0:0.1:10.0 +sim = solve(monte_prob,Tsit5(),num_monte=N,saveat=data_times) +data = Array(sim) +```` + + +```` +2×101×10 Array{Float64,3}: +[:, :, 1] = + 1.0 1.06108 1.14403 1.24917 1.37764 … 0.956979 0.983561 1.0337 +6 + 1.0 0.821084 0.679053 0.566893 0.478813 1.35559 1.10629 0.9063 +7 + +[:, :, 2] = + 1.0 1.01413 1.05394 1.11711 … 1.05324 1.01309 1.00811 1.03162 + 1.5 1.22868 1.00919 0.833191 2.08023 1.70818 1.39972 1.14802 + +[:, :, 3] = + 1.5 1.58801 1.70188 1.84193 2.00901 … 2.0153 2.21084 2.4358 +9 + 1.0 0.864317 0.754624 0.667265 0.599149 0.600943 0.549793 0.5136 +8 + +[:, :, 4] = + 1.5 1.51612 1.5621 1.63555 1.73531 … 1.83822 1.98545 2.15958 + 1.5 1.29176 1.11592 0.969809 0.850159 0.771088 0.691421 0.630025 + +[:, :, 5] = + 0.5 0.531705 0.576474 0.634384 0.706139 … 9.05366 9.4006 8.8391 + 1.0 0.77995 0.610654 0.480565 0.380645 0.809383 1.51708 2.82619 + +[:, :, 6] = + 1.0 1.11027 1.24238 1.39866 1.58195 … 0.753107 0.748814 0.7682 +84 + 0.5 0.411557 0.342883 0.289812 0.249142 1.73879 1.38829 1.1093 +2 + +[:, :, 7] = + 0.5 0.555757 0.623692 0.705084 0.80158 … 8.11213 9.10669 9.9216 +9 + 0.5 0.390449 0.30679 0.24286 0.193966 0.261294 0.455928 0.8787 +92 + +[:, :, 8] = + 2.0 2.11239 2.24921 2.41003 2.59433 … 3.22292 3.47356 3.7301 +1 + 1.0 0.909749 0.838025 0.783532 0.745339 0.739406 0.765524 0.8130 +04 + +[:, :, 9] = + 1.0 0.969326 0.971358 1.00017 … 1.25065 1.1012 1.01733 0.979304 + 2.0 1.63445 1.33389 1.09031 3.02672 2.52063 2.07503 1.69808 + +[:, :, 10] = + 2.0 1.92148 1.88215 1.87711 1.90264 … 2.15079 2.27937 2.43105 + 2.0 1.80195 1.61405 1.4426 1.2907 0.95722 0.884825 0.829478 +```` + + + + + +Here, data[i,j,k] is the same as sim[i,j,k] which is the same as sim[k][i,j] (where sim[k] is the kth solution). So data[i,j,k] is the jth timepoint of the ith variable in the kth trajectory. + +Now let's build a loss function. A loss function is some loss(sol) that spits out a scalar for how far from optimal we are. In the documentation I show that we normally do loss = L2Loss(t,data), but we can bootstrap off of this. Instead lets build an array of N loss functions, each one with the correct piece of data. + +````julia +# Building a loss function +losses = [L2Loss(data_times,data[:,:,i]) for i in 1:N] +```` + + +```` +Error: UndefVarError: L2Loss not defined +```` + + + + + +So losses[i] is a function which computes the loss of a solution against the data of the ith trajectory. So to build our true loss function, we sum the losses: + +````julia +loss(sim) = sum(losses[i](sim[i]) for i in 1:N) +```` + + +```` +loss (generic function with 1 method) +```` + + + + + +As a double check, make sure that loss(sim) outputs zero (since we generated the data from sim). Now we generate data with other parameters: + +````julia +prob = ODEProblem(pf_func,[1.0,1.0],(0.0,10.0),[1.2,0.8]) +function prob_func(prob,i,repeat) + ODEProblem(prob.f,initial_conditions[i],prob.tspan,prob.p) +end +monte_prob = MonteCarloProblem(prob,prob_func=prob_func) +sim = solve(monte_prob,Tsit5(),num_monte=N,saveat=data_times) +loss(sim) +```` + + +```` +Error: UndefVarError: losses not defined +```` + + + + + +and get a non-zero loss. So we now have our problem, our data, and our loss function... we have what we need. + +Put this into build_loss_objective. + +````julia +obj = build_loss_objective(monte_prob,Tsit5(),loss,num_monte=N, + saveat=data_times) +```` + + +```` +Error: UndefVarError: build_loss_objective not defined +```` + + + + + +Notice that I added the kwargs for solve into this. They get passed to an internal solve command, so then the loss is computed on N trajectories at data_times. + +Thus we take this objective function over to any optimization package. I like to do quick things in Optim.jl. Here, since the Lotka-Volterra equation requires positive parameters, I use Fminbox to make sure the parameters stay positive. I start the optimization with [1.3,0.9], and Optim spits out that the true parameters are: + +````julia +lower = zeros(2) +upper = fill(2.0,2) +result = optimize(obj, lower, upper, [1.3,0.9], Fminbox(BFGS())) +```` + + +```` +Error: UndefVarError: BFGS not defined +```` + + + +````julia +result +```` + + +```` +Error: UndefVarError: result not defined +```` + + + + + +Optim finds one but not the other parameter. + +I would run a test on synthetic data for your problem before using it on real data. Maybe play around with different optimization packages, or add regularization. You may also want to decrease the tolerance of the ODE solvers via + +````julia +obj = build_loss_objective(monte_prob,Tsit5(),loss,num_monte=N, + abstol=1e-8,reltol=1e-8, + saveat=data_times) +```` + + +```` +Error: UndefVarError: build_loss_objective not defined +```` + + + +````julia +result = optimize(obj, lower, upper, [1.3,0.9], Fminbox(BFGS())) +```` + + +```` +Error: UndefVarError: BFGS not defined +```` + + + +````julia +result +```` + + +```` +Error: UndefVarError: result not defined +```` + + + + + +if you suspect error is the problem. However, if you're having problems it's most likely not the ODE solver tolerance and mostly because parameter inference is a very hard optimization problem. + + +## Appendix + + This tutorial is part of the DiffEqTutorials.jl repository, found at: + +To locally run this tutorial, do the following commands: +``` +using DiffEqTutorials +DiffEqTutorials.weave_file("model_inference","02-monte_carlo_parameter_estim.jmd") +``` + +Computer Information: +``` +Julia Version 1.4.2 +Commit 44fa15b150* (2020-05-23 18:35 UTC) +Platform Info: + OS: Linux (x86_64-pc-linux-gnu) + CPU: Intel(R) Core(TM) i7-9700K CPU @ 3.60GHz + WORD_SIZE: 64 + LIBM: libopenlibm + LLVM: libLLVM-8.0.1 (ORCJIT, skylake) +Environment: + JULIA_DEPOT_PATH = /builds/JuliaGPU/DiffEqTutorials.jl/.julia + JULIA_CUDA_MEMORY_LIMIT = 2147483648 + JULIA_PROJECT = @. + JULIA_NUM_THREADS = 4 + +``` + +Package Information: + +``` +Status `/builds/JuliaGPU/DiffEqTutorials.jl/tutorials/model_inference/Project.toml` +[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.5.0 +[593b3428-ca2f-500c-ae53-031589ec8ddd] CmdStan 6.0.6 +[ebbdde9d-f333-5424-9be2-dbf1e9acfb5e] DiffEqBayes 2.16.0 +[1130ab10-4a5a-5621-a13d-e4788d82bd4c] DiffEqParamEstim 1.15.0 +[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.15.0 +[31c24e10-a181-5473-b8eb-7969acd0382f] Distributions 0.23.4 +[bbc10e6e-7c05-544b-b16e-64fede858acb] DynamicHMC 2.1.0 +[429524aa-4258-5aef-a3af-852621145aeb] Optim 0.22.0 +[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.41.0 +[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 1.5.2 +[731186ca-8d62-57ce-b412-fbd966d074cd] RecursiveArrayTools 2.5.0 +[f3b207a7-027a-5e70-b257-86293d7955fd] StatsPlots 0.14.6 +[84d833dd-6860-57f9-a1a7-6da5db126cff] TransformVariables 0.3.9 +``` diff --git a/notebook/model_inference/02-monte_carlo_parameter_estim.ipynb b/notebook/model_inference/02-monte_carlo_parameter_estim.ipynb new file mode 100644 index 00000000..de51837b --- /dev/null +++ b/notebook/model_inference/02-monte_carlo_parameter_estim.ipynb @@ -0,0 +1,204 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "First you want to create a problem which solves multiple problems at the same time. This is the Monte Carlo Problem. When the parameter estimation tools say it will take any DEProblem, it really means ANY DEProblem!\n\nSo, let's get a Monte Carlo problem setup that solves with 10 different initial conditions." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "using DifferentialEquations, DiffEqParamEstim, Plots, Optim\n\n# Monte Carlo Problem Set Up for solving set of ODEs with different initial conditions\n\n# Set up Lotka-Volterra system\nfunction pf_func(du,u,p,t)\n du[1] = p[1] * u[1] - p[2] * u[1]*u[2]\n du[2] = -3 * u[2] + u[1]*u[2]\nend\np = [1.5,1.0]\nprob = ODEProblem(pf_func,[1.0,1.0],(0.0,10.0),p)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "Now for a MonteCarloProblem we have to take this problem and tell it what to do N times via the prob_func. So let's generate N=10 different initial conditions, and tell it to run the same problem but with these 10 different initial conditions each time:" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "# Setting up to solve the problem N times (for the N different initial conditions)\nN = 10;\ninitial_conditions = [[1.0,1.0], [1.0,1.5], [1.5,1.0], [1.5,1.5], [0.5,1.0], [1.0,0.5], [0.5,0.5], [2.0,1.0], [1.0,2.0], [2.0,2.0]]\nfunction prob_func(prob,i,repeat)\n ODEProblem(prob.f,initial_conditions[i],prob.tspan,prob.p)\nend\nmonte_prob = MonteCarloProblem(prob,prob_func=prob_func)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "We can check this does what we want by solving it:" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "# Check above does what we want\nsim = solve(monte_prob,Tsit5(),num_monte=N)\nplot(sim)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "num_monte=N means \"run N times\", and each time it runs the problem returned by the prob_func, which is always the same problem but with the ith initial condition.\n\nNow let's generate a dataset from that. Let's get data points at every t=0.1 using saveat, and then convert the solution into an array." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "# Generate a dataset from these runs\ndata_times = 0.0:0.1:10.0\nsim = solve(monte_prob,Tsit5(),num_monte=N,saveat=data_times)\ndata = Array(sim)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "Here, data[i,j,k] is the same as sim[i,j,k] which is the same as sim[k][i,j] (where sim[k] is the kth solution). So data[i,j,k] is the jth timepoint of the ith variable in the kth trajectory.\n\nNow let's build a loss function. A loss function is some loss(sol) that spits out a scalar for how far from optimal we are. In the documentation I show that we normally do loss = L2Loss(t,data), but we can bootstrap off of this. Instead lets build an array of N loss functions, each one with the correct piece of data." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "# Building a loss function\nlosses = [L2Loss(data_times,data[:,:,i]) for i in 1:N]" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "So losses[i] is a function which computes the loss of a solution against the data of the ith trajectory. So to build our true loss function, we sum the losses:" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "loss(sim) = sum(losses[i](sim[i]) for i in 1:N)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "As a double check, make sure that loss(sim) outputs zero (since we generated the data from sim). Now we generate data with other parameters:" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "prob = ODEProblem(pf_func,[1.0,1.0],(0.0,10.0),[1.2,0.8])\nfunction prob_func(prob,i,repeat)\n ODEProblem(prob.f,initial_conditions[i],prob.tspan,prob.p)\nend\nmonte_prob = MonteCarloProblem(prob,prob_func=prob_func)\nsim = solve(monte_prob,Tsit5(),num_monte=N,saveat=data_times)\nloss(sim)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "and get a non-zero loss. So we now have our problem, our data, and our loss function... we have what we need.\n\nPut this into build_loss_objective." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "obj = build_loss_objective(monte_prob,Tsit5(),loss,num_monte=N,\n saveat=data_times)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "Notice that I added the kwargs for solve into this. They get passed to an internal solve command, so then the loss is computed on N trajectories at data_times.\n\nThus we take this objective function over to any optimization package. I like to do quick things in Optim.jl. Here, since the Lotka-Volterra equation requires positive parameters, I use Fminbox to make sure the parameters stay positive. I start the optimization with [1.3,0.9], and Optim spits out that the true parameters are:" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "lower = zeros(2)\nupper = fill(2.0,2)\nresult = optimize(obj, lower, upper, [1.3,0.9], Fminbox(BFGS()))" + ], + "metadata": {}, + "execution_count": null + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "result" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "Optim finds one but not the other parameter.\n\nI would run a test on synthetic data for your problem before using it on real data. Maybe play around with different optimization packages, or add regularization. You may also want to decrease the tolerance of the ODE solvers via" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "obj = build_loss_objective(monte_prob,Tsit5(),loss,num_monte=N,\n abstol=1e-8,reltol=1e-8,\n saveat=data_times)\nresult = optimize(obj, lower, upper, [1.3,0.9], Fminbox(BFGS()))" + ], + "metadata": {}, + "execution_count": null + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "result" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "if you suspect error is the problem. However, if you're having problems it's most likely not the ODE solver tolerance and mostly because parameter inference is a very hard optimization problem." + ], + "metadata": {} + } + ], + "nbformat_minor": 2, + "metadata": { + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.4.2" + }, + "kernelspec": { + "name": "julia-1.4", + "display_name": "Julia 1.4.2", + "language": "julia" + } + }, + "nbformat": 4 +} diff --git a/pdf/model_inference/02-monte_carlo_parameter_estim.pdf b/pdf/model_inference/02-monte_carlo_parameter_estim.pdf new file mode 100644 index 0000000000000000000000000000000000000000..176cb2f0a568b4ed1e156583681a73d8717ab2e2 GIT binary patch literal 53826 zcmbTdQ?Mvq+O4^4+qP}nwr%fa+qP}nwr%cZ+wSjlcXU_9U#Fv@FES%<^6Hs6<`^U2 zNva?sM$1UY3PpN!cl`#%$U?wCU~gmv#lu4{W@+PW>O?PQW9V!uVrpz}VoEP#YG>|j zLBPz&&ceqB^~c%C)X)~neRD)-D)zD)cK4O?3#dqH85d)^P(VAs3$2N|Gd&~Qc}4X1 zRX=%oB{6%rA?1f-d}0_q@k1gZS;Fi6mIy8Oci70VVU+naE$^4c??dP34ZQ)teAHi$ zuP*9HMyjdX(m-qI!Q=gu`}88}5AEIO10a|m_FPAB9Qy?54-gNv?;XV$HIA;NIl!F7 zJK#r8Zi)S~?5i1;GtW}O_&_l;)IrMr>=U%5tpZmJG zzI_2h#7Id~e=sT(^I2MZC=96_Fe%}LH5L@45HQTp{JI1RZ*hoZNg2BPbY1+Qeg%!u z=La%WWi=BL)9w(5Mlxh}eMP1G1<_!~tJTw6sy+k0^5f2KNc(iBa;oDD_KW})NvGTa zkwmE;aMu8q08_eSis2$b)-;Y-Pk2H6WlZ0Fj-YgrES0nrWEn!80EI(+KzdbGU?w;X zAtmMZL80Ny;+gVE>-}t0Ie~j;HfxcIR{=_)l)p4)g-mTSm#s7zdW=eWej$#!Ybr3+ zuPkY{Rh5#9>sp5cW$=GKl#Ho3Ag7D{-eg4_Dm&FVG&flPv*qu{Y};cFC!s(1t# z_Y0B-W+2kl3vm`TJg{1gHdQo5JT43OnmQk%J&bzY$$3veQUNa9s8Z!jimsa;aq?-7 zJM#nc1MoD3P7q4~3k062FMYS@YdlALl>b_bBqgm#;S_H{)7rz#lWJpqBI$s(`glAY z%y50nTSoh3_kADJ4-@YD<{mZah*t(j=IJp8C^s?1cxAKZRGf2~yd5+ZY?Ll6Qgaqy zG6WMcu@n*6e6ln+2hXd1;D%zDy{16P$OUT=<=v5W<$X(zdBdkP%|nBw#e`!q_h9OF zvJW-(Cmf29L>Nb66NSo{a5%^IeV{Oby4qMlU~`!1Zo-#NERWP@Ot?i!%?~Ygf9CXt zLPdwo;X9!=gCo>T@BsR|LZZw?r&~CAxFu4oRld z&}7{m_ZF|S-xh}vFelNbfrVVhnih$>gRa0}%a&rR;!9fho^;a}zz|~oX9CB1GVPj2 zPWOwu1X+ao_r?Q^kKHB5gzcz3F;{2sq0}3`;h7^tb$0G`w$y)vRLnPrY1#pI(jXz! z#Stf%x?`W|XEr9|NpVz+yp=CrHyy7eb&#Q+^t0oI?h?S%K(9UEVQmtBq~0YrsklXW z&XO1=_tb=?iCCviIRKId)Gd05foGCb$(JDjoZ}<~&;2VH#sJ8tB)y=;S~`zRHTy5K zJ4R#g*KW1~8;r^H;E@xYq@0NywmhjX{#;)dlJgk8FFZCM%K5b^w%X?S#9Jsu!hjRC-K12Lk~EU>ub zlhlF#H3V%2xO>7q1PK)`PCcF>53>=GkXpg}I5EWj(_qSL6~cpifQtCl*dfyLwj+OA zVf4tcA|-Y&(IgAZ^pEVET;q~rDIqyKQ%Bi&D#miWDmL|M(A$2YV5Fen$z}02sZ`4u zAKMpFoR-H4hIEAH99GoQS=}^?F{NT9VcGZLx%$BZ_1_A{sub(6LlA+mgEBysgRQ0$ zBrlNS!kOrPt}$t}Xv8Gb6*3d^KKS5&AqTccsj+2Kyo zp43;kS`tCVi*bpMye@Y;)akM9y=_Wy(3hAelChPC-FrZh(w3Am46})u_-VA{3pqA? zwsL2WOhz23I?p)e7A?wRIIxd;eE|IZ5U}}ecV1$A(SAZ30g)BhJZ@~R56)>=Lox1l zi^WdXd&~)wwP_l7CA3`(;!+)%SgIGHhTT9i7hKXK@7pK>qNW)&0g?Hw3h$Hy)~Z5# zC+=}fEleX-Dx0=y-7kJ4c0DTo-lWIO$i&(9L~cS#+=*xn`17H7*bA$RcM!?G zgzX#j&8kE__8AGvH`2)|uP!U_qeUjqHs>lZjW#a+E0=1z9_JsQywdXNmsMKU&{*C>A=zi1t2XDl`B@%&y%@1i!rM@u1s;k*mPAwe zw(Rlus=%5Qu0k#_Et=8f#)(!t(Ic;H#*n}Dm{I%-(WqGm!}$sg6$-l)$?TPfw<`(3 zwMTevN~@u@p=?U?U706N_%gOpTXybMNJmKbw<}QSSfk83 zrza(@0AdWbdPT_fiGcXd55MP3fHN%kW6z6p$i^AVL19}v=C&_(A(o^Wtqd0z<7yhx zHh-j#Ebi&;GFo7i5cafREbQ6fyq+K8?x1PWpSig_KjqE|3GM3Nn-qKnbWk zd*>sOq0|-`EdPRq8){4Mf2b*Q%WtP=s}p(;m2`^Frt4CH&rRt))Q9#*8zeM~Y zokoX|(sCX;gGkT1tF&e^$5O~$;4Q2&lV0C0;O;g0>V&TEl&m1G5<5y*ug#O@wMrOy zG^kOkZGnMjp&zgTbRu0Crf+sC@YqZ1$p3>P2NBn!O&@qBrg${^?Oivh-2i^`u?5H1 zbD97U+11Nno`9C@sJ4%pKWY~xuuSsu&K^*__?e4`V0*7ZD~kV6*}-M7m}@t|45Y&( z<2AGQ)7-VFMP`7pt*@1)?3wl%_YU&(vpDBuCtMuBe!7wc zzYFoILXF=hhXEr639Y@Tqf4XwR!l1x0obfAkyPvakd`syfXc*=R;Rfp*lpE9_skk>juw`5xyI} zR^ua=n{tk2@`tbH?+Gv!Eg_W2CsL3drJD(jY056_s46SS4|i)c!_EBYG)MaOX2s&4 z?)@5?zE&^~7NpbEHF2oE9w^&em>7TaHHnE<+D3%S)8b3$99R3H7M3T^NEk{%WD&q{ zAWvCp$M=JsyY;EP-#dN;+fG6#Q#+IYJ8=B-^uHns6Vrc_RhU?qIR7=Pv}keLe~rW5 zK2oZjK+Rc>d*DUW0v!PrOaG;lY-s%TQjOexr62Tg(8H>6wZ(@Kw2m0szZ7?B@A|S& zn3D2yYN*#Z40)XT^h@jOrse;wM5FHpZ0g=gMJESHMMpt^Wvv1pmH=i!I_`XHxu1C1d1f^VDh+x_YN zdZx7mib@`|hK->TnDiL$ERTa;Ct$5n*Vg%btSi+8B|!0PRMD#^RU8Z*H75+5FAi*s zr$nHe>?127COv9TX5=LpA5LazJtF}VnS(Lau7;Suv_<<m&#XwpvBd8nrB=ycO`j^weop9;o3bj*ujnhdgV9l@d+ts#+}O~sk-gu(m^%y?L4 zzlka#GCG@cpr$v@&0sd9sD>=Yj7Ff}jF$4gFv;dcrggGwFDK8+#Ulc0jFylT^+nRC zt@mMYF#GH()0mbTPa%Fe`asEMoq|etB$*K{U(!=mRAK?c*hH32pd3iK8@#UuQ#yNQ z%b5wxUP{r};;JqDd}t5;r+w9tp^b1Q%$yZ>MDZlu#D;JjwV1+&j_za7@u#eOI*rNo z{d6xByze)aDtLE|%YNj7bERBKVKu~^rq3so5rTfFAEVBV`1dmFLRFyw*Ok=EC#R`k z&GpC|JVSp^v6i4&r3lii+;50&8CwFQ&GJ6J`6ISwJX~wTv8opiehY`y-03+m+Tkc1 zzp1uxP$Y#KO7GXSpLGk1IzczZ90ia5bLg4whgGSJMFFVRUbqbD)%i5*f>oYm*%u7_G(7_eO% zvs3x4e8&C3+;WDGvcxG6<`a{W73$7c=UoE`{^fSIllhC#1ZZjLzTD)p_LVu*F?xRxS5BJfrqI6S`85YKKcFBFDb=O3Wb)Q zy1J4$I8oYgP{yhf|G~qR$9E$u(tSRp6i$r;0!9G=;4+y2)zZMRA$7!Mme2y%~&gzQkwn zw7%>E(K)UuHWntyLIy83JfoUk z;c3Q%&a0Z-*r$k^AonOx$fjJZxuGu#+#sl)tDGy(6*xMQ9M6N{GW9aLx)BK zLS#!pVUe;bD>0weil~^28pedMo;k?qYQ`v_D>iP>fTqFqFhK_&B#>R`po?wXGb^bH zR&RBBl_1n?w9;h2D(cL^<<7vIO;+eV4VTTH!e-7xpl67gOUs*^64AY4C7wr8DM~56 zA`4zexWx6&Y3ezsnA*kUKt2?#6j0m54@KTZoWYEytX1MzSGY_6ewdn6?MvmSu$7ok zM7mILlkCc5iD?1jbf+TsxtxsyIM6`v>cwR?j?opx zdoxI^vgR5Dx4LW+cqk;e3iD7bRql{PUk9iaG(S=zvhZFt9I?M$%sk~QqXQd3{tk11 zLXz~w%tM^{j`C}-dR1Q4{4(Ot#feEd{4Ao9+u*z%7w}pY&3=BQk+ay`Y&pgZg#MCm z3~9qub;kUQ{WNYI=BwhA#fhgSA;1aDKEN}p#n+mI1@4{HYVB=a%J|qn=s3q1P|3H$LEH+mTDf?bnmeFc* zHvXl17~PHQYR%^Qi5CQTrHFnj-xj1k{8YXHjDGL--nF8>94uO7NKP2zwSQN2x5Sq) z?xiQD)SnsJTPvDFSPQYK683kR+&eMxtw#Jb^tTcyE*~)XBi>^22bA{r1#yMl*GMs`rE5KD&_%zbL4(_U%RW*Abu0 zHnJ|o^&!0g6Yh)xC6oMcBKNw`#_i=94ucol^k=vLw0CrB>V*>_m&JS)b3tmOchZiV zB9mRc<*?>Fuw#giso7a&d+qX2RHxJ3a+2c~2({_V!`iIm(65{s1MYvM`vKV>U_;0Q z&P~`a1M+oaps%%o2e4_}j)^f3X4`@_FM8>o-RG7Wg^sX>uOvxr@HnQ?ez<8UzGx_|N8BTepa&Qq7|8|JN+WzjB;bZL60e(Amj z|ML4b92mvuXZLx3y*`}Ur{(!7qW8HIPaki6Kb)Wu(%CyTDCzgs{B|>u;dt}?ZR+i_ zEU0~^hwW=NnD%r(IFC@d{rz}Zo6_%<{Tns%IfcCzbr0Fw0SFVuu)D{h%oi{$e!<)m zx`&h+j}&I12Aja-`@q1CYoziOY;J@R2u%nJF=C-Zqj0xGUpuJl&_ZGAALajs@AtEW z%}(a4`4)L-tp)0!Si#X<$ZV9cs@(2cniUGdqEh(F{R1o9yt0B|0^Un>1H8k{gYL>G zb<2#q5rS~fPQ zB!%EY0$1fw8MZLUrI7rUXq5ZNw<;o^l3M%j6& z{a$nzaY>FRPCQKKH|6(jKIm5Q*WR9c7JU)|tLqCY&u@;vqJg1`8v;v3=ueo2)7uE- zDKcqTdI2BJ;Lyc{0m&C?Ot8m*+yVUo>P2FxEgcQf|x0Ge0`kwG(s*pmxuJxA=;bHhXW;0b=D0}11G)V4PykLfzcc-GrkaY}+T z3=NW23_&{qpWuDK=nz)GcF04p{(L!W};VXtJwXHW15uTzzFv=R!sT5o_83E#f`uc<5Rz&S!c&7~lA`1+8 zuZ9C&-|7Mvh;xXNrQQl}6}^5kdUJBB<095e+FZ7UqzZo;+Zu^b%<=g(DO;wb;OW&W zCSyfRYh}c1oowNk0bhFg!T{@x%J5JEmz)wk18FAf*Xb}Lt^-#Cvnd;sYbsBX&eUmb zlvctEZJHJFo|D{s*Xp=}e4@}r5w+{;Ci4FK5k76krr#pmF;z+6KNx?A>5+ruX&`$j z|7ia%RVt#C%mz?sy#7-w>;2qYzy57IYJX1mb$~6%S3bhL=^63^VjAw>evS`0O9|*E zQ1IY3RMmhrit*JnSh>t69Idwj;8pBjOLa7J)Owk^;^@`%W}9H@Kt#Czq`a`aMMFYM z;Rn2bNSaS=Y2P9Gt`An}Dwj+M!wm+5)K#1n$>`4i0)lm6+SIU43Z6Qfm(oleW=3b{W@9~!rug!oV@0BWlC%dMt2BfFNohfbaLbZ_{nQ(c7+64m6;}02?`GlP9 z<|&o-=-t3M7q=@_3HJTFUMGsU?J7lfbs&1P8KIws^cVWtvd=v8KUQDFvUTV}1GIf-TQSe8re?c{fVzXuTbpGMh(x#>u`g+)rt7rD=B1$b&T+bs7W1 zUC46@z*2BA@LC8b0#H;-a5L$x1S>MshMhhNRV&N$p(wpv7f5%<-!P)ug@tg8R#paz zEpm3K);$_uUh>RGsf8`CyxzUXb1XjSA#fr(?}<+LVhbwdplWUL%)spU71}B zR>_kOdCF7rFySx_?pdnXi$8W;NK?D-YLN&HcF9)UC(`;hFc*dD&|&2?q%mzKLc%^2 z?Xi_gD2}#Pozwc(;8JjeY8t9kZvsb-H5Uq{Ol=*zsIWed%~dWVW!kLMla-{^QnHOM z2O)fykvf6Gk)lIb2U;s+K+B_=!;6O^lSkglBFl;@j z)^5_Oc9po3#iBh~pIDDIA5-cs1XQC5WGHJfNcj@YIJYW>sFypZCk+k$FkvWc5Nr)n z<7Ld9BSE0ngn1|PwiPS{9nUf5;jj^$7FC!fj!X?X#*!ZrISX2>gFovtKX zl{*n|`IQZP7QYm0nshr`wwSqF>_1s-oS6YuiSKj_w5y7ERcCo+#`aIAzla>Xo@w8 zTjG+ksncgfxV(!Y$kEVe2UiC6ba(=))(k?!A}&O<1D~fC8$}i#zRw1to6|jxrmGgi zagcp;heeOkCsg~6JV>GYG_HUzC_eay_B9UL>9CcBWy>`K&D+ zo2z35b4QSD02p2{Y0k^5b3C5hUys8hTr0qV-l9@&(8gUWyr3aT-wHE+dEqusa3S~) z@?|Nz-ego_Wu02YPOxD38dUD1TDn@Uh(^d-8qWzOss%Sr#>e1V$`%V&*UF80ZgV`P zhPGSj%W)t!BJwDDli=L1o(L>0q&1+yt>`w+E9Y=ILR!7mAIlFWnjjhFB$# z7jro5LPkpymrb+2+TUDxF`z^p!(NrXe|A6=1bBM-&kNvp zI5%uxeYkdmr9DHa7Tk{*C(4t6Z9W<=JluzmVIx_J*O}Bk#)PK9$ykw(ZWt~72s@|& z^IGh@H{wW95(oTjAD)%Iyb%#IIzLM4%jX7h+Ws8g01;YU!LsN6G}^ZMtw{GdE)( zF*{#L(weHxDxmhb3Srv4^nJFC5DxJYqMK|>Q2EZ-M72)TLrpd5uziZ*0!~@Bm}e7d z8<%q#$h64^-s63^a;b~B7ec+B%p@jmPrzwU)o~(%2;+=Hy;?w!EhbH8HWAS_0(G zA{adY*?9GwyO$e=%~?Pg=@!m60`lU-&Lm%*1lV_{(2bz>(t!or`$I=)0qu}l2aaYS3D5xy>J&-spGCQx>z>a20*X7AZkVW<;3aL_b5u-5ov<;d5I1I zHJt!g=WoVZL2|QVm4vAF7FIGNU6+xsu{N#y`iB1nPwAZI{GW-7`M(;*ZhXJbUt4#NkA$G%MHz|l^`R4g$EKLOv-5a3w*QtO&<5|%VsB`SPBwWY z%)qOE8owOTE|##m6z`1y&h_s(=!b&@F&e+%abaLFPb|XSZyrJV z7wr%0-05)vrE-exx)-0OipO{MyjH#S#1?)%$5?X>*|SeYS2Is#izvyB2rOX`L@~^w z5aGppaJ1Qqq=sU5J24pw&s1QNV{=I>*G zrJU9ub@UJdL)s|pP{gPjBeK??2CgU%1rtTH%*yCz+()K%uJ@2t{!OGwxNwC|7{`;} zMNpzAA%P5}A{qr^EG1*E9m_LzbdWV6tAqJf@X@_W#N6`@$@pqmeA12g4*a0aYghu&J>7io_G~3tAY#3aRwWTXI=CyEu*0 zN9X_j9#pXvU+KpnZz8kC?kJeY=uLk9+oT{x8i06y9G`fb0mhwKzppjK9Ka*V(Ra^0 zAZEN%lUCnU<4-TQ62`5@E|Z?fXuMrM%S?qjz&gI$I#BZ_uX@|{WbnY^dAZ6HynKLv zorebnOcqTT!fVMh+Ip!LLy!=u!WUJ^p4f*h;QC;keq{_|OF}s!zOZIi;&uMqvqDl+ z&+t(!=wFz=@>G8Fq8A2X*9^Jpb>%6t-;vHL{JuvLFz9514qYM!V@4bTJ^>A3^@*{J zst8TxaLpw10pbIYnnmX;>?1&EWBirRm!BL@>SObK4`QkRJ03XC5vO_l(c$b#+dZg{b7YlyGJ3LmtCs-}=^k zWVMvut()Xm6X`A7C_{>NT2FY9b}|sCBtsY)J%5q;f~Fci?cVl^)|vmb`_G8n-$Z?9 zV3wU@$G_=90|ENV8pU63V69HGF1|%m=^5;9UensJz0$i6^1^)&i z@t4NE6bc2rCI9NrP)HniD|G);dBzjxWwHHJcUV1EFp6%hx!A(X8Md!Ko%U8-HX&$o z0r<~qXH%blJGP0&?jml~CyN}%!=6rp=5ileVg}SS!3Sqe^E$&~XByt_F#2k-+DV|+ zIc}091H<1cR9qGmzjV0UCHe?)^DbF(o#o~{$7O<4%GEo%$ebn|8gtFMTPOKjt}JxL z(tuiDp`+`tPa2@bwYT<$FnrfTYPKm?4ZBQP&Qu=AVXCLD=0KUTsenMlR{HC2bY9dP zb=7WBBi!b&!JC%26u4K(>?&=}P6UtQMz8R=`^rt(jOkQ+ryX61BOaDBhoGMReJX;XsTc zqdQD1o&j0MXhW%hODll_H-n17i%HsEywoTZ(fM2NmmHgR6AeVm!)KnNPqzAJbd+zA z@sDPj-fjY+_+~e=I>hzkbJJ_p@#t5j4VRN&&G77Pq>8GDc3HuWU4+&tzv}rk&hibQ z?dwuX)J&n;Y|5?Bk*jYT(ZHgUgoBAfTtbJ4b~Tb^aeZWBD6AEgx&PVzQsJPWr$wzn z0wq72C9%6HD$R7;e>B^iSdNVqXbIlvki8pq!4n2#d@=Jj8Qcke@O6rY;H7@90pkH! z(-IZJUnyNnkVR{HYzEjKh6+(_f9?4? zK8DffJ?C7O;)8nR%cqmppWmY21w&sVpPWG(B$sqZGH0QcanYU&RSP;~$Ko~1jhW-n ziygx)T9s1hQ6SgaHQ0#lC&L^13CXyoi+tFGqZ-|mGJ*saU2~O1#AH$WeI(UWA1h*& z>bwb#iQRzL+ce5mJvpMnM0vK96~&7FgUyYQc+^%M7h978oiLTZwHpyXt@OF7bE#Q& zvp7do1uc9tI%mPyn5#ce#CH*1(b@?y-M#)X7#5iBDjF$9w{?rNB=ILSBob*Yv%Gp~ zXEP?+_-#It$>jqv#2Q-=TkxjS`}74j2x5x!KK;%Utu2~#)iuRnNIK*rvke5zRebPwM_R+OE4HOuDR4@4UGkrv97}(xw^tGW_$atBZN`Fj%WhF4_i{` zrPOHsf=p%d%EfrUoum=>+Q775jv z*D#o-x2RMz(9jILxZirCg{5-bWmtQIM0zI!H(S`wWn&Ld<~NN*nw*;;4&98!7ZF}= z5wiVJjU%pm%ciaAy&W6K-)*RbAZ9g=a}CT;_u&)%7tn<4DB?eu5%&M4{$*rf=lIw9 zH%6T|?vMlK`d0l7p2R%qqdq`164W*78cxEkohXxJ2krGY3;7>wJsV>V25gL@l{!8) zChqy+rd8VUL#td$S$5rRx3qP8{^|1ZNG?0&v};s+{*f(rnPI_o&-o|2_O-8X=FbuI z@rR#BmN}Ite&WxG5`AZOdZ9*|=S~&rTHA$>E=yP;^^)4H1Wd*zkC^Kwy*>k)5KzV-D9bPQ_n-G9f+f~nmTr~i*08`jCxPn zomO~LExMO37(WCIam1I|bNzl-%?gn_bd0dg zas>9HfVdyaKuYTP?4Zm8m&C$@6>I&*_tYDbyrQZD(KOWHQHrfCah*e&x*^fqMipV2 zX3M`B;-L^~a-5Wuv)n{&F4gW7E9`*LNn7Aj%5nrX8*)At{x3lY0GIzfVWn>U?d@w$0OP%KxND=u^QOv5zR@N# zn2dp|>yr^A9HIe8u~E{g?4N^P2@twspev6A(I5yJfjmKap3r>j;F7wDLq-r?>d3H} zBV1s^@y8dN9bx}iAc%v@u7y?2+ye1wx+GjQ-w**7!tEXopkbB235(r5WvF>BK*Us7 zm`wa^9dsd8Zr~PRN$7x?e`^`mUvXVXGkv2!uY_Y4uZO%^`lhVj6H$Q9PmCPU1=q#+ z9+f4{>7w1O97eRHk0>%8AMGRM3k zo|w@E-pInr`ZS&?(m&J&khKqn$+#VP!a?cdcvR>RUcY2G*k-+2yPbH#G4p+hAg-=c zHcS(;w%c#`yx8{nKw99O+j6eV>pivAhNaWHpVasJ0{3v|z4@Qo!@~LhYLBvuk+X+` zDZR=+H%8(AkAn{^3=IFOM8+lthK2?vCdu%Sf8hU_MSi(|mZh11sGXLboTFr8nTM5< zrJs?c0ita9CY-qw4!sYfS%T+0!GEA}g6~OB`gvqH{)T>k-aoZb(*(<>5+s|oXBiIKc1XB|eYfV2Wg(Kz#1~xwZ3Y;vQeu_o{B_@zpN0-o|B2lH{IHz1gu!JI{M2sR7q#G27x*hZ-3>?V% z?9h%YyH4Ln!OjSNq(+wisK&s+$keFfqd)YTm+=TPBwveKxRb=~gEdxI`|WK{mI>;qn-qCp5)QYm#(m z1h;m&pFKRZGBebZVa@zH*c+kz;aT<~b9JdZIMh^MXYb_q`((-N<>klqa0dh+cbf0{ zKN+(!{{LXi{O{5C|6%<281(?Mb0`k8#Yqkepi2QH+J zL_a-7{pX^k(`zugKR|WzE`j#rkk4W}^tusZ6setns> z$n9k|SB;04OUi5Q<>f{8@+6nHJ~adlFbm=2{U4D0cYl0Y*#7l9(_P%L3SvM3`0JHN zl!i^itYyuPP_Ts7y$;{x-hijA4(!7#I!fQ$3 zc`-UD|L{n!h>WfXn020>n0;>2RX&>0vnRVKxp{tn7@a+KIIF6$cLuFg$~DqNQ$XVF z`XariBO%PV3d>z*i%%BXES}e+9j(Q_0=-5&T?aE5)Z?Zs%7IjIz5`1` zuCdAM0l?&w{>JKYC{wQ8rFN%i(8Zi4=0Wv-fqJ7ru6)=Uz2K!BW6V~HG1Ze)`wI0r zkV-*26zMs%L6-F5sAbIi9JA$dlesS01J?-Ugxg?+u15WU-9Y=i6ZO(@PhOj5v>o=@ zZ4iH@vnnX~vt>c5uITF(yv?zC^grPFuN)aG>%X-_y2;m(MGP1rzh0>CoY04AvNEL# zP&uwp-wpUS=5+m5R}idHD(k!7R}~*R+9-RKWVlL;#J>qg}<83YSjmpclURi?$~iJ_S$1Er=x^7o|mIu z_J_Wm+nG$gEsuG(n}a2qSJy6`s6Srpt8?coSuVPR?C+0g1C^@@mGPY z60=1ufD2ckBFynp4Dd^+ylVuWiShY^o2b>SkhxeMRPIzqg?}p~WzrtcB!xOJV1u1m z|9R74!muWmK*JiTniOf;kQW^!u!WG(nk2y(7L&Z6Yv>6lHz9w8y=1eZ)&w!aQcZkHGxr`IA3_wanQ2=f)B@ITv!AQD)|k;GvC zX1syz-8QTn?gI<=&DSHoTpi2a&hfDxT>aR)s@(3JxCqDb;C@p! za!gZPQ2p&k36=LRAS4jy_JK=5h!9%Fln>Z~3y2>&f1weFAiMbH_dL7(w= zcEAos7Yl};wmIsxu}V2bzpeSOASd9{HIqHfQv+tb<9=>V8=O|5LO7}QuPng^7aNDsbob za_4WTJ7<$p<<;D|b)C%3#VzG35mK~cRyB{ty(uc;RgA$aU5sK2xmcRP71{PV&%CP= zD?z%LRpb+{LWs*?*Myre8;i6oLrk(P(+~>6t31s$fW6pbJ-ZsxB6~@|bK*~{wPn_2 z>DcFZ{YLPYQ9-Jbb+X5=&c4C-KT!G?zx2PS^gong?|&)7$gNML&g@{PDsmEKvQSx$ z5d91MRt9)IR+kVf5~}O2ukoxQU2Ahl9wXshcU7G$nk_rauAJR{ z3o5(UZ72!7$kGZk<BIR`*#le7^lg&6lVev$0h z@ypr8(?Xd2Em)d)MceDK;=8?_U7wj|YHwf8A6l#}qn$PeS(Kxx^AE`g8h&mzt#yKE=i)uvAdds00y=iLoYn0$|QS zyy8@xt)Pxq{;Z@B7nvX~9soH+Az=x6jPe1q>EMZ+qf!Vu)<`;2=J7xSQ6gX^g5&v; z;{+X0OcEE2L4b_F5wrPaOeX%hg>$swe3Ew+@wx+LWCIyOq$t^zNIePxdJ^Sjo}dM} znZce?WTupc}=VfQcoe_ zo%r^wGn$zlkz>|n(ieGUTuc^jq|dU3tk7U6DXRXAVsw>y5p)Zc`FMqaM5y&wq>Io9 z2^o;Qf9}=fR%ROVNOKre7V!i5bHEw2M|Cg3e{>Fz1-~+GiN~`32^?nZay15YU#z7`+Ow&F_rNnq!*opheK=T!W%)o zLzZ(Cf>FB5JmSai#wNLSoAbReL!#6{Cr9VUCBC?{=aF2r;F~a-_yRAzSiAPQJoSmQ zyJ3-J?)0Fm@nQD$#yZ9He(yV(f?*MS<;6h?b@F)spusa^&Wq*9K_XP#7Zf$p73S`V zIi1pLW6i%KV7JJA|FocYFlpNIBKL=Zu(8;(P$!nllX>2#D_5t?J{c!4tQs*)**!D( z=}zCfZLHyPJ&daD+5@}QJYzB|CXz@=)893mv4@})xQA%;DO3ZJix6z{I9!#$L_vmC zB7j2}_5env`W>fYng~@K!TLrJ>Iq6Hv=Nk-cy#z(Bjk~>fc_Z!EF1&58>PG_Fs&W~ za)pY_nZYJPMXWM|syx3Csup@VKVP%x5h@e)fl@`A`C9O-CGZkzN7X@>dB7z^Z!J3@ zf_32%OYg!5hcGwAHof)o!9)9FrGp{yUT(y2@5E{G$>X}ZVXNZp8+Z)jqWC|M{}<=} zzgHV(_J50Gx{3d%O&t4QYU5fzO`>%z5G+gnWQ1RR4ZQAt6`UE2Vd(QYf^$A160$e( z{+5@Gtp4d)2PS?%Uv&(Ka}K_unf_!$8qcf)hUzHTDcG4Zs2^Jsad6SYihm73Qu=6e z{B2j4a^(5~J?CW?rNWqV(&>1w*jAT_CA@tT43Rk(@YQr%TB)UKr>#hLWxHmnPS@9Q5k1I-L9bq^t(E!i^a~RtG-lP|?$xm87cTfnMB~WY>NI zLLmLA@_jT1OvS|pd4bdZ%lZx(c%y2Nm`e^}vY!{?KggOr_dCPF3_hE2suk3mJ8?QK z&drR);`M-uG?#{*M6H4>;gSkqXjIk5t?B|;%K6@xO!D9O{REY?_!;(0An5gcI0#i2 zK&a#sa;qdRnHI>DTITr{7#FmQWglvkQ)E#srID_QQYf9e-drPIDC0yLqACfvua0}s z*t9RR;-;dEVUZD_Za!HzgR9?u0JiEgzy1RP=Kn@BS=kx?O*4~sA`b*mf^Q#Dd^n-A zSS5dF{2^jv|Bn3B&+-V!aM|@c3y}vrU$dW6ofwIR%srhetK!X`IPm5Aa7WD4&8wrc zN0dIhI2iMXO-598weU?qE4HnOMDD4o!%dvr_2Bp>=OAY};|X*9PzFhpzqVN~mcA~f zzXUxl89SUa=-=9=>~nkwH27M>;Xc z*K&N8qORExjUK%c80Gc}b6;HiT=ZA>d_y-wR`-UN#`cl?I;Mvl+;$3uOeXBm>!u8g zFRa<&Opk@m?ddk?uo`pLb5qyo9D9FmJ5(vU*j9}G@aFc>jWx*}cS(dvfT@@yMDkLM zO~9!5L&b~Wo^vMpg7)~w=JQ^VKYuBoHnDcP#-!;w@uCOkgA7NK!1_<;xMcGmyB)r{lrGV|7@#*ORzY#7BjO1>)v+ z3b?`Mq`EngF_>eWTg3-3Ne42;OfV^!NoE6m`GIMXl=dq?<6z1$mG&Eu(Lj1g?Ly`x z?j)3r(qNXik{nfR1mMm8AC$dAbSUA%Ets6xwrx8(v2EM7ZQHhO+qP}nPWtqqSO0gt z+r8eK)bOiGt*Tv{+Y3P9xOt2JuI_~VP;i(K8vpJv6*Tgs#@;|I~ z=nK3a+82RfJ;)bj2pl{TL3pVkPerJ1WyP_n`_ZHO-qb@`_By+CvGMkV$L54JFZ*Xp z>$GQmefQRL=hXe*AL%=BnpB85aYA^dN_V$VH!C44ZGWA5A`=rSBY&p(vZd5dtYnb* zfRPM86f-G3qCX_OS~wQ zuVdu_NfdY=u^KolpkbZem8O*%;$5R2)FMpz#_{ublhNxMI-YE^h#WRMe}Wg zO5=)=36nHv<@{GfgtF=RlvPfD(sYDsqDt*@TJw1o3z@d10gbwbi@HKVhPK66eO-MA zvvrL6q=-mEm3pOdBQ@$_R~dOEi=r{Tv50Y)xB7GUfm@(?qeyJHAtEb|qCQi-`dZ~0;T4S)NDE@*;EGz^Z*{5#hS?WWph=ee zOpK&TZ>rDHk$;GPEbKxl7icAH2T-8*{Ow~8xQ@dERR0$Bvv{6ddtfDHvm;8e3}*dH zxg)$tcaVV(9LyZ}BhM$+m~B@;TDMeursI+#*T1p|f5d@9r{rX=Y~jtbBlLGNOI7Dy z{(KN}NZ~xCBv^Sf&M#C#JOB`2q$|(RLvQ{rV4XJNAnnADl74HFDt_WA0S4xezi7Q3 ziJN24cqFW-j*{nEZ&`jEc|qoN8S3Z}B}Jr1Smxav1m>W@Ms68}9gzq#@Y**QaAtA) zYy=K87HEOM4u7@w^F!hMAyW0zLm&BqN`<=OeyLG=f-#F>L=~%pALlcX4Mg<&JagN0 zYcnU^^7=z&X_kc-;`c^G^Trj?rA7!QkU~jW8F9Cu+aLBPOtm_59Z1FyrPcxy)$Hr& zw{HM1ZvsfwwvZR>AzEPKRzXOVDDcR(Mc$IqQs*DZpZJH6;+}cT%bOjySW64bvUE(8 zFL&}46gE7fqk3{X?I@vT_J}o=hud9)NH~&!H@SwK6}U?|MBxTC!$D$jwxO4H32v5cBQQ|#6g$L) zXM?s25BS3chjI{(`{l;(1>kc%fk(E%+nu>va;-w8^Sjc5Shl15kdEAikS86(^E(=x zp)ZJ=jyF|51Y5HAHp*LYXX;}Kbk(T-DzqHmt{d!R*s|PiZ!EN!+EQ`Owd+AT*{0Q(EjtSWTu$swBgHhIf4T^T+d@K~_qDx}rhly+$t zR4BA7b3-{zl(2T-gqp-4#_NT29-R5GZ7j2j3K{KdGA}bQRGoV(y5p)PWW#%NI*Ihy zu6P|~w<63p-#a=kS+Qm8td!kpBAT4`*(29(`G)#4wGrc(mHGhd5Ca=8KNzOEnJIt^ z!`JbP&Kv$Ur+O!aj*MU`VVj?4e91<HuYLCA1m~T zynBlF4J{@ehxg-rc44}Dky2Y2|LrUL9UR^-upMibs6%iI7WFA zQWp|s0!&`F5Q4emZgGc5Wx`L5A{!`caPsdoX<5)ZF!}ASM;|qKw2TiF$l{ssU5v^E z5q6R&>I}-{5f=lz#laV?i_ZRsf15x|aBYy>asQRcc7QZgql6imnB79;pWjMfGGWE8 z?;qv^h==Q8;h(Wugx#mDcSGyy#)6cSs!fN%DNjk4TZ`&#x)s<4Hhx*6)u+1hI)kJW^cNPVxvy~10I5|hsJ>lzDF66IOt8+|^J zIKy0*hXBHkg07~bg1>QB7ZW>Y8xjQO8yPH2c((T`47sQ+DXnXNup%_7)AcbO6_%{N zk87=&G3M0z+pmx*8o9FFp&`kvXicS?CzVyD6JKJG91I`#w+Jrw``YQgLBoABP6%k} zj&6Zi;pPksLt=$nPOZa-6b?b9Wi%eOF+2Q^5uogW(~Spl_SKnmoq>B~<9y$sE(RVs zhy(kewDGhH1HU!%>i-nGsc((zIev9!v=iL2K}5Rv=hxIjl#Y%qfet%n5r*L%LTp?+ zv{z``xJnzlb_?v^UsX7*Ds^rvD=$areuL$Q(Ukocp!&Z>H8cK4^2<_BZ%o^TrKt#dFk z5P0v4@kjcnzR0$m!h?!{%RBu3!|3gN5r+hnRt0$nb#Zhln1sNj)d z66AY~pWf~qUpu4R4N4t9ZZTw-(s&?YviyEUed3W$KLMqC6r)2b_HPt4ILw_g$lGKg z@(*CVgR`BVSN*F7-8*5i0L?4(J(OJ8pCcr)iQ@`C0QWi#9G=Z-X&e-0aSV(j)UaCl z_TU%oB8SFqn}dR6sxDFw+Iu(8;f+?Za>L~ZKI5V-UcPHA#EOik(Fj#?qCqb~h_8P5 z!En_uh9NZrY=)dPJ}F{S1e+n`f$TI9DKZmyrtsBp^&zYwYy(<`*fd#c{J&wrLxctl z_ptRzY9elkli?^slm;yIA!?Ft2s=ZZ20-=ccLXg+JwvdBap0r^VcA2pgi&(v5`|>3UXKk-5kuZ+0(`dgkzNyikOaw8N)%pZD>OOc;FG= z$ca!h!79e_rV+!C`Q^tPitkR1Re*AWrA-ZAF|bShnO>pOJ+A(OiQIqLEe*`9jPH%BVSAc z>c><}{e;ie2bjz1Vtdzuh-Z>g`l=92Qr}V1)peo=FSZA(oLk1;hwVwH2YD)5Q2o)) zbUL@x#nv19wKrG@CZncSp)dK8lQ$@}{GLF^!Ti`{brq1k$#!kV0~k8MAv{gGVjNYX zU{@JL0H%f2!B=XE#_THJAAk`ufzNU;O8l~cc#btNAS+%0V^4ty@-DCQL0a?bpKR4b z3{k~{Z5QuDr3_@e2GPYx;TEX!CY;2rLBV@_WNW z8|T!HKK@ZH6M`s|xr$Pe8Vgy2O$ruGH4cd7wJ>G5|nfVvJ3A+7o_IeCXMlf+7@$G~A4bo{( zE3S70@;wtrMP3nY%SDnpDiFx%{hTk?Vu|c>Jy?E6e`hLU-EW1WmU5hQ_X4fo^I4?Zz?imXPfX$o+9MZxdCe4`X^$oEX+<)$tY2={WS&U4 z_{k!LH15(GALtnqj#UH@PkbLMF63~}!N#G)s@nb_1z$T!Z{04)Ol z{p?zr3jOG5LyZ|-9&SV+w1((91<+8vNYGBCg zFkIzr0xxoi%xP5}wFD3v9q`r z2IEQ)wY?AY)e+h>7yl1F^w{-;!s!_?CbMQx7n30_3*8Y{C!nlZbv+{S?4}IAW7FGD zGyDV@oM0UxL0=58VY@zgF(v!SeK`z?r#c_xLe zMqHja^#lGto1f9l46sAPDS-+yZz38BfRyY_1NXMG>+i0W668_mhI+F`>pOP~XJ>Ol zyVNc1u9?I=~_$~_?3|jBIJdI14LEw&M?X=9xKA31U9YS>jWTu0vl*-W2(x#jB|@!612Iexri zkfQgOb7P2F;oE*_ANO^%BEZ7guCV}W&%DdO9)Xx%B6mTW#&&&JaJo8&Naay_`FP^bM2K)~7*H7G z?H-lHn|d>QM9{_u!Qu|I>^RkQ{DD&7xZsy?5Xc(H)!9Z)FEwr@ zwT97lQAchd1eX^8GX0j; z_sQK`1NscEHMb9{UJC5{TiZZ{;fhV<{&*-5HGY$5?7L}1+Go_YM; zsn&{C;}H5KpENMiM^SfZ6wBYI@7m9z?n};lLytJ#87AtC>=XScN4Ban+HG#yDkbnNIwnxl@VzY(sa41YXD-w-_3$YLze%=lD?&nI@qF z5;~rs!+%*PT4rqR%2{_>YLQw=bnFHygN#})vv1n7m1w5w%9*pX(AY-q@EU$b=-FtwHQRd-Pi*{ zDY3A8u8TTAxI@VQB4aB#goRu#LJ_NH+bI+tj^SyGBbWJK0DO;-r>Ew`i++A}%q7 zO842NIYM}gfbTbW@AZ_33m`w}m2l5PeU6(b&Re%2I0`bL05j>|%+etNid{$b*_s}9 zmnB;Grwir)F9^mhCHIDZ6wzF`XfDyxzlyVKMb~{aPR$N+zvd_|j7eh}Gmn($6U}`R zD1fJuq&N}DF^yfvGNg!y1kxZ~vDlBSu*^Q8%mUkW{VQZLO`t#Q{g6~>Xk36y^Ga-I z4h`(UMH4rV6O-V+stn_SQGWETxcl{{YWrj!Wf@4q+#>mupfJFZW1MjAL-tpH0kD#{ z>y{LIN|d9bkYOE>dt$6ri~B$#G3 zb{IfO02{@tYw>yG?f;a$zQa)2{`FPUFQ4;x)HxnazEg-k+mkNqq(JT3lg4Nw9Gz5* zA7p1^MuzETJ1(8o;F3{p;1W-IcqmD(a8{R&w~Aq2oy%Gj&_tp}Kzac5T2$#U>g??@w!TwvyexX!I>kfo4CLXG}vVR!kBLgLUY^*JV4mroSnEG)aV= zUvC&n#6yHy8gmdFN$EvBFftE*CBF$YR?P-FdxU)4@^#(~^{OU6=y_skLw=jxMQ1;c z$u3-<-0XVscRuFd5^L0xHE#CD95ufHl@;zASOc}EZC#wH$I|M5WMK0iZ|(1tH0A+R z9P+dg0+{7z(y2E^7GRtaS5B814%{{+NQi8G$Bzd@inBZ&kK;r^x#nMIj~B*wr#YXj z#GAqAn=jcV4L%VS?5BVDAv6y#vUD0k8Clan&DmMNIB=RnlaKUedTFmPnvAA;s`-_S zNc69tq`cSD;SlK0s@1Qt> zQP!63eAY)u+}?N7Oai`BV8-P$G*C51+-Xm#Oy4t^XPszxblTXi<>*rcBRuOA%P4&f zUd(dwmUhM;p<)gfEY1>XHb+~1Rr5w#Xqm_mFYl}*6jYtQfcc#XKV{+{jJdGpKgE4u zKMLPpQ^JBkEgCJhPzwsG9f3<0FXxp}W=>Dp0w+Ts+X2CXnsw1=r~RbDDLqB|Y(ICe zmiOP`{UNc0wq4^9+(*sq$IMrgxDWHuDgA7uAuNAFr)X}P#gZoIgxmW5I}04OZk1u< zD#Bh}CPPtKMLi;8E5*FaRIIaPk4ww-JM3N`2hMGbR~8#`YX}r;L-j<*2MZ&}hHPuu zE9UD}^9c3Ufik_8Gf0~0SK~(kd;$q;(`egq=8R8~z1=S0p#PXahzA;O+S%LDzEwIL zYx!U+@l9MqacyNF9$?RM~~T|+P45X-&06`Q-=Pp!ccv5rn~&zYA)w4 zi$qvPp94x_r}R|`=My%FUg+O=1%ZRz63JU*q(KEnR6pd#8+aqFMb}@SZWJuukbm4qW^Y_W(MO%ta^lcG*#G98 ze{9>NQ}SLc&;1PJ@Z(G&2?e%Qp;w>iyWnb@oNhsI|FRc-i^_DLhdn`*=zgm>ps;8p zod1yQGTWf8;tTUSp77rEf+_oxR;zja($;rN32wK{ALurqvjTH9_J+1+)vfg#mxSQ% zDU0gU=V$r;`09GYG1rG*H`m!`%koi`T93jcXEp(0~pq`M-a3>t|3Vm(&l)=;$V zR$*XG?OS8mg?5Te2+zi;phhY9oaA04n<056jUU8Q>cjbR^Kf1=GLX$WNBKKHpzvr& zhbESU+nF>twnnLf5D{prR~ulMnk01)(Dmjgs{83*ySO*`zP9Yn?-y#@t!t#j=*_P* zwej=4xlYEmO!>KGU7pBf@$)jCR#l~a0q!7#oX^dZAZ9`$?`ZeD`F+|sR*YX&MaM!n zOg*Zxh3;HRjA3e8tt_2WX1)1OhH=^$xkdq^P|+m1bmpED-x8V0kPD6BE>Ym0a1ox$ zqv4|{YpB!d1zg_2f!ds`^FJQ0u+c=Rfq5g@RzlPLM~YJKI(@gNZx;7RM4##qauXhC#D$pov|6y1{Kbf0#!1ndLkh0k*pHT zB%}JIR+|rmHmWu};Af&*wLNS2&@_hnRp$&prz6s~uoJIM1ou+kR~`n@?x**9`I_$H zcn2aU@IJ%vLSMy~u0i((h1_0)zx?YBRH3VvOEDJIne$csO~CmKrWlWm=c7&eqOVgb z-rd3@P`ndfq{4Y>;wGghU}#_kYM*K_B1#}Mu>aLzv*!G9Zb2(MHGuhg0>MFr_Z!?= z`s^g`&#Im@eV9@&Hi4Z%o!GQ`+!T%$XHbxGxoh`+E$ zd!@b-q-2C8Vi|1tC=4mMQ^WLQ&JJ-FEUxW032zxg6%vFhYXHgL@-NIW=xl0Ui&Gj) zXNqt+8JXhy$UoGt7Se>4T6glss$1jV4&%=@>U>f1w^1bk(`G1J8!bX(d_SpTWO_xi z`i5wha5+?aOdqJ)8nA**xC?213@hKtQTiKia%Wc%Eej&RsZ4F8|56jlbggWEM#7$J z%e|yuj_Go>cfZvwr3?A>89>FFWL$Rswv-gb?6Q8})R^myMm6g94zNTqMDG#?s00glfS@r58$E z*cX>9sJd*-vbp$o)5IC|<^C+uz&kEEkpIjd*!s(r;7-YzmL8ee-VO?4bugCL^aY2C zb80PUuHkHjR^K@_z_nQmXjSiqW#U=3xIQx7TGJ)6nNnW|z(^rC_xj!K&r?|Z4Ex33VgqpgW_QgM*Apm{8pEhmP z_)$D}!p7Vd*tx#d#Z8Tkp%RT!2ks~-C*wFQZ}gbfOlNGyy0d_Dddr5_JLfc`Hr?m1 z>?}<}H_%Cm9`6bKB?i{im4*hOaJ|;VXWpgZ2kBeY759t7dpUP#`pV6rbu7;OA+OYD~w)qc58 zFW~pDOJTaklT)So1QFR!GzN;-EB!xO z7ayi_6CGu_+e~Ta39>-g_Xdu3T4z|l?e^sINeDz-+qq(J-@|(DN;x( z${g$IB|yb9GbtX%LvGwJEO*<7zX?_WK;yX}NRg#1E; z@&-04I|p0EOe)7G8Brz5E9n^Kl2V}3N6N5l!ICcHfMD=3$o4|eqT%gs5L9UQj>iWK z?MJ+sZ??b2Kd!0F82mz;y0KbZq=V@1ngfeLf)A}kFgEX8;FO5c1BtsHL%*}cku0Pr z-2#GG)gG0GDjp&tCkl$;gH@{+RZCS-QS0ULmtGd)oX`0R+C~kU5l=}AW%)c(r4Ebb z%cPxi6{QMNA`-%vNg)CgN@BHij%9=(}fMvB2>*1wI$K084{9(0??) z+(rluz~jO^{15TbV=6+yXYmFMgb`0e@YEu13HZxn)(E5y;XSeGn2JQS@VBU;B$~w5 z2}}yZLvAe%pXWu z+CbR=c}$GA4b7PjSIy8p;L0*9+VHJuKc|j~sf}4fn~OV^#edRj$5KS@!UPx@;&_)C(3XV=-cV*SY8c-W(J5Saat}eAh{- z$+&=JgM1)|`(NrL*$)ltwX-&3;ZPg7{8Okz!74SE@S$i}7=t&6q_8-~$nvs722eBh z3`VCHo@waOFDp(0S$n!&kRGwez@7OMq^GZdDL1AwwLSiXmz_D zLXoi8L{5T{MDxjs51M-a)G05oRrH~Ll9r~AjlX}F1qK0`wXQn7*KT4_W0M-`5TJ+mz(nl7Vqj>E;0n~u z(L7E2{W(e3yOdHduUDzOHryER^uh;nMrPU=AQAxHoO>a<+LL3r#Anh2keF-XA)wZm z6xTTrrOV#k1*62|VS7f(?!F!!Q`=#LTom~#q&H=*c`fm&(!ZK%dnHAUwF7Gh?%MSE zI-K}3QU}t8v#HUF$dp1orG!JtX+-v##9>`REZEpw{zuQZc53H9Lq0v}-KZFVcz-;Z z>u2E+%PY@5b%85!OazDPr8zZG*`{PLS!oux%o12e_P2)Mz{l*VW_De3M-bQIX{qy) z)XMgv!UZYRlN*RQca+GjM?@Q-4%LY1`#9avzweyM#~AQE=@&%R*c=65vRugs22>O7 z<6}cSwsqV&8{P9ngFMa_sST~)_AZZe-~~}5-$bcmAN*Ca%v5dfC+m7u4iiinej3_v z*C7{nAHf12M#sB`2^U<=(pkPL_6-24 zE3X^RSA*PE?0!HU>z5!pl-=~e8)*yjSssFS*U!*bvxQ~3@LQ-&Nr$i?cmXQ}1ouEw z!rC)GPae|2nUTO3hEkrIttN4`2th0-J6^Zv^;`L(&7()(Wu)=EF%HTVKxtuhyfHUh z%b5WBO6c$wtJx4i$_EqOadA1t-;q(?G{9lGuRXw_zAh23nGEB})Y^7?)31YJD%IY1 zGn=J8x+zbXU*FsxN3bcJ&~|}+N(?l|<6e6izwok943%>9W_XeZh=$ZHq&K0a0A;Yi zt2(egp)~wY8aQVH&H%`gApH~GI<_`Q8^W7yZ;O``*$Fl)09*_a;YI%ZHNwz(6;T84 zV4f17s5^h^C4FQkIw4{!ukue8(GIwb;anWj7 zQ%4mJOTJ8nYX|r>*DrxhZ0ilA&u_>FA!oDKQgILM@vrLlX(E-q?ata^zB2U4V5!99 zhLQdg*46r1|NW=+4V^zT=f6Kz%!*_np}klkMS2e~^a#HL%f(Cds9b8JR#hAPM{P3R z`~p}~$<3lsRuSkN&!iW|9m@k_Ds@jl>Pj3u{3*UCt_2iqXPT`eGCQBmS)jmms=x@S zL5HP!$6Ag}X(@o#LPs;G(oJb1&S(@RAmlDC{Ra|<4kM>Pb;>UVTRn0o_`Pt_80$Wk zW4lK<0BZil_wG?wepz#82{+E;ghlA(;6s52-u>MDOBU&9y4S97V2gC(6IghdC zAm$Z%zAPLhnq~!9E*xJkfLsM0_?((r^0s?QK{&4CsUi%4ONCnUiL_F_P+SVBa70M) zmTi3hXV%vRkbrIc``0*E+vFsvwrBS>2Sr{)BuR{zQMb|Kd4px95eFORxc_1Gbk_Oh z1p@9#?kuMosRaCYPdd(z+#d;|%3n{PyvwSH5n&&NG(-cb&_}iETl%L{;tMVUktCWV z7j(*vo+Vx4=TI$C)roB+O{$a_G${;IR=NnDCHf=$LtE(k)MZ7-LDrUw?pCijm(j!j z`0G7VWm2#J#fPUyA`;~BlrF3AQ_DpEpAddq&&7v}fx2h9k8owIZ)?79X1-4xOYHhh zt&E`k20en**cUy0Fzhi)adr}SFcqczLC^(?g*kpJA4AmHxMD>SytaX5EZj)$g zbCU4y+N7p_{?l`do9Chk8kDAc-NHb>ZL6zYTI=vuJ1VOF&JJpigz?ye`(x5RPl`kE zY!|nY)LK`cKfnH?)8om8jh#dKrCzH?_2Ni#$^?1wFmC(6z@+B+J{r0SM^Rm(7QzvN zSUj4fbbbk<=Tx4FLajord}G=A61D}cvG{bMPARG!RJrgH;suNaMRSyfOm*H)K^3{> zM2LoTb^cDtiu}ux$ArH_QC692nJa~<1<|p9bpflQX4&-;?L^>1k&B{E5?gcVL-7~6 zFQp%3Q9zlVgiNqpeo_ErVUiLsim)&i;%Qtejd7)8A_9u!cPzPx)2?1?Kak*zUFN(u#8f7yb$BxAY$T9X((co*0kXX>cPb%j>QQqMoz@*Kt~^TD6oCXdquEoj~W--i=Lj2fYs%Qbi;j$DlgtwTEuojD?int zCC@s{4RXbJ=*x0y4}6NKu-;W*YXfK$tYm(Xc;GWjEVY>J6b)TOpz1Y9CCbbaTPWYy zV~+y5lz{ycGBbPnYOgY8Hi4IiR#Mr}Tr1{uk|nYum5e{%4PFi%{VmzMSfOIKIWP+mJ&aqkZVUy3h|b>Rc*yJgYxgof4h7HD zK<2CB=??KjO_K0?He6#paXN}}09_jAsP;e)%};Jpn`3^oV0>gdrZ2OF!I|Iz!%U)j zteWi)!#?aZfefP0xM+|qV1n7flLpp6c6+2ymRpyFm^74d>eLqEmX8-AVpx;#PE;k3>&YIPd z>gYoz#ZOx&jB6#m-q>^1kWU*~M8``+OF}`GPp=L!9r==yJIB%KQi1PnP6Of}&W4*7 zsUgBxZ)6qvjLj|6C^3S|>-cS4Vi~I9m>K65XDEow)frJxRcKnHWQ}p#yY#sUwvnNWa|L1ob)oqCduj-7!m%|q6L<}u$-X*0Iysu44>W?7k6d!UZcpR^ z+@cQq4-hKbAVfHBA!MeRXJw&eBTz(A#g<$B-v2?5Gsh&*A?m0NrZlKBYr_u~j?_oZ ziK;u5`$~0e{=VNIRNWIFECSO5j1+%6l)pWwCzdsZx|cPFf885jOT*$zgL4))0gcu( zuWSaHgeZ{;kdw+IO^Ow@5g_<3r5=PjL-n1F(56P6&I!bb86P~vo%v`!HB*#iGRX3J zdm}nS)t=82M3ZzAoN%&o()>8dfMNU0*=h4R)OZy^%jSWv0d$SV5Wf-MK-$33IB_K? zmE0gJKe7DwdJhvbid=?=O(f3oq8>lnuy}I2Q?DW$NbJit>A?47!Wj<5TN(xyx9TqX zE<9t$K|Gh5+rYzPwu1k+C=lqsCb@d2%G(2p#B3rJw#R}FJ72CY|F2Q_tpu6F(%AW= zQVn)yfo0i?Z0@pm6zM)~o#|dJ=6tAxcM53L(qNc<)S-agWkbAY#tym3h zEw0i(8?*6H1C3@W9?m2t4;bql=aXuMEHFayyux=cSDwQjD>qDH}z(zy3 z%6C+P_3~Z-QQagK*DEw-0HM6GuoK1rc8xvfjAPaDBhfRprPIL}nBMxDl~y?FQK1vJ z`sCqr_ixA1-cWkLW%MVTkU*3k>GhgmS7`41y{DV*E|3=mSnpyn*c&YDEW^d!DH6g@ zw#*aRET=R4cNFa<@vM{1{c5)`uZ2G3JtTi2Fo8Cw(fa*h%ilS?TVBE?>d^cFlrMeA z?^Kz8?)&5DyEL0#9hfLGfP=j68%x1{nKu0Hj(}M0`;!@Avq!0LpFq%tez-Q*kvh%% z99pvF2AZVULlZu_v={)&Ma8FXxHpv~^{eIlno2VfPS+-Ly?T)+J}3@i&p-X@AA^~@ zhx0u*=KDF|Sp?FIe|ud3&g`0rGgMntkI}9+BR4WTDrg&DW!Khf)f)F<-aoq0j)WlvVrbp8T7-eZ*u! zjgN##?bsoIX=TLs8i|a5hRvKP9A;9T?h!!=rus zmcpWld)b#RIy<)lhe>PoP&f5MbMKAN9nZ}yvooMhiNO9sQIR}5cYFjF)&@KmN2=~u zSMsY)YtV|CXE)c$V0myTe%=_HV3Qz2o-A~fwMdI^I()ehnlLLPUTStHW_ONQjMqr~ zK9I*C$cxug?DI43G78%h8Z3;cci!8Qjide^03+!Q&~@t)j%*Y#eN3g0Dkbg7e{y!F z$Vl%wDz>*d^-#wfmXMMpA1^&Ld~G#+d6RS9)-_~$J6#s?#QYm7$XGyJjIxlpe0;5p z5`?^5c&SfzE~dfO356-1xw)yy{40|!I`teNIU4&f{^B%h0*pMO;99B8WU7I~rK7yH zQzL=`GyFx<0yp$!;GN;czLeeuUx?N!v57p?djl`}TX(!R;cVk0FAtveQ@|hYiZ_oz zQ#9$MFID1jQDNj>b)P)G^m=|%XTiLo6w8X zvDG7<>H>XKTO$$2Iv{69%6NY0hj`$RIfLfqCWI^heV>MD@2$oggEcm}IXm5V+a8@V z555V8mI?U^HUXH!X^?{cdXOsO*MIJ>oc8WkAxC8mTbn#NIC%W9ad&dEq2Vn6sg5~6 z%~MC>KpEs22x!UHb@c5c=`r3ARkn2PWk^3kH%3<$w=7zYGI7@i5bAs5w)46kY#Ucs zrd&YBSi&wL$P70x?$rS|m0P>a_!AJ{sZ@t!C5k#n7Y9uGGFxoub}OfQgo8LnwZ1pB zwZ$%zJbH>co{^l4FmacrH&oR1%)JY;2pUzuy$W9y z8HGPgr5$z~pTq)_?mcd&;rcRqi_wE91^StfF!N}RV>L($ySz()@?_5x-zLp@@bRO^ zB(l{a1%!gxD9PRF!E^0TaTHmO^e}FL;01tqKyuA8q%&oJ6U~QZK0Ea`_is?tMe<1; zLIy4;OOPyA1H5=Wc;tB;k|$@oEEGP*9HMG}Q{K0W7(W2+DD6-7mYJR6NuYnTpmIV& z_G+D^&w$Rc&uswn;nKiDlY)!y`c6Z)md<{voSMgu+A^NR<;on@3~_8+XQ~sS$*HGa z)D5>j^Y4v)%7c&+;JUh!50qm=U z7a=IpX+PgwA`82l=77_ybM{3ARb;&=yAC2tlD6v){aG(G#>4C(ZWw-A4(S$--TdTW zj@4P}>D(Is7$Tln)FhPtF0RO4$dB9b-E!?_Gb%0ky*mfRZt zd!g|YpWxW4E)h(j_35PRZ@uW{Hyz(#(-)qcV8uj(&KWW8yuf!i3tPHw5I`Gku#;6; zb|8!?^)Q|ONBwva0E`Bj1Px9N zbxp1v_pj;aOE%sP6JFiis}f9l5@GJZDOXIXxpx7QYN`zR^YCztG0$*I1*H~yD;c2{ zt&*|{^4c59zh2nSbh&gGHRtbx{I_}j)ZPf=q>_{%ROg15p=tqj54W6Zaj-z{@a*ug zFe9zFJfz>>C)zEI^OJGURNqf&*Cv0$Az`K#lPSFKuvY_uQy6|_A*cUZj?+H=82VfM zY_{DO{pKhznG6=WELHpzt0yFv-Ye(B=(EFe9_zCe}XC}As?qszL=B=2Hw+J z#YEnxCDyuH3A3KwT2+~1K`nWTmvafd9nlpn^`~XSdC%r)W>$L)s0n)*S!juC@S}oB z)|X+^wWLmO1#(Qyh}w^rW7sz=4UC1Lt_crRzk$pz#5eBnIyMpjZLcuRaT4D1Uoq}2 zS`UDDP;IHNwnqP+#My8spbqwl1nLLtoN-{P-TE(G_SW`zHBIJ>^1GIrt>Xa)28Y?a%8mO zx#L1mXec*#O>7k|bsi>NS-qZT(sy{-=n+89Qi-d=vHG8UkkV5Vk&>f0yV+1W8`kG9 zPmo;7Yd-^HE;=$)(?$ua}K%XU7`cLmo5cy7vzKjX3R4Ih%lrLna0~Jl{*=pqy~C65!0TRY%OFJY6Q5{hk1uYvX zDMps@wW*SYnmJfHQ7~z^)GU1b35DSGlUdzXCNO-P46`x}*AzRp8&oVX0zNRd@Kav? zlNEgga?Ayn_<5);S`_mPej>)H{22d4uMMV0sJ%gr-?9yCvOw80vdT+F2+DE7PDZuP z_U;ISUDUUL0Bm#a{V|L%&>%dB`~WypVCsA3c==*F`tejm|B7ivgxzf`xckF{UcFsy z9Pmv2eG)_EQeAc`_gD*%4Yh`#k#m=Kvq6nQ`lNfK3GCih(st$nxMSwt@L z@i1@bZ9Putz6l0tU6(SB_@}cdIlO=*27!4Kc!_N{`u4G0_B^kAUV{dxuR6?<*a!tOf~i(1yk{{!^&@tl8g}OZVAa z-*4~{LO0f>w9DC44^doNY5ju`xYUrfAv{bp-i&h}Z+C+}i61_$fDT?qaxJOZm5}X= zi;dN-@H8C23e27eoWHh=!&FYh5ShPd2}vL2d2?yX{y;gPu(^Ifj;do=3*H}cRO|iS z;kCq8F&=a@E&6X5zyB1&+Yw}UEz1d?vX@Wdi8KweDI!f&vs{W$iIVTV!jgV)AXFo<7j zSE>?XpqkpFc;p0b(6U3-fiS`sXaaR3RlumleFdy)Lw`|U&gaBRS0)goDX(LCC7dN3 zr{j@^7izEltT}>5od4}y21~~r4JNi%HQ~64h7{VvYR4&Ny8yOBq8^0k$E#dOA)_x% zp*d216#6u$wO?w=&89XCTj+GhPEvIMK)j1aDLeISf?!uqu&N-lJP2EvbMMOcc@vin zlSp%#Gx3HrV#ZQ6F*wr$(CZQHi{v~AnAPusR_+tz>PU%eAEF%jRcR=KE%s?5A`JyiFm zNsv4wh?wri$(WKRE+I_^Hv#%;$v!CTEWM#uj;OB@Des<+;%mg7!caNkL^`+~^9F%~_a-`&aAg^p=zfNceq?DM z!vh|quMlie8P(?LKBW|c;G!Uxx-_aGWDEL~K5Lh#IU z{RUOXSsa&`_TL0{)ws!Uaaax&u|PNlhAe)lafAQ(O_rtm!@d}-2!LEbZ>Z7NgT})T zgbia&Ox^O3@a{m&?lJYc#w!ro_Z^=G(At>myP{uKS84U+yj-*jX-3PGtd`5E@n$=B ztNXiRgDPKH9Q+xx2c}&IF(9Bu3xT>mjfYOr?|}#3Oe2+q9g$l+Xh_Jo7mg0@uaC$E zdWp0BJcgf9B=cEvOEn>MK%NhD=g0r}?e;p<0he>MwTe{_RHopA?-arK@GA8=gQ(Sr zPNZ*KH)!G_Cnm#VUsmr{B+(_-(4MaE2w|K>A+KW}{akQ|EVyKb!2cYR->TrnI1eO3 zEU{TZu$Qq`kx2cP(IWMl#<2N2s;N1?PT$KBYY#G}(+tuS zh$?9VRY^a{7%}YjFll$>4I$9o2y_(p+MEFuVVI;^T#F8y2O))H565(esM&r&|~kTGW^M3G)7RmkVr_;*xp*<>{b4!elN}lL;hRf?fC(P zgVBoXHDVW+mDp#$s61nh1J9o6cmBNh>!*&trq!m5f|uFQ+^j%_LX`)%&7gnaxR5AB zBRz*0m8aXFpIJUsKM%bNiCiNkHGw=)`XM;5c_1fH)iJQeF;FBSyORoK3K11*f!TTC zp$VjrMo{U*EI=laG8Lk6m`F3>=$=DRDOQ0fK5EDl=teg?BQJw&e2X$3^=1kMx|lN~ zFeQtP4DvuUM5!OkP|U1=M47=pSAbxA3gt+mM3jCik3Pk8G-7D|vFWcj$>XEAUE&mc zTo!sS9i!z=X7DjJ5_@fUQmO=~CQ*Zdc`1u*I(@|q3} zTvi-!2yYDSfL>p!0lgu;ah<{an$#7m4QeZvc2Ldmis5Ag?3x%2h9k%b^eiWM7IWw)3mWqsgC5@PM!NbaizrPv2RobJmAMBZ93=St%3sq% zzbKjLEoMH;uWLwB{F6%dCt*>a$$;*pq!L4;&jODi_(Ie0X&wW__X`ec2 zmHfA1Z)pw(89yD#ClOAfDBa|I!3-%J1lHPbEu#;In>+-%#IwVi)Csiiwgt z({&WVm8IDcDLz`E(2*tzxz+P?|0MMWA@0tR)(7=+TK?+_f;pm;?m0j?25FWQwdM1E z9oLcafGW^1f%{;{-}OL~UO< z!~-H0fjI_u#A`G4(Q=B=i0L16@QBrtnS-(G1ScSbF%!TZ)cA<-#EcBttM#` z%@*WPfQsg24a@SKYQN2O!eX^mcSTdKELzyS;E3Wq`<$Hq2Ga`ysP~-}^gOXPcQsB3 zpt&)@+1Ie{fw((=T0orm4GMZ^+)`W8P_JT7k>$sqM$TbK^`y6hdh-8xb@!Vqd0KQg z2$2&;WiI{D`ETyHx8Z8VzEOZFhVSgnp3Q(`g2yrIs_ZgKZv$Mn=%1Ytwm&sc(YdEa zAzqDM>HNgtaNx-TQMi*B)5@^(X4efQ_=NNajf3bwoc7WD4e@J9o>c@Frk+<}Nc%Hi z4fUJ*2a%42^YiCv$l;3T!+srdw1M)K{4!+z87y3e1^(B+u#yHsx!u9``OC_iJlZAE zPWfLGMu%`IU;zY&6{##ACJY`a#VbimT5S!*8_W927H;lk_wTS_3-`DHp)rE$MVkeH zSG=($72I^NqFRC-nW)YImX%En+}5Yz^UPd?w~`0 zAjV)nsE`8wM^-TZtU&d&k!7)WIP~r9%XJa|rEBhUQW-mZPhhFTdaNijJ4=cC70@O) zc!uvU`Rc4N{O|I4rnw-XW6SNKduaJis#0yFf@J(P>;R?U&fe->SRL`&(6lDnnheeFOp;Cxsb zyzaFN@U7u%zV#i=*Wo_4h%NadGFY&l)Ee+iQ0vB+TiAA6UZbl@hDFkNIavv83GH>w z&yGWTzj>eo!PqA!2)&a+_8&37{BU4oSsc_@J)>NQ>q^$8$S-R3b}XWOYlO}-FH`jp z&5^l(a>=RDa+D66T{q^kU91%1DzE|))RqD z-)l1i@1R6Ksh|fRhBMz#u#vx6p6)~#;?>Bfi)w)irOcjyUWc+!fg^-C5Qtfxrhi3Q zZ)-~$8X?YV`h%VbPi1v>L+V1Rm2^uK_8JLY;>3gct^wFy^&&rB9Vp%#xc^2Igg%_n zX)zi+UQpMi!_;umP7G#|D=M&kXr)lYi+72tiBcLyKlVf+4`xY)%km0?nmju0=k-l% z8{sdCCuO_+b*W23WwP=J`fuqP7$k~wR=+G&>4L~ax0Qb1{J=c^%(^=WnNZ+`^BBoE z-ak z9uZHs%x(?@oj)q@azU)ZRay0i9spGgWHtTt7G>4{C{4QtN~az1aC*Mcb9hadt4Loo`bJxihea zK$RCvKX1wIY6-H`)TMEtpV-ja`o2akxDk1wDq-=SsOMeqVktuqDdq7<8Y$`~M-0D< zCUw_%d>gnfO>mwY_?dMZaygJ{K4}w6)s%IZr46W*z+{$GR{WM6J8YVp?c7Gy?0Pht zT^ia1%i6vq_+MS(59BS_mRvJ%4gmB6-V}D1ag*xaDGgG>RXkL@_Xczz_`lR^sLwuC zHPCX+yNbn71#TMsPsCI{spa;k6EW17SiJ0+J2QyilDM>HZbve!J){Sh<233F40@lt^}*f9*xt5XT|AqQFYLL6zMx|@eQrG8CHaI!JF@Y59-Vz zfmbA!+=H#R=K73qy0SubgtLn-_E#-=OVkXSG1A2vmzS?_L1M_WX7Nac&jXr|s&xm1 z_Lx&~$IraF(nGr*CQnNc(~4Fob&i3^ErJHhX~>c~gM)Ts5oYvZ&~^Vb(e<4kby-*Q25v*Qg@Q@~Cx$0E zp!o&fJu&(kh7^hp{pdT2^BOzG!%F4p6sB3r8bD`8&i-`bqiv_K!23o@2sPc22qig- zA96v~uu#hhPxxGfs+k9kg2nb%yvh4}KIaiTg%=tE@$*by?bmGIXYTp{`FG*;rwZZq zo?dXC{RDQkmVn6$G%-bjR^(F<;$7=W(D`Y@k_$4d-Ax+CJLUccE2#Bl1dAqy8l%hi zFLeYx*>_+6=vUDHL#xwk2SK?ec8&)%5L)c%WTZ50^^(GLN0T`z%L}_~%|+644su0Y zmNxZC&A>jQ-&VsN z!9SX7Bo_~9JHWz`Bgc1u7TL{DqJE{+(PcoEBAe%9{>HpjlY3IHIt&YP_4h_>`>qvCM_d71fYpuL-H@1AYxCa;N%TZ z6%I9thd|BXZp%R3!}%cBHM^Mq^V(mqeSo0O?j-l3$Pm_id2)++b@%unFLoZrwJ!ixE$PYg$$ zA4wV?-^x_+kF&wRhqGo-lpdf`q-@!jtIngax%4zd3N51)i8jJV07}%XX^2(rR2Kn! zpvn3V2|9Yxk9ml-M{*Eq`|93K^{+ie@41=`k=Kf5t|baiBIjuqB@&j|xX*sJS%Uzr z+><0t_Nb&K>qKz>N{*{!-cmQh%^2!I;Y8us-0rw1S8nT5S{K5K#mMn&?O}QGh~MUh z?G?u5AL)QxhHSo}xAi`jqEFckDsOM{kdHCr`#&A;d1nsiS1;s~Ljy-U4`JViFyA;E zmw@TE17dKnn^lM7Na&g)D_Sf>dphHN%5~^V=TQ%XPk>J#kH48tqSflHza?yI+1CI@ zFDvq1ZMs5<6@bfdX_XUG>qMH>4Yq6o(yVNpSc#iFR9r@<1uFahd4qq|=;eaYhZd`c z;z;s$RugW7wsA4%D#CCX^F4{W7~>6HM9na%ppq{KAYDH{5t7OodW$-{xX`MLd$Q7kPUhTTKOGRS z^a8ro_L2)6#LcnCCmZzU#1D?ZiEwiwH@n6KO-ZLh;fmME&PWxaoEf~`}URx$cHbw6A) zOd3|>{8SocL2Xghr0X67j=Zzcc_pd587R}i0%9+R%*CIY6>DH~pBg~^<7P}2Eai}b z)_;7fH;)$udwbheJc2sQq?;nRX4mHB`5(VoNpMVE?Wx?O*0e{*q1(k9A|_*! z9=6gDG5W$^6_evys>A3A(6(zsxhvBBvcoMI*%a*NXd@?k3F;i~+7CAzJbt>#0SoYo z7qHS%@jhZf4Vp zT^TZTKK`SqZwH%q4}0~aF$tqS{6jGP(6`Qdo!4l+t*^LrjKyi3JT3dIdD)l{ja=8D=+pBq^@IdJ^Ck7&j_NROyUc$GE2YRR6smlH<~cNTK|cGn3H7 zF^%3@5#lYunJyyk=C(=id|aDrjzaqazAJ9lB_^M}56ZFVJ+_mqp2r69HC3x0cp6WB zV;G|n$mFiTqVfTor?thUlZnBY?at3m`?(u8d3j_{s%~c68#fkF9vO^Uw(C6&Uvd9G zl`Sc0>gleiMqJjdR@zXUzcxE_0`g7EU}or}{OIy2mN?H~8~U@&OeQw`(e{chmGQ*IgM*qYB$y`!>7@iLu|i?Xutk|7S7t{|BKr5ALW|=LrVH@<0k$| zNqY;Bq6+=j_M-d0V)U;4@ecNNN$~b z+0YtlD0^P=@^a^LcFxJ2@ldz-o1|rS&qB*9l3V^hY~%j?zWLYMnLYjd_VJ$#Fft~D z3v*sg{pzaDndOU5`UKLXtsWnfS+7J=8=oj@ff^ZF2m7h4?-pRO2U zGxQ{_P*nS;ug8S1rAi%@B6${v`ce}oj$;+d(4eZRBA9?G63`J{BpZb!##4-OwL2+O zDH5bmLc&&*ldN?L>5!sP4HFgCijk-xP*tRRxShrV89{3n5jmM0A7GWFOGpwzlxWCE zs8YQjT3S|eNm5nB7pvCuO*B_x%S|+&xl?Mlo}v~X*-U$p*^e?;RC=K*3Q|ENXpknW zQY!aKam1;pKvxmGXpIKm10(~L2xiO;&fOMF&p9iwB?BJ|y(<8!K&1#Wm18M6_v|&E zHLPM;$fAsDj16jzTAtuUts_5z-_&f~ujvoR_Mz!v-7Zufsp1c#j_PXVG<|6K%+lT) z&;+!CEImbG&6q>$Y{0KT`*=-kSv6Heuhp?WD?`3?@|c+B?4HpvCA>yVeT z+WE8K(>CkfaBxB$be7$Pj7@DRKC#M|0=BUXLb?>x#>rhC0_^a__s*2K2U_-%liJ%* z(oyfW!FyxN48X;0pNsf_eez+u_2SvJp4-BP&q#yHD|zLNTh(XO$4q7!bO%>16T#!r zYXCF=fSnugj#QhqLggl;Ca9j}o-8gV1ytPaI+pG3yMG*#eWUO2hu}e}qTIwYx?2Qu zU`qKIm4sS`!kkDumQt_Byjl$b>*%DYzHi{&l%6=$(&VBEUpIAm{o~OGUd52wjG9i6+R}q> zP_xK;bT!o9{X}GA@BY^hq?s)X(>r&5uxL7ReIR(6jiukNp*givB%mr#0eCr(Htr!4t0MHnt`fUU+!L3u z5}n@g0i-~zlG4hi!aSy1g~1BurE#Q=@!NpU=OwFLD|0GchJJP}fUU(7>;)XoQ>s)TLu$~YZN&@7#KulAkty;^x>!G*2y_u|)n@7;(7 zwnI>r&5Ip;NW2T)fIRTHfgW|79X%7f+k6UJ`Xrl$PYU<2ZA=%#qMfl6XN&a8)fFsk zSlBoFbgRELu~zfM?@1R8{yMu%5g5g!frPSQr5&Nw39oIzVxi#TYZBRjtfv#0bm4!z zL4VjLOxP5}sLh(vvuHm4>e9w+n(|iwsMkyCv<|9ui*5#^#7f7-yu!kOHtaO)Vat&i zo{zyMXdQ^QYie~#+Brp$f+AfN-j@VyTl~~k8bR_IS8%BWR$951SatOjlBNl5Y^mw@ z+2LD^&bWun#1$b{mbUEZnA?r8Gc(t>XEbIcsM%ZFi?@}{2H^s?-#B>pa(B56_=8pV z3^1uR%MW(2?`db?trh-5-WNW%KEb_td+V%$8gyb>rBPMvlX|2#xmcbX3jUFWA+D-n zT~>|7TFNiEc-G4T&#Qm7tU~NASvDI=Upode!{&u>&#-Z0ow)OW^{Ui6EfDGvoWK~# z=u9)m%GBhT6lUcp1Q{rFNWJ9zG zZhxE!`cX}lW^T0hX2fWyTGKr52WddFe0Fy|AV>nlYH;j^tHD4eQebs71Nk&m%Hnbb zx3h!DTza{w7@2(MgfzWyDPUy|n1$L&;BOr6)G2CX%F|(+AV}_lv|v;Agv^btVCFj~ zMxA9dhDgPl*h>3HjcpH;HAW65tdi283s{Vb&XE--JBh9~YxVh-)^p+g{9WagZ4(Hs0 znR|z2o$lZao}EM7K_Ua_{5(w6wNC~}O)VSI7QEBGz`WguqZ9Dv%s%I#3z6{~WrJQX zAu{5wj&7zJ-vRG8AFex!na^J5=B;1-P!_=YZk8Lb+rDwm+ z6fbNU06PPd2SZ$1KD5TzM4Lb-jYz*^6_#Lhi@g2nI@P@g8$Vi3Ms^!F~o3)}u#793)ATNS{mWnYd?)~}LYddiRl4UCK zb-vAxUU4GAzTEM&;CSevYeAH3qmucloFY7CL3d*x0ZIs1)h3y5z76C048|ZM5o@;)o%J9uv3C zWnO@lQ%|$fO2WP2C?FG{^8VY^*_ye3&WC%rsv2FIuDvNPT~-y(umKaTWN~V4sN--ng|fngbkJ8g291BPLzmxS2Bs z`z)7dY2SFe^m_={zlQk;Gv18mltfcfOA}gd3zJ<~iYrsCP^xAdEp?XUAr75c5pA67 z9Zgd96a&vBQwW{d03zUnK0l4L~v0|B22jJI!0@ge1JJ-cC zs_Ml;D%5Piuw<)bC2Shhp$x&lvJ81~SivMVfq%(#yNZy*Vdx>|3neT-nlEHaGX5^W zRI_Y|Hg~$emMB{I9T*MepOU2xq|b_P1QC9?3N<+?6`YSZA6-%B&*1>Ds*CrPGrwTtQuIIyn=fy21BTN&kH3)EQz6&b`kfn8`9u>6Ep&aq`2^8NvHTfgFoW!gW#12zNOzB8T~ zfsRRS21GFuxKa|+>jrG@v1`xu8FM23K5kDU;U5HagrYOZF{j z>(0(WFA5ZE{3PRTqKUWgfavyT!7U(KaMKJxBNUslU|M#HK!#70xMq8&Bt+M^u6hdY zeWyifG&%t0lwqJLvn_DRN~D{)06m2HkhY>+PhE|b8<)NrOvtZa6Fr7|mnoHnHn;m` z0AJIzPr&8DP9l`FN`l#4WQyFZSfS-rK?foT{N%pFR1y=z#tu54(0^LW!LGr_K8273 zmg~V*So$U)d5N@J%1aw7i$!kllW$HdjCyJYg}FPa<@zo>SNk<4LBEF@?oZ^b_N0!b zN6!S{G*UU`ac`QOMIvRAGf7{p{yO5tRobDlgYVtA$2C_zF{7^Ia016C49J9wzp%2g zu(7>iYh}jmkc{sOd5eF0f&%~mJUHP0GTW8&KXi}(*Ij^(h2_6t@nimfEdFXPL`FUl)naqH z*VMW+pt#F{(AXRz_f+>lCNUf$4)NX@V)jz7_p*v9UEe=wd~H@1T2D332zE&916m_0Lunb}9kTxd;%g4+7~wJIX$W))(~zJcONW;X zF&S|-1bU=tNLCZ6AzMRu4+J=nWJ8!6;Wp$t5V9fE24goQtjX9AZ$;z;gYQdvAo2y% z?~A`71BwtJihv+@Mj!}@1tS>-XG4$^g-Z}chLah|yfNJC!vzPi8I+1t_vy^B4VVq_ zAxkH@lC_il8ef=d8GFgTwbsXE%bbJ<6N%3gHc6U^5|c;!a71FMvTis%aTh6z?Xu=s z?r@EGO&LNj4D1}-uMIdZaV)IVvB)T`6kAm0V z_)!{`?$v77^yTe04%A{2Zj^M2%$2a#d+~eURPWM%XMYeV!*o7t-J8&4O|NWMF9Na6 z4ueErrL%)K2S&u3yhglH7i-veaZoZ-zzQ0Ly+H&EZ6YMJwIO|^)m$HGZKGZ@C`FYC zVF2K=NC4c_w3>Qu?iI<*5j;fNja4 zK5~aMMjj}1f<{IN4zKOl8frexx~iP(cpDkG*&3u$$yZ`MB2x7P)s|zwHa*^+^24Nt zPN<&L%R5IgZ6#eXT(yEs`MMdk#&^x3U;70tOa)uX8{?B;XSc4&_)X%-j(Lg{Q*p2U zbz_UtJDp#4zG>@b2;HXzN|LY>O$I3BW7%gD`16;ihNJ~yyh>-Y8>{cixawgbF_B7A z*Tz#QlZlub*Dy)dvk$@=iR$10DbYwsrGjQ6Cs>va>Tw|j)Y3)BM#;rT(#n$y*i3G1 zN-_N(;R=gBqE@n!hKg>CTExJKaa1rJck)9-+Nd(IvXYwGN&$^wHao;|y4?m*S5sMU z6z2UR8tdhOKDGN>|GDX#p)v0x2)cvNE~Sa`KEuWe3JT^<_NA;uSr7nufBM$%r{UY7 z(HB`cv77?F1f6SYq|{@nFK{-F|0aYKCWPOV5R86M&4b1dadV&5Ahu>NWptI#N(`Go z>YvFQmDaiXe(Ao$zqz=2*bH%QfK&!X-(p%7y4=SxZ8$9lyF7Uv>OPd#W`(wrkGyOi z=S$Y!c*BMlaA(IvC=(ZQB3KeA1A*(M-Y)ALWUkzz^$WO`@P#6OnMB1t*bN2)CpL0~-2fVe7?UNtpx~|GXi7$ubv}zS-7>k2iCnpV( zyI7+nq$6g-dKYza(t?rYCZ~^&ZQltV!frO-DPv7y#2B??U0LV)h8Nflf@)JtV|w&a zvrKY{%7{rx#`UXH-9vGMMkJzM46Xr4X+M(Kj+cx6?JWIrYM!-0D%ovER_I5_CL9oL zdGdb7L@>+|Rs*Z9;H;+O`~0{YLcg>v>^di@jAAPlx$r^CYKF3eavauW(- zzt0GNgt#bBk#FO7Hsn?vY916%*CIU?2}`}`yrxF+1g!{WzhaNJ@w8IYA#)*VglfaG z?QAS?0Wfu$3bdVdeyF=CMK;o4wWmQ4IYHO*-p7ZaX|f_gYXU(JcwIlth#u*P@rofd zjR0Hr10TF2>(#&JMdIzGtK5y_(h| zNBa`a()&4u-vcGDp*E2g8)XFoY{VEG2i~E;l(D>+ll~o*sGGAo(NAK&mTt*P@WwVP z=Yp!Ni@snQ|b2HLKwsIKOu0uV73dF%N8TDP>!CXKs+jV{bSM+jiVFAi1JD6BG z`6zku%i|Z!`6r&vxH?kiIdQh725^JoZDSt1taZ7QZ;SUJ)@acC`8sV(ZUG}J!dl{6 zH48Z(8UI}dbLk$qF9jYE`tHnWEQ%HN5bkF|#h9!FiG zkE{5bF0?p`Y2={qwHpKtejLy=eFyv9<`#rFXr1jgxgF1Ga5Yf#t_5?OLF>_%>)!86A={$4srfpC zWAqy<<1c2LLE?R7_5VNfB%Fu65R>yJZ(B}s)3bfE!SbgilT;{<&Y8FG2{VykV z?6ElyIXeAgj~8nT{PTu3UQS|YQZPtvAYQ9C#rNfs=_4;|=i)qoVlQ~Ly!C02*b^C{ zT)AVS&t62MjNkQ@`rGh&->KK;#w=Hw@fJd?l}*X%dyO1gI_SIIgl6Jq`Q%nmFchit zi|TZFxQct7=IF?@lN&iC_jL)VWkeU7CAFPJTRtR@ehwf|<@htZ0TCXgyZY0FNn7iS zuAcAj!!km??gxG!S<>n^j7!x=sZTRwc1v38+3F=s)G#+=hNQ=|!?-Yh=+MJ-csKt{ zEDom|o!x}PaIW4bmAX#oN9G^ls~(@SmgfgzUscdgIF^FSQl%FnUwpB8=9aR3DB^&V z!B_iQL7gmx@~Knn4by*9y+oXh9F=Mk10<`0Kb<=tF^ah8|CRyz^8rM7c}49d$8@}n zuX}(&!^~R@>8;^+VY$*8mZP}=G3nLDc^&36355lq+4exdxQpH$Z{v@OVyGR|skNK| z^<7x1Xsvu>Eu$jan0nl>^Zbz{ws*ZqEhse( zwMJ_bS)M#rRCMz6RgsqV^^J<9`RmiNcR{N$oOM3CJy*6sG?=*xwo7JCClM#Rk+_K% z3I+8jB}D@3=o)#K?vXLr_Y&rAm+KY$x81g|AXG04bKdGEe&wi$ec!caUfsZV zo+8Orn!RZmMgvBU{7?fv1OEX&r2&dds$e3{lxiKAk8wYOkh;eh{`%HXK3AZm1pKpB zvYc{Dd9jx}cxLgnd1=sFyng~Q7rxAZ$_A;G{XBY8%!a*f&dr3%1FF>Tg8}K_;A4& zjvQq_0~LTBma945au^n>yP?OCMTm*4=hFeS!R;vrk+473jO9aUVz;xu!{q^f8e9Gn zxA>gcp&TXUEu#LFF~#y;O6iH3aq8q=!?x1!+3~_u@>7A8S3XUpM-*(SJLSxEA-1Mq zDTNH}J4;nuvv-$_I>}Nh7};q$J(o;s8kEk`Cf6>WvW~Wl$_mc_9jYm2L5Tjocb;Et z6}0ltK^ZT<(Dcfq?O%j7*i);Nip0>R%-_HVUqXtiUPJx1PqBOPDx|{O7ZVXC#unY# zmPXu|JOogYb9uk9bVpxUFdUKPhwd+VT#aRG{@epnidXeTbzkSDzfMv8F}qK~ z!N==pvTQ_~{bx8#RWvQECnT*WCz{u?kw(OppZ*7%xBn<0Q76ytPLdkeX7za4=BNsK zvuq2v$DZJrwZ)20uY1Dz_Lzr$UoKF9v2I8XPC z3j&9yZs;DI2atK;g!oZ>nK7U%&h?yl&&$<$vrp_eoWXa&$sx+63s9PRvM@zOH66jB z|DdsFceI~vUr;`-;EXmRleJMfx@~BZqm6?wai!+ zi)S48u_CSaiEDIBD6a>^umuL-Zm_mFInS_HnfM8(%M* zgZBMM>*x`-*{QV}`!3JpMe=ky9ZwNCHMyTdJhyi(+T;Aajk@A{%1xR?dM~%kQDo<4Py@|UowObr2ghy8C1+=?dDECC;7j;TG zUK8xI{-LvbJO07)1n+SRIbQu5UEHESdXt@2jQ*7+bO({mU%^Vj^GN_6-xWua=rSPZ z$|{;#$@+wp)8%SoX*F%`;@D{8HoJ)iNy4pRtAt@a)HB{^CxMwrBlSAybticbV}J3! z2P4bM=?&J)I?wN()|QG3qc=!$kd> zwDJ9E=OyK^2o;Vv(`2lKg z<+P-nIF2>DV^^U8|9Q~Wa8zhD6eKz>4bU9u`$ zWJ;3kZ>%)R;8cafHxh4ga9Q{=#bB3{r&4&3R?E~V+{~$!78|>ah1QA9!Pf?o(UR9{ z%jKyUi?$+|b#cq&rF(Szz~h8G?tklCleY=1xivM>g_lX$MH~0=5>sp>bOE5M>|$<< z>6i4kX7}3m%yS=2i5{P8^+Pgtjg&zZ(wiE zcoAJ~%im;~1RkjHG>kyzaDWB$7vU<%2F|$CZqZO0IYBMdcXF$VIx0F!=NVQwJ@LF| ztd}=ni(vZ>8mOIQ;ITl6IV_6kPP|4)D4|4MSOEOj;M1 zA)oDBDId!8bsJ4iI(+#;&Zi7#GpF>A^oLeVTF0iWh^3C1BC1$+4}zREVj(RNT2UoC2RR<(-Tw`>6}2#1 z-dGq+Y@Mkcm49##J}xZZ+NgHQ(NqtsRw|~}Sh~8N5u5k08|nh6Gh5%WD5WmtJ#;I? zT>jD9m=;lw(Ern06J$PVo>0z2TqRT6WN{ELcuZ;|kNR-^Nl=vpMvgWm-GY`P;&i$k zs5c!ZJmn(18?!~Q{#1EK_;ibO9*O}w4RImS1qn`WT+&UF0DEzx;;L52v*v@H4H&Rg zUCojw@2YEg&JBKhS8^PdCK?e5~MtfqsSl#!5u z`0t(_ZDU(gY0Y))>H_8$@oOI7zURWA`F_+!?%t2H-peuRcO7|N03S<6Ji%C%O@GL_RK@6f%*=rYvs?v?$*u5oF#C`IQ^k&%seHquVI{g+DC}=HQqM2=!lr&p!DA&U5PY9=2%~AHI7aUZa z;NvOA3FfNhmgl__NuoSV*6O#!&v)dn6C?`;c9C|BM!dRmU8xT!e_ znS>D58rG8d&sN*kv|ab&)Ht`CI#WC7z*k|HQI1_iv=mJ4^AI3S-@{E1`3F*FT`Q00 zr1s!WpeEJ@<>@xD-CW!#5h_GkQa#`a1ofoYvJ{lIJnlOC1R`dTC<%)^j!9G zk4~|^>mlghDmJ>fqnQfOmh7!yY3!N(ijF3hCfq3?A67@_s|go+J%HlgG#`!q0+Wd9W+cMGsCzuL~uxHp-Am+HZg&$N22Z3GX-fduv}RP3210 z(bo2MzP5Oen5zUdR_rgLsOYu>rhKm<^K(}e?Vwd5oF-=0I{s2=8DwNtJn?W7GEU)! zg=%FXCc}3xwm>a+!kA4wW8j~?io5U-e!W=Z=Eg3L`W2M2kx$=tIMvFkRE{7=qfT$a>j0 zd)-J|#68AuAp=-jbb(nqjrGNW)AUI$;`Q?oyI)xg%7g>#;7KjUyHh5T1bL!7VVJUL z&HenP0eF7Ey?FPMgJ6JvMM2W`g$+G5{rqbDg1kFmyI`?d1PdVaNARJN5}BEepRm6^ z(sna)=5qSPsdux>{-EYP0a^F_n+rEHTaSuB@CHY+_vdPI$Rm=tRpaj*;>aIX4gIg5 ziIk|M2G+*8pNJ8TShl<;C_#j#7JnzlWu}U3ZACa}>PHn4EYk}rX8iHYIJLNiXgm4} z#p6nqN*-0;fOb!J+@;Q5b)=>K8x{Er^<+g#y%` zY0u~H4*-CS3DnvDT%3`OfPvt@@&r6Q^eUbXru6bgR?056^l}7j%=BWGPR=d_EF28~ zo@`u9o#@4E3|&k`OpWbLOzCAz?aW;)2$NlH#U8 zRLB!hF;i}jlNlkON5yDCfEn1T#x%7lU3xs<#0G>Ju59p~M$QsfF4WnvO9)%8%5f|G z>10k0+6?xokK3lGy}+4)Hph6aBBkGLuo+e3Hr-8(=?r>5Gl6eGSa?i1{ zTq7O+r5WgVUxzK!$0Fa+fWm-aLKq(KP@spLc_=>25D*^ZGeJ0cG=Km;bs^rG0q&N= z*!bnLZeEhHD%=36rzc#Fpv>=2)djed+;_1gP;M@7)XeGdv#A>C{QDV6$HY%<;bP~# zEL-)sqO9Da`Rva;V&lP*6l3uf)?=)RbPtStY}XrXVqMRc`TV|(n#%m~fw8DP{p0%S zS$VpAqo>f8C)PJwj+k&449m~L83Hza(dn=%;JCN{w%KEJ%De6Zr>W<7Rt(~T( z47#e`6xeH${O3Hy{=B4~pp8p{t~#Uv*Bzc6l2g7wiVCnDB)=a5^s06Gu^>0uhEaOS zvR^HyLvQMchqI$pUiW5Y2fcqpVOCd%t`Ck5y8C?8FbjU%gM41>FY@w{l_Kh+4-YN) zZ|)5}EYhb30dD=-TOVOO#x!{5p_3;!1Jw_g(#D|ofUC#Aaa1pKJnf_X$&-~&TFN-y zZd{NXOD2I2;^M2nUtbPhutWja-@&F~m&}DyMKv|jkVygt4zjw^^_~51U$(3x|66wg zjF`UcW~|vJO;N|@ZI{H^$w(j?NpDU)Uic``)JU&oQ8l!vcULtgMg4n#KrKtEJ}+1BZ!J5fg)j@^~@xbrO7h)c>z=GYf}; zQR8@h)O=Ycq)^7#mBbh`wrK`~kV2M`tqz&750hyavXrvRzBEkDl$3RZ36iu8*pZDgy`Ca_ph#2$Yy4KI=6fvH~GX>3L%~s@Os7`&y zbY(+;>)3^buEjAWl3x8U0m81}>Zk-LRo1xH#048ja>o3a%5jrKJQfS`g_X4m zj5^y^&vr&v2R>-#cVjlbWQmsbxD%+|-={s~Gvs(f)^>bc8Ud0Ob^U5-<~%)mn*sbr6l7Wt=B zR`*|EP3EEe@5`o95kbTqMz)8?{L0Qo=Tw`UBVTwx>Xb^sJUJ!idH}L_ZSI`%bZKdI%h~(o(Y^?@vFSFTzWN+c^6|>ED$C3(iG)IE+-FTYlzQ4_Lb`XRxT?)9p=h!< z_NwhCS}GID)a*y+TkvYC79DF)e_-%8F_-u$ugl=5=@wYraqkg;xZog0l`Pr1 zbG;w0D>)Y05(J*fH34Sz$ux&uplP(iApuTwp=~FN5;9E7OPxN!XtdmY5Swc7#|~Xx zeNg)0&}lzW%*!mEYgrP)IlLwGtrtLSZ9vLZD+7BfQW&X>Iuz5_GYhA8J?}mzB19^uBW+va~r(fowU*Em{le(T)l@v)DRCEX`ww6E0 ze0|g(5CgOJ>`xGC3sXafaH^YPl1A3)h--x9p6jJjPD&A94TWz9nikhj6~=|WEb(la zC>1yhLlu}MbqSqz3fNICSm(0j%GOQ@j$@0v@TOqFn6~E$tu|B&r_FH6SAPmAxH5c) zraB?g6GVRmud@=13~BeeUoQclpyR=DQTu1tDqR3P2ANEu(rbEq?0yJ{aJ*GZIk158m1WY$=q)Zolclo(bn_4F!fpa}Ify z#u26zTbCnGejE%~htqjj!%@I{-7H>!jgm2izopLJ)Je8$b>aKm()NcXJCmuWHo3J) z9a9nZpMk9JwNte0uQi~!3evkvJ3_K=l5IWqQddf5W+ntjzwJ1FTQ#oMeh$v(g=RaW zRP_ey@dB>_=6g{t8N(r%Ge6+bKAq<#km4BnHdA#O5#u}C6SIc(6fQ*k)2Qa<8&|{S8_pBOJBN! zhMF$+HLcQGc7FNWgC==ElR(`q+r@XvGF+Yudq78DMmGtA#w)vWftR=Bxm0HGtdJb3 z&&_9i%nq;&UXl*#qbvQvDKg&c15=+ynL2iUc0nQXGReCZ4R1X*eM!_g$k#R4@*@P( zCp?Olbik7&j?nr$ltnfE`mT7_PO zd6TPn>)hfSj%Oi=lb7k>fWdj=4)PnXJ*Y<~-0CqQtYCci!1OWovz~s(D313F z6hr-2KY#iz^<+-H(xS~nCLx6&BE2bZe#-2TN;L>^_d^7f0YYnm(|KbIqN{55!gGhJ9-jmVJ+%g^3QH-PH5wnHQoWT5v3pW>8K ze_4UjysYIb+XpA1I@PB(1Wjb-)}=r#h_-7jguQ~XvzL{MbTWGF%+EFypR(#M`t)kL z$W43Kwv^5;!~|^Tk2f-JW;X6Ky!nrxv=iwz9|BTt7Wq}-Z6&Je?YydLo4yo|>2WMo eUi|%%fkfBa#NgXl93Nam1E$HRpkQoyiSM6EDtQ$E literal 0 HcmV?d00001 diff --git a/script/model_inference/02-monte_carlo_parameter_estim.jl b/script/model_inference/02-monte_carlo_parameter_estim.jl new file mode 100644 index 00000000..4148c5ca --- /dev/null +++ b/script/model_inference/02-monte_carlo_parameter_estim.jl @@ -0,0 +1,74 @@ + +using DifferentialEquations, DiffEqParamEstim, Plots, Optim + +# Monte Carlo Problem Set Up for solving set of ODEs with different initial conditions + +# Set up Lotka-Volterra system +function pf_func(du,u,p,t) + du[1] = p[1] * u[1] - p[2] * u[1]*u[2] + du[2] = -3 * u[2] + u[1]*u[2] +end +p = [1.5,1.0] +prob = ODEProblem(pf_func,[1.0,1.0],(0.0,10.0),p) + + +# Setting up to solve the problem N times (for the N different initial conditions) +N = 10; +initial_conditions = [[1.0,1.0], [1.0,1.5], [1.5,1.0], [1.5,1.5], [0.5,1.0], [1.0,0.5], [0.5,0.5], [2.0,1.0], [1.0,2.0], [2.0,2.0]] +function prob_func(prob,i,repeat) + ODEProblem(prob.f,initial_conditions[i],prob.tspan,prob.p) +end +monte_prob = MonteCarloProblem(prob,prob_func=prob_func) + + +# Check above does what we want +sim = solve(monte_prob,Tsit5(),num_monte=N) +plot(sim) + + +# Generate a dataset from these runs +data_times = 0.0:0.1:10.0 +sim = solve(monte_prob,Tsit5(),num_monte=N,saveat=data_times) +data = Array(sim) + + +# Building a loss function +losses = [L2Loss(data_times,data[:,:,i]) for i in 1:N] + + +loss(sim) = sum(losses[i](sim[i]) for i in 1:N) + + +prob = ODEProblem(pf_func,[1.0,1.0],(0.0,10.0),[1.2,0.8]) +function prob_func(prob,i,repeat) + ODEProblem(prob.f,initial_conditions[i],prob.tspan,prob.p) +end +monte_prob = MonteCarloProblem(prob,prob_func=prob_func) +sim = solve(monte_prob,Tsit5(),num_monte=N,saveat=data_times) +loss(sim) + + +obj = build_loss_objective(monte_prob,Tsit5(),loss,num_monte=N, + saveat=data_times) + + +lower = zeros(2) +upper = fill(2.0,2) +result = optimize(obj, lower, upper, [1.3,0.9], Fminbox(BFGS())) + + +result + + +obj = build_loss_objective(monte_prob,Tsit5(),loss,num_monte=N, + abstol=1e-8,reltol=1e-8, + saveat=data_times) +result = optimize(obj, lower, upper, [1.3,0.9], Fminbox(BFGS())) + + +result + + +using DiffEqTutorials +DiffEqTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file]) +