Skip to content

Instantly share code, notes, and snippets.

@bwengals
Last active April 30, 2019 15:25
Show Gist options
  • Save bwengals/0cb96e9a756f09c66fcbf9db4ab3f235 to your computer and use it in GitHub Desktop.
Save bwengals/0cb96e9a756f09c66fcbf9db4ab3f235 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2017-03-17T22:14:37.816696",
"start_time": "2017-03-17T22:14:35.777690"
},
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/bill/anaconda3/lib/python3.5/site-packages/matplotlib/__init__.py:913: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.\n",
" warnings.warn(self.msg_depr % (key, alt_key))\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.9.0dev4.dev-3c0be3d94102ac6864b2e5ab52ae96d07c6375c6\n"
]
}
],
"source": [
"import theano\n",
"import theano.tensor as tt\n",
"import theano.tensor.slinalg\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib as mpl\n",
"mpl.rcParams['axes.color_cycle'] = ['b', 'k', 'g']\n",
"%matplotlib inline\n",
"\n",
"import pymc3 as pm\n",
"import numpy as np\n",
"\n",
"print(theano.__version__)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Set up 3 simple models, using one of tau, cov, chol as input\n",
"\n",
"The `s` random variable is there to make sure the gradients go through the covariance"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"ExecuteTime": {
"end_time": "2017-03-17T22:14:37.884808",
"start_time": "2017-03-17T22:14:37.879070"
},
"collapsed": false
},
"outputs": [],
"source": [
"def make_data(n):\n",
" A = np.random.randn(n, n)\n",
" cov = np.dot(A.T, A) + 1e-8 * np.eye(n)\n",
" chol = np.linalg.cholesky(cov)\n",
" tau = np.linalg.inv(cov)\n",
" mu = np.random.randn(n)\n",
" y = np.random.randn(n)\n",
" return y, mu, cov, chol, tau"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2017-03-17T22:14:38.381845",
"start_time": "2017-03-17T22:14:37.947110"
},
"collapsed": false
},
"outputs": [],
"source": [
"logp_timing = []\n",
"dlogp_timing = []\n",
"\n",
"y_, mu_, cov_, chol_, tau_ = make_data(300)\n",
"\n",
"y = theano.shared(y_)\n",
"mu = theano.shared(mu_)\n",
"tau = theano.shared(tau_)\n",
"cov = theano.shared(cov_)\n",
"chol = theano.shared(chol_)\n",
"\n",
"with pm.Model() as model1:\n",
" s = pm.HalfNormal(\"s\", sd=100, testval=10)\n",
" like = pm.MvNormal(\"like\", mu=mu, cov=s*cov, observed=y)\n",
"\n",
"with pm.Model() as model2:\n",
" s = pm.HalfNormal(\"s\", sd=100, testval=10)\n",
" like = pm.MvNormal(\"like\", mu=mu, chol=s*chol, observed=y)\n",
"\n",
"with pm.Model() as model3:\n",
" s = pm.HalfNormal(\"s\", sd=100, testval=10)\n",
" like = pm.MvNormal(\"like\", mu=mu, tau=s*tau, observed=y)\n",
"\n",
"models = {\"cov\": model1, \"chol\": model2, \"tau\": model3}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Look at Ops usage and timings by call type"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"ExecuteTime": {
"end_time": "2017-03-17T22:15:05.576449",
"start_time": "2017-03-17T22:14:38.444871"
},
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"MvNormal with cov log posterior\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Function profiling\n",
"==================\n",
" Message: /home/bill/pymc3/pymc3/model.py:600\n",
" Time in 1000 calls to Function.__call__: 1.081794e+00s\n",
" Time in Function.fn.__call__: 1.050849e+00s (97.139%)\n",
" Time in thunks: 1.036325e+00s (95.797%)\n",
" Total compile time: 4.657686e-01s\n",
" Number of Apply nodes: 19\n",
" Theano Optimizer time: 3.939168e-01s\n",
" Theano validate time: 1.563072e-03s\n",
" Theano Linker time (includes C, CUDA code generation/compiling): 3.501344e-02s\n",
" Import time 1.019335e-02s\n",
" Node make_thunk time 3.361845e-02s\n",
" Node CGemv{inplace}(AllocEmpty{dtype='float64'}.0, TensorConstant{1.0}, InplaceDimShuffle{x,0}.0, Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}.0, TensorConstant{0.0}) time 1.204896e-02s\n",
" Node Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 + (i3 * sqr(i0))), i4) + log(Abs(exp(i5))) + (i6 * i7) + (i8 * i9) + (i10 * i11))}}[(0, 0)](Elemwise{exp,no_inplace}.0, TensorConstant{0}, TensorConstant{-4.830961538632819}, TensorConstant{-5e-05}, TensorConstant{-inf}, Reshape{0}.0, TensorConstant{1.8378770664093453}, Shape_i{0}.0, TensorConstant{2.0}, Sum{acc_dtype=float64}.0, TensorConstant{-0.5}, InplaceDimShuffle{}.0) time 3.001690e-03s\n",
" Node Elemwise{exp,no_inplace}(s_log_) time 2.281904e-03s\n",
" Node InplaceDimShuffle{x,x}(s) time 1.635790e-03s\n",
" Node Elemwise{Log}[(0, 0)](ExtractDiag{view=False}.0) time 1.595020e-03s\n",
"\n",
"Time in all call to theano.grad() 3.124356e-02s\n",
"Time since theano import 4.328s\n",
"Class\n",
"---\n",
"<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>\n",
" 74.0% 74.0% 0.767s 7.67e-04s Py 1000 1 theano.tensor.slinalg.Cholesky\n",
" 16.6% 90.6% 0.172s 1.72e-04s Py 1000 1 theano.tensor.slinalg.Solve\n",
" 7.5% 98.1% 0.078s 1.56e-05s C 5000 5 theano.tensor.elemwise.Elemwise\n",
" 1.2% 99.3% 0.013s 1.26e-05s Py 1000 1 theano.tensor.nlinalg.ExtractDiag\n",
" 0.2% 99.5% 0.002s 5.87e-07s C 3000 3 theano.tensor.elemwise.DimShuffle\n",
" 0.2% 99.6% 0.002s 1.60e-06s C 1000 1 theano.tensor.blas_c.CGemv\n",
" 0.1% 99.7% 0.001s 6.56e-07s C 2000 2 theano.tensor.basic.Reshape\n",
" 0.1% 99.8% 0.001s 1.08e-06s C 1000 1 theano.tensor.elemwise.Sum\n",
" 0.1% 99.9% 0.001s 7.18e-07s C 1000 1 theano.tensor.basic.AllocEmpty\n",
" 0.0% 99.9% 0.000s 3.67e-07s C 1000 1 theano.compile.ops.Shape_i\n",
" 0.0% 100.0% 0.000s 3.19e-07s C 1000 1 theano.compile.ops.ViewOp\n",
" 0.0% 100.0% 0.000s 2.18e-07s C 1000 1 theano.compile.ops.Rebroadcast\n",
" ... (remaining 0 Classes account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Ops\n",
"---\n",
"<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>\n",
" 74.0% 74.0% 0.767s 7.67e-04s Py 1000 1 Cholesky{lower=True, destructive=False}\n",
" 16.6% 90.6% 0.172s 1.72e-04s Py 1000 1 Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}\n",
" 6.1% 96.7% 0.063s 6.34e-05s C 1000 1 Elemwise{mul,no_inplace}\n",
" 1.2% 97.9% 0.013s 1.26e-05s Py 1000 1 ExtractDiag{view=False}\n",
" 0.9% 98.8% 0.010s 9.69e-06s C 1000 1 Elemwise{Log}[(0, 0)]\n",
" 0.2% 99.1% 0.002s 2.27e-06s C 1000 1 Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 + (i3 * sqr(i0))), i4) + log(Abs(exp(i5))) + (i6 * i7) + (i8 * i9) + (i10 * i11))}}[(0, 0)]\n",
" 0.2% 99.2% 0.002s 1.60e-06s C 1000 1 CGemv{inplace}\n",
" 0.2% 99.4% 0.002s 1.56e-06s C 1000 1 Elemwise{sub,no_inplace}\n",
" 0.1% 99.5% 0.001s 1.08e-06s C 1000 1 Sum{acc_dtype=float64}\n",
" 0.1% 99.6% 0.001s 9.19e-07s C 1000 1 Reshape{1}\n",
" 0.1% 99.6% 0.001s 8.66e-07s C 1000 1 Elemwise{exp,no_inplace}\n",
" 0.1% 99.7% 0.001s 7.81e-07s C 1000 1 InplaceDimShuffle{x,0}\n",
" 0.1% 99.8% 0.001s 7.18e-07s C 1000 1 AllocEmpty{dtype='float64'}\n",
" 0.1% 99.8% 0.001s 5.41e-07s C 1000 1 InplaceDimShuffle{x,x}\n",
" 0.0% 99.9% 0.000s 4.39e-07s C 1000 1 InplaceDimShuffle{}\n",
" 0.0% 99.9% 0.000s 3.94e-07s C 1000 1 Reshape{0}\n",
" 0.0% 99.9% 0.000s 3.67e-07s C 1000 1 Shape_i{0}\n",
" 0.0% 100.0% 0.000s 3.19e-07s C 1000 1 ViewOp\n",
" 0.0% 100.0% 0.000s 2.18e-07s C 1000 1 Rebroadcast{1}\n",
" ... (remaining 0 Ops account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Apply\n",
"------\n",
"<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>\n",
" 74.0% 74.0% 0.767s 7.67e-04s 1000 10 Cholesky{lower=True, destructive=False}(Elemwise{mul,no_inplace}.0)\n",
" 16.6% 90.6% 0.172s 1.72e-04s 1000 11 Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}(Cholesky{lower=True, destructive=False}.0, Elemwise{sub,no_inplace}.0)\n",
" 6.1% 96.7% 0.063s 6.34e-05s 1000 9 Elemwise{mul,no_inplace}(InplaceDimShuffle{x,x}.0, <TensorType(float64, matrix)>)\n",
" 1.2% 97.9% 0.013s 1.26e-05s 1000 12 ExtractDiag{view=False}(Cholesky{lower=True, destructive=False}.0)\n",
" 0.9% 98.8% 0.010s 9.69e-06s 1000 14 Elemwise{Log}[(0, 0)](ExtractDiag{view=False}.0)\n",
" 0.2% 99.1% 0.002s 2.27e-06s 1000 18 Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 + (i3 * sqr(i0))), i4) + log(Abs(exp(i5))) + (i6 * i7) + (i8 * i9) + (i10 * i11))}}[(0, 0)](Elemwise{exp,no_inplace}.0, TensorConstant{0}, TensorConstant{-4.830961538632819}, TensorConstant{-5e-05}, TensorConstant{-inf}, Reshape{0}.0, TensorConstant{1.8378770664093453}, Shape_i{0}.0, TensorConstant{2.0}, Sum{acc_dtype=float64}.0, TensorConstant{-0.5}, InplaceDimShuffle{}.0)\n",
" 0.2% 99.2% 0.002s 1.60e-06s 1000 15 CGemv{inplace}(AllocEmpty{dtype='float64'}.0, TensorConstant{1.0}, InplaceDimShuffle{x,0}.0, Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}.0, TensorConstant{0.0})\n",
" 0.2% 99.4% 0.002s 1.56e-06s 1000 2 Elemwise{sub,no_inplace}(<TensorType(float64, vector)>, <TensorType(float64, vector)>)\n",
" 0.1% 99.5% 0.001s 1.08e-06s 1000 16 Sum{acc_dtype=float64}(Elemwise{Log}[(0, 0)].0)\n",
" 0.1% 99.6% 0.001s 9.19e-07s 1000 4 Reshape{1}(s_log_, TensorConstant{(1,) of -1})\n",
" 0.1% 99.6% 0.001s 8.66e-07s 1000 0 Elemwise{exp,no_inplace}(s_log_)\n",
" 0.1% 99.7% 0.001s 7.81e-07s 1000 13 InplaceDimShuffle{x,0}(Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}.0)\n",
" 0.1% 99.8% 0.001s 7.18e-07s 1000 3 AllocEmpty{dtype='float64'}(TensorConstant{1})\n",
" 0.1% 99.8% 0.001s 5.41e-07s 1000 7 InplaceDimShuffle{x,x}(s)\n",
" 0.0% 99.9% 0.000s 4.39e-07s 1000 17 InplaceDimShuffle{}(CGemv{inplace}.0)\n",
" 0.0% 99.9% 0.000s 3.94e-07s 1000 8 Reshape{0}(Rebroadcast{1}.0, TensorConstant{[]})\n",
" 0.0% 99.9% 0.000s 3.67e-07s 1000 1 Shape_i{0}(<TensorType(float64, matrix)>)\n",
" 0.0% 100.0% 0.000s 3.19e-07s 1000 5 ViewOp(Elemwise{exp,no_inplace}.0)\n",
" 0.0% 100.0% 0.000s 2.18e-07s 1000 6 Rebroadcast{1}(Reshape{1}.0)\n",
" ... (remaining 0 Apply instances account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Here are tips to potentially make your code run faster\n",
" (if you think of new ones, suggest them on the mailing list).\n",
" Test them first, as they are not guaranteed to always provide a speedup.\n",
" - Try the Theano flag floatX=float32\n",
" - Try installing amdlibm and set the Theano flag lib.amdlibm=True. This speeds up only some Elemwise operation.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"MvNormal with cov grad of log posterior\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Function profiling\n",
"==================\n",
" Message: /home/bill/pymc3/pymc3/model.py:600\n",
" Time in 1000 calls to Function.__call__: 8.509416e+00s\n",
" Time in Function.fn.__call__: 8.454453e+00s (99.354%)\n",
" Time in thunks: 8.187661e+00s (96.219%)\n",
" Total compile time: 3.527660e-01s\n",
" Number of Apply nodes: 55\n",
" Theano Optimizer time: 2.537792e-01s\n",
" Theano validate time: 4.625559e-03s\n",
" Theano Linker time (includes C, CUDA code generation/compiling): 5.474687e-02s\n",
" Import time 1.768398e-02s\n",
" Node make_thunk time 5.185175e-02s\n",
" Node Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 * i0 * i0), i1) + (i3 * i0) + (i4 * i0))}}[(0, 0)](Elemwise{exp}.0, TensorConstant{0}, TensorConstant{-0.0001}, Reshape{0}.0, Sum{acc_dtype=float64}.0) time 2.070189e-03s\n",
" Node Elemwise{Composite{((((i0 + i1) * i2) - i3) * i4)}}[(0, 2)](Solve{A_structure='upper_triangular', lower=False, overwrite_A=False, overwrite_b=False}.0, InplaceDimShuffle{1,0}.0, Tri{dtype='float64'}.0, Diag.0, <TensorType(float64, matrix)>) time 1.996040e-03s\n",
" Node Elemwise{Mul}[(0, 1)](TensorConstant{(1,) of 0.5}, ExtractDiag{view=False}.0) time 1.859903e-03s\n",
" Node Shape_i{1}(<TensorType(float64, matrix)>) time 1.859188e-03s\n",
" Node Elemwise{Composite{((i0 * i1) + (i0 * i1))}}(TensorConstant{(1,) of -0.5}, Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}.0) time 1.821280e-03s\n",
"\n",
"Time in all call to theano.grad() 5.143404e-01s\n",
"Time since theano import 13.994s\n",
"Class\n",
"---\n",
"<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>\n",
" 29.5% 29.5% 2.418s 6.04e-04s Py 4000 4 theano.tensor.slinalg.Solve\n",
" 24.0% 53.5% 1.966s 9.83e-04s C 2000 2 theano.tensor.blas.Dot22\n",
" 15.1% 68.6% 1.237s 1.24e-03s Py 1000 1 theano.tensor.slinalg.Cholesky\n",
" 11.7% 80.3% 0.955s 3.18e-04s Py 3000 3 theano.tensor.basic.Tri\n",
" 11.5% 91.8% 0.944s 8.58e-05s C 11000 11 theano.tensor.elemwise.Elemwise\n",
" 3.0% 94.8% 0.243s 2.43e-04s C 1000 1 theano.tensor.blas.Dot22Scalar\n",
" 1.4% 96.2% 0.117s 5.87e-05s Py 2000 2 theano.tensor.basic.Diag\n",
" 1.0% 97.3% 0.084s 8.35e-05s C 1000 1 theano.tensor.elemwise.Sum\n",
" 0.8% 98.0% 0.062s 2.05e-05s Py 3000 3 theano.tensor.nlinalg.ExtractDiag\n",
" 0.7% 98.7% 0.060s 5.96e-05s Py 1000 1 theano.tensor.nlinalg.AllocDiag\n",
" 0.5% 99.2% 0.041s 4.14e-05s C 1000 1 theano.tensor.subtensor.IncSubtensor\n",
" 0.5% 99.7% 0.038s 3.77e-05s C 1000 1 theano.tensor.basic.Alloc\n",
" 0.1% 99.8% 0.011s 1.12e-06s C 10000 10 theano.tensor.elemwise.DimShuffle\n",
" 0.1% 99.9% 0.005s 1.10e-06s C 5000 5 theano.tensor.basic.Reshape\n",
" 0.0% 99.9% 0.002s 9.65e-07s C 2000 2 theano.compile.ops.Shape_i\n",
" 0.0% 100.0% 0.002s 5.62e-07s C 3000 3 theano.compile.ops.Rebroadcast\n",
" 0.0% 100.0% 0.001s 1.33e-06s C 1000 1 theano.tensor.opt.MakeVector\n",
" 0.0% 100.0% 0.001s 1.06e-06s C 1000 1 theano.tensor.elemwise.CAReduce\n",
" 0.0% 100.0% 0.001s 9.07e-07s C 1000 1 theano.tensor.basic.ScalarFromTensor\n",
" 0.0% 100.0% 0.001s 7.26e-07s C 1000 1 theano.compile.ops.ViewOp\n",
" ... (remaining 0 Classes account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Ops\n",
"---\n",
"<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>\n",
" 26.9% 26.9% 2.205s 7.35e-04s Py 3000 3 Solve{A_structure='upper_triangular', lower=False, overwrite_A=False, overwrite_b=False}\n",
" 24.0% 50.9% 1.966s 9.83e-04s C 2000 2 Dot22\n",
" 15.1% 66.0% 1.237s 1.24e-03s Py 1000 1 Cholesky{lower=True, destructive=False}\n",
" 11.7% 77.7% 0.955s 3.18e-04s Py 3000 3 Tri{dtype='float64'}\n",
" 4.1% 81.8% 0.335s 3.35e-04s C 1000 1 Elemwise{Composite{((((i0 + i1) * i2) - i3) * i4)}}[(0, 2)]\n",
" 3.8% 85.6% 0.309s 3.09e-04s C 1000 1 Elemwise{Composite{((i0 * i1) - i2)}}[(0, 0)]\n",
" 3.0% 88.5% 0.243s 2.43e-04s C 1000 1 Dot22Scalar\n",
" 2.6% 91.1% 0.213s 2.13e-04s Py 1000 1 Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}\n",
" 2.1% 93.3% 0.173s 1.73e-04s C 1000 1 Elemwise{mul,no_inplace}\n",
" 1.4% 94.7% 0.117s 5.87e-05s Py 2000 2 Diag\n",
" 1.3% 96.0% 0.108s 1.08e-04s C 1000 1 Elemwise{Composite{(i0 + (i1 * i2))}}[(0, 0)]\n",
" 1.0% 97.0% 0.084s 8.35e-05s C 1000 1 Sum{acc_dtype=float64}\n",
" 0.8% 97.8% 0.062s 2.05e-05s Py 3000 3 ExtractDiag{view=False}\n",
" 0.7% 98.5% 0.060s 5.96e-05s Py 1000 1 AllocDiag\n",
" 0.5% 99.0% 0.041s 4.14e-05s C 1000 1 IncSubtensor{InplaceSet;:int64:, :int64:}\n",
" 0.5% 99.5% 0.038s 3.77e-05s C 1000 1 Alloc\n",
" 0.1% 99.6% 0.007s 1.11e-06s C 6000 6 InplaceDimShuffle{1,0}\n",
" 0.1% 99.6% 0.006s 6.34e-06s C 1000 1 Elemwise{TrueDiv}[(0, 1)]\n",
" 0.0% 99.7% 0.003s 1.57e-06s C 2000 2 Reshape{1}\n",
" 0.0% 99.7% 0.003s 3.14e-06s C 1000 1 Elemwise{sub,no_inplace}\n",
" ... (remaining 18 Ops account for 0.29%(0.02s) of the runtime)\n",
"\n",
"Apply\n",
"------\n",
"<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>\n",
" 15.1% 15.1% 1.237s 1.24e-03s 1000 19 Cholesky{lower=True, destructive=False}(Elemwise{mul,no_inplace}.0)\n",
" 13.0% 28.1% 1.065s 1.07e-03s 1000 44 Solve{A_structure='upper_triangular', lower=False, overwrite_A=False, overwrite_b=False}(InplaceDimShuffle{1,0}.0, Elemwise{Composite{((i0 * i1) - i2)}}[(0, 0)].0)\n",
" 12.3% 40.4% 1.006s 1.01e-03s 1000 36 Dot22(InplaceDimShuffle{1,0}.0, Elemwise{Composite{(i0 + (i1 * i2))}}[(0, 0)].0)\n",
" 11.7% 52.1% 0.960s 9.60e-04s 1000 39 Dot22(InplaceDimShuffle{1,0}.0, Cholesky{lower=True, destructive=False}.0)\n",
" 11.1% 63.3% 0.912s 9.12e-04s 1000 46 Solve{A_structure='upper_triangular', lower=False, overwrite_A=False, overwrite_b=False}(InplaceDimShuffle{1,0}.0, InplaceDimShuffle{1,0}.0)\n",
" 5.3% 68.6% 0.434s 4.34e-04s 1000 8 Tri{dtype='float64'}(Shape_i{0}.0, Shape_i{1}.0, TensorConstant{0})\n",
" 4.1% 72.6% 0.335s 3.35e-04s 1000 50 Elemwise{Composite{((((i0 + i1) * i2) - i3) * i4)}}[(0, 2)](Solve{A_structure='upper_triangular', lower=False, overwrite_A=False, overwrite_b=False}.0, InplaceDimShuffle{1,0}.0, Tri{dtype='float64'}.0, Diag.0, <TensorType(float64, matrix)>)\n",
" 3.8% 76.4% 0.309s 3.09e-04s 1000 43 Elemwise{Composite{((i0 * i1) - i2)}}[(0, 0)](Dot22.0, InplaceDimShuffle{1,0}.0, InplaceDimShuffle{1,0}.0)\n",
" 3.3% 79.7% 0.267s 2.67e-04s 1000 7 Tri{dtype='float64'}(Shape_i{1}.0, Shape_i{1}.0, TensorConstant{0})\n",
" 3.1% 82.8% 0.254s 2.54e-04s 1000 11 Tri{dtype='float64'}(Shape_i{0}.0, Shape_i{0}.0, TensorConstant{0})\n",
" 3.0% 85.8% 0.243s 2.43e-04s 1000 34 Dot22Scalar(InplaceDimShuffle{0,x}.0, InplaceDimShuffle{x,0}.0, TensorConstant{-1.0})\n",
" 2.8% 88.5% 0.228s 2.28e-04s 1000 29 Solve{A_structure='upper_triangular', lower=False, overwrite_A=False, overwrite_b=False}(InplaceDimShuffle{1,0}.0, Elemwise{Composite{((i0 * i1) + (i0 * i1))}}.0)\n",
" 2.6% 91.1% 0.213s 2.13e-04s 1000 21 Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}(Cholesky{lower=True, destructive=False}.0, Elemwise{sub,no_inplace}.0)\n",
" 2.1% 93.3% 0.173s 1.73e-04s 1000 16 Elemwise{mul,no_inplace}(InplaceDimShuffle{x,x}.0, <TensorType(float64, matrix)>)\n",
" 1.3% 94.6% 0.108s 1.08e-04s 1000 35 Elemwise{Composite{(i0 + (i1 * i2))}}[(0, 0)](IncSubtensor{InplaceSet;:int64:, :int64:}.0, Dot22Scalar.0, Tri{dtype='float64'}.0)\n",
" 1.0% 95.6% 0.084s 8.35e-05s 1000 51 Sum{acc_dtype=float64}(Elemwise{Composite{((((i0 + i1) * i2) - i3) * i4)}}[(0, 2)].0)\n",
" 0.8% 96.4% 0.064s 6.38e-05s 1000 41 Diag(Elemwise{Mul}[(0, 1)].0)\n",
" 0.7% 97.1% 0.060s 5.96e-05s 1000 30 AllocDiag(Elemwise{TrueDiv}[(0, 1)].0)\n",
" 0.7% 97.8% 0.054s 5.35e-05s 1000 49 Diag(ExtractDiag{view=False}.0)\n",
" 0.5% 98.3% 0.041s 4.14e-05s 1000 33 IncSubtensor{InplaceSet;:int64:, :int64:}(Alloc.0, AllocDiag.0, ScalarFromTensor.0, ScalarFromTensor.0)\n",
" ... (remaining 35 Apply instances account for 1.74%(0.14s) of the runtime)\n",
"\n",
"Here are tips to potentially make your code run faster\n",
" (if you think of new ones, suggest them on the mailing list).\n",
" Test them first, as they are not guaranteed to always provide a speedup.\n",
" - Try the Theano flag floatX=float32\n",
" - Try installing amdlibm and set the Theano flag lib.amdlibm=True. This speeds up only some Elemwise operation.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
"\n",
"\n",
"MvNormal with tau log posterior\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Function profiling\n",
"==================\n",
" Message: /home/bill/pymc3/pymc3/model.py:600\n",
" Time in 1000 calls to Function.__call__: 8.873688e+00s\n",
" Time in Function.fn.__call__: 8.838151e+00s (99.600%)\n",
" Time in thunks: 8.824074e+00s (99.441%)\n",
" Total compile time: 1.624010e-01s\n",
" Number of Apply nodes: 17\n",
" Theano Optimizer time: 1.092184e-01s\n",
" Theano validate time: 1.314640e-03s\n",
" Theano Linker time (includes C, CUDA code generation/compiling): 1.570606e-02s\n",
" Import time 3.223419e-03s\n",
" Node make_thunk time 1.455784e-02s\n",
" Node Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 + (i3 * sqr(i0))), i4) + log(Abs(exp(i5))) + (i6 * (((i7 * i8) - i9) + i10)))}}[(0, 0)](Elemwise{exp,no_inplace}.0, TensorConstant{0}, TensorConstant{-4.830961538632819}, TensorConstant{-5e-05}, TensorConstant{-inf}, Reshape{0}.0, TensorConstant{-0.5}, TensorConstant{1.8378770664093453}, Shape_i{0}.0, LogDet.0, Sum{acc_dtype=float64}.0) time 2.853870e-03s\n",
" Node AllocEmpty{dtype='float64'}(Shape_i{1}.0) time 1.860619e-03s\n",
" Node Elemwise{Mul}[(0, 0)](CGemv{inplace}.0, Elemwise{sub,no_inplace}.0) time 1.777649e-03s\n",
" Node InplaceDimShuffle{x,x}(s) time 1.252174e-03s\n",
" Node CGemv{inplace}(AllocEmpty{dtype='float64'}.0, TensorConstant{1.0}, InplaceDimShuffle{1,0}.0, Elemwise{sub,no_inplace}.0, TensorConstant{0.0}) time 7.073879e-04s\n",
"\n",
"Time in all call to theano.grad() 5.143404e-01s\n",
"Time since theano import 23.351s\n",
"Class\n",
"---\n",
"<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>\n",
" 98.2% 98.2% 8.661s 8.66e-03s Py 1000 1 pymc3.math.LogDet\n",
" 1.0% 99.2% 0.090s 9.00e-05s C 1000 1 theano.tensor.blas_c.CGemv\n",
" 0.8% 99.9% 0.067s 1.33e-05s C 5000 5 theano.tensor.elemwise.Elemwise\n",
" 0.0% 99.9% 0.001s 6.87e-07s C 2000 2 theano.tensor.basic.Reshape\n",
" 0.0% 100.0% 0.001s 6.30e-07s C 2000 2 theano.tensor.elemwise.DimShuffle\n",
" 0.0% 100.0% 0.001s 1.06e-06s C 1000 1 theano.tensor.basic.AllocEmpty\n",
" 0.0% 100.0% 0.001s 9.56e-07s C 1000 1 theano.tensor.elemwise.Sum\n",
" 0.0% 100.0% 0.001s 4.74e-07s C 2000 2 theano.compile.ops.Shape_i\n",
" 0.0% 100.0% 0.000s 2.43e-07s C 1000 1 theano.compile.ops.ViewOp\n",
" 0.0% 100.0% 0.000s 1.99e-07s C 1000 1 theano.compile.ops.Rebroadcast\n",
" ... (remaining 0 Classes account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Ops\n",
"---\n",
"<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>\n",
" 98.2% 98.2% 8.661s 8.66e-03s Py 1000 1 LogDet\n",
" 1.0% 99.2% 0.090s 9.00e-05s C 1000 1 CGemv{inplace}\n",
" 0.7% 99.9% 0.061s 6.12e-05s C 1000 1 Elemwise{mul,no_inplace}\n",
" 0.0% 99.9% 0.002s 2.09e-06s C 1000 1 Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 + (i3 * sqr(i0))), i4) + log(Abs(exp(i5))) + (i6 * (((i7 * i8) - i9) + i10)))}}[(0, 0)]\n",
" 0.0% 99.9% 0.001s 1.30e-06s C 1000 1 Elemwise{sub,no_inplace}\n",
" 0.0% 99.9% 0.001s 1.06e-06s C 1000 1 AllocEmpty{dtype='float64'}\n",
" 0.0% 99.9% 0.001s 1.03e-06s C 1000 1 Elemwise{exp,no_inplace}\n",
" 0.0% 99.9% 0.001s 9.85e-07s C 1000 1 Elemwise{Mul}[(0, 0)]\n",
" 0.0% 100.0% 0.001s 9.56e-07s C 1000 1 Sum{acc_dtype=float64}\n",
" 0.0% 100.0% 0.001s 9.45e-07s C 1000 1 Reshape{1}\n",
" 0.0% 100.0% 0.001s 6.53e-07s C 1000 1 InplaceDimShuffle{1,0}\n",
" 0.0% 100.0% 0.001s 6.07e-07s C 1000 1 InplaceDimShuffle{x,x}\n",
" 0.0% 100.0% 0.001s 5.56e-07s C 1000 1 Shape_i{1}\n",
" 0.0% 100.0% 0.000s 4.29e-07s C 1000 1 Reshape{0}\n",
" 0.0% 100.0% 0.000s 3.91e-07s C 1000 1 Shape_i{0}\n",
" 0.0% 100.0% 0.000s 2.43e-07s C 1000 1 ViewOp\n",
" 0.0% 100.0% 0.000s 1.99e-07s C 1000 1 Rebroadcast{1}\n",
" ... (remaining 0 Ops account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Apply\n",
"------\n",
"<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>\n",
" 98.2% 98.2% 8.661s 8.66e-03s 1000 12 LogDet(Elemwise{mul,no_inplace}.0)\n",
" 1.0% 99.2% 0.090s 9.00e-05s 1000 13 CGemv{inplace}(AllocEmpty{dtype='float64'}.0, TensorConstant{1.0}, InplaceDimShuffle{1,0}.0, Elemwise{sub,no_inplace}.0, TensorConstant{0.0})\n",
" 0.7% 99.9% 0.061s 6.12e-05s 1000 10 Elemwise{mul,no_inplace}(InplaceDimShuffle{x,x}.0, <TensorType(float64, matrix)>)\n",
" 0.0% 99.9% 0.002s 2.09e-06s 1000 16 Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 + (i3 * sqr(i0))), i4) + log(Abs(exp(i5))) + (i6 * (((i7 * i8) - i9) + i10)))}}[(0, 0)](Elemwise{exp,no_inplace}.0, TensorConstant{0}, TensorConstant{-4.830961538632819}, TensorConstant{-5e-05}, TensorConstant{-inf}, Reshape{0}.0, TensorConstant{-0.5}, TensorConstant{1.8378770664093453}, Shape_i{0}.0, LogDet.0, Sum{acc_dtype=float64}.0)\n",
" 0.0% 99.9% 0.001s 1.30e-06s 1000 3 Elemwise{sub,no_inplace}(<TensorType(float64, vector)>, <TensorType(float64, vector)>)\n",
" 0.0% 99.9% 0.001s 1.06e-06s 1000 6 AllocEmpty{dtype='float64'}(Shape_i{1}.0)\n",
" 0.0% 99.9% 0.001s 1.03e-06s 1000 0 Elemwise{exp,no_inplace}(s_log_)\n",
" 0.0% 99.9% 0.001s 9.85e-07s 1000 14 Elemwise{Mul}[(0, 0)](CGemv{inplace}.0, Elemwise{sub,no_inplace}.0)\n",
" 0.0% 100.0% 0.001s 9.56e-07s 1000 15 Sum{acc_dtype=float64}(Elemwise{Mul}[(0, 0)].0)\n",
" 0.0% 100.0% 0.001s 9.45e-07s 1000 4 Reshape{1}(s_log_, TensorConstant{(1,) of -1})\n",
" 0.0% 100.0% 0.001s 6.53e-07s 1000 11 InplaceDimShuffle{1,0}(Elemwise{mul,no_inplace}.0)\n",
" 0.0% 100.0% 0.001s 6.07e-07s 1000 8 InplaceDimShuffle{x,x}(s)\n",
" 0.0% 100.0% 0.001s 5.56e-07s 1000 1 Shape_i{1}(<TensorType(float64, matrix)>)\n",
" 0.0% 100.0% 0.000s 4.29e-07s 1000 9 Reshape{0}(Rebroadcast{1}.0, TensorConstant{[]})\n",
" 0.0% 100.0% 0.000s 3.91e-07s 1000 2 Shape_i{0}(<TensorType(float64, matrix)>)\n",
" 0.0% 100.0% 0.000s 2.43e-07s 1000 5 ViewOp(Elemwise{exp,no_inplace}.0)\n",
" 0.0% 100.0% 0.000s 1.99e-07s 1000 7 Rebroadcast{1}(Reshape{1}.0)\n",
" ... (remaining 0 Apply instances account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Here are tips to potentially make your code run faster\n",
" (if you think of new ones, suggest them on the mailing list).\n",
" Test them first, as they are not guaranteed to always provide a speedup.\n",
" - Try the Theano flag floatX=float32\n",
" - Try installing amdlibm and set the Theano flag lib.amdlibm=True. This speeds up only some Elemwise operation.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"MvNormal with tau grad of log posterior\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Function profiling\n",
"==================\n",
" Message: /home/bill/pymc3/pymc3/model.py:600\n",
" Time in 1000 calls to Function.__call__: 2.421551e+00s\n",
" Time in Function.fn.__call__: 2.381437e+00s (98.343%)\n",
" Time in thunks: 2.362702e+00s (97.570%)\n",
" Total compile time: 2.096241e-01s\n",
" Number of Apply nodes: 24\n",
" Theano Optimizer time: 1.514192e-01s\n",
" Theano validate time: 1.968861e-03s\n",
" Theano Linker time (includes C, CUDA code generation/compiling): 1.817369e-02s\n",
" Import time 3.171444e-03s\n",
" Node make_thunk time 1.684165e-02s\n",
" Node Elemwise{mul,no_inplace}(TensorConstant{(1, 1) of -0.5}, InplaceDimShuffle{x,0}.0) time 2.072096e-03s\n",
" Node Gemm{inplace}(InplaceDimShuffle{1,0}.0, TensorConstant{1.0}, InplaceDimShuffle{0,x}.0, Elemwise{mul,no_inplace}.0, TensorConstant{0.5}) time 2.008915e-03s\n",
" Node Elemwise{Mul}[(0, 0)](Gemm{inplace}.0, <TensorType(float64, matrix)>) time 1.682043e-03s\n",
" Node Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 * i0 * i0), i1) + (i3 * i0) + (i4 * i0))}}[(0, 0)](Elemwise{exp}.0, TensorConstant{0}, TensorConstant{-0.0001}, Reshape{0}.0, Sum{acc_dtype=float64}.0) time 9.992123e-04s\n",
" Node Reshape{0}(Rebroadcast{1}.0, TensorConstant{[]}) time 6.287098e-04s\n",
"\n",
"Time in all call to theano.grad() 7.757266e-01s\n",
"Time since theano import 26.471s\n",
"Class\n",
"---\n",
"<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>\n",
" 83.0% 83.0% 1.961s 1.96e-03s Py 1000 1 theano.tensor.nlinalg.MatrixInverse\n",
" 10.6% 93.6% 0.251s 3.58e-05s C 7000 7 theano.tensor.elemwise.Elemwise\n",
" 4.2% 97.8% 0.100s 1.00e-04s C 1000 1 theano.tensor.elemwise.Sum\n",
" 1.8% 99.6% 0.043s 4.29e-05s C 1000 1 theano.tensor.blas.Gemm\n",
" 0.1% 99.8% 0.003s 6.97e-07s C 5000 5 theano.tensor.elemwise.DimShuffle\n",
" 0.1% 99.9% 0.003s 6.84e-07s C 5000 5 theano.tensor.basic.Reshape\n",
" 0.0% 100.0% 0.001s 3.64e-07s C 3000 3 theano.compile.ops.Rebroadcast\n",
" 0.0% 100.0% 0.000s 3.98e-07s C 1000 1 theano.compile.ops.ViewOp\n",
" ... (remaining 0 Classes account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Ops\n",
"---\n",
"<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>\n",
" 83.0% 83.0% 1.961s 1.96e-03s Py 1000 1 MatrixInverse\n",
" 7.7% 90.7% 0.182s 1.82e-04s C 1000 1 Elemwise{Mul}[(0, 0)]\n",
" 4.2% 94.9% 0.100s 1.00e-04s C 1000 1 Sum{acc_dtype=float64}\n",
" 2.7% 97.6% 0.063s 3.16e-05s C 2000 2 Elemwise{mul,no_inplace}\n",
" 1.8% 99.4% 0.043s 4.29e-05s C 1000 1 Gemm{inplace}\n",
" 0.1% 99.5% 0.002s 2.08e-06s C 1000 1 Elemwise{sub,no_inplace}\n",
" 0.1% 99.6% 0.002s 9.90e-07s C 2000 2 Reshape{1}\n",
" 0.1% 99.6% 0.001s 1.50e-06s C 1000 1 Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 * i0 * i0), i1) + (i3 * i0) + (i4 * i0))}}[(0, 0)]\n",
" 0.1% 99.7% 0.001s 4.80e-07s C 3000 3 Reshape{0}\n",
" 0.1% 99.8% 0.001s 1.28e-06s C 1000 1 Elemwise{exp}\n",
" 0.0% 99.8% 0.001s 1.01e-06s C 1000 1 InplaceDimShuffle{1,0}\n",
" 0.0% 99.8% 0.001s 4.28e-07s C 2000 2 Rebroadcast{1}\n",
" 0.0% 99.9% 0.001s 7.70e-07s C 1000 1 InplaceDimShuffle{0,x}\n",
" 0.0% 99.9% 0.001s 7.09e-07s C 1000 1 InplaceDimShuffle{x,x}\n",
" 0.0% 99.9% 0.001s 6.81e-07s C 1000 1 Elemwise{Composite{(sgn(exp(i0)) / Abs(exp(i0)))}}\n",
" 0.0% 100.0% 0.001s 5.50e-07s C 1000 1 InplaceDimShuffle{x,0}\n",
" 0.0% 100.0% 0.000s 4.43e-07s C 1000 1 InplaceDimShuffle{x}\n",
" 0.0% 100.0% 0.000s 3.98e-07s C 1000 1 ViewOp\n",
" 0.0% 100.0% 0.000s 2.35e-07s C 1000 1 Rebroadcast{0}\n",
" ... (remaining 0 Ops account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Apply\n",
"------\n",
"<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>\n",
" 83.0% 83.0% 1.961s 1.96e-03s 1000 12 MatrixInverse(Elemwise{mul,no_inplace}.0)\n",
" 7.7% 90.7% 0.182s 1.82e-04s 1000 18 Elemwise{Mul}[(0, 0)](Gemm{inplace}.0, <TensorType(float64, matrix)>)\n",
" 4.2% 94.9% 0.100s 1.00e-04s 1000 20 Sum{acc_dtype=float64}(Elemwise{Mul}[(0, 0)].0)\n",
" 2.6% 97.5% 0.061s 6.14e-05s 1000 10 Elemwise{mul,no_inplace}(InplaceDimShuffle{x,x}.0, <TensorType(float64, matrix)>)\n",
" 1.8% 99.3% 0.043s 4.29e-05s 1000 16 Gemm{inplace}(InplaceDimShuffle{1,0}.0, TensorConstant{1.0}, InplaceDimShuffle{0,x}.0, Elemwise{mul,no_inplace}.0, TensorConstant{0.5})\n",
" 0.1% 99.4% 0.002s 2.08e-06s 1000 2 Elemwise{sub,no_inplace}(<TensorType(float64, vector)>, <TensorType(float64, vector)>)\n",
" 0.1% 99.5% 0.002s 1.79e-06s 1000 9 Elemwise{mul,no_inplace}(TensorConstant{(1, 1) of -0.5}, InplaceDimShuffle{x,0}.0)\n",
" 0.1% 99.6% 0.001s 1.50e-06s 1000 21 Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 * i0 * i0), i1) + (i3 * i0) + (i4 * i0))}}[(0, 0)](Elemwise{exp}.0, TensorConstant{0}, TensorConstant{-0.0001}, Reshape{0}.0, Sum{acc_dtype=float64}.0)\n",
" 0.1% 99.6% 0.001s 1.28e-06s 1000 0 Elemwise{exp}(s_log_)\n",
" 0.1% 99.7% 0.001s 1.25e-06s 1000 22 Reshape{1}(Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 * i0 * i0), i1) + (i3 * i0) + (i4 * i0))}}[(0, 0)].0, TensorConstant{(1,) of -1})\n",
" 0.0% 99.7% 0.001s 1.01e-06s 1000 14 InplaceDimShuffle{1,0}(MatrixInverse.0)\n",
" 0.0% 99.7% 0.001s 7.70e-07s 1000 6 InplaceDimShuffle{0,x}(Elemwise{sub,no_inplace}.0)\n",
" 0.0% 99.8% 0.001s 7.51e-07s 1000 13 Reshape{0}(Elemwise{Composite{(sgn(exp(i0)) / Abs(exp(i0)))}}.0, TensorConstant{[]})\n",
" 0.0% 99.8% 0.001s 7.26e-07s 1000 1 Reshape{1}(s_log_, TensorConstant{(1,) of -1})\n",
" 0.0% 99.8% 0.001s 7.09e-07s 1000 7 InplaceDimShuffle{x,x}(s)\n",
" 0.0% 99.9% 0.001s 6.81e-07s 1000 11 Elemwise{Composite{(sgn(exp(i0)) / Abs(exp(i0)))}}(Reshape{0}.0)\n",
" 0.0% 99.9% 0.001s 5.84e-07s 1000 4 Rebroadcast{1}(Reshape{1}.0)\n",
" 0.0% 99.9% 0.001s 5.50e-07s 1000 5 InplaceDimShuffle{x,0}(Elemwise{sub,no_inplace}.0)\n",
" 0.0% 99.9% 0.000s 4.77e-07s 1000 8 Reshape{0}(Rebroadcast{1}.0, TensorConstant{[]})\n",
" 0.0% 100.0% 0.000s 4.43e-07s 1000 15 InplaceDimShuffle{x}(Reshape{0}.0)\n",
" ... (remaining 4 Apply instances account for 0.05%(0.00s) of the runtime)\n",
"\n",
"Here are tips to potentially make your code run faster\n",
" (if you think of new ones, suggest them on the mailing list).\n",
" Test them first, as they are not guaranteed to always provide a speedup.\n",
" - Try the Theano flag floatX=float32\n",
" - Try installing amdlibm and set the Theano flag lib.amdlibm=True. This speeds up only some Elemwise operation.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
"\n",
"\n",
"MvNormal with chol log posterior\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Function profiling\n",
"==================\n",
" Message: /home/bill/pymc3/pymc3/model.py:600\n",
" Time in 1000 calls to Function.__call__: 3.402998e-01s\n",
" Time in Function.fn.__call__: 3.153903e-01s (92.680%)\n",
" Time in thunks: 3.047023e-01s (89.539%)\n",
" Total compile time: 1.776102e-01s\n",
" Number of Apply nodes: 18\n",
" Theano Optimizer time: 1.267796e-01s\n",
" Theano validate time: 1.578093e-03s\n",
" Theano Linker time (includes C, CUDA code generation/compiling): 1.334190e-02s\n",
" Import time 0.000000e+00s\n",
" Node make_thunk time 1.204419e-02s\n",
" Node Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 + (i3 * sqr(i0))), i4) + log(Abs(exp(i5))) + (i6 * i7) + (i8 * i9) + (i10 * i11))}}[(0, 0)](Elemwise{exp,no_inplace}.0, TensorConstant{0}, TensorConstant{-4.830961538632819}, TensorConstant{-5e-05}, TensorConstant{-inf}, Reshape{0}.0, TensorConstant{1.8378770664093453}, Shape_i{0}.0, TensorConstant{2.0}, Sum{acc_dtype=float64}.0, TensorConstant{-0.5}, InplaceDimShuffle{}.0) time 1.914740e-03s\n",
" Node CGemv{inplace}(AllocEmpty{dtype='float64'}.0, TensorConstant{1.0}, InplaceDimShuffle{x,0}.0, Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}.0, TensorConstant{0.0}) time 8.981228e-04s\n",
" Node Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}(Elemwise{mul,no_inplace}.0, Elemwise{sub,no_inplace}.0) time 7.765293e-04s\n",
" Node ViewOp(Elemwise{exp,no_inplace}.0) time 6.928444e-04s\n",
" Node Elemwise{exp,no_inplace}(s_log_) time 6.840229e-04s\n",
"\n",
"Time in all call to theano.grad() 7.757266e-01s\n",
"Time since theano import 27.199s\n",
"Class\n",
"---\n",
"<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>\n",
" 72.1% 72.1% 0.220s 2.20e-04s Py 1000 1 theano.tensor.slinalg.Solve\n",
" 22.4% 94.4% 0.068s 1.36e-05s C 5000 5 theano.tensor.elemwise.Elemwise\n",
" 3.7% 98.1% 0.011s 1.13e-05s Py 1000 1 theano.tensor.nlinalg.ExtractDiag\n",
" 0.5% 98.6% 0.002s 5.08e-07s C 3000 3 theano.tensor.elemwise.DimShuffle\n",
" 0.4% 99.0% 0.001s 1.18e-06s C 1000 1 theano.tensor.blas_c.CGemv\n",
" 0.3% 99.4% 0.001s 4.85e-07s C 2000 2 theano.tensor.basic.Reshape\n",
" 0.2% 99.6% 0.001s 7.59e-07s C 1000 1 theano.tensor.elemwise.Sum\n",
" 0.1% 99.8% 0.000s 4.56e-07s C 1000 1 theano.tensor.basic.AllocEmpty\n",
" 0.1% 99.9% 0.000s 3.41e-07s C 1000 1 theano.compile.ops.Shape_i\n",
" 0.1% 99.9% 0.000s 2.18e-07s C 1000 1 theano.compile.ops.Rebroadcast\n",
" 0.1% 100.0% 0.000s 1.86e-07s C 1000 1 theano.compile.ops.ViewOp\n",
" ... (remaining 0 Classes account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Ops\n",
"---\n",
"<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>\n",
" 72.1% 72.1% 0.220s 2.20e-04s Py 1000 1 Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}\n",
" 18.7% 90.8% 0.057s 5.71e-05s C 1000 1 Elemwise{mul,no_inplace}\n",
" 3.7% 94.5% 0.011s 1.13e-05s Py 1000 1 ExtractDiag{view=False}\n",
" 2.6% 97.1% 0.008s 7.89e-06s C 1000 1 Elemwise{Log}[(0, 0)]\n",
" 0.4% 97.6% 0.001s 1.33e-06s C 1000 1 Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 + (i3 * sqr(i0))), i4) + log(Abs(exp(i5))) + (i6 * i7) + (i8 * i9) + (i10 * i11))}}[(0, 0)]\n",
" 0.4% 98.0% 0.001s 1.18e-06s C 1000 1 CGemv{inplace}\n",
" 0.3% 98.3% 0.001s 1.02e-06s C 1000 1 Elemwise{sub,no_inplace}\n",
" 0.2% 98.5% 0.001s 7.59e-07s C 1000 1 Sum{acc_dtype=float64}\n",
" 0.2% 98.8% 0.001s 7.28e-07s C 1000 1 Elemwise{exp,no_inplace}\n",
" 0.2% 99.0% 0.001s 6.78e-07s C 1000 1 Reshape{1}\n",
" 0.2% 99.2% 0.001s 6.56e-07s C 1000 1 InplaceDimShuffle{x,0}\n",
" 0.2% 99.4% 0.001s 5.45e-07s C 1000 1 InplaceDimShuffle{x,x}\n",
" 0.1% 99.6% 0.000s 4.56e-07s C 1000 1 AllocEmpty{dtype='float64'}\n",
" 0.1% 99.7% 0.000s 3.41e-07s C 1000 1 Shape_i{0}\n",
" 0.1% 99.8% 0.000s 3.22e-07s C 1000 1 InplaceDimShuffle{}\n",
" 0.1% 99.9% 0.000s 2.92e-07s C 1000 1 Reshape{0}\n",
" 0.1% 99.9% 0.000s 2.18e-07s C 1000 1 Rebroadcast{1}\n",
" 0.1% 100.0% 0.000s 1.86e-07s C 1000 1 ViewOp\n",
" ... (remaining 0 Ops account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Apply\n",
"------\n",
"<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>\n",
" 72.1% 72.1% 0.220s 2.20e-04s 1000 10 Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}(Elemwise{mul,no_inplace}.0, Elemwise{sub,no_inplace}.0)\n",
" 18.7% 90.8% 0.057s 5.71e-05s 1000 9 Elemwise{mul,no_inplace}(InplaceDimShuffle{x,x}.0, <TensorType(float64, matrix)>)\n",
" 3.7% 94.5% 0.011s 1.13e-05s 1000 11 ExtractDiag{view=False}(Elemwise{mul,no_inplace}.0)\n",
" 2.6% 97.1% 0.008s 7.89e-06s 1000 13 Elemwise{Log}[(0, 0)](ExtractDiag{view=False}.0)\n",
" 0.4% 97.6% 0.001s 1.33e-06s 1000 17 Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 + (i3 * sqr(i0))), i4) + log(Abs(exp(i5))) + (i6 * i7) + (i8 * i9) + (i10 * i11))}}[(0, 0)](Elemwise{exp,no_inplace}.0, TensorConstant{0}, TensorConstant{-4.830961538632819}, TensorConstant{-5e-05}, TensorConstant{-inf}, Reshape{0}.0, TensorConstant{1.8378770664093453}, Shape_i{0}.0, TensorConstant{2.0}, Sum{acc_dtype=float64}.0, TensorConstant{-0.5}, InplaceDimShuffle{}.0)\n",
" 0.4% 98.0% 0.001s 1.18e-06s 1000 14 CGemv{inplace}(AllocEmpty{dtype='float64'}.0, TensorConstant{1.0}, InplaceDimShuffle{x,0}.0, Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}.0, TensorConstant{0.0})\n",
" 0.3% 98.3% 0.001s 1.02e-06s 1000 2 Elemwise{sub,no_inplace}(<TensorType(float64, vector)>, <TensorType(float64, vector)>)\n",
" 0.2% 98.5% 0.001s 7.59e-07s 1000 15 Sum{acc_dtype=float64}(Elemwise{Log}[(0, 0)].0)\n",
" 0.2% 98.8% 0.001s 7.28e-07s 1000 0 Elemwise{exp,no_inplace}(s_log_)\n",
" 0.2% 99.0% 0.001s 6.78e-07s 1000 4 Reshape{1}(s_log_, TensorConstant{(1,) of -1})\n",
" 0.2% 99.2% 0.001s 6.56e-07s 1000 12 InplaceDimShuffle{x,0}(Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}.0)\n",
" 0.2% 99.4% 0.001s 5.45e-07s 1000 7 InplaceDimShuffle{x,x}(s)\n",
" 0.1% 99.6% 0.000s 4.56e-07s 1000 3 AllocEmpty{dtype='float64'}(TensorConstant{1})\n",
" 0.1% 99.7% 0.000s 3.41e-07s 1000 1 Shape_i{0}(<TensorType(float64, matrix)>)\n",
" 0.1% 99.8% 0.000s 3.22e-07s 1000 16 InplaceDimShuffle{}(CGemv{inplace}.0)\n",
" 0.1% 99.9% 0.000s 2.92e-07s 1000 8 Reshape{0}(Rebroadcast{1}.0, TensorConstant{[]})\n",
" 0.1% 99.9% 0.000s 2.18e-07s 1000 6 Rebroadcast{1}(Reshape{1}.0)\n",
" 0.1% 100.0% 0.000s 1.86e-07s 1000 5 ViewOp(Elemwise{exp,no_inplace}.0)\n",
" ... (remaining 0 Apply instances account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Here are tips to potentially make your code run faster\n",
" (if you think of new ones, suggest them on the mailing list).\n",
" Test them first, as they are not guaranteed to always provide a speedup.\n",
" - Try the Theano flag floatX=float32\n",
" - Try installing amdlibm and set the Theano flag lib.amdlibm=True. This speeds up only some Elemwise operation.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"MvNormal with chol grad of log posterior\n",
"\n",
"\n",
"\n",
"\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Function profiling\n",
"==================\n",
" Message: /home/bill/pymc3/pymc3/model.py:600\n",
" Time in 1000 calls to Function.__call__: 1.586063e+00s\n",
" Time in Function.fn.__call__: 1.550441e+00s (97.754%)\n",
" Time in thunks: 1.522541e+00s (95.995%)\n",
" Total compile time: 2.752922e-01s\n",
" Number of Apply nodes: 36\n",
" Theano Optimizer time: 2.067325e-01s\n",
" Theano validate time: 2.783298e-03s\n",
" Theano Linker time (includes C, CUDA code generation/compiling): 2.502918e-02s\n",
" Import time 1.274109e-03s\n",
" Node make_thunk time 2.299953e-02s\n",
" Node Elemwise{Composite{((i0 + (i1 * i2)) * i3)}}[(0, 0)](IncSubtensor{InplaceSet;:int64:, :int64:}.0, Dot22Scalar.0, Tri{dtype='float64'}.0, <TensorType(float64, matrix)>) time 2.413034e-03s\n",
" Node Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 * i0 * i0), i1) + (i3 * i0) + (i4 * i0))}}[(0, 0)](Elemwise{exp}.0, TensorConstant{0}, TensorConstant{-0.0001}, Reshape{0}.0, Sum{acc_dtype=float64}.0) time 1.215696e-03s\n",
" Node Tri{dtype='float64'}(Shape_i{0}.0, Shape_i{1}.0, TensorConstant{0}) time 1.049280e-03s\n",
" Node Solve{A_structure='upper_triangular', lower=False, overwrite_A=False, overwrite_b=False}(InplaceDimShuffle{1,0}.0, Elemwise{Composite{((i0 * i1) + (i0 * i1))}}.0) time 9.133816e-04s\n",
" Node Elemwise{Composite{((i0 * i1) + (i0 * i1))}}(TensorConstant{(1,) of -0.5}, Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}.0) time 7.977486e-04s\n",
"\n",
"Time in all call to theano.grad() 1.014499e+00s\n",
"Time since theano import 29.681s\n",
"Class\n",
"---\n",
"<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>\n",
" 28.1% 28.1% 0.427s 4.27e-04s C 1000 1 theano.tensor.blas.Dot22Scalar\n",
" 26.4% 54.5% 0.402s 2.01e-04s Py 2000 2 theano.tensor.slinalg.Solve\n",
" 15.6% 70.1% 0.237s 2.37e-04s Py 1000 1 theano.tensor.basic.Tri\n",
" 11.5% 81.6% 0.175s 2.18e-05s C 8000 8 theano.tensor.elemwise.Elemwise\n",
" 5.6% 87.1% 0.085s 8.45e-05s C 1000 1 theano.tensor.elemwise.Sum\n",
" 4.7% 91.8% 0.072s 7.21e-05s C 1000 1 theano.tensor.basic.Alloc\n",
" 3.3% 95.1% 0.050s 4.97e-05s Py 1000 1 theano.tensor.nlinalg.AllocDiag\n",
" 3.2% 98.4% 0.049s 4.93e-05s C 1000 1 theano.tensor.subtensor.IncSubtensor\n",
" 0.9% 99.3% 0.014s 1.44e-05s Py 1000 1 theano.tensor.nlinalg.ExtractDiag\n",
" 0.2% 99.5% 0.003s 6.67e-07s C 5000 5 theano.tensor.elemwise.DimShuffle\n",
" 0.2% 99.7% 0.003s 6.58e-07s C 5000 5 theano.tensor.basic.Reshape\n",
" 0.1% 99.8% 0.001s 3.29e-07s C 3000 3 theano.compile.ops.Rebroadcast\n",
" 0.1% 99.9% 0.001s 8.97e-07s C 1000 1 theano.tensor.opt.MakeVector\n",
" 0.1% 99.9% 0.001s 4.27e-07s C 2000 2 theano.compile.ops.Shape_i\n",
" 0.0% 99.9% 0.001s 5.84e-07s C 1000 1 theano.tensor.elemwise.CAReduce\n",
" 0.0% 100.0% 0.000s 4.69e-07s C 1000 1 theano.tensor.basic.ScalarFromTensor\n",
" 0.0% 100.0% 0.000s 3.12e-07s C 1000 1 theano.compile.ops.ViewOp\n",
" ... (remaining 0 Classes account for 0.00%(0.00s) of the runtime)\n",
"\n",
"Ops\n",
"---\n",
"<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>\n",
" 28.1% 28.1% 0.427s 4.27e-04s C 1000 1 Dot22Scalar\n",
" 16.1% 44.1% 0.245s 2.45e-04s Py 1000 1 Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}\n",
" 15.6% 59.7% 0.237s 2.37e-04s Py 1000 1 Tri{dtype='float64'}\n",
" 10.4% 70.1% 0.158s 1.58e-04s Py 1000 1 Solve{A_structure='upper_triangular', lower=False, overwrite_A=False, overwrite_b=False}\n",
" 7.3% 77.4% 0.112s 1.12e-04s C 1000 1 Elemwise{Composite{((i0 + (i1 * i2)) * i3)}}[(0, 0)]\n",
" 5.6% 83.0% 0.085s 8.45e-05s C 1000 1 Sum{acc_dtype=float64}\n",
" 4.7% 87.7% 0.072s 7.21e-05s C 1000 1 Alloc\n",
" 3.6% 91.3% 0.055s 5.51e-05s C 1000 1 Elemwise{mul,no_inplace}\n",
" 3.3% 94.6% 0.050s 4.97e-05s Py 1000 1 AllocDiag\n",
" 3.2% 97.8% 0.049s 4.93e-05s C 1000 1 IncSubtensor{InplaceSet;:int64:, :int64:}\n",
" 0.9% 98.8% 0.014s 1.44e-05s Py 1000 1 ExtractDiag{view=False}\n",
" 0.1% 98.9% 0.002s 9.64e-07s C 2000 2 Reshape{1}\n",
" 0.1% 99.0% 0.002s 1.71e-06s C 1000 1 Elemwise{Composite{((i0 * i1) + (i0 * i1))}}\n",
" 0.1% 99.1% 0.002s 1.69e-06s C 1000 1 Elemwise{sub,no_inplace}\n",
" 0.1% 99.2% 0.002s 1.58e-06s C 1000 1 Elemwise{TrueDiv}[(0, 1)]\n",
" 0.1% 99.3% 0.001s 4.54e-07s C 3000 3 Reshape{0}\n",
" 0.1% 99.4% 0.001s 1.31e-06s C 1000 1 Elemwise{exp}\n",
" 0.1% 99.5% 0.001s 1.05e-06s C 1000 1 Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 * i0 * i0), i1) + (i3 * i0) + (i4 * i0))}}[(0, 0)]\n",
" 0.1% 99.5% 0.001s 8.97e-07s C 1000 1 MakeVector{dtype='int64'}\n",
" 0.1% 99.6% 0.001s 7.82e-07s C 1000 1 InplaceDimShuffle{x,x}\n",
" ... (remaining 12 Ops account for 0.42%(0.01s) of the runtime)\n",
"\n",
"Apply\n",
"------\n",
"<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>\n",
" 28.1% 28.1% 0.427s 4.27e-04s 1000 30 Dot22Scalar(InplaceDimShuffle{0,x}.0, InplaceDimShuffle{x,0}.0, TensorConstant{-1.0})\n",
" 16.1% 44.1% 0.245s 2.45e-04s 1000 16 Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}(Elemwise{mul,no_inplace}.0, Elemwise{sub,no_inplace}.0)\n",
" 15.6% 59.7% 0.237s 2.37e-04s 1000 7 Tri{dtype='float64'}(Shape_i{0}.0, Shape_i{1}.0, TensorConstant{0})\n",
" 10.4% 70.1% 0.158s 1.58e-04s 1000 24 Solve{A_structure='upper_triangular', lower=False, overwrite_A=False, overwrite_b=False}(InplaceDimShuffle{1,0}.0, Elemwise{Composite{((i0 * i1) + (i0 * i1))}}.0)\n",
" 7.3% 77.4% 0.112s 1.12e-04s 1000 31 Elemwise{Composite{((i0 + (i1 * i2)) * i3)}}[(0, 0)](IncSubtensor{InplaceSet;:int64:, :int64:}.0, Dot22Scalar.0, Tri{dtype='float64'}.0, <TensorType(float64, matrix)>)\n",
" 5.6% 83.0% 0.085s 8.45e-05s 1000 32 Sum{acc_dtype=float64}(Elemwise{Composite{((i0 + (i1 * i2)) * i3)}}[(0, 0)].0)\n",
" 4.7% 87.7% 0.072s 7.21e-05s 1000 9 Alloc(TensorConstant{(1, 1) of 0.0}, Shape_i{0}.0, Shape_i{1}.0)\n",
" 3.6% 91.3% 0.055s 5.51e-05s 1000 13 Elemwise{mul,no_inplace}(InplaceDimShuffle{x,x}.0, <TensorType(float64, matrix)>)\n",
" 3.3% 94.6% 0.050s 4.97e-05s 1000 25 AllocDiag(Elemwise{TrueDiv}[(0, 1)].0)\n",
" 3.2% 97.8% 0.049s 4.93e-05s 1000 28 IncSubtensor{InplaceSet;:int64:, :int64:}(Alloc.0, AllocDiag.0, ScalarFromTensor.0, ScalarFromTensor.0)\n",
" 0.9% 98.8% 0.014s 1.44e-05s 1000 18 ExtractDiag{view=False}(Elemwise{mul,no_inplace}.0)\n",
" 0.1% 98.9% 0.002s 1.71e-06s 1000 21 Elemwise{Composite{((i0 * i1) + (i0 * i1))}}(TensorConstant{(1,) of -0.5}, Solve{A_structure='lower_triangular', lower=True, overwrite_A=False, overwrite_b=False}.0)\n",
" 0.1% 99.0% 0.002s 1.69e-06s 1000 4 Elemwise{sub,no_inplace}(<TensorType(float64, vector)>, <TensorType(float64, vector)>)\n",
" 0.1% 99.1% 0.002s 1.58e-06s 1000 22 Elemwise{TrueDiv}[(0, 1)](TensorConstant{(1,) of 2.0}, ExtractDiag{view=False}.0)\n",
" 0.1% 99.2% 0.001s 1.31e-06s 1000 0 Elemwise{exp}(s_log_)\n",
" 0.1% 99.3% 0.001s 1.22e-06s 1000 34 Reshape{1}(Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 * i0 * i0), i1) + (i3 * i0) + (i4 * i0))}}[(0, 0)].0, TensorConstant{(1,) of -1})\n",
" 0.1% 99.3% 0.001s 1.05e-06s 1000 33 Elemwise{Composite{(Switch(Identity(GE(i0, i1)), (i2 * i0 * i0), i1) + (i3 * i0) + (i4 * i0))}}[(0, 0)](Elemwise{exp}.0, TensorConstant{0}, TensorConstant{-0.0001}, Reshape{0}.0, Sum{acc_dtype=float64}.0)\n",
" 0.1% 99.4% 0.001s 8.97e-07s 1000 8 MakeVector{dtype='int64'}(Shape_i{0}.0, Shape_i{1}.0)\n",
" 0.1% 99.4% 0.001s 7.82e-07s 1000 10 InplaceDimShuffle{x,x}(s)\n",
" 0.1% 99.5% 0.001s 7.77e-07s 1000 27 InplaceDimShuffle{0,x}(Solve{A_structure='upper_triangular', lower=False, overwrite_A=False, overwrite_b=False}.0)\n",
" ... (remaining 16 Apply instances account for 0.50%(0.01s) of the runtime)\n",
"\n",
"Here are tips to potentially make your code run faster\n",
" (if you think of new ones, suggest them on the mailing list).\n",
" Test them first, as they are not guaranteed to always provide a speedup.\n",
" - Try the Theano flag floatX=float32\n",
" - Try installing amdlibm and set the Theano flag lib.amdlibm=True. This speeds up only some Elemwise operation.\n"
]
}
],
"source": [
"for key, model in models.items():\n",
" print(\"MvNormal with {} log posterior\".format(key))\n",
" lp = model.profile(model.logpt, n=1000).summary()\n",
" print(\"MvNormal with {} grad of log posterior\".format(key))\n",
" dlp = model.profile(pm.theanof.gradient(model.logpt, [model.s_log_]), n=1000).summary()\n",
" print(\"\\n\\n\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From the ops used it looks like the gradients are working properly. It will be good to come back to this after the performance tests.\n",
"\n",
"# Performance\n",
"\n",
"These tests use the same models as before, but increase the size of the covariance (or precision, or Cholesky factor) from $10 \\times 10$ to $1000 \\times 1000$."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"ExecuteTime": {
"end_time": "2017-03-17T22:28:19.900253",
"start_time": "2017-03-17T22:15:05.638817"
},
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-------- cov matrix size: 10 x 10 --------\n",
"cov\n",
"Avg logp: 0.0006245107650756836 \n",
"Avg dlogp: 0.0006097595691680908 \n",
" ---\n",
"tau\n",
"Avg logp: 0.004484299778938294 \n",
"Avg dlogp: 7.514762878417969e-05 \n",
" ---\n",
"chol\n",
"Avg logp: 0.00023196661472320556 \n",
"Avg dlogp: 0.0002770888805389404 \n",
" ---\n",
"-------- cov matrix size: 25 x 25 --------\n",
"cov\n",
"Avg logp: 0.00047759620348612465 \n",
"Avg dlogp: 0.0006503653526306152 \n",
" ---\n",
"tau\n",
"Avg logp: 0.003033899943033854 \n",
"Avg dlogp: 9.696483612060547e-05 \n",
" ---\n",
"chol\n",
"Avg logp: 0.00019541136423746746 \n",
"Avg dlogp: 0.00027946734428405763 \n",
" ---\n",
"-------- cov matrix size: 50 x 50 --------\n",
"cov\n",
"Avg logp: 0.00040795266628265383 \n",
"Avg dlogp: 0.0009344248771667481 \n",
" ---\n",
"tau\n",
"Avg logp: 0.0023460907936096193 \n",
"Avg dlogp: 0.00014599895477294923 \n",
" ---\n",
"chol\n",
"Avg logp: 0.00017871296405792237 \n",
"Avg dlogp: 0.0003032219409942627 \n",
" ---\n",
"-------- cov matrix size: 100 x 100 --------\n",
"cov\n",
"Avg logp: 0.00038023056983947754 \n",
"Avg dlogp: 0.0012469961643218994 \n",
" ---\n",
"tau\n",
"Avg logp: 0.0020779722213745116 \n",
"Avg dlogp: 0.0003375685214996338 \n",
" ---\n",
"chol\n",
"Avg logp: 0.00017246441841125488 \n",
"Avg dlogp: 0.00036454248428344725 \n",
" ---\n",
"-------- cov matrix size: 250 x 250 --------\n",
"cov\n",
"Avg logp: 0.0004370990991592407 \n",
"Avg dlogp: 0.004349910974502563 \n",
" ---\n",
"tau\n",
"Avg logp: 0.0028498863379160565 \n",
"Avg dlogp: 0.0017523510456085205 \n",
" ---\n",
"chol\n",
"Avg logp: 0.00018805797894795736 \n",
"Avg dlogp: 0.0008606476783752442 \n",
" ---\n",
"-------- cov matrix size: 500 x 500 --------\n",
"cov\n",
"Avg logp: 0.0007895946162087576 \n",
"Avg dlogp: 0.030420325756072997 \n",
" ---\n",
"tau\n",
"Avg logp: 0.006698529652186803 \n",
"Avg dlogp: 0.010900997877120972 \n",
" ---\n",
"chol\n",
"Avg logp: 0.00029938908985682896 \n",
"Avg dlogp: 0.004142542362213134 \n",
" ---\n",
"-------- cov matrix size: 750 x 750 --------\n",
"cov\n",
"Avg logp: 0.002144187867641449 \n",
"Avg dlogp: 0.08897358822822571 \n",
" ---\n",
"tau\n",
"Avg logp: 0.016014666199684145 \n",
"Avg dlogp: 0.03246113181114197 \n",
" ---\n",
"chol\n",
"Avg logp: 0.0010221514403820039 \n",
"Avg dlogp: 0.013307257890701293 \n",
" ---\n",
"-------- cov matrix size: 1000 x 1000 --------\n",
"cov\n",
"Avg logp: 0.003685814486609565 \n",
"Avg dlogp: 0.15237263059616088 \n",
" ---\n",
"tau\n",
"Avg logp: 0.03578151358498467 \n",
"Avg dlogp: 0.05463961386680603 \n",
" ---\n",
"chol\n",
"Avg logp: 0.0015350072118971083 \n",
"Avg dlogp: 0.01749673891067505 \n",
" ---\n"
]
}
],
"source": [
"logp_results = {\"chol\": [], \"cov\": [], \"tau\": []}\n",
"dlogp_results = {\"chol\": [], \"cov\": [], \"tau\": []}\n",
"nvec = [10, 25, 50, 100, 250, 500, 750, 1000]\n",
"for n in nvec:\n",
" print(\"-------- cov matrix size: {} x {} --------\".format(n,n))\n",
" y_, mu_, cov_, chol_, tau_ = make_data(n)\n",
" y.set_value(y_)\n",
" mu.set_value(mu_)\n",
" cov.set_value(cov_)\n",
" tau.set_value(tau_)\n",
" chol.set_value(chol_)\n",
" \n",
" for key, model in models.items():\n",
" print(key)\n",
" lp = model.profile(model.logpt, n=1000)\n",
" dlp = model.profile(pm.theanof.gradient(model.logpt, [model.s_log_]), n=1000)\n",
" \n",
" avg_logp = lp.fct_call_time / lp.fct_callcount\n",
" avg_dlogp = dlp.fct_call_time / dlp.fct_callcount\n",
" \n",
" logp_results[key].append(avg_logp)\n",
" dlogp_results[key].append(avg_dlogp)\n",
" \n",
" print(\"Avg logp: {} \\nAvg dlogp: {} \\n ---\".format(avg_logp, avg_dlogp))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Results\n",
"\n",
"The plot below summarizes average runtimes for evaluating the log probability and its gradient on matrices of increasing size."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"ExecuteTime": {
"end_time": "2017-03-17T22:28:20.419422",
"start_time": "2017-03-17T22:28:19.964039"
},
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1AAAAGDCAYAAAAlLvGZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VPXVx/HPISA7KosrhIiiCAhaI7ZWtOICWhRrUUGK\nomjax1KoghtRKgjFBVQQraKiQiO4PLaC4lqLPLiV4I6gRQwQXNgRiOzn+ePe4BCTMCSZ3GTm+369\n5pW565w7k9yTc3+/+V1zd0RERERERGTPakQdgIiIiIiISHWhAkpERERERCROKqBERERERETipAJK\nREREREQkTiqgRERERERE4qQCSkREREREJE4qoEQqmJm9ZGaXVdJrjTSzVWb2bQXtz83siIrYV3Vg\nZkPN7JGo4xCR5Fedc0OimdnjZjYyfN7ZzD6POiYAM/ulmf3XzDaa2fnFLM8zszOiiC0KZtbHzF6N\nOo6qwHQfKCmNmc0COgIHufuWiMOpcszsVuAId/9dBK+dDnwOtHT3FRW0Twdau/uiithfIpnZ40C+\nu98cdSwiqUa5oXTJlhsSLVHnczPLA65099fLuP2/gOnuPi4R+69MUf5OJiO1QEmJzCwD6Aw4cF6C\nXqNmIvYb52ubmVXnv4F0YHVZEmSU73tVUZ73QO+fpDLlhiov0tyQZOfHlsD8qIOoCsqZM6v739RP\nubseehT7AIYBbwF3Ay/EzD8R+BZIi5n3G+Dj8HkN4EbgS2A18DTQOFyWQZB0+wNLgdnh/GfCfa4H\nZgPtYvbdBJgBfA/MBUYCc2KWtwFeA9YQXHW7qJRjmgWMCo/rB+AIIA84I2adW4G/F4n3sjDeVUB2\nuKwbsBXYBmwEPop5jSvD5/3C17oHWAcsBk4K5y8DVgCXxbx2bWBM+FrfAQ8CdYs5jjPC+HeGr/14\nOP88gpP9ujCOo2O2yQNuAD4GtgA1i9mvE1yhAtgXmAysBJYANwM1wmVpwNjw/fgKGBBu+5N9xrz2\nTcBnwFrgMaBOzPKrgEXhZzgdOCScb+F7tyL8/D8B2gNZ4fu+NTz+GeH6hwD/G8b8FTCwyOf6LPD3\ncF9Xxn7WFfH+6aFHKjxQboiNN1Vyw1nhe7geeAB4s4RjWR1+DocDb4TTq4AcYL+Y/R0HvA9sAJ4C\npgEjw2W/ImiNKlx3T+f1pwly1YbwGDPDZVPC9+GH8L24voTPvqT882WR7WsXs20e4e9I+BndC3wd\nPu6N3Qa4HvgmXHYlMfm2hN/H0cB/CH6/nyf8W4nj87wBWB6+H58Dp1Py7+S+wKNhXMvDzy6tlM+1\nH7v/jZ1E8Le3Pvx5Uml/U1Gfuyr0PBh1AHpU3Ud4QrkaOD78ozswZtmXwJkx088AN4bPBwHvAs3D\nE8pDwNRwWUZ40pgM1CdMAMAVQMOYE9CHMfueFj7qAW0JksuccFn9cPpyoCbBSXkV0LaEY5pFkIDa\nhevXIr4k+TBQl6DLypbCkxVF/gGPeY3YxLI9jC8tPAEtBe4Pj/Ws8CTXIFz/HoITeOPw/ZgBjC7h\nWH7F7knmSGATcGZ4XNeHn+E+4fI84EOgBcUk3nCd2AJqMsFJu2H4PnwB9A+X/YGgGGoO7A+8zp4L\nqE/D125McEItTJZdws/sZ+F7ch8//vPUFZgH7EdQTB0NHBwue7xwH+F0jXDdYcA+QCuCf0q6xnxW\n24Dzw3XrFvmsy/3+6aFHKjxQboiNN+lzA9CU4J/4C8L3ZlD4uRc9lj+Fy+sSFKBnhsfSjKD4vTdc\nfx+Ci3LXhPH0DPf3kwKK+M7rm4FzwvdxNPBuTOy7fYbFHFuJ+SfO7XctB0YQ/H4fEB7z28Bt4bJu\nBBcC2hH8vv6dPRdQywkuGNYnKCD3mKuAowh+7wuLwAzg8FJ+J/9B8HdYP4z7P8DvS/lc+/Hj31hj\ngguifcPlvcPpJiX9TUV97qrQ82DUAehRNR/AyeEJrWk4vRC4Jmb5SGBS+Lxh+MfcMpxeAJwes+7B\n4b5q8mPSaVXKa+8XrrNveELcBhxV5LUL/4AvBv6vyPYPAX8pYd+zgBFF5u12gqT4JNk8Zvl/gF5F\n1y3yGrGJ5b8xy44J9xf7D8dq4FiCAmFT4ckuXPYL4KsSjuVX7J4kbwGejpmuQXAC/lXMcV6xh8/d\nCRJfGsHVqrYxy34PzAqfv0F4kg2nz2DPBdQfYqbPAb4Mnz8K3BmzrEH4mWcQJLcvgJ8Ttn7FrPc4\nuxdQJwJLi6xzE/BYzGc1u8jy2M+63O+fHnok+wPlhpTLDcClwDsx00bwT3rssSwtaftwnfOBD8Ln\npxC0wljM8rcpvoCK57z+esyytsAPJX2GxcRVYv6Jc/tdywkuHpwTs6wrkBc+n0RMwUuQZ/dUQN1e\n5Li2Evzel/h5hvtdQZCTaxXZ563s3uPiQIKiv27MvN7Av0v6XNm9gOoL/KfI8neAfiX9TSXTI5n6\nqUrFugx41d1XhdNPhvPuiZl+28z+h+Cq1PvuviRc1hL4h5ntjNnfDoI/1kLLCp+YWRpBM++FBFdt\nCrdrSnDFo2bs+kWetwRONLN1MfNqEjTdl2RZKctKEjuSUQHBSTZe38U8/wHA3YvOa0Bw7PWAeWZW\nuMwITpjxOITgqh7ha+w0s2XAoTHrxHvsTQmubC2JmbckZl+HUPJnUpLYdZaE+yjc1/uFC9x9o5mt\nBg519zfMbALBVdmWZvYcMMTdvy9m/y2BQ4r8LqQB/xdnnBX5/okkK+WG3aVCbtjtfO/ubmb5RdbZ\nbXszOxAYR/BduYYE/+Svjdnfcg//yw7F5ppY8ZzXi34GdcysprtvL+WYCpWYfwiKo72x2/vMT/Nc\nbsyysuTMWgS/+yV+nu4+y8z+TFAstTOzV4Br3f3rYvbfMtznNzG/VzWIP7cXPd7COFMiZ6qAkp8w\ns7rARUBazBCotYH9zKyju3/k7p+Z2RLgbOASgqRZaBnB1ay3itl3Rvg09sR5CdCD4IpJHsHVxbUE\nCWIlQRNyc4KWCAi6GcS+1pvufuZeHKIXmd5EkJwKHVSOfZXHKoKE2c7dl5dh+68JrmICwZc2Cd6r\n2H3FG+8qgqtwLQm66kHwxeTCfX1D8JkUiv1MShK7TnoYb2HcLWPirk/w3YblAO4+HhhvZgcQ9HW/\njuAKXNFjWUZwRbZ1KTGUdvwV+f6JJB3lhpTNDbud78PtmxdZp+j2fw3nHePua8IhwCfE7O9QM7OY\nIiqdoAWnqHjO66XZ0+dQav7ZS4X7Khx0IjbPVUTO3Ebwu1Dq5+nuTwJPmlkjglbXOwhai4rLmVsI\nWpNLKjb3lDNbFpmXDrwc5/bVWnKNiCEV5XyCq4JtCboPHEvw3ZP/I2jKL/QkQV/oUwj6uRd6EBhl\nZi0BzKyZmfUo5fUaEvwRryZIVn8tXODuO4DngFvNrJ6ZtSkSwwvAkWbW18xqhY8TzOzovTjeD4Fe\n4baZBP2x4/UdkFERo8u4+06C/vT3hMUCZnaomXWNcxdPA782s9PNrBYwmOB9fbsMsewI9zfKzBqG\nn+W1BP22C19rUBjffgRfWt2TP5pZczNrDGQTfHEYYCpwuZkda2a1CT7/99w9L/wsTwyPZxNBX/fC\nq9DfEfSHL/QfYIOZ3WBmdc0szczam9kJcR52hb1/IklKuSF+yZQbXgSOMbPzw5HY/siei8mGBIMV\nrDezQwkufBV6h6D4HRi+txcAnUrYT3nP60XzRFEl5p849190XzeHv9dNCb63FZszLzezo82sHsFF\nwD35nZm1DdcfATwbk5uL/TzN7Cgz6xIey2Z+HFAEivxOuvs3wKvAWDNrZGY1zOxwMzs1zuOdSfA3\ndomZ1TSziwnODS/EuX21pgJKinMZQf/ipe7+beGD4OpRH/txKMupwKnAGzHdOSBotp8OvGpmGwi+\nVHliKa83maDZdzlBa8e7RZYPILjy+C1B94upBCcL3H0DwZdtexFcDfmW4GpL7b043lsIRgxaCwxn\n9yume1L4z8FqM3u/1DXjcwPBl0HfNbPvCQZnOCqeDd39c+B3BF+CXQWcC5zr7lvLGMufCIqWxcAc\ngvdlUrjsYYIT78fABwQn0u0E/1yV5Mlwm8UEVxpHhnG/TvAZ/C/BVbrDCT5PgEbha60l+B1ZDdwV\nLnsUaGtm68zsn2Fi6U7wT91X4XvwCMHvzh4l4P0TSTbKDfFLmtwQfoYXAncSnIPbEnRHK+3+X8MJ\nBmZYT1CAPRezv60E3Tv7EYx8d3Hs8iKvXa7zOsGgEjeHeWJIMfsvLf/srZEE78vHBCPGvs+Pee4l\nYDzwb8LPMdymtPdwCsF3fb8F6gADw32V9nnWBm4P539LMDDETeH+ivudvJRg8InCEXKfJfhu4h65\n+2qCz2Ywwe/F9UD3In/zSUs30pVqx8zuILh542VRxyIBMzsbeNDdizbnFy7Po5rcbFBEqiflhsoR\ntmDkA33c/d9Rx1MdhS2hnxIMc/6T7nMW3Kj67+7+SGXHJvFRC5RUeWbWxsw6WKATwX1C/hF1XKks\n7EpxTthsfyjwF/SZiEglUm6oPGbW1cz2C7uGDSX4HlrRFkEphZn9xsxqm9n+BK2hM+Ic6EKqIBVQ\nUh00JGje30TwvZmxBPcnkugYQReNtQRd+BYQ9PcWSQpm1s3MPjezRWZ2YzHLTzGz981su5n1LLIs\n3cxeNbMFZvaZ/ThAglQs5YbK8wuCrteFXcbOd/cfog2p2vk9wRDjXxJ0d/+faMOR8lAXPhERkRgW\nDJ/9BcGNKvOBuUBvd/8sZp0Mgu/oDQGmu/uzMctmAaPc/TUzawDsdPeCSjsAERFJKA1jLiIisrtO\nwCJ3XwxgZtMIhtPeVUAVjtJlu9/TCDNrS3BD6dfC9TZWUswiIlJJ1IVPRERkd4ey+w0g89n95pCl\nORJYZ2bPmdkHZnZX2KIlIiJJIiVaoJo2beoZGRlRhyEiktLmzZu3yt2bRR1HgtUEOgPHAUsJvpvT\nj2DY/d2YWRaQBVC/fv3j27RpU3lRiojIT8Sbp1KigMrIyCA3NzfqMEREUpqZLYk6hjgtB1rETDcP\n58UjH/gwpvvfP4GfU0wB5e4TgYkAmZmZrjwlIhKtePOUuvCJiIjsbi7Q2swOM7N9CG6sOX0vtt3P\nzAqvYHYh5rtTIiJS/amAEhERiRHem2UA8ArBEP1Pu/t8MxthZucBmNkJZpYPXAg8ZGbzw213EIzM\n9y8z+4RgyP+HozgOERFJjJTowiciIrI33H0mMLPIvGExz+cSdO0rbtvXgA4JDVBERCKTsgXUtm3b\nyM/PZ/PmzVGHUuHq1KlD8+bNqVWrVtShiIhIGSVzngLlKhGpvlK2gMrPz6dhw4ZkZGRgZlGHU2Hc\nndWrV5Ofn89hhx0WdTgiIlJGyZqnQLlKRKq3lP0O1ObNm2nSpEnSJSUzo0mTJkl7xVJEJFUka54C\n5SoRqd5StoACkjIpQfIel4hIqknm83kyH5uIJLeULqBERERERET2hgqoOOXkQEYG1KgR/MzJiToi\nERGRHylPiYhUDhVQccjJgawsWLIE3IOfWVnlT06TJ0+mQ4cOdOzYkb59+5KXl0eXLl3o0KEDp59+\nOkuXLmX9+vW0bNmSnTt3ArBp0yZatGjBtm3bKuDIRET2LCcnh4yMDGrUqEFGRgY5+s+8ylGeEpFU\nVtkXkFJ2FL6ifvWrn87r3h2GDIHsbCgo2H1ZQQEMGAB9+hS//axZpb/e/PnzGTlyJG+//TZNmzZl\nzZo1XHbZZbsekyZNYuDAgfzzn//k2GOP5c033+S0007jhRdeoGvXrhr2VUQqRU5ODllZWRSEJ8El\nS5aQlZUFQJ/CE6BUCuUpEZGfKryAVHgOLLyABD+e/yqaWqDisHRp8fPXrSv7Pt944w0uvPBCmjZt\nCkDjxo155513uOSSSwDo27cvc+bMAeDiiy/mqaeeAmDatGlcfPHFZX9hEZG9kJ2dvat4KlRQUEB2\ndnZEEUlxlKdEJFWVdAEpkWlKLVCh0q7EpacH1WxRLVvGt315nXfeeQwdOpQ1a9Ywb948unTpkrgX\nExGJsbSE/8xLmi+JozwlIvJTJaWjRKYptUDFYdQoqFdv93n16gXzy6pLly4888wzrF69GoA1a9Zw\n0kknMW3aNCDoNtO5c2cAGjRowAknnMCgQYPo3r07aWlpZX9hEZG9kJ6evlfzJRrKUyKSqlq0KH5+\nItOUWqDiUNh/Mjs7qGbT04OkVJ5+le3atSM7O5tTTz2VtLQ0jjvuOO677z4uv/xy7rrrLpo1a8Zj\njz22a/2LL76YCy+8kFmJvIQoIlLEzTffzKBBg3brxlevXj1Glec/c6lwylMikqo6dfppa1N5LyDt\nibl74vZeRWRmZnpubu5u8xYsWMDRRx8dUUSJl+zHJyKVo0ePHnz11VesX7+eZcuWkZ6ezqhRo8o0\ngISZzXP3zASEWe2lYp6C1DhGEUmcTz+Fn/0Mjj8evvmm/BeQ4s1TCW2BMrNuwDggDXjE3W8vsrw2\nMBk4HlgNXOzueWbWCZhYuBpwq7v/I9wmD9gA7AC2KxmLiCTGBx98wPTp0xk+fDjDhg2LOhwREZFd\ntm+Hfv1gv/1g+nRo1qzyXjthBZSZpQH3A2cC+cBcM5vu7p/FrNYfWOvuR5hZL+AO4GLgUyDT3beb\n2cHAR2Y2w923h9ud5u6rEhW7iIjAiBEj2HfffRk4cGDUoYiIiOxm6lSYNw+eeqpyiydI7CASnYBF\n7r7Y3bcC04AeRdbpATwRPn8WON3MzN0LYoqlOkDy9zMUEalCPvroI/75z3/y5z//mf322y/qcERE\nRHbTpw/MnAkXXlj5r53IAupQYFnMdH44r9h1woJpPdAEwMxONLP5wCfAH2IKKgdeNbN5ZpZV0oub\nWZaZ5ZpZ7sqVKyvkgEREUsXYsWNp1KgRgwYNijoUERGRXXbsgG+/hRo14OyzwazyY6iyw5i7+3vu\n3g44AbjJzOqEi052958BZwN/NLNTSth+ortnuntms8pu1xMRqeYmTJjAjBkz2H///aMORUREZJd7\n74U2beCrr6KLIZEF1HIgdmT25uG8Ytcxs5rAvgSDSezi7guAjUD7cHp5+HMF8A+CroIiIlJB3J1G\njRpxyinFXp8SERGJxBdfwM03w6mnQkZGdHEksoCaC7Q2s8PMbB+gFzC9yDrTgcvC5z2BN9zdw21q\nAphZS6ANkGdm9c2sYTi/PnAWwYAT1c66det44IEHog5DRGQ38+fP55hjjuHDDz+MOhSJmPKUiFQl\nO3dC//5Qpw48+GA0XfcKJayACr+zNAB4BVgAPO3u881shJmdF672KNDEzBYB1wI3hvNPJhh570OC\nVqarw1H3DgTmmNlHwH+AF9395UQdQ6ycnBwyMjKoUaMGGRkZ5OTklGt/SkwiUhXddtttLFmyhBYl\n3dpdqizlKRFJZhMmwJw5QRe+gw+ONpaE3gfK3WcCM4vMGxbzfDPwk7Ez3H0KMKWY+YuBjhUfaely\ncnLIysqioKAAgCVLlpCVFYxfUZabSQLceOONfPnllxx77LGcdtppfPzxx6xdu5Zt27YxcuRIevTo\nQV5eHt27d+fTT4NGtjFjxrBx40ZuvfXWCjkuEZFYn332GU8//TQ33HADTZo0iToc2QvKUyKS7N59\nNxg04tJLo44kwQVUdfKrX/3qJ/O6d+/OkCFDyM7O3pWUChUUFDBgwIBdiano9rNmzSr19W6//XY+\n/fRTPvzwQ7Zv305BQQGNGjVi1apV/PznP+e8884rdXsRkYo2cuRI6tWrx+DBg6MORYqhPCUiqSwn\nB374Idque4VUQMVh6dKlxc5ft25dhezf3Rk6dCizZ8+mRo0aLF++nO+++65C9i0iEo/PP/+cadOm\ncd1119G0adOow5G9pDwlIslq+nRo1w4OPxzq1Ys6moAKqFBpV+LS09NZsmTJT+a3bNkyru33JCcn\nh5UrVzJv3jxq1apFRkYGmzdvpmbNmuzcuXPXeps3by7za4iIlKZly5aMHz+eiy66KOpQpATKUyKS\navLy4JJL4PTT4fnno47mR1X2PlBVyahRo6hXpOStV68eo0aNKvM+GzZsyIYNGwBYv349BxxwALVq\n1eLf//73riR44IEHsmLFClavXs2WLVt44YUXyn4QIiKlqFOnDgMGDOCAAw6IOhQpA+UpEUk27nDV\nVUGXvfHjo45md2qBikNh//Hs7GyWLl1Keno6o0aNKvMXcwGaNGnCL3/5S9q3b88JJ5zAwoULOeaY\nY8jMzKRNmzYA1KpVi2HDhtGpUycOPfTQXfNFRCpSdnY2rVu3pl+/flGHImWkPCUiyeaRR+D11+GB\nByCmMb1KMHePOoaEy8zM9Nzc3N3mLViwgKOPPjqiiBIv2Y9PRCrGl19+yVFHHcXAgQO5++67E/pa\nZjbP3TMT+iLVVCrmKUiNYxSRvbdsWfC9p8zMoIiqUUl95uLNU+rCJyKSwkaNGkWtWrW47rrrog5F\nREQEgPr14aKLglaoyiqe9oa68ImIpKjFixczefJk/vjHP3Jw1HclFBERCTVuHBRPVVUVrOlERKQy\n/PWvf6VmzZrccMMNUYdS5ZhZNzP73MwWmdmNxSw/xczeN7PtZtazmOWNzCzfzCZUTsQiItXf11/D\nmWfCwoVRR1I6tUCJiKSoLl260Lp1aw455JCoQ6lSzCwNuB84E8gH5prZdHf/LGa1pUA/YEgJu7kN\nmJ3IOEVEkok7/M//wJw5kJYWdTSlUwElIpKiLrnkkqhDqKo6AYvcfTGAmU0DegC7Cih3zwuX7Sy6\nsZkdDxwIvAxo0AwRkThMnRrcNHfsWGjdOupoSqcufCIiKWbp0qWMGTOGTZs2RR1KVXUosCxmOj+c\nt0dmVgMYS8ktUyIiUsS338Kf/gS/+AUMGhR1NHumAqoK6devH88++2zc6+fl5dG+ffsERiQiyej2\n229n6NChrF69OupQktHVwEx3z9/TimaWZWa5Zpa7cuXKSgit/JSnRCQRRoyATZtg0qSq330PVEDF\nLeeTHDLuzaDG8Bpk3JtBzic5UYckIrLX8vPzefTRR7niiitIT0+POpyqajnQIma6eTgvHr8ABphZ\nHjAGuNTMbi9uRXef6O6Z7p7ZrFmz8sQLKE+JSPV1xx0wYwZUl3txq4CKQ84nOWTNyGLJ+iU4zpL1\nS8iakVXu5DR58mQ6dOhAx44d6du3LwCzZ8/mpJNOolWrVruu8rk71113He3bt+eYY47hqaeeKvcx\niUhquv3223F3brrppqhDqcrmAq3N7DAz2wfoBUyPZ0N37+Pu6e6eQdCNb7K7/2QUv4qmPCUi1dG6\ndbB5MzRsGIy+V11oEInQrx7/1U/mdT+yO0NOGkL2v7Ip2Faw27KCbQUMmDmAPsf0KXb7Wf1mlfp6\n8+fPZ+TIkbz99ts0bdqUNWvWcO211/LNN98wZ84cFi5cyHnnnUfPnj157rnn+PDDD/noo49YtWoV\nJ5xwAqecckp5DldEUtDy5ct5+OGH6devHy1btow6nCrL3beb2QDgFSANmOTu881sBJDr7tPN7ATg\nH8D+wLlmNtzd2yUyLuUpEUk2v/89fPEFzJ0LNatRVaIWqDgsXb+02PnrNq8r8z7feOMNLrzwQpo2\nbQpA48aNATj//POpUaMGbdu25bvvvgNgzpw59O7dm7S0NA488EBOPfVU5s6dW+bXFpHUtHLlSo49\n9li1PsXB3We6+5Hufri7jwrnDXP36eHzue7e3N3ru3uT4oond3/c3QdURrzKUyJS3Tz3HDz9NPTs\nWb2KJ1AL1C6lXYlL3zedJeuX/GR+y31/vIK7pyt58apdu/au5+5eIfsUEQE49thjee+996IOQ8pI\neUpEksXq1cE9n447Dq6/Pupo9p5aoOIw6vRR1KtVb7d59WrVY9Tpo8q8zy5duvDMM8/sGgVrzZo1\nJa7buXNnnnrqKXbs2MHKlSuZPXs2nTp1KvNri0jqeemll0o9z0j1pjwlItXJoEGwZg089hjUqhV1\nNHtPLVBxKOw/nv2vbJauX0r6vumMOn3Urvll0a5dO7Kzszn11FNJS0vjuOOOK3Hd3/zmN7zzzjt0\n7NgRM+POO+/koIMOIi8vr8yvLyKp49tvv+WCCy6gb9++TJw4MepwJAGUp0Skuvj+e5g3D7KzoWPH\nqKMpG0uF5vfMzEzPzc3dbd6CBQs4+uijI4oo8ZL9+EQkfoMHD+bee+9l4cKFtI7w9u5mNs/dMyML\noApLxTwFqXGMIvJTmzdDjRqwzz5RR7K7ePOUuvCJiCSxFStW8Le//Y0+ffpEWjyJiIj8/e+wcSPU\nqVP1iqe9oQJKRCSJjRkzhi1btnDzzTdHHYqIiKSwl16Cvn1hwoSoIyk/FVAiIknK3Zk/fz69e/fm\nyCOPjDocERFJUevXQ1YWtG0L11wTdTTll9KDSLg7ZhZ1GBUuFb7XJiJ7Zma8+OKLbN68OepQpIyS\nNU+BcpVIKrnuOvj6a/jf/4WYOyFUWynbAlWnTh1Wr16ddCdwd2f16tXUqVMn6lBEJEJr167l66+/\nBtD5oJpK1jwFylUiqeT11+Hhh2HwYEiWuxukbAtU8+bNyc/PZ+XKlVGHUuHq1KlD8+bNow5DRCJ0\n1113MW7cOJYsWULTpk2jDkfKIJnzFChXiaSKgw6Ciy+G4cOjjqTipGwBVatWLQ477LCowxARqXBr\n1qzhvvvuo3v37iqeqjHlKRFJBu3bw7RpUUdRsVK2C5+ISLK655572LhxI7fcckvUoYiISIqaNQsu\nvTQYQCLZpGwLlIhIMlq7di3jx4+nZ8+etG/fPupwREQkBW3aBP37gxnUqhV1NBVPBZSISBJ58cUX\n2bBhg1qfREQkMtnZsHgxvPkm1KsXdTQVT134RESSyO9+9zv++9//0qFDh6hDERGRFPTWWzB+PAwY\nAKecEnU0iZHQAsrMupnZ52a2yMxuLGZ5bTN7Klz+npllhPM7mdmH4eMjM/tNvPsUEUlVGzduBODw\nww+POBIREUlF7kHhlJEBo0dHHU3iJKyAMrM04H7gbKAt0NvM2hZZrT+w1t2PAO4B7gjnfwpkuvux\nQDfgITN4bGzbAAAgAElEQVSrGec+RURSzvr162nVqhUTJkyIOhQREUlRZvDMM8Goew0aRB1N4iSy\nBaoTsMjdF7v7VmAa0KPIOj2AJ8LnzwKnm5m5e4G7bw/n1wEK7yIYzz5FRFLOfffdx8qVKznppJOi\nDkVERFLQ6tVBC9QRRyTPDXNLksgC6lBgWcx0fjiv2HXCgmk90ATAzE40s/nAJ8AfwuXx7FNEJKV8\n//333H333Zx77rn87Gc/izocERFJMZs3w8knwx//GHUklaPKDiLh7u+5ezvgBOAmM6uzN9ubWZaZ\n5ZpZbrLexV1EBOD+++9n7dq1DBs2LOpQREQkBQ0fDgsXQo8U6ReWyAJqOdAiZrp5OK/YdcysJrAv\nsDp2BXdfAGwE2se5z8LtJrp7prtnNmvWrByHISJSdW3dupV7772Xc845h8zMzKjDERGRFJObC3fd\nBVdcAV27Rh1N5UjkfaDmAq3N7DCCIqcXcEmRdaYDlwHvAD2BN9zdw22Wuft2M2sJtAHygHVx7FNE\nJGXss88+zJ49G3ff88oiIiIVaMsWuPxyOPBAGDs26mgqT8IKqLD4GQC8AqQBk9x9vpmNAHLdfTrw\nKDDFzBYBawgKIoCTgRvNbBuwE7ja3VcBFLfPRB2DiEhV5u6YGUcddVTUoYiISAr69FNYuhSefBL2\n2y/qaCpPIlugcPeZwMwi84bFPN8MXFjMdlOAKfHuU0QkFY0dO5a3336bqVOnUrt27ajDERGRFHP8\n8ZCXB/vvH3UklavKDiIhIiIlKygo4M4772TTpk0qnkREpFJt2wZTpwbDlqda8QQqoEREqqUHH3yQ\nlStXauQ9ERGpdKNHwyWXwOzZUUcSDRVQIiLVTGHr0+mnn84vf/nLqMMREZEU8sknMHIk9OoFp54a\ndTTRSOh3oEREpOI9/PDDfPfddzzzzDNRhyIiIilk+/Zg1L399oP77os6muiogBIRqWYuuOACdu7c\nSefOnaMORUREUshdd8G8efD009C0adTRREdd+EREqpkWLVpwzTXXRB2GiIikmKOOgj/8AS78yRja\nqUUFlIhINbF582Z69+7NBx98EHUoSc/MupnZ52a2yMxuLGb5KWb2vpltN7OeMfOPNbN3zGy+mX1s\nZhdXbuQiIolzwQXwt79FHUX0VECJiFQTjz76KNOmTWPt2rVRh5LUzCwNuB84G2gL9DaztkVWWwr0\nA54sMr8AuNTd2wHdgHvNLIVuLykiyei++4Luezt3Rh1J1aACSkSkGtiyZQu33347J598MqeddlrU\n4SS7TsAid1/s7luBaUCP2BXcPc/dPwZ2Fpn/hbv/N3z+NbACaFY5YYuIVLwvvoDrr4e33gKzqKOp\nGjSIhIhINTBp0iTy8/OZNGkSpgyWaIcCy2Km84ET93YnZtYJ2Af4soTlWUAWQHp6+t5HKSKSYDt2\nwBVXQN26Qdc9pZ+AWqBERKq4rVu3Mnr0aH7xi19wxhlnRB2OxMHMDgamAJe7e7GdXtx9ortnuntm\ns2ZqpBKRqmfChKDl6d574eCDo46m6lALlIhIFbd161b69OnDGWecodanyrEcaBEz3TycFxczawS8\nCGS7+7sVHJuISKVYtQqGDoVzzoG+faOOpmpRASUiUsU1aNCA0aNHRx1GKpkLtDazwwgKp17AJfFs\naGb7AP8AJrv7s4kLUUQksZo2heeeg3bt1HWvKHXhExGpwmbMmMHMmTNx96hDSRnuvh0YALwCLACe\ndvf5ZjbCzM4DMLMTzCwfuBB4yMzmh5tfBJwC9DOzD8PHsREchohIma1fH/zs2hWaN482lqpILVAi\nIlXUtm3bGDhwIE2bNuXss8+OOpyU4u4zgZlF5g2LeT6XoGtf0e3+Dvw94QGKiCRIXh4cdxyMH6+u\neyVRASUiUkVNmTKFvLw87rvvPn33SUREEs4drrwStm+HU0+NOpqqSwWUiEgVtH37dkaNGsXxxx/P\nr3/966jDERGRFPDww/CvfwVDluvuCiVTASUiUgXl5OSwePFinn/+ebU+iYhIwi1dCkOGwGmnQVZW\n1NFUbRpEQkSkCtq+fTunn3465557btShiIhICnjtteDnI49ADVUIpdLbIyJSBfXv35/XXntNrU8i\nIlIp+veHL7+EVq2ijqTqUwElIlKF7Nixg3/+85/s2LFDxZOIiCTc11/DnDnB82bNoo2lulABJSJS\nhTz11FP85je/4aWXXoo6FBERSXLu8Ic/QLdusGZN1NFUHxpEQkSkitixYwe33XYb7du355xzzok6\nHBERSXJPPgkzZsDYsdC4cdTRVB8qoEREqohnnnmGhQsX8vTTT1ND3+AVEZEE+vZbGDgQfvELGDQo\n6miqF2VoEZEqYOfOndx22220bduW3/72t1GHIyIiScwdrr4aNm2CSZMgLS3qiKoXtUCJiFQBeXl5\nrF+/njFjxqj1SUREEsodOnSAk0+GNm2ijqb6UQElIlIFtGrVikWLFlGrVq2oQxERkSRXowbcemvU\nUVRfuswpIhKxr776iq1bt1KnTh3S1I9CREQS6M9/hpdfjjqK6k0FlIhIhHbu3Mn5559P9+7dow5F\nRESS3HPPwbhxMG9e1JFUb+rCJyISoenTp/Pxxx8zZcqUqEMREZEktno1/M//wHHHwfXXRx1N9aYC\nSkQkIu7OiBEjOOKII+jVq1fU4YiISBIbNCi4We6rr4K+bls+KqBERCIyY8YMPvjgAx5//HFq1tTp\nWEREEuPf/4acHPjLX6Bjx6ijqf6UsUVEIjJlyhRatWpFnz59og5FRESS2CmnwCOPQN++UUeSHBI6\niISZdTOzz81skZndWMzy2mb2VLj8PTPLCOefaWbzzOyT8GeXmG1mhfv8MHwckMhjEBFJlKlTp/Lq\nq6+q9UlERBJm48bgRrn9+8M++0QdTXJIWAFlZmnA/cDZQFugt5m1LbJaf2Ctux8B3APcEc5fBZzr\n7scAlwFFv13dx92PDR8rEnUMIiKJ4O4UFBRQs2ZNDj/88KjDERGRJPXSS9CqFXz0UdSRJJdEtkB1\nAha5+2J33wpMA3oUWacH8ET4/FngdDMzd//A3b8O588H6ppZ7QTGKiJSaV5++WUyMjL45JNPog5F\nRESS1Pr1cNVV0KwZtGkTdTTJJZEF1KHAspjp/HBeseu4+3ZgPdCkyDq/Bd539y0x8x4Lu+/dYmZW\n3IubWZaZ5ZpZ7sqVK8tzHCIiFcbdGT58OPXq1eOoo46KOhwREUlSQ4bAN9/AY49BbTVDVKgqfSNd\nM2tH0K3v9zGz+4Rd+zqHj2K/DufuE909090zmzVrlvhgRUTi8Oqrr/Lee+8xdOhQ9lFndBERSYDX\nXgsGjRg8GDp1ijqa5JPIAmo50CJmunk4r9h1zKwmsC+wOpxuDvwDuNTdvyzcwN2Xhz83AE8SdBUU\nEanyClufWrRoQb9+/aIOR0REktSTT8KRR8Lw4VFHkpwSWUDNBVqb2WFmtg/QC5heZJ3pBINEAPQE\n3nB3N7P9gBeBG939rcKVzaymmTUNn9cCugOfJvAYREQqzJw5c3jnnXe46aab1PokIiIJM2kSvPkm\n1K0bdSTJKWFj57r7djMbALwCpAGT3H2+mY0Act19OvAoMMXMFgFrCIosgAHAEcAwMxsWzjsL2AS8\nEhZPacDrwMOJOgYRkYp08sknM2PGDM4888yoQxERkST0/vtwwAHQvDkcdFDU0SQvc/eoY0i4zMxM\nz83NjToMEUlh7k4JY96kDDOb5+6ZUcdRFSlPiUh5bdoEHTrA/vvD3LmQ4imnTOLNU1V6EAkRkWTR\no0cP7rnnnqjDEBGRJDV0KCxeDHffreIp0VRAiYgk2KxZs5gxYwZpaWlRhyIiIklozhy47z4YMABO\nOSXqaJKfCigRkQQbPnw4Bx10EFdddVXUoYiISJIpKIArroCMDBg9OupoUoMKKBGRBJo9ezazZs3i\n+uuvp66GQ6o2zKybmX1uZovM7MZilp9iZu+b2XYz61lk2WVm9t/wcVnRbUVEKtLmzdC+fXDfpwYN\noo4mNSRsFD4REYERI0Zw4IEH8vvf/37PK0uVYGZpwP3AmUA+MNfMprv7ZzGrLQX6AUOKbNsY+AuQ\nCTgwL9x2bWXELiKpp3FjeO65qKNILWqBEhFJoOzsbCZMmEC9evWiDkXi1wlY5O6L3X0rMA3oEbuC\nu+e5+8fAziLbdgVec/c1YdH0GtCtMoIWkdSyeTNcfjn8979RR5J61AIlIpJAp512WtQhyN47FFgW\nM50PnFiObQ8tbkUzywKyANLT0/c+ShFJacOHw+OPQ+/e0Lp11NGkFrVAiYgkwHvvvcfAgQNZu1Y9\nt6R47j7R3TPdPbNZs2ZRhyMi1cjcuXDnndC/P5x1VtTRpB4VUCIiCXDrrbcydepUatWqFXUosveW\nAy1ippuH8xK9rYjIHm3ZEnTdO/hgGDs26mhSkwooEZEK9p///IeXX36ZwYMH00BDIlVHc4HWZnaY\nme0D9AKmx7ntK8BZZra/me0PnBXOExGpEGPHwvz5MHEi7Ltv1NGkJn0HSkSkgg0fPpzGjRvzxz/+\nMepQpAzcfbuZDSAofNKASe4+38xGALnuPt3MTgD+AewPnGtmw929nbuvMbPbCIowgBHuviaSAxGR\npPSHP8ABB8A550QdSepSASUiUoHmzp3LzJkzGTVqFA0bNow6HCkjd58JzCwyb1jM87kE3fOK23YS\nMCmhAYpIytm2LfjZuDFceWW0saQ6deETEalADRo0oHfv3gwYMCDqUEREJImMHg0nnggbNkQdiaiA\nEhGpQEcffTRPPvkkjRo1ijoUERFJEp98AiNHQps2oM4N0VMBJSJSQe677z7+qzsaiohIBdq2Dfr1\ng/33h/Hjo45GQAWUiEiF+Oijjxg4cCBPPvlk1KGIiEgSuesueP99uP9+aNo06mgEVECJiFSIESNG\nsO+++zJo0KCoQxERkSSxbRtMmQI9ewYPqRo0Cp+ISDl9/PHHPPfccwwbNoz99tsv6nBERCRJ1KoF\nc+fC1q1RRyKx1AIlIlJOt912G40aNeLPf/5z1KGIiEiSmD0bfvgBGjQIhi6XqkMFlIhIOezYsYO6\ndetyzTXXsP/++0cdjoiIJIEvvoCuXeGGG6KORIqjLnwiIuWQlpbG5MmTcfeoQxERkSSwYwdccQXU\nrQtDh0YdjRRHLVAiImW0ePFiPvzwQwDMLOJoREQkGUyYAG+9BePGwUEHRR2NFEcFlIhIGd188810\n7tyZDbotvIiIVIAvv4SbboJf/xp+97uoo5GSqIASESmDhQsXMm3aNK6++moa6rbwIiJSATZtguOO\ng4ceAnVsqLr0HSgRkTIYOXIkdevWZciQIVGHIiIiSaJDh6D7nlRtaoESEdlLX3zxBVOnTuXqq6+m\nWbNmUYcjIiLVVE4OZGRAjRrQqBE88kjUEUk81AIlIrKX3n33XRo1aqTWJxERKbOcHMjKgoKCYHrD\nBhg4MBh9r0+faGOT0qkFSkRkL1166aXk5+dz4IEHRh2KiIhUU9nZPxZPhX74IZgvVZsKKBGRvbBo\n0SIA6tevH3EkIiJSnS1dunfzpepQASUiEqcvv/ySNm3aMGHChKhDERGRaq5Fi+Lnp6dXbhyy91RA\niYjEafTo0dSsWZMLLrgg6lBERKSa++tfoU6d3efVqwejRkUTj8RPBZSISBzy8vJ44oknuOqqqzjk\nkEOiDkdERKqxFSuCgSIeeQRatgzu+dSyJUycqAEkqoOEFlBm1s3MPjezRWZ2YzHLa5vZU+Hy98ws\nI5x/ppnNM7NPwp9dYrY5Ppy/yMzGm+k2YyKSeH/961+pUaMGN9xwQ9ShiIhINbZqFRx7LAwbFhRL\neXmwc2fwU8VT9ZCwAsrM0oD7gbOBtkBvM2tbZLX+wFp3PwK4B7gjnL8KONfdjwEuA6bEbPM34Cqg\ndfjolqhjEBEB2LhxI88++yxXXnklzZs3jzocERGpptzhyith9Wr47W+jjkbKKq77QJlZfeAHd99p\nZkcCbYCX3H1bKZt1Aha5++JwH9OAHsBnMev0AG4Nnz8LTDAzc/cPYtaZD9Q1s9pAY6CRu78b7nMy\ncD7wUjzHISJSFg0aNOCLL77A3aMORfZSGfOXiEhCTJwIzz8Pd98NHTtGHY2UVbwtULOBOmZ2KPAq\n0Bd4fA/bHAosi5nOD+cVu467bwfWA02KrPNb4H133xKun7+HfYqIVJjNmzfj7jRt2pRmzZpFHY7s\nvbLkLxGRCrdgAVxzDZx1FgwaFHU0Uh7xFlDm7gXABcAD7n4h0C5xYYUvataOoFvf78uwbZaZ5ZpZ\n7sqVKys+OBFJCYMHD6Zz587s2LEj6lCkbCLJXyIiRS1cCAccAI8/DjU0jFu1FncBZWa/APoAL4bz\n0vawzXIgdoT75uG8Ytcxs5rAvsDqcLo58A/gUnf/Mmb92C8gFLdPANx9ortnunumrhqLSFksX76c\nRx55hLZt25KWtqdTnlRRZclfIiIV7je/gS++gIMPjjoSKa94C6g/AzcB/3D3+WbWCvj3HraZC7Q2\ns8PMbB+gFzC9yDrTCQaJAOgJvOHubmb7ESS6G939rcKV3f0b4Hsz+3k4+t6lwPNxHoOIyF654447\n2LlzJ0OHDo06FCm7suQvEZEK8/rr8MQTwQAS++wTdTRSEeIaRMLd3wTejJleDAzcwzbbzWwA8ArB\n1b5JYfIaAeS6+3TgUWCKmS0C1hAUWQADgCOAYWY2LJx3lruvAK4m6L9el2DwCA0gISIV7uuvv2bi\nxIlcdtllZGRkRB2OlFFZ8peISEVZuRL69oXGjaFXL6hdO+qIpCKUWkCZ2QygxGGn3P280rZ395nA\nzCLzhsU83wxcWMx2I4GRJewzF2hf2uuKiJTXuHHj2L59u1qfqqny5i8RkfJyh/79Yc0aeOUVFU/J\nZE8tUGPCnxcABwF/D6d7A98lKigRkajdfPPNdO7cmVatWkUdipSN8peIROrBB2HGDLjnHujQIepo\npCKVWkCFXR8ws7HunhmzaIaZ5SY0MhGRCDVs2JDu3btHHYaUUXnzl5l1A8YRdEF/xN1vL7K8NjAZ\nOJ5g8KOL3T3PzGoBjwA/I8ixk919dEUck4hUH8uWwbXXQteuMFCdhpNOvINI1A+/eAuAmR0G1E9M\nSCIi0fnuu+/IzMzk7bffjjoUqRh7nb/MLA24HzgbaAv0NrO2RVbrD6x19yOAewhuuQFBt/Ta7n4M\nQXH1ezPLqIDjEJFqpHlzePhhDVmerOIaRAK4BphlZosBA1pShnsziYhUdWPGjOGDDz7QTXOTR1ny\nVydgUTjgBGY2DegBfBazTg/g1vD5s8CEcHRYJyjaahIMdrQV+L5iDkVEqoNVq6BpU/jd76KORBIl\n3lH4Xjaz1kCbcNZCd9+SuLBERCrfihUreOCBB7jkkkto3bp11OFIBShj/joUWBYznQ+cWNI64aiz\n64EmBMVUD+AboB5wjbuvKe5FzCwLyAJIT0+P+5hEpOp6+WXo2TP4efLJUUcjiRJvCxQEXREywm06\nmhnuPjkhUYmIRODuu+/mhx9+IDs7O+pQpGJVZv7qBOwADgH2B/7PzF4vbM2K5e4TgYkAmZmZJY4Y\nKCLVw4oV0K8fHHYYHH981NFIIsVVQJnZFOBw4EOCxABBNwUVUCKSFFatWsWECRPo1asXbdq02fMG\nUi2UMX8tB1rETDcP5xW3Tn7YXW9fgsEkLgFedvdtwAozewvIBH5SQIlI8igcsnzdOnjtNahbN+qI\nJJHibYHKBNq6u66QiUhSatCgAXfeeSennXZa1KFIxSpL/poLtA4HnFhOcJP3S4qsMx24DHgH6Am8\n4e5uZkuBLgQ3ia8P/By4t5zHICJV3AMPwAsvwL33wjHHRB2NJFq8BdSnBPfR+CaBsYiIRKZOnTpc\nffXVUYchFW+v81f4naYBwCsEw5hPcvf5ZjYCyHX36cCjBEXSImANQZEFweh9j5nZfIJBKx5z948r\n7nBEpCr65BPo1k1DlqeKeAuopsBnZvYfYNeXb3UndxFJBuPHj6dOnTpkZWVFHYpUvDLlL3efCcws\nMm9YzPPNBEOWF91uY3HzRSS5PfggbNkCZlFHIpUh3gLq1kQGISISlTVr1nDzzTfTtWtXFVDJ6dao\nAxCR5HX33XDWWdC+PdSuHXU0UlniurVXeEf3hUDD8LGg8C7vIiLV2bhx49iwYQPDhg3b88pS7Sh/\niUiivPwyDB4Mjz4adSRS2eIqoMzsIuA/BN0SLgLeM7OeiQxMRCTR1q1bx7hx47jgggs4Rt/6TUrK\nXyKSCIVDlrdvD6NHRx2NVLZ4u/BlAye4+woAM2sGvE5ww0ARkWpp3LhxrF+/Xq1PyU35S0QqlDtc\nfnkwZPnrr0OdOlFHJJUt3gKqRmHyCa0mztYrEZGqql27dlxzzTV07Ngx6lAkcZS/RKRCPfEEzJwJ\n48cHLVCSeuItoF42s1eAqeH0xcBLiQlJRKRy9OzZk5491ZsrySl/iUiFuvBC2LABBgyIOhKJSryD\nSFwHPAR0CB8T3f36RAYmIpIo33//PXfffTebNm2KOhRJMOUvEakomzfDpk1Qvz786U8asjyVxTuI\nxGHATHe/1t2vJbiil5HIwEREEmXChAkMHjyYBQsWRB2KJJjyl4hUlBtugBNOCIooSW3x9gN/BtgZ\nM70jnCciUq1s2LCBsWPH8utf/5rMzMyow5HEU/4SkXJ76aXgO09nnhm0QElqi7eAqunuWwsnwuf7\nJCYkEZHEuf/++1mzZo1G3ksdyl8iUi7ffRcMWX7MMXDHHVFHI1VBvAXUSjM7r3DCzHoAqxITkohI\nYmzcuJGxY8fSrVs3OnXqFHU4UjmUv0SkzNzhiitg/Xp48kkNWS6BeEfh+wOQY2b3Aw7kA5cmLCoR\nkQT4+uuvSU9P5y9/+UvUoUjlUf4SkTJbtQqWLIExYzRkufworgLK3b8Efm5mDcLpjQmNSkQkAY48\n8khyc3MxDZ2UMpS/RKQ8mjWD3FyoXTvqSKQqiXcUvgPN7FHgGXffaGZtzax/gmMTEakwb731FmvW\nrFHxlGKUv0SkLH74AW65JbjfU506GrJcdhfvd6AeB14BDgmnvwD+nIiAREQqWkFBARdccAGXX355\n1KFI5Xsc5S8R2UvXXw8jR8J770UdiVRF8RZQTd39acKhYN19O8FQsCIiVVZOTg4ZGRnUr1+fFStW\n0LFjx6hDksqn/CUie+XFF2HCBBg0CM44I+popCqKt4DaZGZNCL6Ai5n9HFifsKhERMopJyeHrKws\nlixZsmve2LFjycnJiTAqiYDyl4jE7dtv4fLLoUMHuP32qKORqireAupaYDpwuJm9BUwG/pSwqERE\nyik7O5uCgoLd5hUUFJCdnR1RRBIR5S8RidvVVwffe9KQ5VKaUgsoMzvBzA5y9/eBU4GhwBbgVYKh\nYEVEqqSlS5fu1XxJLspfIlIWo0bBE09Au3ZRRyJV2Z5aoB4CCu/gfhKQDdwPrAUmJjAuEZFyad68\nebHz09PTKzkSiYjyl4jEbd264OfRR8NFF0Ubi1R9eyqg0tx9Tfj8YmCiu/+vu98CHJHY0ERE9p67\nM3ToUHbu3EndunV3W1avXj1GjRoVUWRSyZS/RCQuP/wAv/wlXHNN1JFIdbHHAsrMCm+2ezrwRsyy\nuG7CKyJSWXbu3MmgQYMYPXo055xzDg899BAtW7bEzGjZsiUTJ06kT58+UYcplUP5S0Tict118Nln\ncPbZUUci1cWekshU4E0zWwX8APwfgJkdgUYxEpEqZMeOHVx11VU89thjXHvttYwZMwYzo2/fvlGH\nJtFQ/hKRPXrhBbj//qD16ayzoo5GqotSW6DcfRQwmOBGhCe7u8dst8dRjMysm5l9bmaLzOzGYpbX\nNrOnwuXvmVlGOL+Jmf3bzDaa2YQi28wK9/lh+DggngMVkeR23XXX8dhjj/GXv/xlV/Ekqau8+UtE\nkt833wRDlnfsCKNHRx2NVCd77Mbg7u8WM++LPW1nZmkEX9g9k2DEo7lmNt3dP4tZrT+w1t2PMLNe\nwB0EfdU3A7cA7cNHUX3cPXdPMYhI6hgwYABHHHEEV199ddShSBVR1vwlIqnh44/BLBiyvHbtqKOR\n6iTe+0CVRSdgkbsvdvetwDSgR5F1egBPhM+fBU43M3P3Te4+h6CQEhEp1saNG7nnnnvYuXMnrVq1\nUvEkIiJx69oV8vKgbduoI5HqJpEF1KHAspjp/HBeseu4+3aCfulN4tj3Y2H3vVtM/XREUtK6devo\n2rUrQ4YMITdXDdIiIhKfjz6CRx8Fd6hXL+popDpKZAGVKH3c/Rigc/go9hviZpZlZrlmlrty5cpK\nDVBEEmvVqlV06dKFuXPn8vTTT9OpU6eoQxIRkWrghx/gkkvgllvg+++jjkaqq0QWUMuBFjHTzcN5\nxa4TDje7L7C6tJ26+/Lw5wbgSYKugsWtN9HdM909s1mzZmU6ABGper755htOPfVUFixYwPPPP89v\nf/vbqEMSEZFqYsiQYMjyJ56AffeNOhqprhJZQM0FWpvZYWa2D9ALmF5knenAZeHznsAbMSMl/YSZ\n1TSzpuHzWkB34NMKj1xEqqx58+axfPlyXnrpJc7WTTtERCROM2bAAw/A4MFw5plRRyPVWcJuJuju\n281sAPAKkAZM+n/27juuqvoN4PjnyxQEUXGiAu69d+YemWZmWalkOQpHOzUrzJGZpeXK1KhsKPaz\nUsvcIzUrcZa5tyBiDlT2vHx/fxymQDmAA9zn/XrdF3jOufc+9xb33Od8n+/z1VofUUq9A+zTWq8G\nvgCWKKVOA9cxkiwAlFLngRKAg1LqEaAHEARsTEmebIEtwGd59RqEEAVHbGwsTk5OPPTQQ5w7d45S\npUqZHZIowpRSPYG5GOeaz7XW79+y3xH4BmiOUTnxpNb6fMq+RsCnGOewZKCl1lqaIglhohs3YNgw\naAWa3yAAACAASURBVNIEpk0zOxpR2OXpauxa63XAulu2TczwexzweA739c7hYZvnVnxCiMLh0KFD\n9OrVi0WLFtG7d29JnkSeupdlOFLK0ZcCg7XWB5VS7kBiPr8EIcQtSpWCOXOgeXNpWS7uXWFsIiGE\nsCJ79+6lY8eOaK2pXr262eEI63DXy3BgVEv8rbU+CKC1DtNaW/IpbiFENlKbRfj4QJ065sYiigZJ\noIQQBdavv/5K165dKVmyJDt37qSOnPlE/riXZThqAVoptVEpdUAp9XpOTyLdYoXIe3/9BZ6esH69\n2ZGIokQSKCFEgXTo0CF69uxJpUqV2LlzJ1WrVjU7JCFuhx1wP+CT8rOfUqprdgdKt1gh8lZMjNGy\n3NkZWrY0OxpRlEgCJYQokOrXr8+4cePYsWMHlSrdevFfiDx1L8twhAC/aq2vaa1jMOYBN8vziIUQ\nWYwdC8eOwTffQJkyZkcjihJJoIQQBcoPP/xAUFAQNjY2TJkyhXLlypkdkrA+97IMx0agoVLKOSWx\n6ggcRQiRr376CRYuNJKobt3MjkYUNZJACSEKDH9/f5544gmmTJlidijCiqXMaUpdhuMY8F3qMhxK\nqYdTDvsCcE9ZhuM14I2U+94AZmEkYX8BB7TWa/P7NQhh7X77DZo2lZblIm/kaRtzIYS4XbNmzWLM\nmDH06tWLTz75xOxwhJW7x2U4lmK0MhdCmGTmTIiKAgcHsyMRRZGMQAkhTKW15p133mHMmDH079+f\nVatW4eTkZHZYQgghCqElS+Dvv43fXVzMjUUUXZJACSFMFRUVxbJly3jmmWf49ttvcZDLhUIIIe7C\nn3/C8OFStifyniRQQghTJCcnk5iYiKurK7///juLFy/Gzk6qioUQQty56GgYOBDKloUFC8yORuS3\ngIAAvL29sbGxwdvbm4CAgDx9Pvm2IoTId0lJSQwfPpyEhAQCAgJwd3c3OyQhhBCF2GuvwcmTsHkz\nyCnFugQEBODr60tMTAwAQUFB+Pr6AuDj45MnzykjUEKIfJWQkMCAAQP45ptvqF+/Pkops0MSQghR\niK1fD/7+RsvyrtkuWy2KsrfeeisteUoVExODn59fnj2njEAJIfJNbGwsjz32GOvXr2fWrFm8+uqr\nZockhBCikOvY0Zj3NHas2ZEIM1y4cCHb7cHBwXn2nDICJYTIF1pr+vfvz4YNG/D395fkSQghxD1J\nTjbmPjk7w1tvSctya6G1ZsWKFbRq1YqrV6/i6emZ7XE5bc8NkkAJIfKFUooXXniBpUuX8txzz5kd\njhBCiEJu1ixo1gyuXTM7EpFftm7dSuvWrenfvz9RUVFcvHiRadOm4ezsnOk4Z2dnpuVhO0ZJoIQQ\neerKlSusXLkSgAcffJBBgwaZHJEQQojC7sABY9SpQQNpGmENoqOj6datG926dePy5ct8+eWXHDp0\niCZNmuDj44O/vz9eXl4opfDy8sLf3z/PGkiAzIESQuShixcv0q1bN0JCQmjfvj1ly5Y1OyQhhBCF\nXGrL8nLl4LPPQHoRFV1hYWG4u7tTvHhxKlasyJw5cxg5ciSOjo6ZjvPx8cnThOlWkkAJIfLEuXPn\n6Nq1K9euXWPdunWSPAkhhMgVr74Kp07B1q1QurTZ0Yi8EBwczJQpU/j22285evQo3t7eLFmyxOyw\n0kgCJYTIdceOHaNbt27ExcWxdetWWrZsaXZIQgghioCoKAgMhNdfh86dzY5G5LarV68yffp0Pvnk\nEwCef/55XF1dTY4qK0mghBC5bvny5VgsFrZv307Dhg3NDkcIIUQR4eICe/aAjcziL3KuXbtGjRo1\niIqKYsiQIUyaNClPO+ndC/nfTwiRa5KSkgCYNGkSf/75pyRPQgghckVyMsyYARERUKyYtCwvKuLi\n4tiwYQMAZcqU4Z133uHw4cN88cUXBTZ5AkmghBC5ZNu2bdStW5fTp0+jlKJixYpmhySEEKKI+PBD\nGD8efv7Z7EhEbkhKSuLLL7+kVq1a9OrVi7NnzwLw8ssvU7duXZOj+2+SQAkh7tm6devo1asXDg4O\nFC9e3OxwhBBCFCH794OfHzz2GMhKGIWb1pqVK1fSqFEjhg0bRoUKFdi8eTPVqlUzO7Q7IgmUEOKe\nrFixgkceeYR69eqxY8cOGXkSQgiRa6KjjaSpQgXw95eW5YVdcHAwTz75JFprVqxYwe7du+natavZ\nYd0xaSIhhLhrP//8M0888QRt27Zl7dq1uLm5mR2SEEKIIuT1142W5b/8Ii3LC6t9+/axdu1aJk2a\nhJeXFzt27KBVq1bY2RXeNERGoIQQd61du3aMGjWKjRs3SvIkhBAi1734IixYAJ06mR2JuFPHjx/n\n8ccfp2XLlsyfP58rV64AcN999xXq5AkkgRJC3IXly5cTFxdH6dKlmT9/vsx7EkIIkauio0FrqFMH\nRo40OxpxJy5fvsyzzz5L/fr12bBhA5MnT+bMmTOUK1fO7NByjSRQQojbprVm4sSJDBgwgIULF5od\njhBCiCLIYoE+feC558yORNwJrTUAycnJrFy5kpdeeomzZ88yadIkSpQoYXJ0uatwj58JIfKN1pox\nY8Ywe/Zshg8fzksvvWR2SEIIIYqgDz+EbdvgqafMjkTcjsjISGbPnk1gYCBr166lYsWKXLhwoUhX\np8gIlBDiP1ksFkaMGMHs2bN56aWX8Pf3x9bW1uywhBBCFDH79sGECdC/PwwdanY04t/Ex8czd+5c\nqlevzqRJkyhWrBjR0dEARTp5AkmghBC34fTp0yxbtoy33nqLOXPmYGMjHx1CCCFyV1SUtCwvLPbv\n30+tWrV45ZVXaNCgAYGBgaxcuRIXFxezQ8sXUsInhMiRxWLB1taW2rVrc+TIEby8vMwOSQghRBF1\n8CBcuQI//QSlSpkdjbiV1porV65Qvnx5atSoQe3atfn888/p1q0bysqyXbmMLITIVkxMDL169WLe\nvHkAkjwJIYTIU+3aQVAQdOxodiTiVtu2baNNmzZ069YNi8WCm5sbmzZtonv37laXPEEeJ1BKqZ5K\nqRNKqdNKqTey2e+olFqesn+3Uso7Zbu7UmqbUipKKTX/lvs0V0odSrnPPGWN/9WEyGMRERH07NmT\nLVu24OrqanY4QgghirALF+Dzz4225bKkYMGyb98+evToQZcuXQgNDeXVV181O6QCIc8SKKWULfAJ\n8CBQDxiolKp3y2HDgRta6xrAbOCDlO1xwNvA2GweeiHwHFAz5dYz96MXwnqFhYXRtWtXdu3axbff\nfstQmcUrhBAij1gsMHgwvPIKhIaaHY3IaPXq1bRs2ZIDBw7w0UcfcerUKYYNGyZNpMjbEahWwGmt\n9VmtdQLwP6DvLcf0Bb5O+f0HoKtSSmmto7XWv2EkUmmUUhWBElrrQG00m/8GeCQPX4MQViUmJoZO\nnTpx6NAhVq1axRNPPGF2SEIIIYqwGTNgxw74+GOoVMnsaERISAh//PEHAN27d2f69OmcPXuW1157\njWLFipkcXcGRlwlUJeBChn+HpGzL9hitdRIQDrj/x2OG/MdjAqCU8lVK7VNK7bt69eodhi6EdXJ2\ndmbgwIGsXbuWhx56yOxwhBBCFGF79sDEifD44zBkiNnRWLewsDDGjh1LjRo1GDZsGFprnJyceOON\nN4rcIri5ocg2kdBa+2utW2itW5QtW9bscIQo0E6fPs2BAwcAeOutt+jatavJEQlhrrudw5thv2fK\nPN7sStGFsHpxcUbL8ooV4dNPpWW5WaKiopg6dSrVqlVj9uzZDBw4kI0bN1plY4g7kZdtzC8CVTL8\nu3LKtuyOCVFK2QFuQNh/PGbl/3hMIcQdOHLkCN26dcPNzY0jR45IbbOwehnm8HbHqHTYq5RarbU+\nmuGwtDm8SqkBGHN4n8ywfxawPr9iFqKwKVYM3n3XKNuTluXm+eGHH5g4cSL9+vXj3XffpV69W9sV\niOzkZQK1F6iplKqKkeQMAAbdcsxq4BlgF9Af+CVlblO2tNaXlFIRSqk2wG7gaeDjvAheCGuwf/9+\nHnjgARwcHFi5cqUkT0IY0ubwAiilUufwZkyg+gKTU37/AZifModXK6UeAc4B0fkXshCFR2wsODnB\ngAFmR2J9LBYLS5cuxdbWlqeeeoqnnnqKBg0a0KJFC7NDK1TyrIQvZU7TC8BG4Bjwndb6iFLqHaXU\nwymHfQG4K6VOA68BaWUSSqnzGFfwhiilQjJ08BsNfA6cBs4gV/iEuCu///47Xbp0wcXFhZ07d8pV\nJyHS3fUcXqWUCzAemPJfTyJzdYU1Cg6GatXghx/MjsS6aK358ccfadSoEUOGDGHZsmUA2NnZSfJ0\nF/JyBAqt9Tpg3S3bJmb4PQ54PIf7euewfR/QIPeiFMI6TZ8+nQoVKrBlyxaqVKny33cQQtyOycBs\nrXXUf80h0Fr7A/4ALVq0yLH6QoiiwmKBp56CqCho2tTsaKxHYGAgr7zyCrt376ZWrVp8//33PPbY\nY2aHVajlaQIlhCh4tNYopfj222+JiYmhfPnyZockREFzL3N4WwP9lVIzgJJAslIqTms9HyGs3Pvv\nw86d8NVXUL262dEUfann+9DQUEJCQvjss88YMmQIdnby9f9eFdkufEKIrJYvX06XLl2Ijo7G1dVV\nkichspc2h1cp5YAxh3f1LcekzuGFDHN4tdbttdbeKVUUc4D3JHkSAnbvhkmT4Mkn4emnzY6maDt5\n8iRPPvkkU6dOBaBfv36cOnWKZ599VpKnXCIJlBBWYvHixQwaNAiLxYLFYjE7HCEKrHudwyuEyGrD\nBqPj3qJF0rI8r4SEhODr60u9evVYu3Ytjo6OACilcHJyMjm6okX9S9O7IqNFixZ63759ZochhGnm\nzZvHyy+/TI8ePVi1ahXOzs5mhySskFJqv9ZaZitnQ85TwhrcuCEty/PKp59+yiuvvILFYmHUqFH4\n+flRrlw5s8MqdG73PCUjUEIUcXPnzuXll1+mX79+rF69WpInIYQQ+WbdOvjzT+N3SZ5yV1RUFOHh\n4QDUrl2bJ598kpMnTzJ37lxJnvKYJFBCFHEdO3Zk1KhRfPfdd2nD+UIIIUReCwqCQYPg1VfBCgqe\n8k1CQgLz58+nevXqTJo0CYBOnTrx1Vdf4e3tbW5wVkISKCGKoOTkZNasWQNAkyZNWLBggUwcFUII\nkW8sFhg8GJKTYfFimfeUGywWC0uWLKF27dq8+OKL1K1blwGyGrEpJIESooixWCw8++yz9OnTh61b\nt5odjhBCCCs0fbrRsvyTT4yFc8W9e/HFF3n66acpVaoUGzZsYNu2bbRp08bssKySXJIWoghJTExk\n8ODBLF++nEmTJtGlSxezQxJCCGFl9u2DyZNh4EBj4Vxx93bs2IG3tzdeXl6MHDmSTp060b9/f2xs\nZAzETPLuC1FExMXF8dhjj7F8+XJmzpzJ5MmTUVIzIYQQIp/Vrw9vvgkLFkjp3t06cOAAPXv2pFOn\nTnz00UcANGrUiCeeeEKSpwJA/gsIUURs3LiRNWvWsGDBAsaOHWt2OEIIIaxQXBw4OcHUqVCypNnR\nFHwBAQF4e3tjY2ODt7c3H330EQMGDKB58+bs2bOHGTNm8MEHH5gdpriFlPAJUchprVFK0bdvXw4f\nPky9evXMDkkIIYQV+t//YMIE2LoVvLzMjqbgCwgIwNfXl5iYGACCgoIYP348NjY2+Pn5MXbsWEpK\nFlogyQiUEIXYtWvX6NixI7/99huAJE9CCCFMERQEI0dCuXJQqZLZ0RQOfn5+aclTKovFQvny5Xn3\n3XcleSrAZARKiEIqNDSU7t27c/bsWaKioswORwghhJVKSjKaRSQnQ0AAyKoZ/+3mzZsEBQVlu+/i\nxYv5HI24UzICJUQhFBQURIcOHQgODmb9+vX07NnT7JCEEEJYmYAA8PYGe3v47Tfw8YGqVc2OquCb\nMWPGvy546+npmX/BiLsiCZQQhUxQUBDt27cnLCyMLVu20KlTJ7NDEkIIYWUCAsDX1yjdS/XNN8Z2\nkVV8fHza7+fPn6djx4689957ODs7ZzrO2dmZadOm5Xd44g5JAiVEIVOhQgU6duzI9u3bad26tdnh\nCCGEsEJ+fnDL9B1iYoztIl18fDzz58/H29s7bb7yxx9/zE8//cSbb76Jv78/Xl5eKKXw8vLC398f\nHx8fk6MW/0WqVIUoJA4cOICnpydlypRhyZIlZocjhBDCSl2/nnnkKaPg4PyNpaBKTEzkq6++YurU\nqVy4cIEOHTqkjTbZ2tqmHefj4yMJUyEkI1BCFAK//vornTp1YuTIkWaHIoQQwkppDUuWQJ06OR8j\n03cgOTmZFi1a4Ovri4eHB5s3b2b79u00a9bM7NBELpEESogCKuPiep06dcLFxYW5c+eaHZYQQggr\ndOkSdOsGTz8N1arBe+/BLdN3cHYGa52+k5yczIYNG9BaY2Njw+jRo1mzZg27du2iW7duKKXMDlHk\nIkmghCiAUhfXCwoKQmuN1pqbN2+yfft2s0MTQghhhdzcjNK9hQvhjz/gzTfB399YMFcp46e/v9GJ\nz5porVm9ejVNmzblwQcfZMuWLQCMGDGC3r17S+JUREkCJUQBo7XOdnG92NhY/GR2rhBCiHyyfTs8\n9BDExRmjS/v3G4vl2qR8e/TxgfPnjfWfzp+3ruRJa82mTZto3bo1ffv2JSYmhqVLl9KlSxezQxP5\nQBIoIQqII0eO8Oqrr1KjRg2Cc5iFm9N2IYQQIrdcuwZDhkDnznD0aHrDCBv51pgmOjqagQMHcvny\nZb744guOHTuGj49PpgYRouiSPwUhTBQdHc2XX37JfffdR4MGDfjkk09o3rw5Hh4e2R4vi+sJIYTI\nK1rDl19C7drGek5vvgmHDxv/FrB7925GjBiBxWLBxcWFLVu2cPLkSYYNG4adnTS2tiaSQAmRz7TW\nJCQkALBnzx6GDRvG9evX+fDDD7l48SLfffcdH3zwgSyuJ4QQIl8lJcFHH0HduvDXX9k3irBGf/75\nJ3369KFNmzasWrWK06dPA9C0aVMcHR1Njk6YQRIoIfJJeHg4CxYsoFmzZrz++usAdOrUiT/++INj\nx44xZswYypYtCxjrQsjiekIIIfJabKyRKIWHg709bNkCv/4K9eubHZn5rly5wuOPP06zZs347bff\neO+99zh79iy1ZUjO6sl4oxB57I8//sDf35/vvvuO2NhYGjduTNOmTQFQStG2bdts7yeL6wkhhMhL\nmzfDqFFw5gxUrmy0KK9QweyozBcXF0exYsUoUaIER48eZeLEibz66quULFnS7NBEASEJlBB5IDw8\nHDc3NwBmzZrFpk2bePrpp3n22Wdp3ry5tDUVQghhmsuX4bXXYNkyqFnTGHXq2tXsqMwXFBTE1KlT\n2b59O0eOHKFYsWL8/fff0hhCZCEJlBC5JDk5ma1bt/L555/z448/cujQIWrVqsWcOXMoWbIkLi4u\nZocohBBCMGyYkTRNnGg0iihWzOyIzBUaGsp7772Hv78/SilGjRpFfHw8jo6OkjyJbEkCJcQ9un79\nOgsXLuSLL77g3LlzlC5dmtGjR+Pk5ARA5cqVTY5QCCGEtTt0CMqVg/LlYdYso+NenTpmR2W+ffv2\n0b59e5KSkhg+fDh+fn5UqVLF7LBEASdNJIS4C0lJSVy6dAkwaqUnT55M1apV+fbbb7l48SKzZ8+W\nD2AhhBCmi4mBN96AZs1g0iRjW+3a1p083bhxg99++w2AJk2a8NJLL3HixAkWLVok525xW2QESog7\ncPbsWRYvXsyXX35JvXr12Lx5Mx4eHoSEhFC+fHmzwxNCCCHSrF8Po0fD+fNG2Z61r4QRERHBnDlz\nmDVrFo6Ojly4cAEHBwc++OADs0MThUyejkAppXoqpU4opU4rpd7IZr+jUmp5yv7dSinvDPveTNl+\nQin1QIbt55VSh5RSfyml9uVl/EKkWrduHd27d6d69epMnz6dpk2b8sILL6Ttl+RJCCFEQTJzJvTq\nZcxv2rEDvvgC3N3Njsoc0dHRzJgxg6pVqzJp0iQ6d+7Mli1bcHBwMDs0UUjl2QiUUsoW+AToDoQA\ne5VSq7XWRzMcNhy4obWuoZQaAHwAPKmUqgcMAOoDHsAWpVQtrbUl5X6dtdbX8ip2IQCOHTtGjRo1\nsLe3Z/fu3Zw8eZIpU6YwdOhQGeIXQghR4FgsEBEBpUpB//6QkABjx4K1r/W6detWxo8fz4MPPsg7\n77xDixYtzA5JFHJ5OQLVCjittT6rtU4A/gf0veWYvsDXKb//AHRVRn/nvsD/tNbxWutzwOmUxxMi\nT8XExPD1119z//33U69ePdauXQvA+PHjOXv2LBMnTpTkSQgrcLcVFEqp7kqp/SmVEvuVUl3yO3Zh\nnf78E9q2haeeMhpEVK0Kfn7WmTwlJibi7+/PvHnzAOjTpw979+5l3bp1kjyJXJGXCVQl4EKGf4ek\nbMv2GK11EhAOuP/HfTWwKeXE5JsHcQsrdPPmTUaPHk3FihUZMmQIV69eZcaMGbRr1w4AZ2dnaWUq\nhJXIUEHxIFAPGJhSGZFRWgUFMBujggLgGtBHa90QeAZYkj9RC2sVFWWs6dSiBQQFGQmUtUpKSuLr\nr7+mdu3ajBgxgrVr16K1RikliZPIVYWxC9/9WutmGCe255VSHbI7SCnlq5Tap5Tad/Xq1fyNUBQK\n4eHh7NtnTKMrXrw4GzZs4OGHH2bHjh0cP36ccePGUbZsWZOjFEKY4K4rKLTWf2qtQ1O2HwGclFJW\nOAYg8sOePVCvHsyeDc89B8ePw8CBYI1rtW/dupUGDRowZMgQSpUqxdq1a9mwYYMsXC/yRF4mUBeB\njLVOlVO2ZXuMUsoOcAPC/u2+WuvUn1eAVeRQ2qe19tdat9Bat7jbL8EBAQF4e3tjY2ODt7c3AQEB\nd/U4ouDQWvPHH38wdOhQPDw86NevHxaLBXt7e06ePMmSJUvo0KGDfOAKYd3upYIio8eAA1rr+Oye\nRC70ibultfHT0xO8veH332HRImPukzXRWhMbGwuAnZ0ddnZ2rFy5kn379tGrVy85l4s8k5cJ1F6g\nplKqqlLKAaMpxOpbjlmNUeIA0B/4RWutU7YPSKkxrwrUBPYopYorpVwBlFLFgR7A4bwIPiAgAF9f\nX4KCgtBaExQUhK+vryRRhdiqVato0KAB7dq144cffsDHx4cVK1ZgY2P8GdjZSVd/IUTuUErVxyjr\nG5HTMblxoU9Yl6QkmDMHHnwQkpOhQgX49Ve47z6zI8tfWmvWr19Pq1atGD9+PAAdO3bk77//pl+/\nfpI4iTyXZwlUyhW5F4CNwDHgO631EaXUO0qph1MO+wJwV0qdBl4D3ki57xHgO+AosAF4PqUDX3ng\nN6XUQWAPsFZrvSEv4vfz8yMmJibTtpiYGPz8/PLi6UQeSE5O5pdffuHy5cuAManU1dWVzz//nEuX\nLuHv70+rVq3kg1YIcat7qaBAKVUZo0Liaa31mTyPVliFffugdWt49VWwsTG67Vmjbdu2cf/999Or\nVy+uXr1Ky5Yt0/alXhAV1ifgUADec7yxmWKD9xxvAg7l7YCH0qnjwEVYixYtdOpcl9tlY2NDdu+N\nUork5OTcCk3kgUuXLvHll1/yxRdfcPbsWaZPn84bb7yRNpFUCGEOpdR+rXWBn8mdkhCdBLpiJEp7\ngUEpF/dSj3keaKi1HpmyDMejWusnlFIlgR3AFK31ytt9zrs5TwnrEBkJEybA/PlQrhzMm2e0KLfG\n09nbb7/Nu+++i4eHBxMmTGD48OGylpMg4FAAvj/7EpOYPvDhbO+Mfx9/fBr63NFj3e55SmqWcuDp\n6UlQUFC220XBlJCQwBNPPMGaNWuwWCx06tSJd955h0cffRRAkichxG3RWicppVIrKGyBxakVFMA+\nrfVqjAqKJSkVFNcxytTBqLyoAUxUSk1M2dYjZd6uEHcsORlWrIBRo2DaNHBzMzui/HXgwAHc3d3x\n8vLisccew93dnREjRuDk5GR2aMJEZ2+c5fT104RGhvLqhlczJU8AMYkx+G31u+ME6nbJCFQOUudA\nZSzjc3Z2xt/fHx+fvPmPIe7c+fPnCQwMZMAA47vLE088QdWqVXn22WepWbOmydEJITIqLCNQZpAR\nKJHR+fMwaxZ89BHY2xujUK6uZkeVvw4fPsykSZNYuXIlvr6+fPrpp2aHJPKQ1pobcTco7VQagO3n\nt7Prwi5CI0MJjQrlUuQlbsTd4OjooyileHrV0yz5+99XiVAokifdWdWYjEDdo9Qkyc/PL20kavbs\n2ZI8FQAJCQmsXr2azz77jM2bN+Pg4MCDDz6Im5sb3333ndnhCSGEEHclMdFoEjF5slGiN3gwtGxp\nXcnTyZMnmTJlCt9++y2urq5MnjyZV155xeywxF3SWnM99jqXoi4RGhlK92rdUUqx7NAyVh5baSRI\nkaFcirpEgiWB+AnxONg6sOLoCubvnU/JYiXxcPWgoktFapSuQWJyIg62Doy9byy+zX3xcPWg01ed\nuBBxIctze7rlXdWYJFD/wsfHBx8fH/7++28aN25MdHS02SFZvdWrV/Pss89y9epVqlSpwsSJExk2\nbBhu1lbTIEQ+CggAPz8IDjbaJk+bBnItSYjctWsXjBgBhw5B377GXCdrnDUwdepUfvzxR9544w3G\njh1L6dKlzQ5JZCN1xCgtAYo0EqRX2ryCk70T83bPY3bgbC5FXiLekr6Sw9VxVynjXIbg8GCOXj2K\nh6sH7b3a4+HigYerB5ZkC9jCu13eZUb3GTjZZ1+q2ah8o7Tfp3ebnu0cqGldp+XZ65cSvtsUEBDA\ngw8+KH/I+SwmJoYVK1ZQr149mjdvzuHDh5k4cSLPPfccPXr0wNbW1uwQhSjSAgLA1xcyNiV1dgZ/\n/ztPoqSEL2dSwmfdLBaoWxdiY+Hjj+GRR8yOKP+EhIQwbdo0Ro4cSePGjQkNDcXOzo5y5cqZannp\nNAAAIABJREFUHZpVuxR5icNXDmcaIQqNDGVh74WULV6Wd399l7e3vZ3lfmdeOkO1UtX4/sj3rD65\nGg8XDyq6VsTD1UiQWnq0xNEu99cWDzgUgN9WP4LDg/F082Ra12l3Nf/pds9TkkCJAungwYN89tln\nLF26lPDwcF566SXmzp1rdlhCWIWbN+GHH+D4cViwwPhSdysvL2Oexp2QBCpncp6yPlrDqlXQs6dx\nUeLoUahSxXrK9S5fvsz777/PwoULSU5OZv78+fj6+podVpGVOmLk4uCCg60Dhy4fYu2ptVkSpA0+\nG6hdpjZzA+fyysb00kk3Rzc8XD34acBP1HSvyZ6Le/jjwh9p5XUerkai5GzvbOKrvHcyByoPLF++\nnGPHjjF58mSzQymytNb06NGDLVu24OjoyGOPPcZzzz1Hx44dzQ5NiCIjORmCgowE6dgx4+fx4zBg\nAIweDdHR8Nxz4OgI8fHZP0ZwcP7GLERRcvas8be2cSPMng2vvAL16pkdVf6ZPHkyM2fOJC4ujmee\neYa3336bqlWrmh2WKe515CRjKd2lyEs092hOaafS7Azaydzdc9OSo9RSul3Dd9Gmchv2X9rPm1vf\nxM3RLW2E6H7P+7GzMVKDR+s+SrOKzXJMjFpVakWrSq1y9b0oTCSBugOBgYF88sknjBo1ivLly5sd\nTpGgtWb37t2sXr2aadOmoZSie/fu9OnTh6eeekpKJoW4B7GxcPJkeoJUpw48+aSRIFWrln5cmTJG\n+ZBzyvnRwwPOnDFGmapXN5KtW1nj3Awh7lVCgtFZ7513jO568+YZiZQ1iI6Opnjx4gBcv36dRx55\nhEmTJlGrVi2TIzPPresXBYUH4fuzMQo3qMGgTIlRaiLUv15/apSuwaYzmxixZkSWOUYbfDbwQI0H\niEyI5MjVI3i4etCuSru0ErrKJSoD8GT9J3m83uMUdyiebWxV3KpQxa1KtvuElPDdkePHj1O3bl1K\nlixJeHg4np6eTJs2TTrz3YXr16+zZMkSPv/8cw4fPkzx4sU5ePAg1atXNzs0IQoVreHaNSNBsrGB\ndu2MbXXrGslT6ke8UsaoUmon4KVLoWpVqF3bSKByInOg8oeU8FmHQYPg22/hscdg7lyoVMnsiPJe\ndHQ0H3/8MTNnzuTHH3+kffv2srB9Cq85XgSHZx3O93LzYkHvBfRe1jvLvu8f/57+9frz1z9/8dGu\nj9LL51J+NirfCLdi0ljrbkkJXx7Yv38/NjY23Lx5E4CgoKC0el1Jom7f1q1b6d27N/Hx8bRs2ZJP\nP/2UgQMH4mothd9C3AWLBW7cSE92xoyB3buNErzr141tnTvDL78YydJDD0GJEsaoU506ULMmZFx3\n8qmnbu95Uz/apAufEHfn+nXj4kbJkvDaazBwIPTpY3ZUeS8uLo5FixYxffp0rly5Qu/evdOqSqw1\neboSbaynXa54OQ5cOpBt8gQQHB5Mo/KNmP3A7Ezziyq6VEwbMWpSoQlL+v37Okgi78gI1B3w9vZO\nWxMqIy8vL87f6WzqIiwgIAA/Pz+Cg4Px9PTk9ddfJzIyEi8vLwYMGEBkZCQTJkxg2LBhNG7c2Oxw\nhSiQ1qyBPXvSy+9OnoT69WH/fmP/ww9DeHh6glS3rjGHoiCX1skIVM5kBKro0doY6R0zBh59FBYt\nMjui/GOxWKhfvz4nTpyga9euTJ06lbZt25odVr6yJFs4cOkAgSGBBF4MJDAkkLM3zjKp4yQmd5pM\nRHwElWZVIiohKst9vdy8OP/K+fwPWkgXvoxy68RkY2NDdu+XUork5Dtb6bioCggIwNfXl5iM9T4p\nZCVxIdJdvQpHjqQnSMeOGQ0btm839j/4IGzaZMxVSk2QGjcu3CM/kkDlTBKoouXkSRg1yhgRbtPG\nKJ1t1Oi/71fYZLxgWqVKFR599FFmzZqFUorFixdTtWpVOnfubHaYeU5rTUhECIEhgRR3KE6vmr2I\nSYzB7X03kpKTqOhSkbZV2tK2clseqP4ADcs3BLLOgQJj/SL/Pv531YJb3Dsp4csDnp6e2Y5AVahQ\nwYRozKW1JiwsjLNnz3L27FnOnTtHeHg4//vf/7JNnjw8PCR5ElYnKQnOnUtPkIKDYf58Y9+LL8Ly\n5cbvzs5GktSggXHVWin46iuj5Mcx95fLEELkoaVLYfhwo2R24UJjDqGNjdlR5b5bL5gGBwczZ84c\nnJyceO+99xg2bJjJEea9j3d/zPag7QSGBBIaGQpAt2rd6FWzF872zqwdtJa6ZepSuUTlbMsWU5Ok\n3Fi/SOQvGYG6AzmNrtjZ2TFnzhxGjx5dpOp6Y2NjOX/+fFqCdPHiRaZPnw7AkCFD+PrrrzMdX6lS\nJUJDQ2WUTlidqCg4cSK9FbitLUyeDNOnG123UlWoAKdOgYuLMX8ptQSvcuWi+QXrVjIClTMZgSr8\nkpLAzs64WDJtGnz4ofE3X1RVqVKFkJCQLNtzuthcWGmtOXPjjFGKFxJITGIMi/suBuD+xffzT9Q/\ntKncJu3WqHwjHGwdTI5a3C0p4csgN09M2c3vWbNmDevXr6dnz5489NBDzJw5M21/Qe7Sl5ycTGho\naFqCdO7cOSZMmICdnR2vv/46M2fOzHS8s7Mzly9fxsXFhbVr13L69GmqVatG1apV8fb2xsXFReaJ\niSJLa/jnH6OJg709/PQTfPKJkTRduJB+3JkzRtndzz/D77+nl9/Vrm2MKFkzSaByJglU4XXtGowd\naywbkDqqXJRcu3aN3bt3ExgYSGBgIC1atGD69OlFdlpDdEJ0WqOGCb9MYNG+RYTFhgHg4uBCe8/2\nrB20FqUUcUlxFLMrZma4IpdJCV8e8fHxyZIQjRo1ioULF/Lyyy+zcePGtA+UgtClLzw8PC1BOnv2\nLCNHjsTFxYWZM2cyYcIEEjJcHldKMXz4cKpUqUKXLl0oUaJEWoJUrVo1ypUrlzbC1rt31taaANOm\nTcsySufs7My0adPy9oUKkctOnIBVqzIvNhsRAQcOQNOmxlpKN25Ax45GgpTazCG1iUOfPtbRaUsI\na6W1UWo7bpwxmjxunLFIdWEeTU5MTOTSpUt4pnyQtWrVir179wJga2tLo0aNKFWqFJDzSJNnQe5k\ncwtLsoVj146ljS4FhgRyMuwkN8bfoLhDcSq6VKRv7b5po0v1ytbD1sY27f6SPFkxrXWRvzVv3lzn\nBw8PDw1kuZUuXVr/8ssvOigoSFssln99jKVLl2ovLy+tlNJeXl566dKl/3p8fHy8PnXqlN60aZNe\ntGiRvnbtmtZa6y+//FKXLl06SywHDx7UWmu9efNm/frrr+uFCxfqjRs36lOnTun4+PhceR/u9DUI\nkVeWLtXay0trpYyfGf9XjInRes8erb/5Rus339S6Xz+t69bV+tdfjf0//KA1aF2pktZdu2r9/PNa\nf/yx1pcumfFKigZgny4A54SCeMuv85TIHWfOaN2xo/EZcf/9Wh8+bHZEdyckJER///33esyYMbpd\nu3a6WLFiukGDBmn7p0yZoj/44AO9Y8cOHRUVlem+S5cu1c7Ozpm+Yzg7Oxfoc/7V6Kt6zYk1OiIu\nQmut9bs73tVMRjMZXfqD0rpXQC/9zvZ39PWY6yZHKsxyu+cpKeHLRTkNZ2f0/vvvM378eC5fvsx7\n771H9erV0267d+9m9OjRWUZvPvzwQ5o2bcrZs2fp0qULFSpUYPXq1bz00ktcuHAh01D5tm3b6NSp\nEzt37mTZsmVUq1Yt7Va1alVKWnsNkbAa2S0Aa2dnrGc0eTJs2QLdu6dvr1HDGEkaPx5atzbKcRIT\njbWURO6QEr6cSQlf4RIcDO3bw8SJMHRo4Rh1iomJ4cCBAxw5coQRI0YA8Mgjj/DTTz/h6OhI8+bN\nad26Nffddx/9+/e/rce8dVpDQZu2EBoZyspjK9NGl87cOAPAxqc20qN6D45ePcr+0P20qdyGGqVr\nFKl57OLuyByoDPLrxJTT/B8PDw++/vprzpw5Q9u2bWnUqBF79+6lS5cuREVl7f//b3788Uf69u3L\nnj17mDdvXqbkqFq1anh4eGBra/vfDyREEZGQAOfPG3OPTp+Gli2NtsEeHnDpUtbj3d2NOQs3bxot\nw+vWNeYs2dvnd+TWRxKonEkCVfBt3gzff2+0JFfKuMBS0D83du7cyfLlywkMDOTgwYMkJSUBEBYW\nRunSpdm3bx9aaxo3boyDQ+FufJDaRjwwJJBH6jzC/Z73szNoJx2+6pDWRrxNpTa0rtyaFh4tcLZ3\nNjtkUQBJApVBfp2YsuvS5+zsjL+/f7ZXZLTWXL16lTNnznDmzBkGDx6c42P//PPPaclSsWJScyus\nS3S0kSCdOQNVqkCLFkZDh/vug6AgY95BqgkTYOpU44pwdh9vSmU+XuQfSaByJglUwXX5Mrz2Gixb\nBjVrws6dUL682VFlFh4ezp49e9IaPXz22Wd4eHgwe/ZsJk6cSKtWrWjTpg1t2rShdevWlCtXzuyQ\nc8WV6Cs8v+55dl3YxcXIiwA42Dow+4HZjG45mvikeK5EX8mxjbgQt5IEKoP8PDHdy3C2dLAT1kpr\nuH7dSJCcnKBhQ6OErnt3Y9s//6QfO3q00f0uKQmGDDFGj2rUgOrVjZ/lyhlJkre3kVzdysvLGLES\n+U8SqJxJAlXwJCfD558bZb0xMfDmm/DGG2D2NUyLxYLFYsHBwYFt27bxwgsvcOzYMWNehlLUq1eP\nb775hmbNmhEbG4uDg0OhrkzRWnP2xtm00aVdIbvoWrUrH3T/gARLAo0XNaZphaZpjR4al2+Mo50s\noCfujnThM0l2Xfpul3SwE0VZcrJRUhcXZyQ7WsNTTxnd7s6cMUrqAAYNMuYvFStmzD/q1Ss9OUr9\nCca8paVLc36+adOyzoFydja2izsTcChAFnoUVufGDSNpatIEFi0yliIww5UrV9JGlnbv3s2ePXv4\n9NNPGTRoEGXKlKFq1aoMHDiQNm3a0LJlS9zc3NLu6+TkZE7Q9yAiPoJ/ov6hlnsttNZUn1edczfP\nAVDcvjitKrWiRmnjROBg68Cx54+ZGa6wUpJAFSCpiVdBnpApxL9JSjISoTJljH9PmgQHD6aX38XG\nGgnR2rXGKFFoqHFs69bpyVH9+sZ9lYJ16+4+ltQ/Gz8/Y8K3p6eRPMmf051ZcnAJI9aMIDYpFoCg\n8CB8f05ZnkGSKFHExMQYo07PP2/Ml9yzxxjlzq/qr4SEBA4ePEjx4sWpV68eJ06coE6dOgDY2dnR\nuHFjnnnmGWqkXElq2LAha9asyZ/g8sjxa8f5Pfh3Y4TpYiBHrhyhacWm7Pfdj1KKEc1HUMqpFG0q\nt6F+2fqZ2ogLYRYp4RNC3JGkJGP0B2DxYti/P72BQ1AQtGplLCALxhyliIjMI0iNG0O7dubFb02u\nx17nYsRFIhMiiYyPJCI+gsiESAY2GIiTvRPrTq1j9YnVWfZvGbyFUk6lmLRtEu/8+k62j+3l5sX5\nV87fUTxSwpczOU+Zb/16o0T4/HmjYUS3bnn/nBaLhZUrV6aNMO3fv5/4+Hh8fX359NNPSU5OZs6c\nObRq1YpmzZrh7Fy4Gx9ci7nG7pDdnLp+ilfavAJAr4BerD+9nlLFSqWV4bWr0o6u1bqaHK2wRlLC\nJ4S4Z5s2wb596SNIp09DqVJw6JCx/5tv4K+/jOSoRQt48klo1Cj9/n/8YU7chVGCJSEtialcojL2\ntvacDDvJX//8ZSQ28ZFEJhj7J3SYQMliJQn4OwD/A/5p+1MToFMvnqJyicp8vPtjJu+YnOW5ulXr\nhqebJ0euHGHV8VW4OrhSwrEEro6ueLh6kJRsdOrq6N0Rfs0+3uDw4Dx8N4TIWwEB6aPTHh5QqZIx\n2lSnjtGds2PH3H/O6Oho9u/fT2BgIE5OTrz44ovY2NgwevRooqKiaN68OS+88AJt2rThvvvuA4zl\nUV577bXcDyaX3E5574bTG1h2aBm7QnZx+vppAOxt7Hm22bO4OLgwvet05vScQ83SNaXRgyg0JIES\nwoodOgR//pmeHJ05Y5Sw/P23sX/hQvjxR6hQwRg96tYN6tVLv//69cZcpYJ6zsvruTuJlkTC48Mz\nJS+R8ZF08OpAcYfi7A7ZzfrT6zPtj4iP4IuHv6Cia0UW7l3IxO0TiYyPJN4Sn/a4p188TfXS1fnx\n+I+M3zI+03M62TkxqsUoShYriUZjo2yoUqKKkQA5uOLq6EoxO2OWe/96/alfrn6mBKmEYwk8XD0A\nGNduHOPajcvx9XWp2gUvNy+CwrN24/B088yNt1CIfHfrGnEXLxq3/v2NeZWOudx/YNKkSfz888/8\n/fffWCwWAHr06MGLL76IUopdu3bh6elZ6NqIBxwKwPdnX2ISjTcyKDyI4T8N59tD3xIRH8GKJ1ZQ\ntnhZDv5zkM1nN9O2cluea/YcbSq3oXnF5hR3KA5A4wqNzXwZQtwVKeETogj75x84ciRzgnTxIuza\nZSQ9Q4bA118bLb89PdNL7RYsMLb98w+4uBi3wubWkzuAs70zn/T6hMGNBmNrY8vlqMucDDuZKbmJ\njI/kmSbPUMa5DBtPb2TxX4szjf5Exkey8amN1HSvyaxdsxizaUyW5z486jD1y9Vn3u55vLzhZYrb\nF09LXlwdXPnu8e+oVqoam89sZuWxlZmSG1cHV/rV7UfJYiW5HHWZazHXcHV0TUuO7Gzy97pXTu+j\nfx//O05GpYQvZ3KeyjsnThijSufPG2XGK1YY68fd6l46dN64cSNTG/HLly9z4MABAIYPH05wcHBa\nG/FWrVpRtmzZu305+S4+KZ7QyFCuxlzlWsy1tNtHuz4iNDI02/u0rdyWzx/+nHpl65FoScTOxk5G\nl0ShICV8QhRyGUtMcmqAEB8P585lLrE7e9b4guDoCO+/D3PnGsc6OBiToatXN9ZVcnGBt9+Gt94y\nWn5nd/GzQoU8f5n35HrsdU6GnSQsJoyw2LC0n8+3fB6/rX6ZvvQDxCTGMPSnoTSt0JTGFRqz4tgK\nnl/3fJbH7eTdiTLOZbgWc42D/xxMS26qlqxKCccSONgab1aXql2Y13NeWgKUOtJTtVRVAEa1GMXz\nLZ/PcdJz9+rd6V69e46vr7xLecq7mLvgTGqSJF34REF15gxs3JieIKX+3LTJKCnevh1GjjTmbnp6\nZp88gfFZezuSkpI4evQoDRs2RCnFuHHj+PDDDwFQStGgQQNat25NYmIi9vb2fPHFF7nxMnNFbGJs\nWgLk6eaJu7M7p6+fZsnBJcb22PQEafYDs+lStQubzmzi4f89nOWxFNknRArFH8PT67ftbQv4asNC\n3AVJoP7F7XyBFSIv3FpiEhQEw4bBzz8brb3ffx9Kl4bp02HKlPT7ubgYI0hhYUZd/3PPQd++RtJU\nqRLcuhRI9er595pyEpMYkyUB6uzdmbLFy7Lrwi4W7V+Utv167HXCYsJY57OOVpVa8ePxHxm+enim\nx7NRNvSu2ftf5+iULW5c/e1Tqw+13WtnGv1xdXTFxcEYcvNp5INPo5z/6JtUaEKTCk1y3F9Uvjj4\nNPSRhEnkO62NkfLz543PvlsTpB9+MOYq7dtndM1zcDBGkby8oE+f9PWanngCeveGihWNz8Cc1ojz\nzKEqNSwsjN9++y1tdGnv3r1ER0dz7tw5vL296dy5M6VKlaJNmza0aNGCEiVK5M0bkoOI+AgOXDqQ\naXToWsw1fBr60LJSS3Zd2MWAFQO4FnMt00WlZY8uY2DDgYREhDD116m4O7tTxrkMZZzLUK1UtbRS\n4GYVm/FV36/S9qXeGi1qlO3nrJT3CmsgCVQOsvsC6+trdBdbt06SKmFISDDadsfGpt9iYoyrnq6u\nxojQ778b2zIeM3IkVK4M27YZ64vcuv+ffzKvX5T6XMuXG22/X3rJSKD69TOSoNTSu7JlM89Hql8/\nvS14XrMkW7gRd4OwmDA8XD1wdXTl9PXTrD6xOj1BSkmSZnafSXOP5vzv8P8YuGJglsfaPHgz3ap1\n40r0Fbaf3467kzvuzu7GFVMnd0oWKwlA92rdWTtobdp+dyd33Iq5YaNs8HTzzHbujpebV9ocoCpu\nVajiViVv3xghRBZaGzcbG6Os+PvvsyZIX35pXAA6ccL4zCtWzEh+vLyMpjWlSxuP1bu3sSRC+fLG\n492qVCnjlmraNBg6NIDERD8gGPDE3n4a06b5EB8fz59//klgYCCPPvoonp6erFq1iueeew47Ozua\nNm3KsGHDaNOmDaVSHrRXr1706tXrnt6PZJ3MjdgbmRKgBuUaUL10dc7dOMeUHVOyJEhze87lmSbP\ncPTqUTp/3TnT47k5utGqUitaVmpJ2eJl6ezdOUsC1NKjJQDtPduT+HZijiPllUpU4pkmz2TZ/l7X\n97It753WVRbbE0WfJFA58PPL+gU2Jsb4sps6bSw1qYL0JOpeR62KwqiXWa8hOTk9ASlRwrgaee0a\nnDqVOTmJjYUHHoBy5Yw1ilasSE98UvdPmQK1asGaNTB5ctb9mzYZaxctXQrDh2eN5Y8/oG1b2LnT\nGDnKyMbGOOFXrgzXrxsxODkZi7w6ORlfCv76K/vXqBRcvZr+78aNjVtu0lqnNSe4HnudA5cOZBkh\nGt50OA3LN2TbuW34rvElLCaMm3E30Rh/HD8P/JmHaj3EsavHGLNpDLbKltJOpdOSnMTkRACaVmjK\n9K7TMyVA7s7uVCtVDYC+dfrSt07fHGP9twRoWtdpcnIXIhv5cZ7SGiwWo2zu8mXjszI1OUpNkObN\ng2eeMRKoV1+F4sXTE6T77jM+IwE6dDAuKpUrl33DmjufpxlAcv2h0CkR3IDwIBK3DmbSpEkMG3aB\nhJQav/Lly+Pp6UmfPn34/fffadq06W0tTKu1JjIhEhtlg4uDCzfjbvLT8Z8yJ0Cx1xjSeAh96/Tl\n78t/0/TTpiTr5EyPM6/nPF5s/SLxlnh+OfdLptGhMs5lqOleE4D6Zevzy9Pp+92d3dPKjAFqlK7B\nV498lWO8d7uukpT3CmsmTSRyYGMDukEAdPUDt2AI94St0+BQ1g+GKlWMk8GyZZlHrcD4Uuzvf3sn\np1tHve70/gVBdq/ByQlmzICePdOTkDp1wM3NeN927sya4AwZYpxI//gDPvkkawLz6afQtKlx1dLX\n19gWn97EjF9/hfbtjZP24MFZ40xtUxsQYOx3csp8CwiAli1h61aYNStzguPkBC++aIz6HD9uHHPr\n/latoGRJYw2ksLDMj21v/99d67y9IahE1v//vCJ87miSsyXZkinxSf3ZwasDNUrX4NDlQ0zeMZmw\nmJTyuJRjlj66lP71+rP5zGZ6LO2R6TFdHVxZ0m8Jfev05eA/B3n/9/eNxCdDEtTJuxOVSlQiLimO\nuKQ43BzdTJlAnNdd+MSdkSYSObvbJhLZJTOQc4Jzr+eZgAB4euZokjv7g5sFwm2x2ebL43UW4O6e\nOUF67z1j5OjECeMzv0SJ9ATJ2xsGDjQuNCUkGJ+V7u6319FTa010dDRKKYoXL05MTAy//fYbERER\nhIeHExERQUREBJ06daJjx46EhoYydOjQtH3H7Y+jH9KQcd5nAtiut2VM9zG0adOG1q1b4+FhjFRb\nki1cjLyYZQSopUdL2lZpy8WIiwxeNTjTvsTkRGZ0m8G4duM4c/0MNT42Fr61s7FLS3TG3TeOpxs/\nTVhMGHN3z80yQlS1ZFVKOZXK5h0QQuSV2z1PSQKVgzKdAwi7zxccMpxlEpzhz2eg9rosSZWdnTEC\nklw/65dep9M+9Opl1F7b2Rm3ceOgQQOjhfTixca2xYshwivr/UsE+fDss8bjWyzGzxdeME5Ie/ca\nraYz7rNYYMIE4/F/+81IXm7dP2MGNGtmfPl/++2s+z/91BhhWb8eXn456/7//Q/uv99ocT1sWPr2\n6OjbSzy3bIGuXY2StAEDsr7/W7dCly5G3ftrr2VNcN5/3xh52bfPWIsodXtqEtOvn3H18uJFo1X3\nrfsrVTLKQZKTjRN2QWsONHphAAsv+oJ9hv//Ep3wqfAerz7ansolKlPepTwhESEs/nNxlhK5CR0m\n8HDth9kZtJMOX3XI8viLH17M0KZDOXDpAINXDU5LfkoXM0aJBjUcRJMKTbgZd5O/L/+dvt+pdKYr\nm0LcicKUQCmlegJzAVvgc631+7fsdwS+AZoDYcCTWuvzKfveBIYDFuAlrfXG/3q+uzlPZZfMsNUX\npUB3yZzgzBi8gHbtjJK4K+VHQ9fM93E5t4D+/TNfrBowwBhhv3oVmjc3tl2rOBr6LMySfPDzKEqF\nLMDbOz1J6tfPGD2yWCAyEtzcNHFxcWmJTokSJahQoQLR0dH88MMPmRKg8PBwevXqRd++fbl06RI9\nevRIS4wiIiJITk5m6tSpTJgwgeDgYLxqeIFNyn+tlJ9vjX2LaROmEfpPKD2f6YmTixPOLs5sL7cd\nshuxiodBLQalJUFDmwzlhVYvcDHiIpVnV85yuF97P97t8i5hMWE8svwRI/FxSk+AOnp3pIVHCxIt\niQSHB1PGuQwlHEtINzohCrAC0YUvL05A//WYuaabHyTdUsPnEAOtFoFKSTpLBkEfX5yc4ZUuPkxf\nEwB9MiRdKftjf4bjx31ISjJOJElJxuR+MK7SBQQY2yO8sr9/xM/g7++DjY0xMmZra6xXUaeOURqx\nebOxLXWfjY1xNQ+ME96FC1n3pyxFgZ2dUTZx6/7UdTBKlTJOnLfuT6099/Q0rlqmbp+9JfvXAPDN\nOJ+0BCa17KxnT82howk4Olmwd7TgWCwZewcLxR2cACd69kqkcYdLWJItWLQl7Wf54uUBd+o1juGp\n8YfTtifrZCzJFihRE6iMa5kIEqvuIE5buJacsj/WgiW6GbWK1eJm/HVWHluZ5fE7eXeiSYUmXI66\nzGcHPsvy+H1q9+G+KvdxMeIi7//2ftp9k3UyFm3Bp6EPXat1JehmEOM2j8uyf2TzkfQ4dJdOAAAS\nEElEQVSp3Ycz188wbPWwLM8fHBmcOXkCsI8lIOxVAj6Dhb0XMrLFSK7FXGPS9kmUcCyRaQTI0db4\nD1i7TG0+fvDjLCVyFVyM9nrNKjbjyOgjOf4ZlCxWkg5eWRMwIYoypZQt8AnQHQgB9iqlVmutj2Y4\nbDhwQ2tdQyk1APgAeFIpVQ8YANQHPIAtSqlaWmtLbsc5Yv5okntnSGZKWqDvQrTCOEOmbEvuvZCx\nP16FWeOh9nTosBLsM9ynz0Kidl1kzZEXcHCMx8EhAXvHRE7frAa0xMYujgqdPsXWLpFrZRZlTp7A\n+PeDi2huf4qYuBjOxcdwNDoW2yMP0aHDh1y8Ekzt0bVJsCSQTHJaktO3Xl9+nPojV25eYciWIWnb\nlb3C1t6Wf47+Q9++fbHYWTjz8BmwAW2jsVf2JKtkTnmeAsDRzRH8sr4/ia2NUmHnks4cuu/Qf7+h\nDrA7ZDdlnMtQ0aUipYoZoz9li5flsz6fZRkhSt3v7uzOzqE7c3xYe1t7qpcuAB17hBC5Js8SqLw4\nAaXc578eM1dcT8qhg5e6ZcTOIYa4B59mSYk3oN8lsLFk2c8jQ7H1mIEdKu3K0wW314GBNOxwhuof\nGEMw+y78Df9v797DrKruM45/3zkzA3IRgREFUUDFC4oCRUTUxKhRtPESa4porQmKkUgN1rTqo0ZJ\nQqPmovi02PAgto8aNTVRKZpovcW2iYo+qKCigqLiDW/x2nD99Y+9znBm5sxwBmY4w5z38zznYe+1\n9t5n7bX3nB9r77XXzq0usv4k9t55Zv3zJRHBO3XfA05hr3FLqbv0G/XPreTzl/e8iHGcyuBRL7Nm\n8tdZTZC/2xgEL3a5hAP4G/rv+xKvHf+1BusCPMNljOAMeu/+Io+NHd9k+4esms4wvkW3XZZw55Aj\nifz2v/4uVK1vug8nnc63X59c34C4pu4aptZN5fU/L2a/X+3XpJrzfb+XvL+E/f61+fxlHy7jwDkH\nNpv/2p9eKzr86szxM9mj7x689elbTP7PyUXXH7HjCFZ+vpLLHr6sPj2nXPbi0l47M27ncXy86mNu\nWXQLuaocOeXIVWX5h+5yKACr1q1i0cpF9Xn59fPP5UhCiC7VXRoss/DthU3KlHfXhLsY2X8kAMP7\nDWf1paubHe2tX/d+TB0ztdltmVlRY4ClEfEKgKTbgBOAwlhzAnBFmr4D+GdlP/AnALdFxCrgVUlL\n0/b+2NaF/Hzs7KaNmWJRvRY48o5UzGbyvzyP95nXIPmBD0fzYxbwxZr3WLDbtJYL0y14gAeyhlnP\nLGnAJ1kXuG7du7F+v/V0UReqVU11VTU1VTUM3nMwAHV96xgxdgRda7rStbYrtdW11FTVcOLeJ2b5\nveuYMGYCNVU12SeX/XvY4MMA6N2jN1cecWV9ev7f/XbIYkeP2h7cNeGu+vSTbz2ZT9Z90mQX+tb0\nZel5S5tWT66Ws0ad1fL+m1lFac87UO0RgChhm22iuRG8ignWc9SuRzH36bnFF8itYch2Qxo0QPLD\nJOeqcvTr3i9rfDRuPNWvv5q+3fqiggZYt5puANRU1dQPGVqYv22XbBjVLtVd2Kturyw//WcdoO82\nfQHYpnobRg8YXZ+f307+DkX32u71dyAKt59/cL9nbU/G7z6+Pn/OwjnF90FR/z6cKlUxqv8oAHbs\nsSMzDp/RpIFx6KCsATJw24HMOW5OgwZKTrn69QdvN5j5E+c3acAM7ZM9XLt7n91ZMHlBk+3n92+P\nvnvw+rTXm2w/X7/79NuH1Zeurt9uY8O2H8aHF35YfJ/T9l8494Vm83ftvSuPfPORJumDrx3c7Ahy\nhYMq5Kpy5OovNZtZG9kJeKNgfgXQ+EpN/TIRsVbSx0DflP5Yo3V3KvYlks4GzgbYpbkxtFvSqxU3\ntQIOf/dwHtrhIYq+vidgUm5S1nipqaG2ppZxw8cB0G+7fly929XU1tRy/pN/T/RY32T1qk9zvDP9\nbWpztdTkaqjN1ZJT9ttUt20dq6avarJOXs+uPVl4bvMXjbpWd+XGE25sNr82V8uFh1zYbH51VXWD\n381ZJ8xi0p2TWB0bYm6tapl5/Mxmt2FmVqjdnoGSdDIwPiLOSvOnAwdGxNSCZRanZVak+WVkQeoK\n4LGIuDml3wD8Nq3W4jYLtl0YmP7itWIvfWjBLYtuaTKCl1B9I6jQoF6DWD5teYv/6V0+bflGv3Nz\n1+8IOsM+dATFzr9uNd2YfdxsD4JgW62t5Rmo9ohfEdHM7Z/MpjwDlbugmvXbltaIqvokx7qfraX6\nH6pZ16PpOrnPcqz9ydqNbuc713+H61c0fQZqysApzJoyq9Sil50HmDGzYkqNU0XemNA5RMTsiBgd\nEaO33377Vq9/2vDTmH3cbAb1GoQQg3oN4pzR59TfmcgrHBZ5xhEzWszfmM1dvyPoDPvQERQ7/9x4\nMtti3gQKx8cfmNKKLiOpmmxA7A9KXLdNfHv3s7MBHAqtTZ9Cq9OywNm7ng1rGuWvSeklmDVlFlMG\nTiH3WQ4ia3htbY0nyH5jl09bzvrL17N82nL/tppZq7RnF77WBKAVrQhAWyQwQfYD2/hH9eBdDm72\nqtXmvhOhM7xToTPsQ0dR7Pwzsy1iATBU0hCyGHMKcGqjZeYBZ5A923Qy8FBEhKR5wC8l/ZzsGd6h\nwBPtUchZU2bB9TD7ldms676O3Oe5+oZQ47R8A6e5dVrTAJo1ZRaz2LoaTGZmbak9u/BVAy8BR5AF\noAXAqRHxXMEy5wLDI+KcNIjESRHx15L2AX5J9tzTAOBBsiCkjW2zmE19v4aZmbWdraULH4CkY4Fr\nycaGmxsRMyT9AHgyIuZJ6grcBIwEPgROKXg+9xJgEtm9oGkR8duiX1LAccrMrPzKPox5eqh2KnAf\nGwLQc4UBCLgBuCkNEvEh2VU+0nK/IhscYi1wbn4I2GLbbK99MDOzyhQR9wL3Nkr7fsH0n4FvNLPu\nDMD9ls3MOql2fQ9UewSgYts0MzMzMzPbEjrtIBJmZmZmZmZtzQ0oMzMzMzOzErkBZWZmZmZmViI3\noMzMzMzMzErkBpSZmZmZmVmJ3IAyMzMzMzMrkRtQZmZmZmZmJXIDyszMzMzMrESKiHKXod1Jeg94\nrZWr1QHvt0NxtiauA9cBuA7AdZC3ufUwKCK2b6vCdCaOU5vMdeA6ANdBnuthC8WpimhAbQpJT0bE\n6HKXo5xcB64DcB2A6yDP9dCx+Hi4DsB1AK6DPNfDlqsDd+EzMzMzMzMrkRtQZmZmZmZmJXIDqnmz\ny12ADsB14DoA1wG4DvJcDx2Lj4frAFwH4DrIcz1soTrwM1BmZmZmZmYl8h0oMzMzMzOzErkB1Yik\n8ZJelLRU0kXlLk97kbSzpIclPS/pOUnfTel9JP2XpJfTv71TuiRdl+rlWUmjyrsHbUdSTtJCSfPT\n/BBJj6d9vV1SbUrvkuaXpvzB5Sx3W5K0naQ7JC2R9IKkgyrtXJB0fvpbWCzpVkldO/u5IGmupJWS\nFhektfq4SzojLf+ypDPKsS+VxHGqsn6bwHEKHKfAcaogrexxyg2oApJywL8AxwDDgImShpW3VO1m\nLXBBRAwDxgLnpn29CHgwIoYCD6Z5yOpkaPqcDVy/5Yvcbr4LvFAwfxVwTUTsDnwEnJnSzwQ+SunX\npOU6i5nA7yJiL2B/svqomHNB0k7AecDoiNgXyAGn0PnPhX8DxjdKa9Vxl9QHuBw4EBgDXJ4PZtb2\nHKccp5LO/ttUjOOU41Re+eNURPiTPsBBwH0F8xcDF5e7XFto3+8Gvgq8CPRPaf2BF9P0L4CJBcvX\nL7c1f4CB6Y/vcGA+ILIXsFU3PieA+4CD0nR1Wk7l3oc2qINewKuN96WSzgVgJ+ANoE86tvOBoyvh\nXAAGA4s39bgDE4FfFKQ3WM6fNj9ejlMV9NuU9sNxynHKcaqDxSnfgWoof3LmrUhpnVq6rTsSeBzY\nISLeTlnvADuk6c5aN9cC/wisT/N9gT9FxNo0X7if9XWQ8j9Oy2/thgDvATemLiJzJHWngs6FiHgT\n+CnwOvA22bF9iso7F6D1x73TnQ8dXEXWt+OU4xSOU45TG5Q9TrkBVeEk9QB+DUyLiE8K8yJrpnfa\nYRolfQ1YGRFPlbssZVYNjAKuj4iRwOdsuB0OVMS50Bs4gSxIDwC607TLQMXp7Mfdtg6OU45TOE45\nTjWjXMfdDaiG3gR2LpgfmNI6JUk1ZEHploj4TUp+V1L/lN8fWJnSO2PdHAwcL2k5cBtZ94iZwHaS\nqtMyhftZXwcpvxfwwZYscDtZAayIiMfT/B1kgaqSzoUjgVcj4r2IWAP8huz8qLRzAVp/3Dvj+dCR\nVVR9O045TiWOU45Thcoep9yAamgBMDSNaFJL9nDevDKXqV1IEnAD8EJE/Lwgax6QH53kDLI+5/n0\nv00jnIwFPi64fbpVioiLI2JgRAwmO9YPRcRpwMPAyWmxxnWQr5uT0/Jb/dWuiHgHeEPSninpCOB5\nKuhcIOsSMVZSt/S3ka+DijoXktYe9/uAoyT1TldIj0pp1j4cpyrot8lxKuM4BThOFSp/nCr3g2Ed\n7QMcC7wELAMuKXd52nE/DyG75fks8HT6HEvWP/ZB4GXgAaBPWl5kIz8tAxaRjQJT9v1ow/o4DJif\npncFngCWAv8BdEnpXdP80pS/a7nL3Yb7PwJ4Mp0PdwG9K+1cAKYDS4DFwE1Al85+LgC3kvWlX0N2\nhffMTTnuwKRUF0uBb5V7vzr7x3Gqsn6bCurDccpxynGqg8QppY2amZmZmZnZRrgLn5mZmZmZWYnc\ngDIzMzMzMyuRG1BmZmZmZmYlcgPKzMzMzMysRG5AmZmZmZmZlcgNKKtokkLSzQXz1ZLekzR/I+uN\nkHRsC/mjJV23GeVaLqkuTf9hU7ezOSTNkTSsHN9tZmYZx6kWy+A4ZWVRvfFFzDq1z4F9JW0TEf8H\nfJXS3k49AhgN3Ns4Q1J1RDxJ9r6KzRYR49piO5vwvWeV43vNzKwBx6nmv9dxysrCd6DMsuDyl2l6\nItlL2wCQNEbSHyUtlPQHSXtKqgV+AEyQ9LSkCZKukHSTpP8FbpJ0WP7qoKSZkr6fpo+W9KikBn97\nkvpKul/Sc5LmkL0MLp/3Wfr3MEm/l3S3pFckXSnpNElPSFokabe03PaSfi1pQfocnNKvkDRX0iNp\n/fNSendJ90h6RtJiSRNS+iOSRqfpiek7Fku6qrBskmakdR+TtEMbHhczM8s4TjlOWQfiBpQZ3Aac\nIqkrsB/weEHeEuDQiBgJfB/4p4hYnaZvj4gREXF7WnYYcGRETGy0/YvJgthXgOvI3oC9vtEylwP/\nExH7AHcCuzRT1v2Bc4C9gdOBPSJiDDAH+Lu0zEzgmog4APirlJe3F3A0MAa4XFINMB54KyL2j4h9\ngd8VfqGkAcBVwOFkVzQPkHRiyu4OPBYR+wOPApObKbeZmW06xynHKetA3IXPKl5EPCtpMNlVvcZd\nHXoB/y5pKBBATQubmpe6VzTe/heSJpP9cJ8fEcuKrPsl4KS0/D2SPmrmOxZExNsAkpYB96f0RcBX\n0vSRwDCp/uLgtpJ6pOl7ImIVsErSSmCHtO7P0hW7+RHx342+8wDgkYh4L33vLam8dwGrgXw//KfI\nupaYmVkbcpxynLKOxXegzDLzgJ9S0C0i+SHwcLridRzQtYVtfN5C3nDgA2DA5hQSWFUwvb5gfj0b\nLohUAWPTVccREbFTRHxWZP11QHVEvASMIgtQP8p34yjRmoiIwu21Yl0zMyud45TjlHUQbkCZZeYC\n0yNiUaP0Xmx4WPebBemfAj1L2bCkQcAFwEjgGEkHFlnsUeDUtPwxQO+SS97U/WzoJoGkERsp3wDg\ni4i4GfgJWZAq9ATwZUl1knJkV0B/vxnlMzOz1nOccpyyDsINKDMgIlZERLHhXK8GfixpIQ2vWj1M\n1v3g6fzDrMUo659wA/C9iHgLOBOYk/qxF5oOfEnSc2RdJF7fjN05Dxgt6VlJz5P1RW/JcOAJSU+T\n9XH/UWFm6opxEdk+PwM8FRF3b0b5zMyslRynHKes49CGu5pmZmZmZmbWEt+BMjMzMzMzK5EbUGZm\nZmZmZiVyA8rMzMzMzKxEbkCZmZmZmZmVyA0oMzMzMzOzErkBZWZmZmZmViI3oMzMzMzMzErkBpSZ\nmZmZmVmJ/h/tImNysAE0kQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f2a1d823198>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize=(14,6))\n",
"ax = fig.add_subplot(121)\n",
"for k, v in logp_results.items():\n",
" ax.plot(nvec, v, \"--o\", label=k)\n",
"ax.set_title(\"Average runtime for log posterior\")\n",
"ax.set_ylabel(\"Seconds\")\n",
"ax.set_xlabel(\"Matrix dimension\")\n",
"ax.legend();\n",
"\n",
"ax = fig.add_subplot(122)\n",
"for k, v in dlogp_results.items(): \n",
" ax.plot(nvec, v, \"--o\", label=k)\n",
"\n",
"ax.set_title(\"Average runtime for gradient of log posterior\")\n",
"ax.set_ylabel(\"Seconds\")\n",
"ax.set_xlabel(\"Matrix dimension\")\n",
"ax.legend();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Zoom in on the 1000 x 1000 case"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"ExecuteTime": {
"end_time": "2017-03-17T22:28:20.498765",
"start_time": "2017-03-17T22:28:20.481808"
},
"collapsed": false
},
"outputs": [],
"source": [
"logp_results = pd.DataFrame(logp_results)\n",
"dlogp_results = pd.DataFrame(dlogp_results)\n",
"\n",
"chol_lp = logp_results.loc[[7]][\"chol\"].values[0]\n",
"chol_lpd = dlogp_results.loc[[7]][\"chol\"].values[0]\n",
"\n",
"rel_increase_tau = logp_results.loc[[7]][\"tau\"].values[0] / chol_lp\n",
"rel_increase_tau_d = dlogp_results.loc[[7]][\"tau\"].values[0] / chol_lpd\n",
"\n",
"rel_increase_cov = logp_results.loc[[7]][\"cov\"].values[0] / chol_lp\n",
"rel_increase_cov_d = dlogp_results.loc[[7]][\"cov\"].values[0] / chol_lpd\n",
"\n",
"lp_vec = np.array([1.0, rel_increase_cov, rel_increase_tau])\n",
"dlp_vec = np.array([1.0, rel_increase_cov_d, rel_increase_tau_d])"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"ExecuteTime": {
"end_time": "2017-03-17T22:35:57.415350",
"start_time": "2017-03-17T22:35:56.987832"
},
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAF1CAYAAADBS4p1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XecbWV97/HPl2ahWAC9DqAHM2pCuII6MZZoUOyiGFsg\ngdhPNBbwWsGCPTEqETXtKIgGxEKJPYFLUOJVUUGkHYzbgqLoURABy9B+94+1BvYZpuzZs/fsKZ/3\n67VfM3u157dnnr1+61nrWc9KVSFJkiRJa8EWow5AkiRJkpaKDSBJkiRJa4YNIEmSJElrhg0gSZIk\nSWuGDSBJkiRJa4YNIEmSJElrhg0gDUWSY5O8pc9135DkuDnmX5hkn+nLJrlrkmuSbNlX0Fo1kjwk\nybdHHYe0WiT5QZJHtL8fnuQDo44JIMkLkvys3ffvOG3euiSVZKtRxdeLxcbZrjs+y7y/THLqTMsm\n+Zckr+svaq0mST6f5BmjjmMp2QBaAboTj6Cq/rCqvjDD9B9W1XZVdQNAki8kee4gykyyTZIT2/9F\nTTXAuuYnyduTXN6+3p4kXfP3TnJ2kt+0P/fudd21FHOPn2vWZD+lqv67qu61FPFIo5bkgCRnJfl1\nkk3t738zrO9kVb2tqha9bx3Agf/WwJHAo9p9/+WLjWm1qarjq+pRs8x7flW9GSDJPkkuHVS5SZ6e\n5Mtt/vjCDPP7zi9zrbvWYu7hM815QnlKVT22qj60FDEtFzaAtGD9JqtV4EvAQcBPZ5i3HngSsBdw\nb+AJwF9D0xABPgkcB9wB+BDwyXb6nOuu0ZgXZQ3XT61BSV4GHAW8A/hfwJ2B5wMPBraZZZ3VcpX8\nzsCtgQtHHchc1ug+6Qrg3cDfTZ+xmPzSw7prLeZFaRtua7MtUFW+lvkL+AHwiFnmPQ/o0HxxPwWM\ndc17FPBt4FfAPwFfBJ47y3beAJwIfAy4GjgH2GtaDK8CzgMmga2APwC+AFxJk4Ce2LX8scC/AKe1\n2/sicLeu+UcBPwKuAs4GHrLAWB7Rtexx7e/rgGpjeytwA/A74BrgfcA/Au+a9rk/Bbx0gf+PS4F9\npk37MrC+6/1zgK92/R9+DKRr/g+Bx8y37gxlvwo4C9iqff+C9m9/62Uc8z5t+a8ENgGX0SSKxwH/\nQ1N3D+9a/v7AV9p6dVn7v9umnXdm+z/+dft//fOu7b+KpqH3b1PTura5G3Ay8HPgcuB9XfOeDWwE\nfgn8J209BQL8QxvzVcD5wJ6j3h/48jX1Am7XfheeMs9yxwL/DHyuXf4RwOOBb7Z1+0fAG6atczBw\nSft9eQ2z7Hfb9w9o9wlXAt/q3tfQ5Ig3A/+PZn9+KrBTO++H7ff5mvb1wBlivxXNQelP2te722n3\nbD/L1Pr/NcO669r5U/vLMZp9/hU0efN5XcvehubA9Jft/uCV3fuQGbZdwEuA7wG/oGmAbtHOe2b7\nef+h/fu9heaE82vbv+km4MPA7abFub79jJcBL+8qa9Z9Yo+xfGnasuNd9eItwLbAb4Ebu/4XY8Bv\ngB271r0vzT506wXU0ecCX5g2re/8Mt+608q5I01ueEL7frv2//5XyzXmru/MW9ptXgN8GtgROJ7m\n+/p1YF3X8jMeTwGPAa4Frmu3862u7b+Vpo7+Fhhvpz23a5vPo/keXA1cBNy36zt0UlsPvg+8ZFo9\n/UYbx8+AI5diP9jva222+laJJA8H/hZ4OnAXmh3rR9t5O9E0Ig6j+eJ8G3jQPJvcH/gEzU7jI8C/\nt10MphxIkzRvT3Nw+GmaZHYn4MXA8Um6ux39JU3i2wk4l+bLO+XrwN5dZX0iya0XEMucquo1wH8D\nL6qma8SLaJLbgVNnO9q/0SPa7S/WH9Ik/infaqdNzTuv2j1E67xp82dbd7p30DRAX5vkHsDbgIOq\n6nfLOGZozkzfGtgFeD3wfporU/cDHgK8Lsnu7bI3AC+lqTcPBPYF/gagqh7aLrNX+3/9WNf27wjc\njeYg4ibt2e7P0Hw/1rUxTH1P9gcOB54M7ExTZ05oV30U8FCaA63b0XzP7GKj5eSBNI2BT/aw7F/Q\nHPRsT3Nl+NfAX9Hszx8PvCDJkwCS7EHTYDqY5oBnR2DXmTaaZBfgszQHbHcEXg6clGTnaWU/iyZX\nbNMuA833C+D27ff5KzMU8RqaBtbeNGfX7w+8tqr+h5v3Obevqof38Df4KM0B8RjwVOBtbR4FOIJm\n/3B34JE0+6f5/BkwQdMw2J/mZMqUP6ZpkNyZ5u/+zPb1sLaM7WgaMt0eBtyDZt/zqtzc9X3WfWKP\nscypqn4NPBb4Sft/2K6qfkJzUPz0rkUPBj5aVdf1uu1ZLCa/zLfuTarqCpq/w/uT3ImmQXpuVX14\nucbc5QCav/cuwO/RNIA/SPMd20hTX6fMeDxVVf9Bc4zwsfZ/ulfXOgfT5MrtaXLjTZI8jeYkx18B\nOwBPBC5vj50+3X62XWjq4aFJHt2uehRwVFXt0Mb88Tk+38jZAFrZ/hI4pqrOqapJmsbOA5Osozm7\nfmFVnVxV1wPvYeZuUN3OrqoT253bkTQHrA/omv+eqvpRVf22nb4d8HdVdW1V/RfNQeaBXct/tqrO\nbGN7TRvbbgBVdVxVXV5V11fVu2iSeHfjab5YFqyqvkZzNWzfdtIBNGd5fraY7ba2a7c95VfAdm0f\n4OnzpuZv38O6m6mqG2l2Si+hOZP591X1zeUcc+s64K3t//OjNIn8qKq6uqoupDnDtBdAVZ1dVV9t\n68YPgH8F/nSez3IjcERVTbb1s9v9aQ54XlFVv66q31XVl9p5zwf+tqo2tt+TtwF7J7lbG/P2wO/T\nnLnbWFWXzROHtJR2An7R1l0A2nsYrkzy2yQP7Vr2k1X1/6rqxvY78IWqOr99fx5Nw3/qe/ZU4DNd\n++/X0XzHZnIQ8Lmq+ly7rdNozgI/rmuZD1bV/7TfzY/THKz16i+BN1XVpqr6OfBGmoO3BWlzz4OB\nV7Wf/1zgAzT7U2gO9N9WVb+sqktpcuZ83l5VV1TVD2muTHXnv59U1Xvb/dhv289xZFV9r6quocnX\nB0zrHvfGdh91Ps3B7oHQ8z5xrlj69SHahmB7IulAmivsi7WY/DLfupupqlNpTqaeTlMn++2qvWQx\ntz5YVd+tql8Bnwe+W1X/t/2ufwK4z9SCPRxPzeTYqrqwXWd6g/a5NMcWX69Gp6ouAf4I2Lmq3tQe\n932P5mTmAe161wHjSXaqqmuq6qvzxDBSNoBWtjG6Wu7tTvVympb5GM0l0al5RXPmay7dy9/IzWfK\nbjF/avvtclMuacueaXvX0HQ7GANI8vIkG5P8KsmVNGfYd1pALP26aYfe/hzEzhyay8s7dL3fAbim\n/btPnzc1/+oe1r2FNgGeQXO28h9XQszA5dUOTkFzyR2aS+R0TdsOIMk9k3wmyU+TXEXTKOmuGzP5\nec1+FWw34JLug8QudwOOag8Yr6SpowF2aRv1U10nNyXZkGT630QapcuBnboPoqvqQVV1+3Zed47v\n3n+T5I+TnJHk50l+RXMyYOp7Nj1//JrZr37eDXja1Heo/R79CU2vhCndJ99+Q/td79Fmea79vZ9c\nMAZcUVVXd03rzlmbfeZpv8+me5npcU1ff6bPsRXNFaI5t9fjPnGuWPr1SWCP9ur8I4FftScSF2sx\n+WW+dWeyAdiT5qC/36v4Sx3z9Pw4Y76Eno6nZjJX/d4N+O4M0+8GjE37rh/OzXX4OTQ9Ji5O8vUk\n+80Tw0jZAFrZfkJTIQFIsi1NV4Uf0/QT3rVrXpilC0OX3bqW36Jd/idd87sPbn8C7Dbt5rm7tmXP\ntL3taC7P/iTJQ2j6Vz8duEObrH9Fc+DZayy9mOlg/Dhg/yR70dzD9O8L3OZsLqS9gtHai5tvzL0Q\nuPe0qyP3njZ/tnVvIcnjabpBnE7TJW7Zx7xA/wxcDNyjmkvph7N53ZjJbA0vaHb0d53lRuQfAX9d\nVbfvet2mqr4MUFXvqar7AXvQ7NhfsdAPIw3RV2i6xO7fw7LTvyMfobmKvFtV3Y7mns2p79llbL4P\nvi1NbpnJj4B/m/Yd2raqbnEjeQ8xzWSzPEeTZxaaC6a2c8ck3Wfdu3PWZjmTrs8/h+5lpsc1/bPN\n9DmuZ/MD29m218s+ca5YenGL/0V7UunjNCcLD2ZwJwwXk1/mW3cz7ZWrDTT3XP1N5hlBdDnEvBA9\nHE/N9h2bL2f+3izTvz/tu759VT0OoKq+U1UH0nR1fTtwYntcuizZAFo5tk5y667XVjRdFp6VZnjF\nW9GcFTqrvULwWeB/J3lSu+wLae6TmMv9kjy5Xf5QmsQ62yXMs2jO5L0yydZphlh+Au29Fa3HJfmT\nNCOdvJnmhsAf0Vz2vZ7mJrqtkryeW54dWUgss/kZTV/rm7RdG75OsyM/qbq6S6V5dtGxs20sya26\n7lPapv0/TO1kPgz8nyS7JBkDXkZzgyk0/ahvAF7SbuNF7fT/6mHd6THsRNNt47nAM4AnJHncTMsu\nl5j7sD3NTZTXJPl9moEeut3i/zqPr9Ec3Pxdkm3bv8GD23n/AhyW5A8BktwuTf9nkvxRe5Z8a5r7\nJX7H7N2ApCVXVVfSdAn7pyRPTbJ9ki3SDLE734HH9jRXRH6X5P409+lMORHYr2v//SZmP144jmY/\n9OgkW7bfr32SzHfCDZoccCNzf59PoLnnced2//f6tswFaXPPl4G/bWO8N80Z66ltfZxmX3CHNPc1\nvWiWTXV7Rbv8bsAhNAP3zPU5Xppk9/aE4NS9Gd1Xpl+X5Lbt/uhZXdubb5+40Fhm8jNgxyS3mzb9\nwzT3Lj2RrgZQbh7CfN1MG5uqCzRXubZo/+ZT9/F+gf7zy3zrTnc4zcH+s2lOGH44s4yCuIxiXoj5\njqd+BqzLwkZ6+wDw8iT3S2M8TbfwrwFXJ3lVktu0f689k/wRQJKDkuzc9tq5st3W8s2ZtQxGYvA1\n94tm9J2a9npLO+/5NJcqr6C5B2fXrvUeQzPK1tQocF8BDp6ljDew+chr36Qd9aMrhkdMW+cPaUZ3\n+xXNPRx/1jXvWG4eBe4amtG7dm/nbQkcQ7NDv4zm7MVN219ILMwyClz7/oHt5/8lzf1LU+sf1C73\nsGmf53S6RgXq8f+wrp0X4O/b/8MV7e/dI77ch2Z0lt/SjGp3n655c647LYaTgX/pev9YmjN9Oy7j\nmPdh8xHZtuqOo532JZrBHKC5Mfritt78N83BV/coRs9v682VNGe9Ntv+LGXeleZq3+U0oyR114eD\naUZ4mxoN65h2+r40N6pe065zPLDdqPcHvnxNf9HcX/I1mpNSP6c5QbWem0dPPJY2Z3St81SarlJX\n0+SO97H5yG7PoBmpqpdR4P6YJhdc0Zb/WeCu7bwvsPnoUs+c9n1+U7vOlcADZvhst6a5H+ey9vUe\n2lEvmbbPn2HdzebTXOH5TBvnd4Hndy27Lc0B/pU0N5m/lua+i9n+5sXNI69dDrwL2HKmz9hO24Km\n8faj9vMeR3PGvjvOqVHgfgq8smvd+faJPcfCDKPAdc07pl3/SjYfUfY7wBenfZ6HtHVixhHh2nKn\n555ju+b3nV/mWndaDPejyf9Tn3dLmpHPXrNcY57lO/OWaXE8Auh0faa5jqd2pMmvvwTOmWn7s5T5\nfJrBs64BLpiKl6Zr5Qk0dfSXNCemp8o6jmaEw2torm49aSn3gwt9pQ1aq1zb+r8U+MuqOmOG+W+g\n2Un0MvLNipbmxuDjaIY7rnbaNjQjm9y7Fj/CjSRpBUvyAuCAqvrTWeYXTZe0ztJGtvSS/Bfwkar6\nQNe019Lce/mvo4tM6t9afDjXmpFmaMKzaM44vILmDMWyHpVj2NrL2YcAH6iu1n9VXUtzT5AkaY1J\nchearnhfoRmK+mXccpjqNaft3jQ1tPZNquoto4lIGgzvAVrdHkhzmf8XNPfnPKluOUTwmpHkD2gu\n7d+FZphQSZKgeT7Rv9J0CfwvmhHQ/mmkEY1Ykg8B/xc4tDYfPU9a8ewCJ0mSJGnN8AqQJEmSpDXD\nBpAkSZKkNWNFDIKw00471bp160YdRt+uvfZattlmm1GHoRGyDghWfj04++yzf1FVO486juXIPKXV\nwHqglV4Hes1TK6IBtG7dOr7xjW+MOoy+dTodxsf7ffiwVgPrgGDl14Mkl4w6hqWS5BDgeTSjZ76/\nquYcOMU8pdXAeqCVXgd6zVN2gZMkqUuSPWkaP/cH9gL2S7JyjwgkSZsZWgMoyW5JzkhyUZIL27Np\n3fNflqSS7DSsGCRJ6sMfAGdV1W+q6nrgi8CTRxyTJGlAhnkF6HrgZVW1B/AA4IVJ9oCmcQQ8Cvjh\nEMuXJKkfFwAPSbJjktsCjwN2G3FMkqQBGdo9QFV1GXBZ+/vVSTYCuwAXAf8AvJLmQWOSJC0bVbUx\nyduBU4FfA+cCN0xfLsl6YD3A2NgYnU5nSeMcpMnJyRUdvwbDeqC1UgeWZBCEJOuA+wBnJdkf+HFV\nfSvJXOuYWLRqWAcE1oOVpKqOBo4GSPI24NIZltkAbACYmJiolXzj8Eq/8VmDYT3QWqkDQ28AJdkO\nOAk4lKZb3OE03d/mZGLRamIdEFgPVpIkd6qqTUnuSnP/zwNGHZMkaTCG2gBKsjVN4+f4qjo5yf8G\ndgemrv7sCpyT5P5V9dNhxiJJ0gKclGRH4DrghVV15agDkiQNxtAaQGlaOEcDG6vqSICqOh+4U9cy\nPwAmquoXw4pDkqSFqqqHjDoGSdJwDHMUuAcDBwMPT3Ju+3rcEMuTJEmSpDkNcxS4L9E8QXuuZdYN\nq3xJkiRJmm6YV4AkSZIkaVmxASRJkiRpzbABJEmSJGnNWJIHoUrSqD32oe8cdQgj9/kzXz7qECTN\nwn2U+ygtnaFdAUqyW5IzklyU5MIkh7TT35zkvHZUuFOTjA0rBkmSJEnqNswucNcDL6uqPWieoP3C\nJHsA76iqe1fV3sBngNcPMQZJkiRJusnQGkBVdVlVndP+fjWwEdilqq7qWmxboIYVgyRJkiR1W5J7\ngJKsA+4DnNW+fyvwV8CvgIfNss56YD3A2NgYnU5nKUIdisnJyRUdvxbPOqDlwDooSdISNICSbAec\nBBw6dfWnql4DvCbJYcCLgCOmr1dVG4ANABMTEzU+Pj7sUIem0+mwkuPX4lkHtBxYByVJGvIw2Em2\npmn8HF9VJ8+wyPHAU4YZgyRJkiRNGeYocAGOBjZW1ZFd0+/Rtdj+wMXDikGSJEmSug2zC9yDgYOB\n85Oc2047HHhOknsBNwKXAM8fYgySJEmSdJOhNYCq6ktAZpj1uWGVKUmSJElzGeo9QJIkSZK0nNgA\nkiRJkrRm2ACSJEmStGbYAJIkSZK0ZgxzGOzdkpyR5KIkFyY5pJ3+jiQXJzkvySlJbj+sGCRJkiSp\n2zCvAF0PvKyq9gAeALwwyR7AacCeVXVv4H+Aw4YYgyRJkiTdZGgNoKq6rKrOaX+/GtgI7FJVp1bV\n9e1iXwV2HVYMkiRJktRtSe4BSrIOuA9w1rRZzwY+vxQxSJIkSdLQHoQ6Jcl2wEnAoVV1Vdf019B0\nkzt+lvXWA+sBxsbG6HQ6ww51aCYnJ1d0/Fo864CWA+ugJElDbgAl2Zqm8XN8VZ3cNf2ZwH7AvlVV\nM61bVRuADQATExM1Pj4+zFCHqtPpsJLj1+JZB7QcWAclSRpiAyhJgKOBjVV1ZNf0xwCvBP60qn4z\nrPIlSZIkabphXgF6MHAwcH6Sc9tphwPvAW4FnNa0kfhqVT1/iHFIkiRJEjDEBlBVfQnIDLM+N6wy\nJUkahCQvBZ4LFHA+8Kyq+t1oo5IkDcKSjAInSdJKkWQX4CXARFXtCWwJHDDaqCRJg2IDSJKkW9oK\nuE2SrYDbAj8ZcTySpAEZ+jDYkiStJFX14yTvBH4I/BY4tapOnb6cj2uQBss6OHprZV9gA0iSpC5J\n7gDsD+wOXAl8IslBVXVc93I+rkEaLOvg6K2VfYFd4CRJ2twjgO9X1c+r6jrgZOBBI45JkjQgQ2sA\nJdktyRlJLkpyYZJD2ulPa9/fmGRiWOVLktSnHwIPSHLb9pl2+wIbRxyTJGlAhtkF7nrgZVV1TpLt\ngbOTnAZcADwZ+Nchli1JUl+q6qwkJwLn0OSyb9J2dZMkrXzDfA7QZcBl7e9XJ9kI7FJVpwG0D0GV\nJGnZqaojgCNGHYckafCWZBCEJOuA+wBnLWAdR9fRqmEd0HJgHZQkaQkaQEm2A04CDq2qq3pdz9F1\ntJpYB7QcWAclSRryKHBJtqZp/BxfVScPsyxJkiRJms8wR4ELcDSwsaqOHFY5kiRJktSrYXaBezBw\nMHB+knPbaYcDtwLeC+wMfDbJuVX16CHGIUmSJEnAcEeB+xIw21BvpwyrXEmSJEmazVDvAZIkSZKk\n5cQGkCRJkqQ1wwaQJEmSpDXDBpAkSZKkNWOYw2DvluSMJBcluTDJIe30OyY5Lcl32p93GFYMkiRJ\nktRtmFeArgdeVlV7AA8AXphkD+DVwOlVdQ/g9Pa9JEmSJA3d0BpAVXVZVZ3T/n41sBHYBdgf+FC7\n2IeAJw0rBkmSJEnqNswHod4kyTrgPsBZwJ2r6rJ21k+BO8+yznpgPcDY2BidTmf4gQ7J5OTkio5f\ni2cd0HJgHZQkaQkaQEm2A04CDq2qq5Kbn41aVZWkZlqvqjYAGwAmJiZqfHx82KEOTafTYSXHr8Wz\nDmg5sA5KkjTkUeCSbE3T+Dm+qk5uJ/8syV3a+XcBNg0zBkmSJEmaMsxR4AIcDWysqiO7Zn0KeEb7\n+zOATw4rBkmSJEnqNswucA8GDgbOT3JuO+1w4O+Ajyd5DnAJ8PQhxiBJkiRJNxlaA6iqvgRkltn7\nDqtcSZIkSZrNUO8BkiRJkqTlxAaQJEmSpDXDBpAkSZKkNcMGkCRJkqQ1Y5jDYB+TZFOSC7qm7ZXk\nK0nOT/LpJDsMq3xJkiRJmm6YV4COBR4zbdoHgFdX1f8GTgFeMcTyJUmSJGkzQ2sAVdWZwBXTJt8T\nOLP9/TTgKcMqX5IkSZKmW+p7gC4E9m9/fxqw2xKXL0mSJGkNG9qDUGfxbOA9SV4HfAq4drYFk6wH\n1gOMjY3R6XSWJsIhmJycXNHxa/GsA1oOrIO9SXIv4GNdk+4OvL6q3j2ikCRJA7SkDaCquhh4FECS\newKPn2PZDcAGgImJiRofH1+SGIeh0+mwkuPX4lkHtBxYB3tTVd8G9gZIsiXwY5r7ViVJq8CSdoFL\ncqf25xbAa4F/WcryJUlaoH2B71bVJaMORJI0GEO7ApTkBGAfYKcklwJHANsleWG7yMnAB4dVviRJ\nA3AAcMJMM+yqLQ2WdXD01sq+YGgNoKo6cJZZRw2rTEmSBiXJNsATgcNmmm9XbWmwrIOjt1b2BUs9\nCpwkSSvFY4Fzqupnow5EkjQ4NoAkSZrZgczS/U2StHLZAJIkaZok2wKPpLlfVZK0iiz1c4AkSVr2\nqurXwI6jjkOSNHhDuwKU5Jgkm5Jc0DVt7yRfTXJukm8kuf+wypckSZKk6YbZBe5Y4DHTpv098Maq\n2ht4fftekiRJkpbEvA2g9inYC1ZVZwJXTJ8M7ND+fjvgJ/1sW5KkXvSbwyRJq1cv9wB9J8lJwAer\n6qJFlnco8J9J3knT+HrQIrcnSdJcBpnDJEmrQC8NoL1onoT9gSRbAMcAH62qq/oo7wXAS6vqpCRP\nB44GHjHTgj5hW6uJdUDLwRqtg4PMYZKkVWDeBlBVXQ28H3h/kj8FPgL8Q5ITgTdX1UIy6jOAQ9rf\nPwF8YI5yfcK2Vg3rgJaDtVgHB5zDJEmrQE/3ACV5YpJTgHcD7wLuDnwa+NwCy/sJ8Kft7w8HvrPA\n9SVJ6tmAc5gkaRXo6R4g4AzgHVX15a7pJyZ56GwrJTkB2AfYKcmlwBHA84CjkmwF/I62i5skSUPS\nVw6TJK1evTSAHlpVl3ZPSLJ7VX2/ql4y20pVdeAss+63kAAlSVqEvnKYJGn16uU5QB9LMjV0NUn2\noOk6IEnScmcOkyRtppcG0NuATyfZLsn9aAYvOGi4YUmSNBDmMEnSZnoZBe6zSbYGTgW2B/6sqv5n\n6JFJkrRI5jBJ0nSzNoCSvBeorkm3A74LvCgJ9p2WJC1X5jBJ0mzmugL0jWnvzx5mIJIkDZA5TJI0\no1kbQFX1oenTktwB2K2qzptvw0mOAfYDNlXVnu20jwH3ahe5PXBlVe3dT+CSJM1msTlMkrR69fIg\n1C8k2SHJHYFzaJ6mfWQP2z4WeEz3hKr686rau230nASc3EfMkiT1ZBE5TJK0SvUyCtztquoq4MnA\nh6vqj4FHzLdSVZ0JXDHTvCQBng6csIBYJUlaqL5ymCRp9erlQahbJbkLTYPlNQMq9yHAz6rqO7Mt\nkGQ9sB5gbGyMTqczoKKX3uTk5IqOX4tnHdBysEbr4DBymCRpBeulAfQm4D+B/1dVX09yd2DWhkuP\nDmSeqz9VtQHYADAxMVHj4+OLLHJ0Op0OKzl+LZ51QMvBGq2Dw8hhkqQVrJfnAH2C5sFxU++/Bzyl\n3wKTbEXTFeF+/W5DkqReDDqHSZJWvl4GQbhnktOTXNC+v3eS1y6izEcAF1fVpYvYhiRJ8xpCDpMk\nrXC9DILwfuAw4DqAdvjQA+ZbKckJwFeAeyW5NMlz2lkH4OAHkqSl0VcOkyStXr3cA3TbqvpaM3Db\nTa6fb6WqOnCW6c/sLTRJkhatrxwmSVq9erkC9IskvwcUQJKnApcNNSpJkgbDHCZJ2kwvV4BeSDMa\n2+8n+THwfeCgoUYlSdJgmMMkSZvpZRS47wGPSLItsEVVXT38sCRJWjxzmCRpulkbQEn+zyzTAaiq\nI4cUkyRJi7LYHJbk9sAHgD1pus89u6q+MuAwJUkjMNcVoO0Xs+EkxwD7AZuqas+u6S+m6ZJwA/DZ\nqnrlYsqRJGkGi8phwFHAf1TVU5NsA9x2ADFJkpaBWRtAVfXGRW77WOB9wIenJiR5GLA/sFdVTSa5\n0yLLkCTpFhaTw5LcDngo8Mx2W9cC1w4mMknSqPXyINRdk5ySZFP7OinJrvOtV1VnAldMm/wC4O+q\narJdZlPk0x66AAAgAElEQVRfUUuS1IM+c9juwM+BDyb5ZpIPtPcQSZJWgV5Ggfsg8BHgae37g9pp\nj+yjvHsCD0nyVuB3wMur6ut9bEeSpF70k8O2Au4LvLiqzkpyFPBq4HXdCyVZD6wHGBsbo9PpDDj0\npTM5Obmi49fqYB0cvbWyL+ilAbRzVX2w6/2xSQ5dRHl3BB4A/BHw8SR3r6qavqCJRauJdUDLwRqt\ng/3ksEuBS6vqrPb9iTQNoM1U1QaaIbaZmJio8fHxQcQ7Ep1Oh5Ucv1YH6+DorZV9QS8NoMuTHASc\n0L4/ELi8z/IuBU5uGzxfS3IjsBNNV4PNmFi0mlgHtBys0Tq44BxWVT9N8qMk96qqbwP7AhcNOU5J\n0hKZ9x4g4NnA04Gf0jw9+6nAs/os79+BhwEkuSewDfCLPrclSdJ8+s1hLwaOT3IesDfwtqFFKEla\nUr08CPUS4IkL3XCSE4B9gJ2SXAocARwDHJPkApoRdZ4xU/c3SZIGod8cVlXnAhODj0iSNGrzNoCS\n7Aw8D1jXvXxVPXuu9arqwFlmHbSA+CRJ6lu/OUyStHr1cg/QJ4H/Bv4vzcNLJUlaKcxhkqTN9NIA\num1VvWrokUiSNHjmMEnSZnoZBOEzSR439EgkSRo8c5gkaTOzXgFKcjVQQIDDk0wC17Xvq6p2WJoQ\nJUlaGHOYJGk2szaAqmr7xWw4yTHAfsCmqtqznfYGmptRp577c3hVfW4x5UiSNN1ic5gkafWatQtc\nkkcneeoM05+S5JE9bPtY4DEzTP+Hqtq7fdn4kSQN3ABymCRplZrrHqDXA1+cYfoXgTfNt+GqOhO4\nos+4JElajEXlMEnS6jVXA+hWVfXz6ROr6hfAtoso80VJzktyTJI7LGI7kiTNZlg5TJK0ws01DPYO\nSbaqquu7JybZGrhNn+X9M/BmmhtT3wy8C5jxYXRJ1gPrAcbGxuh0On0WOXqTk5MrOn4tnnVAy8Ea\nq4PDyGGSpFVgrgbQycD7k7yoqn4NkGQ74Kh23oJV1c+mfk/yfuAzcyy7AdgAMDExUePj4/0UuSx0\nOh1WcvxaPOuAloM1VgcHnsMkSavDXF3gXgv8DLgkydlJzga+TzOC22v7KSzJXbre/hlwQT/bkSRp\nHgPPYZKk1WGuYbCvB16d5I3A1GnDTlX9tpcNJzkB2AfYKcmlwBHAPkn2pukC9wPgr/sPXZKkmS02\nh0mSVq+5usAB0CaL8xe64ao6cIbJRy90O5Ik9avfHCZJWr3m6gInSZIkSauKDSBJkiRJa8a8XeAA\nkuwC3K17+fZBp5IkLWvmMElSt3kbQEneDvw5cBFwQzu5AJOHJGlZM4dJkqbr5QrQk4B7VdXksIOR\nJGnAzGGSpM30cg/Q94CtF7rhJMck2ZTkFs/6SfKyJJVkp4VuV5KkBegrh0mSVq9ergD9Bjg3yenA\nTWfQquol86x3LPA+4MPdE5PsBjwK+OGCIpUkaeH6zWGSpFWqlwbQp9rXglTVmUnWzTDrH4BXAp9c\n6DYlSVqgvnKYJGn16uVBqB8aVGFJ9gd+XFXfSjLfsuuB9QBjY2N0Op1BhbHkJicnV3T8WjzrgJaD\ntVgHB5nDJEmrQy+jwN0D+FtgD+DWU9Or6u4LKSjJbYHDabq/zauqNgAbACYmJmp8fHwhxS0rnU6H\nlRy/Fs86oOVgLdbBQeUwSdLq0csgCB8E/hm4HngYzT09x/VR1u8BuwPfSvIDYFfgnCT/q49tSZLU\ni0HlMEnSKtFLA+g2VXU6kKq6pKreADx+oQVV1flVdaeqWldV64BLgftW1U8Xui1Jkno0kBwmSVo9\nemkATSbZAvhOkhcl+TNgu/lWSnIC8BXgXkkuTfKcRcYqSdJC9ZXDJEmrVy+jwB0C3BZ4CfBm4OHA\nM+ZbqaoOnGf+uh7KliRpMfrKYW1X7auBG4Drq2piiDFKkpZQL6PAfb399RrgWcMNR5KkwVlkDntY\nVf1iwCFJkkasl1Hg7gm8Arhb9/JV9fAhxiVJ0qKZwyRJ0/XSBe4TwL8A76fpCiBJ0krRbw4r4NQk\nBfxr+2iGzfi8OmmwrIOjt1b2Bb00gK6vqn8eeiSSJA1evznsT6rqx0nuBJyW5OKqOrN7AZ9XJw2W\ndXD01sq+YNZR4JLcMckdgU8n+Zskd5ma1k6fU5JjkmxKckHXtDcnOS/JuUlOTTI2oM8hSdJNFpvD\nqurH7c9NwCnA/YccsiRpicx1Behsmi4Aad+/omteAfM9RftY4H00D52b8o6qeh1AkpcArweev4B4\nJUnqRd85LMm2wBZVdXX7+6OANw0rUEnS0pq1AVRVuwMkuXVV/a57XpJbz7fhqjozybpp067qerst\nTRKSJGmgFpnD7gyckgSaPPmRqvqPoQQqSVpyvdwD9GXgvj1M60mStwJ/BfwKeFg/25AkqUcLzmFV\n9T1gr2EGJUkanVkbQEn+F7ALcJsk9+HmbgQ70DxUri9V9RrgNUkOA14EHDFL+Y6uo1XDOqDlYC3V\nwWHlMEnSyjfXFaBHA88EdgXexc3J4yrg8AGUfTzwOWZpADm6jlYT64CWgzVWB4edwyRJK9Rc9wB9\nCPhQkqdU1UmDKCzJParqO+3b/YGLB7FdSZK6DSOHSZJWh3nvAeo3cSQ5AdgH2CnJpTRXeh6X5F7A\njcAlOAKcJGmIbPxIkqbrZRCEvlTVgTNMPnpY5UmSJEnSfGZ9ECpAki2SPGipgpEkaVDMYZKkmczZ\nAKqqG4F/XKJYJEkaGHOYJGkmczaAWqcneUraJ8JJkrSCmMMkSZvppQH018AngGuTXJXk6iRXDTku\nSZIGwRwmSdpML6PAbd/PhpMcA+wHbKqqPdtp7wCeAFwLfBd4VlVd2c/2JUmaT785TJK0es17BSiN\ng5K8rn2/W5L797DtY4HHTJt2GrBnVd0b+B/gsAXGK0lSzxaRwyRJq1QvXeD+CXgg8Bft+2vo4abS\nqjoTuGLatFOr6vr27VdpntAtSdKw9JXDJEmrVy/PAfrjqrpvkm8CVNUvk2wzgLKfDXxsANuRJGk2\nw8phkqQVqpcG0HVJtgQKIMnOwI2LKTTJa4DrgePnWGY9sB5gbGyMTqezmCJHanJyckXHr8WzDmg5\nWKN1cOA5TJK0svXSAHoPcApwpyRvBZ4KvK7fApM8k2ZwhH2rqmZbrqo2ABsAJiYmanx8vN8iR67T\n6bCS49fiWQe0HKzROjjQHCZJWvl6GQXu+CRnA/sCAZ5UVRv7KSzJY4BXAn9aVb/pZxuSJPVqkDlM\nkrQ6zNsASvJvVXUwcPEM0+Za7wRgH2CnJJcCR9CM+nYr4LT2mXRfrarn9x++JEmz6zeHSZJWr166\nwP1h95u2L/X95lupqg6cYfLRPcYlSdIg9JXDJEmr16zDYCc5LMnVwL3bp2df1b7fBHxyySKUJGmB\nzGGSpNnM2gCqqr9tn6D9jqraoX1tX1U7VpUPMJUkLVvmMEnSbHp5EOpnkmwL0D5N+8gkdxtyXJIk\nDYI5TJK0mV4aQP8M/CbJXsDLgO8CHx5qVJIkDYY5TJK0mV4aQNe3z+vZH3hfVf0jsP1ww5IkaSDM\nYZKkzfTSALo6yWHAQcBnk2wBbD3fSkmOSbIpyQVd056W5MIkNyaZ6D9sSZJ60lcOkyStXr00gP4c\nmASeU1U/BXYF3tHDescCj5k27QLgycCZC4hRkqR+9ZvDJEmr1LzPAWoTxpFd739ID/2nq+rMJOum\nTdsI0D4EVZKkoeo3h0mSVq95G0BJHgC8F/gDYBtgS+CaqrrdMANLsh5YDzA2Nkan0xlmcUM1OTm5\nouPX4lkHtBysxTo4qhwmSVq+5m0AAe8DDgA+AUwAfwXcc5hBAVTVBmADwMTERI2Pjw+7yKHpdDqs\n5Pi1eNYBLQdrtA72ncOSbAl8A/hxVe03tAglSUuql3uAqKoOsGVV3VBVH+SW9/ZIkrQsLSKHHQJs\nHF5kkqRR6OUK0G+SbAOcm+TvgcvoseEkSdKI9ZXDkuwKPB54K/B/hhuiJGkp9dIAOpgmWbwIeCmw\nG/CU+VZKcgKwD7BTkkuBI4AraPpi70wzHOm5VfXo/kKXJGlefeUw4N3AK5njmUHeqyoNlnVw9NbK\nvmDWBlCS06tqX+BvqupVwO+AN/a64ao6cJZZpywsREmSFmYxOSzJfsCmqjo7yT6zLee9qtJgWQdH\nb63sC+a6AnSXJA8Cnpjko8BmY1dX1TlDjUySpP4tJoc9uF3vccCtgR2SHFdVBw0vXEnSUpmrAfR6\n4HU0D417F5snjwIePsS4JElajL5zWFUdBhwG0F4BermNH0laPWZtAFXVicCJSV5XVW9ewpgkSVoU\nc5gkaTbzDoJg4pAkrVSLzWFV9QXgCwMJRpK0LAxtOOskxyTZlOSCrml3THJaku+0P+8wrPIlSZIk\nabphPs/nWG75sLlXA6dX1T2A09v3kiRJkrQkenkY3O8luVX7+z5JXpLk9vOtV1Vn0jz3p9v+wIfa\n3z8EPGmB8UqS1LN+c5gkafXq5QrQScANScZpnnewG/CRPsu7c1Vd1v7+U+DOfW5HkqReDDKHSZJW\ngXkHQQBurKrrk/wZ8N6qem+Sby624KqqJDXbfJ+wrdXEOqDlYI3WwaHkMEnSytVLA+i6JAcCzwCe\n0E7bus/yfpbkLlV1WZK7AJtmW9AnbGs1sQ5oOVijdXCQOUyStAr00gXuWcADgbdW1feT7A78W5/l\nfYomCdH+/GSf25EkqReDzGGSpFWgl+cAXZTkVcBd2/ffB94+33pJTgD2AXZKcilwBPB3wMeTPAe4\nBHh6/6FLkjS3fnOYJGn1mrcBlOQJwDuBbYDdk+wNvKmqnjjXelV14Cyz9l1wlJIk9aHfHCZJWr16\n6QL3BuD+wJUAVXUucPchxiRJ0qC8AXOYJKlLLw2g66rqV9Om3TiMYCRJGjBzmCRpM72MAndhkr8A\ntkxyD+AlwJeHG5YkSQNhDpMkbaaXK0AvBv4QmAROAK4CDh1mUJIkDYg5TJK0mV5GgfsN8Jr2NRBJ\nDgGeBwR4f1W9e1DbliRpyjBymCRpZetlFLgJ4HBgXffyVXXvfgpMsidN4+f+wLXAfyT5TFWtyUeU\nS5KGZ9A5TJK08vVyD9DxwCuA8xnMjaN/AJzVnpUjyReBJwN/P4BtS5LUbdA5TJK0wvXSAPp5VX1q\ngGVeALw1yY7Ab4HHAd8Y4PYlSZoy6BwmSVrhemkAHZHkA8DpNDeRAlBVJ/dTYFVtTPJ24FTg18C5\nwA3Tl0uyHlgPMDY2RqezcnvITU5Oruj4tXjWAS0Ha7QODjSHSZJWvl4aQM8Cfh/Ympu7DxTQd/Ko\nqqOBowGSvA24dIZlNgAbACYmJmp8fLzf4kau0+mwkuPX4lkHtBys0To48BwmSVrZemkA/VFV3WuQ\nhSa5U1VtSnJXmvt/HjDI7UuS1Bp4DpMkrWy9NIC+nGSPqrpogOWe1N4DdB3wwqq6coDbliRpyjBy\nmCRpBeulAfQA4Nwk36fpPx2gFjOEaFU9pN91JUlagIHnMEnSytZLA+gxQ49CkqThMIdJkjYzawMo\nyQ5VdRVw9RLGI0nSopnDJEmzmesK0EeA/YCzaUbMSde8Au4+xLgkSVoMc5gkaUazNoCqar/25+5L\nF44kSYtnDpO00jz2oe8cdQgj9/kzX74k5Wwx3wJJTu9lmiRJy02/OSzJrZN8Lcm3klyY5I3DiVCS\ntNTmugfo1sBtgZ2S3IGbuw/sAOyymEKTvBR4Lk03hPOBZ1XV7xazTUmSpgwgh00CD6+qa5JsDXwp\nyeer6qvDiViStFTmugfor4FDgTGaPtRTyeMq4H39FphkF+AlwB5V9dskHwcOAI7td5uSJE2zqBxW\nVQVc077dun3V4MOUJC21ue4BOgo4KsmLq+q9Qyj3NkmuozlD95MBb1+StIYNIocl2ZKm8TQO/GNV\nnTXIGCVJozHvc4AG3fipqh8neSfwQ+C3wKlVder05ZKsB9YDjI2N0el0BhnGkpqcnFzR8WvxrANa\nDtZiHVxMDquqG4C9k9weOCXJnlV1wdR885Q0WNZBLVUd6OVBqAPV9sXeH9gduBL4RJKDquq47uWq\nagOwAWBiYqLGx8eXOtSB6XQ6rOT4tXjWAS0H1sH+VNWVSc6geajqBV3TzVPSAFkHtVR1YNZR4JI8\nuP15qwGX+Qjg+1X186q6DjgZeNCAy5AkrWGLzWFJdm6v/JDkNsAjgYsHF6EkaVTmGgb7Pe3Prwy4\nzB8CD0hy2yQB9gU2DrgMSdLattgcdhfgjCTnAV8HTquqzwwkMknSSM3VBe66JBuAXZK8Z/rMqnpJ\nPwVW1VlJTgTOAa4HvknbhUCSpAFZVA6rqvOA+wwrOEnS6MzVANqPprvao2lGwRmYqjoCOGKQ25Qk\nqcvQcpgkaWWbaxjsXwAfTbKxqr61hDFJkrQo5jBJ0mzmugdoyuVJTkmyqX2dlGTXoUcmSdLimcMk\nSZvppQH0QeBTNE/THgM+3U6TJGm5M4dJkjbTSwPoTlX1waq6vn0dC+w85LgkSRoEc5gkaTO9NIB+\nkeSgJFu2r4OAy/stMMm9kpzb9boqyaH9bk+SpDkMNIdJkla+XhpAzwaeDvwUuAx4KvCsfgusqm9X\n1d5VtTdwP+A3wCn9bk+SpDkMNIdJkla+uYbBBqCqLgGeOKTy9wW+25YhSdJADTmHSZJWoF6uAA3T\nAcAJI45BkiRJ0hox7xWgYUmyDc1ZucNmmb8eWA8wNjZGp9NZwugGa3JyckXHr8WzDmg5sA5KkjTC\nBhDwWOCcqvrZTDOragOwAWBiYqLGx8eXMraB6nQ6rOT4tXjWAS0H1kFJknroApfktV2/32qAZR+I\n3d8kSUM0xBwmSVqhZm0AJXlVkgfSjJgz5SuDKDTJtsAjgZMHsT1JkroNM4dJkla2ubrAXQw8Dbh7\nkv9u3++Y5F5V9e3FFFpVvwZ2XMw2JEmaw9BymCRpZZurC9yVwOFAB9gHOKqd/uokXx5yXJIkLYY5\nTJI0o7muAD0aeD3we8CRwHnAr6vKB8hJkpY7c5gkaUazXgGqqsOral/gB8C/AVsCOyf5UpJPL1F8\nkiQtmDlMkjSbXobB/s+q+gbwjSQvqKo/SbLTsAOTJGkAzGGSpM3MOwx2Vb2y6+0z22m/WEyhSW6f\n5MQkFyfZ2I7UI0nSQA0jh0mSVrYFPQi1qr41oHKPAv6jqp6aZBvgtgPariRJMxpgDpMkrWALagAN\nQpLbAQ/l5jNx1wLXLnUckiRJktaeebvADcHuwM+BDyb5ZpIPtA9GlSRJkqShWvIrQG2Z9wVeXFVn\nJTkKeDXwuu6FkqwH1gOMjY3R6XSWPNBBmZycXNHxa/GsA1oOrIOSJI2mAXQpcGlVndW+P5GmAbSZ\nqtoAbACYmJio8fHxpYtwwDqdDis5fi2edUDLgXVQkqQRdIGrqp8CP0pyr3bSvsBFSx2HJEmSpLVn\nFFeAAF4MHN+OAPc9wCdzS5IkSRq6kTSAqupcYGIUZUuSJElau0YxCpwkSctWkt2SnJHkoiQXJjlk\n1DFJkgZnVF3gJElarq4HXlZV5yTZHjg7yWlV5f2qkrQKeAVIkqQuVXVZVZ3T/n41sBHYZbRRSZIG\nxStAkiTNIsk64D7AWTPM83l10gBZB7VUdcAGkCRJM0iyHXAScGhVXTV9vs+rkwbLOqilqgMjaQAl\n+QFwNXADcH1VOSKcJGnZSLI1TePn+Ko6edTxSJIGZ5RXgB5WVb8YYfmSJN1CkgBHAxur6shRxyNJ\nGiwHQZAkaXMPBg4GHp7k3Pb1uFEHJUkajFFdASrg1CQF/Gvbj3oz3lyq1cQ6oOXAOtibqvoSkFHH\nIUkajlE1gP6kqn6c5E7AaUkurqozuxfw5lKtJtYBLQfWQUmSRtQFrqp+3P7cBJwC3H8UcUiSJEla\nW5a8AZRk2/bJ2iTZFngUcMFSxyFJkiRp7RlFF7g7A6c0g+ywFfCRqvqPEcQhSZIkaY1Z8gZQVX0P\n2Gupy5UkSZIkh8GWJEmStGbYAJIkSZK0ZtgAkiRJkrRm2ACSJEmStGaMrAGUZMsk30zymVHFIEmS\nJGltGeUVoEOAjSMsX5IkSdIaM5IGUJJdgccDHxhF+ZIkSZLWplE8CBXg3cArge1nWyDJemA9wNjY\nGJ1OZ4lCG7zJyckVHb8Wzzqg5cA6KEnSCBpASfYDNlXV2Un2mW25qtoAbACYmJio8fHxJYpw8Dqd\nDis5fi2edUDLgXVQkqTRdIF7MPDEJD8APgo8PMlxI4hDkiRJ0hqz5A2gqjqsqnatqnXAAcB/VdVB\nSx2HJEmSpLXH5wBJkiRJWjNGNQgCAFX1BeALo4xBkiRJ0trhFSBJkiRJa4YNIEmSJElrhg0gSZIk\nSWuGDSBJkiRJa8aSN4CS3DrJ15J8K8mFSd641DFIkiRJWptGMQrcJPDwqromydbAl5J8vqq+OoJY\nJEmSJK0hS94AqqoCrmnfbt2+aqnjkCRJkrT2jOQeoCRbJjkX2AScVlVnjSIOSZIkSWvLSB6EWlU3\nAHsnuT1wSpI9q+qC7mWSrAfWA4yNjdHpdEYQ6WBMTk6u6Pi1eNYBLQfWQUmSRtQAmlJVVyY5A3gM\ncMG0eRuADQATExM1Pj4+gggHo9PpsJLj1+JZB7QcWAd7l+QYYD9gU1XtOep4JEmDM4pR4HZur/yQ\n5DbAI4GLlzoOSZLmcCzNyTlJ0ioziitAdwE+lGRLmgbYx6vqMyOIQ5KkGVXVmUnWjToOSdLgjWIU\nuPOA+yx1uZIkDdIg71V98bP/fVBhrVjvPeZJow5BI+Z9ilqqOjDSe4AkSVqpVtO9qsuBfz9ZB7RU\ndWAkw2BLkiRJ0ijYAJIkSZK0ZtgFTmvCYx/6zlGHMHKfP/Plow5BWjGSnADsA+yU5FLgiKo6erRR\nSZIGwQaQJEnTVNWBo45BkjQco3gO0G5JzkhyUZILkxyy1DHo/7d35zGWlOUex78/URh0xjEqLozI\nEAYwxAXDYFRcQARRYxQ1QaIoajLRRAlq3BLvTYzJFfXGm1yvRiEgLgg6KIojIqBBQEZZhoEZltFW\nVBDXGBeUGbbn/lFvh0PT3fTQPV101/eTVLrOW3Wq3j7vc+o5b62SJEnSMPVxBOgu4H1VtSHJMuCq\nJBdU1fU91EWSJEnSgMz7EaCq+n1VbWjj/wRuAFbMdz0kSZIkDU+v1wC1p2w/G/jZJNPm7AFzfdu2\nbduCrr8WB2NQxoAkST12gJIsBb4JnFBV/5g4fTE9YG5sbMyHe6l3xqCMAUmSenoOUJJH0HV+Tq+q\nb/VRB0mSJEnD08dd4AKcAtxQVZ+e7/VLkiRJGq4+jgAdDBwLvCTJxja8ood6SJIkSRqYeb8GqKou\nBTLf65UkSZKkXq4BkiRJkqQ+2AGSJEmSNBh2gCRJkiQNhh0gSZIkSYPR13OATk3ypySb+1i/JEmS\npGHq6wjQacCRPa1bkiRJ0kD10gGqqouBv/axbkmSJEnDNe/PAZqpJGuANQC77747Y2NjD3pZ737b\nt+eqWgvWZ059Td9VUM9m8x3S4mAMSJL0EO4AVdVJwEkAq1evrlWrVvVco4XNz0/GgIwBSZK8C5wk\nSZKkAbEDJEmSJGkw+roN9hnAemC/JLckeXsf9ZAkSZI0LL1cA1RVx/SxXkmSJEnD5ilwkiRJkgbD\nDpAkSZKkwbADJEmSJGkw7ABJkiRJGoy+7gJ3ZJItScaSfKiPOkiSNBXzlCQtXvPeAUqyE/BZ4OXA\n/sAxSfaf73pIkjQZ85QkLW59HAF6DjBWVb+qqjuAM4FX91APSZImY56SpEWsjw7QCuDmkde3tDJJ\nkh4KzFOStIj18iDUmUiyBljTXt6WZEuf9ZmlxwN/6bMCyfv7XL2MAXV6jYM5iIE956Iei4V5am65\njXpIWOjbKM3eQo+BGeWpPjpAvwP2GHn9lFZ2H1V1EnDSfFVqR0pyZVWt7rse6o8xIDAOFhDzlAbJ\nONBQYqCPU+CuAPZJsleSnYE3AOf0UA9JkiZjnpKkRWzejwBV1V1J3gX8ANgJOLWqrpvvekiSNBnz\nlCQtbr1cA1RV5wLn9rHuniyKUyQ0K8aAwDhYMMxTGijjQIOIgVRV33WQJEmSpHnRxzVAkiRJktQL\nO0ATJHlSkjOT/DLJVUnOTbImybrtXM5FSbb7LhpJTkvy+u19nx6cKdp73x24vst21LI1d5LcnWRj\nks1J1iZ55Bwsc3WS/51m+u5JzprterT4maeGxTylyZinZscO0IgkAc4GLqqqvavqQODDwBP7rZl2\nhPls7yQPB6iq58/1srVD3F5VB1TV04E7gHeMTkxnu7afVXVlVR0/zfRbq8oflZqWeWpYzFOahnlq\nFuwA3dehwJ1V9fnxgqq6BrgEWJrkrCQ3Jjm9bZRIcliSq5NsSnJqkl0mLjTJEUnWJ9nQeulLW/mJ\nSa5Pcm2S/57kfR9re9oOS/LtkfLDk5y9A/7/oZmqvS9N8qm2V2VTkqMB2h64V47PO74XNMnKJJe0\n9t2Q5Plt+iGt/Bzg+lZ2W/u7NMkP2/ybkry6la9MckOSk5Ncl+T8JLu2aauSXJjkmva+vVv5+5Nc\n0eLoo/PyyQ3LJcCq1jZbknwZ2AzsMc13+6Akl7W2ujzJshYP69r0F7c9dxvb9mNZW/7mNn1Jki+2\n2Lg6yaGt/Lgk30pyXpJfJPlkT5+J+mOeGhbzlGbCPLW9qsqhDcDxwP9MUn4I8He6h+E9DFgPvABY\nAtwM7Nvm+zJwQhu/CFhN90Tdi4FHtfIPAv8JPA7Ywr03onhM+3sa8HrgU8DngbThRmC3Ns/XgFf1\n/Xkt9GGa9n4dcAHd7W+fCPwWeDJwFPClNs/Ore13BR4JLGnl+wBXjsTNv4C9RpZ9W/v7cODRbfzx\nwFhr55XAXcABbdo3gDe18Z8BR7XxJW29R9DdsSUtNtcBL+r7s13ow4R2+g7wztY29wDPHWm3yb7b\nO08zJr0AAAP7SURBVAO/Ag5q5Y9uyzkEWNfKvgsc3MaXtukrgc2t7H10t14GeFqLwSXAcW3Zy9vr\n3wB79P15OcxrbJqnBjRM097mqYEPmKdmNXgEaOYur6pbquoeYCNdEOwH3FRVP2/zfAl40YT3PRfY\nH/hJko3AW4A96RLVVuCUJK8F/j3ynv8AllfVO6oBvgK8KcljgOcB398R/6SA7kfDGVV1d1X9Efgx\ncBDdZ35o23v6cuDiqrodeARwcpJNwFq69h53eVXdNMk6AvxXkmuBC4EV3HtKw01VtbGNXwWsTLIM\nWFFVZwNU1daq+jddYjkCuBrYQLcR2mdOPoVh27V9X6+k26if0sp/U1U/beNTfbf3A35fVVcAVNU/\nququCcv/CfDpJMfT/aicOP0FwFfb+2+kSyDj5/z/sKr+XlVb6fbY7jkn/7EWA/PUcJinZJ6ahV6e\nA/QQdh3dXq3JbBsZv5uZf3YBLqiqY+43IXkOcFhb57uAl7RJVwAHJnlsVf21lX2Rrje+FVg7SSBq\n+03X3vdTVVuTXAS8DDgaOLNNeg/wR+BZdHu3to687V9TLO6NwG7AgVV1Z5Jf0+0pgfvH2q7TVCvA\nx6vqCzP9PzQjt1fVAaMF7Wyi0fac9Lud5BkPtPCqOjHJ94BX0CWml3HfuJnOg90WaXEwTw2LeUpT\nMU/NgkeA7utHwC5J1owXJHkm8MIp5t9Ct9djVXt9LN1emFE/BQ4enyfJo5Ls287BXF7dw/beQ7dR\nGncecCLwvbZHhaq6FbgV+AhdktHsTdXefwOOTrJTkt3o9pZe3mb5OvBWupg4r5Utp9uTcg9dDOw0\ng3UvB/7UksqhPMDekar6J3BLkte0eu6S7o4vPwDeNnJO74okT5jB+jV7k3636bYLT05yUCtflnZx\n8bgke1fVpqr6BN0PyadNWPYldD8+aMt8aluuZJ4aFvOUZsM8NQU7QCPaIfyjgJemu93kdcDHgT9M\nMf9Wuo3M2nZY+R6686FH5/kz3fmQZ7TDyOvpgmgZsK6VXQq8d8L71gInA+ekXVwInA7cXFU3zMG/\nO3jTtPfXgGuBa+iSzweqajwGzgdeDFxYVXe0ss8Bb0lyDV3bTrU3bdTpwOoWN2+mO3f+gRwLHN9i\n5jLgSVV1fqvv+rass+hiSzvYVN/tFhdHA59pMXEB9+41HXdCuouXrwXu5P6nCn0OeFhr068Dx1XV\nNjR45qlhMU9pNsxTUxu/sFELQJL/A66uqlMecGZJkuaZeUrSQmAHaIFIchXdHpvDF1IPW5I0DOYp\nSQuFHSBJkiRJg+E1QJIkSZIGww6QJEmSpMGwAyRJkiRpMOwASZIkSRoMO0CSJEmSBsMOkCRJkqTB\n+H+WUQ+plv46lwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f29fb973908>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize=(14,6))\n",
"ax = fig.add_subplot(121)\n",
"\n",
"ax.bar([0,1,2], lp_vec, tick_label=[\"Cholesky\", \"Covariance\", \"Precision\"], width=0.5, zorder=3, color=\"darkslateblue\")\n",
"ax.grid(ls=\"solid\", c=\"lightgrey\", zorder=0)\n",
"ax.set_yticks(np.arange(0,25));\n",
"ax.set_title(\"Log probability, 1000 x 1000 matrices\")\n",
"ax.set_ylabel(\"# of times faster than Cholesky\");\n",
"\n",
"ax = fig.add_subplot(122)\n",
"ax.bar([0,1,2], dlp_vec, tick_label=[\"Cholesky\", \"Covariance\", \"Precision\"], width=0.5, zorder=3, color=\"darkslateblue\")\n",
"ax.grid(ls=\"solid\", c=\"lightgrey\", zorder=0)\n",
"ax.set_yticks(np.arange(0,10));\n",
"ax.set_title(\"Gradient of log probability, 1000 x 1000 matrices\");\n",
"ax.set_ylabel(\"# of times faster than Cholesky\");\n",
"\n"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [conda root]",
"language": "python",
"name": "conda-root-py"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
@sammosummo
Copy link

Am I correct in thinking that s is not equivalent in the three models? Is there are way to convert between them?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment