From 19e9e3312f233c359e15903f0e742881863b51fa Mon Sep 17 00:00:00 2001 From: David Gustavsson Date: Wed, 7 Apr 2021 23:37:08 +0200 Subject: [PATCH 1/4] Add plot recipe for fit result --- .gitignore | 1 + Project.toml | 2 ++ src/LsqFit.jl | 2 ++ src/plot.jl | 75 +++++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 80 insertions(+) create mode 100644 src/plot.jl diff --git a/.gitignore b/.gitignore index 90debbd..2009d1e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ benchmarks/graphs/* *~ *.kate-swp +Manifest.toml diff --git a/Project.toml b/Project.toml index aaa2ef0..5b7a0fc 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NLSolversBase = "d41bc354-129a-5804-8e4c-c37616107c6c" OptimBase = "87e2bd06-a317-5318-96d9-3ecbac512eee" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [compat] @@ -16,6 +17,7 @@ Distributions = "0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24" ForwardDiff = "0.10" NLSolversBase = "7.5" OptimBase = "2.0" +RecipesBase = "1" StatsBase = "0.32, 0.33" julia = "1.1" diff --git a/src/LsqFit.jl b/src/LsqFit.jl index e1ecc91..0d139d8 100644 --- a/src/LsqFit.jl +++ b/src/LsqFit.jl @@ -15,6 +15,7 @@ module LsqFit using OptimBase using LinearAlgebra using ForwardDiff + using RecipesBase import NLSolversBase: value, jacobian import StatsBase import StatsBase: coef, dof, nobs, rss, stderror, weights, residuals @@ -24,5 +25,6 @@ module LsqFit include("geodesic.jl") include("levenberg_marquardt.jl") include("curve_fit.jl") + include("plot.jl") end diff --git a/src/plot.jl b/src/plot.jl new file mode 100644 index 0000000..478caab --- /dev/null +++ b/src/plot.jl @@ -0,0 +1,75 @@ +@recipe function f(model::Function, fit::LsqFitResult; significance=0.05, purpose=:neither) + @series begin + seriestype --> :line + label --> "Fit" + x->model(x, fit.param) + end + if purpose in (:confidence, :both) + @series begin + seriestype := :line + seriescolor := :black + linestyle := :dash + label := ["Confidence interval" nothing] + + [ + x -> model(x, fit.param) - margin_error(model, x, fit, significance), + x -> model(x, fit.param) + margin_error(model, x, fit, significance) + ] + end + end + if purpose in (:prediction, :both) + @series begin + seriestype := :line + seriescolor := :black + linestyle := :dot + label := ["Prediction interval" nothing] + + [ + x -> model(x, fit.param) - margin_error(model, x, fit, significance; purpose=:prediction) + x -> model(x, fit.param) + margin_error(model, x, fit, significance; purpose=:prediction) + ] + end + end +end + + +""" +```julia +_band(direction, model, x, p_conf) +``` + +Evaluates `model(x, p)` for each combination `p` in `pconf` (e.g. `[low, low, low]`, `[low, low, high]`...), +and returns the `direction`st return value. + +# Example +```julia +_band(minimum, model, x, [(0.0, 1.0), (2.0, 3.0)]) == minimum([model(x, [0.0, 2.0]), model(x, [0.0, 3.0]), model(x, [1.0, 2.0]), model(x, [1.0, 3.0])]) +``` +""" +function _band(direction, model, x, p_conf) + l = length(p_conf) + return direction(begin + i = digits(k; base=2, pad=l) .+ 1 + p = getindex.(p_conf, i) + y = model(x, p) + end + for k in 0:2^l-1 + ) +end + +""" +```julia +margin_error(model, x, fit, significance; purpose) +``` +Find the width at `x` of the confidence or prediction interval. +""" +function margin_error(model::Function, x, fit::LsqFitResult, alpha=0.05; purpose=:confidence) + g = p -> ForwardDiff.gradient(p -> model(x, p), fit.param) + c = g(fit.param)' * estimate_covar(fit) * g(fit.param) + if purpose === :prediction + c = c + 1 + end + dist = TDist(dof(fit)) + critical_values = quantile(dist, 1 - alpha/2) + return sqrt(c*rss(fit)/dof(fit))*critical_values +end From e255ff0f9f281c2228a3240133655a061b232e97 Mon Sep 17 00:00:00 2001 From: David Gustavsson Date: Thu, 8 Apr 2021 08:57:52 +0200 Subject: [PATCH 2/4] Remove unneeded function --- src/plot.jl | 24 ------------------------ 1 file changed, 24 deletions(-) diff --git a/src/plot.jl b/src/plot.jl index 478caab..f4d76b4 100644 --- a/src/plot.jl +++ b/src/plot.jl @@ -33,30 +33,6 @@ end -""" -```julia -_band(direction, model, x, p_conf) -``` - -Evaluates `model(x, p)` for each combination `p` in `pconf` (e.g. `[low, low, low]`, `[low, low, high]`...), -and returns the `direction`st return value. - -# Example -```julia -_band(minimum, model, x, [(0.0, 1.0), (2.0, 3.0)]) == minimum([model(x, [0.0, 2.0]), model(x, [0.0, 3.0]), model(x, [1.0, 2.0]), model(x, [1.0, 3.0])]) -``` -""" -function _band(direction, model, x, p_conf) - l = length(p_conf) - return direction(begin - i = digits(k; base=2, pad=l) .+ 1 - p = getindex.(p_conf, i) - y = model(x, p) - end - for k in 0:2^l-1 - ) -end - """ ```julia margin_error(model, x, fit, significance; purpose) From 430f3255519359d403f42c73df13518b05cafdf7 Mon Sep 17 00:00:00 2001 From: David Gustavsson Date: Sat, 3 Sep 2022 16:33:45 +0200 Subject: [PATCH 3/4] Add test --- .gitignore | 1 + Project.toml | 3 ++- test/plotting.jl | 35 +++++++++++++++++++++++++++++++++++ test/runtests.jl | 3 ++- 4 files changed, 40 insertions(+), 2 deletions(-) create mode 100644 test/plotting.jl diff --git a/.gitignore b/.gitignore index 2009d1e..a8bd3a4 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ benchmarks/graphs/* *~ *.kate-swp Manifest.toml +test/*.png diff --git a/Project.toml b/Project.toml index 5b7a0fc..585c560 100644 --- a/Project.toml +++ b/Project.toml @@ -23,6 +23,7 @@ julia = "1.1" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" [targets] -test = ["Test"] +test = ["Test", "Plots"] diff --git a/test/plotting.jl b/test/plotting.jl new file mode 100644 index 0000000..23d8f1b --- /dev/null +++ b/test/plotting.jl @@ -0,0 +1,35 @@ +let + @testset "Plotting" begin + model(x, p) = @. p[1] * x^2 + p[2] * x + p[3] + p[4] * exp(-(x - p[5])^2 / p[6]^2) + p_true = [0.5, 0.7, 0.0, 4.5, 0.3, 0.5] + xdata = -4.0:0.5:4.0 + ydata = model(xdata, p_true) + 0.7 * randn(size(xdata)) + + f = curve_fit(model, xdata, ydata, fill(1.0, 6)) + + p = plot( + map((:neither, :prediction, :confidence, :both)) do purpose + plot( + xdata, + ydata; + seriestype=:scatter, + label="Data", + legend_foreground_color=nothing, + legend_background_color=nothing, + legend=:topleft, + ) + plot!( + x -> model(x, p_true); + seriestype=:line, + label="Ground truth", + linestyle=:dot, + ) + plot!(model, f; purpose=purpose, title="$purpose") + end...; + layout=(2, 2), + ) + @test p isa Plots.Plot + savefig(p, "plots.png") + println("Plots saved to `test/plots.png`") + end +end diff --git a/test/runtests.jl b/test/runtests.jl index ce4b553..11638a6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,9 +3,10 @@ # using LsqFit, Test, LinearAlgebra, Random using OptimBase +using Plots import NLSolversBase: OnceDifferentiable -my_tests = ["curve_fit.jl", "levenberg_marquardt.jl", "curve_fit_inplace.jl", "geodesic.jl"] +my_tests = ["curve_fit.jl", "levenberg_marquardt.jl", "curve_fit_inplace.jl", "geodesic.jl", "plotting.jl"] println("Running tests:") From 8ff8e30bbc34f9e302a71e8445fc140259fe35c8 Mon Sep 17 00:00:00 2001 From: David Gustavsson Date: Sun, 4 Sep 2022 13:24:15 +0200 Subject: [PATCH 4/4] Document plot recipe --- .gitignore | 1 + docs/src/img/plots.png | Bin 0 -> 40730 bytes docs/src/tutorial.md | 27 +++++++++++++++++++++++++++ 3 files changed, 28 insertions(+) create mode 100644 docs/src/img/plots.png diff --git a/.gitignore b/.gitignore index a8bd3a4..bc8c91c 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ benchmarks/graphs/* *.kate-swp Manifest.toml test/*.png +docs/build/ diff --git a/docs/src/img/plots.png b/docs/src/img/plots.png new file mode 100644 index 0000000000000000000000000000000000000000..59ff55a0a207bd033dd9d127ee988093a45d9b9d GIT binary patch literal 40730 zcmd43WmJ`G+cmuCMM^57C?F*u9nyldG)M{39U@3bcOxY&-60|fh=O#tbO=%s(j5Zd z!F@l^kN594#yiHl_8yzPw~Mu|bzbK=k2&Wtj}@w{D2|yPvXF=zvM3!&YMrpoZl4guk!qW5^COcad_Q|W_Hge0V$FONRr6Lq z>dNs_m!pfb^Op99@8(;*X@clI6+VRG$mh>#=TN8u2*o{`B|-?r3*}n=fBqOYrcq}k zeKlKWNS+_yq^hbK8y9z|BG^+wr)ew0t*}qo{Gy!%yMuX~#ugia@`;Lx|qSaIS-WOG`7o|MWpA3J|`8d{*(DQ^f@xQMUs#?Sy zt1la5sLn(-r?P`No#}g9iFN|B>L4C@i$#mwxE=4+tL>NdZjW@oDl1i8f4b$|LOvVo#k{>lVk#d*(n_F67#27xFiqm5G^A3ssYJH?gI_e$|509|ppWgM+ zlB$n*Sz#fDuMWoaVdZhSAjY^jD z&-W>c)UAw-TjC!*dgOH2(^*j(^RFM(=LuD#%q;bZcixo@5r!27&cjTt&gUq7;^`C? zH*pvdKDuP(h(L-#_RL(dE(Bt!#MjIh1D=3TW+A=28dB-M(vM{vL5_S0nHJsW$*eb1 z%|t1lz#zpXLWc2h7a?-5rnip11m-($xBu(R=%1-3|p_FQ6k1HYOdn zXY2p``P0FUUy zccwXs5*F6hFO-!V<8g6tNtS?jovBMa&nIl$2C;W7-dMkEb*?NE|;2M~nuM970HJ@DLLPMT3|M zg{2K0p%vy!+bEATLXp8RV5Uc8CSV|uRESbq1VS9uB)jH&c0BsmSZrA|1=&1RL^&am z-__i6QyM~ziRo}7Bcrb@NCa#_diwl;_?+x)f_p-p@td2Q)77@H1(U7{;y%WcNHq^x za8FK7>PoM}=}8lGp+_1WZ~a2U#EX{8b3fj)ohUbK3Bqq_X^D%AlYV#i_Wk<}CS9?j zSN}*5Zbut$=bJGQTCk*HWXvXRUCtoL!kHNX|VgTzkV&Zuw*@LGd#mW=^5NeUh6UW`&Cpq;vCUs zCS38+Wv(wemtHbaOSk=jyMWKn#H)`#MqHjim>}+3_5L$y+)bnC&?oe$)*SP!XvxT2 z)aXzd_#wXaNSq0gJo}SRfy9a@DQRceX*e!Bn`VOK=H;sA&Z%>uVrL957vec zRd6}$?vm9#G@2u=c!fi)&_qurmg1)n`?UB&pQ>*f=|*lq8v0!YY(eY^S?pKV<$Nf2 zuq_{XrCyb`zT%@x?4ctLjrmaU(F4vBE-r4ZgP~zoV4x(nx#`)_=Jz*is%agKbmCn275R~98?8QBpZaE8Z*(-(rO|QXGUaS@Oa=xs`!x; zb=^Y~S=8cG5}w|yZ~XHX;$Ixl8PtBZk$#XzN|R7AKhG#ODp4zVC!Y1Y`}FiwR~qFx zlAoU+;1f~)cK?rT96d@Sz*~Fdot4c^%GC<(kSezOcqC|sN;E%xJyBCtt*ET5ud5@)X#+43e6yOVQHgb9>mIumDdOM7@m#(8admk) z5+C2m$?5R$(DUEPYh7Jk*aJ}b`CiRlKLxGvc+R`Qs0|09Rc_G4H|Jq&emqNeraG)i zLn~YHG?>6`>hVwNpg>KGjp4gTrXqE{#`=Fgv+Yx>b&Y+dTGt_VdX$;PGt?JZ`p|#< zPe7o`hZo`*X0Sxd%kk;yTfcwX!s4#)M;8XJAl z2GRtv5cKr)P;+UCAHA%bXl+GDd^=i=XEFa&@EJSM`60bC)~n$)IT_6)x}Ksp7^!)v z*4Ptr$QKVMUSB2p#XX!i#n&(@kBx?Z5c|g`{XVmus^kRX0mYJ5_J`NS)X2!l?5qei zH8l@Um3F0>q@*Mk!CfSU;7!Dz-sDP$WjRC=EXWtj;SV1^@Y+l$7bvDFrU|6+JN}uf zvVNAp0<||;=q;y9*W2#*4`+TS_NDS+VPOgK^Z%-QTj{(x1-t=jGuz9$oxQzm#&x(G zx)p=(SYo9g`c2caIEjo^uB~u~Q3uY;1|$DW_f6A#{4w)E6=6 za@_B#80k=>WG~7cFRP)UAtMt~U+?MWc9z6pEaGfMI-N<)fv2Xf{^sx4BWC6jv%VC) zulDovEx}x7y~<4SaL93x$ln@jLxh@*?Khv;h}=mXpeGemB}!FpHG9}J6i{eKP-Yeo z{_t@s@jLq82I60n{QR*E*)Z%16IF|rwx&a)qHf^?mX(!#0w%b-AeqRj7Zw(FlkH{i z^7YyT1_$G!`+YX*}5E_PtuxQBy)VUs%r25BoO+OgP&d%=b z?WN$c8sx5bTzS{O*pAZJ8S~)k^86LzrpVmqr{2nY^D~ z7iX`1H-CLkWVTQ$Be=~O!0hVg=1?#Uj5ZM?OY4Tbd;YKP?l9&Zr_|^&1DetgFK<y z1x8X4*+&t)`w`uN47b(_@A3*{d^D4D)h0P)$x0UH_{yvO0@LnKfR9gaul&}x@rJ9* zzdtj5!=D~PH6Xir^Bp+{&BKRlgIQ1a|MV@kdpJCY!VXSAi!laS1DCj9Oq$@TAhU&=7p{_OU~FRine^fP^fNl03SjY&+6!M@0( z(ddIOH1YwNszva&%i41!VUZ}OXMaEO*-n|*+J5YF_VxAEZE9&To~?VUyQ%MaNT%6H ztzsG!T5hDqW&bGtf-iHf<0 zOsyIf;RivNvJk%)Ya2YW`=;h=-;-HnpQs5|6bW>bC?<{Z6Vn@^pkA=bdhR1R@S(T2k{qHAyd^KYqU!eIUC*Qu_N;DFf{-NbG zUE__ru4Y%vf83a++mlzH`|D1tE-6n+3DTHE5EB#s{P{D5*H%tt1Q;URM&P;MmC}#* zpZ<}jab%0LwRwjA{DCG9L+T_Ez4%w>O`X3KA5jc?lQ{NQ20j5>fJ016j7|Fejk%%W z{PJ?!N{u+e2R^Ve7vWjZBuP&Y1BY!|SBz8VQ{UX%m&4S#-T;nJS<}6Rek@#}?vhK2?$mf3PcG`@H| zzRBDpJAM)|UB*9q(CHcI7xC-z;qk}dq705lZQ+Q`n9)wCV~mbn4`JB$IlgPXA(t)0{I(Z~<1B&0I+$qOQihhf0&eyUwzvN-Zr7=P7 z?(V`-7et(#o<7p~vA3Prn^;ho$*#lXQz9I|=Y~_jOGT0)`VTM+CDZQLw-+^b%OZuI zg^!B|*@^bHWeEgVQGYQ7IB@&kLX)_ogT20}=Yw3!kjwk#vV9NnwYJuCx9!$L2QB-I z-N$lwUE6wXWPaRZ(H&D)CcJa{%ViG5Cn_rHhyd%s+b`5V$>ak(!upMaU*`HZGkUu( zj8YBJi2A0f$SX%E$2V569DfKff4k%+skZ$fMVCx>lNE21w)6#gmDTb9;`icWOP@@? z@KV%DNQ{>O&62b#k>D22uQS?ZTZL9U(morRltpTKam=O)8>lOmDG73U7k-)tBGN9G ziAPm*FY$@{92CB}oFNdduC91@9}^Xjqty5q_&#eyQSTdVIzLl}g4C?^>bjLK_j__(;x5HT|9-<@rLEJ|Yf z98&}xmAVDz#tkrq`uRQ`Evf~F;GK1}FIn7Bl@NcKeg+gs zL_`D#5NPz@g7PL76A!qf6W`i&43Bh6_clFgnsr9)^+Vq71s(-A2LrK~KcGjvJtOko zYFS3&2I;%PqOIu?EQFSh4(p5Rh;4SFSL{UVX+gg6N`W60RrO{%LtZu~5G0`zaT0E2 zNnT!F67-Lb2^rSCv&cnt6?w9Bwh^t4P%U}`TGpc@UQ zYLD7s6HIdkgmu!~J0R-cp$~NZ?TZ&^s%@)%zA?tM47(Pv3CqVd=5MCCQ$HQYxEpf( z3#nEpwXm>Wk-GS2^U|DRd417d*x`pqN8nQo=POiJRgJqp8ZXxHtx?}vdpnz=mKJd5mb~}n zKYoOKDz5@ln|u7=3%9#0YFb*I?0Cl$F~#(sLns|$5~2x*UTC#^Xk7afGp|Bg<+^JW zFE^U?G`N@-r5jj2ukF+u5ozAdpmHUJuRnQ%mB~;?&ejSt%5g(x3pykURg2Pv-Y(3~ ztBg|2pPD$2n7Vlg8;KdC^v49bj+AJ!`4x(PHvlhAi>11zX0Fb98?NwX8y^k;LVqpEGEy=?j#F3*^`x!*INpPqhv zS+B!@qD(5O@?I)usB28_#Y41`?11n3*q2sWmX-sgWIyy?YH4XLtl)LH(x@U`8tBh8D7OifAIJE6RB_pyPU-8~{A4ySc>EUe(DkIswjk-$(Lmx0OH^Tj6K zAIHQ$J-jdM_C797P1vAuwPdjQUbBvw-qo?qF!fVEKR=+Pc>k<-I4C09U1_Lgf0CEo zALjlzWDy`jaxgj(h4Vf8Zy{6q04&cyXVqQ7?e%r{_>33{8bTy3ap*3s zI0*^K542Mohg#h{)%u2+kr8#8m^W)b+nSoh7)eu-lNHj1PJxEqpMIxW)Z5iH-wg`x z+nK4!Wrkw425N#S6oTQpxBJzbQ*7=M==uTpk0#2>SQK6ia`#~vB+1CNOwe!yIbvc? zKNrd02pt?4_ypT-n{7q0X&{+v<+>bIS66@e5)l@*VCH!8+Zi;%zkmO_?*_%VMzBO{xfm{9NLO-f2)jFH*c+&u3;%WIQ`-MrV)PAl6OP5m@t zsTZr5@L;b+fv|OTC?`YM{RouGqZidSQ(t87SCNN`rn9TPVmC|+NWzHOr9v73`2xcc zVZZac3CPqwWz0(Znl*2d$veUrcCM*ZoSYy(0TSWp*#KGu6HuNINOusEkm%{@O=Dp8 zeN)Kr>V@)fP)=)+SpD+lOON0nrhCG(A|N5)k=3xdA8ur@L3Z{W#l7nLAjJ*gnB5q; zJb$bke*|$nPv?X0Ym~j@+fA?iVD##MP?3V);{0fnGr2V2886ZkFTIZ$mj?~4e|VDV zNcYw?{Jz)ui0Ly|efjG=CApqd(dDMjGlO5>oZ=6ANjMSm0T~KnV%v`IqXTj%1L#l4 znFH_`d3l|G=F0A3*!+!C54LL6n%}G{@ka9@-ino2joS(Sx6f#~>iMpMkJ`DK=r`UJ z^0wCIRF%J`t>qoLfj#O5yC*;2>{E=E>{F5!@9fn4<0`aRb$nM5NBn)t=%z$xCTJ@% z3W_tvRD1!nXa730<4R&xkx1Y`WD1JN;k0&{LE$6@%m*(s* zGAyZuZ?62b^YlZ!7(E@somQ?f|;63w#H0u0@Q9C)rVG}P`RvV^}p0?MV=rO4z z1{PLNd~w#THO#Sqj!NP43^U5BkM>GYcntKL(@PAnk;@mm~WosZ#ay_X6NIwjqv6H#4mzod}p zwCASigOY|rx_&u-;DX<^gTD{5ts8D~apTJYZ!MqKnrlCJLhp_4o*=0Gi9;nGNJS($ z>MN^Z&hgN$(LIwPPNn7~oEg2W+`S#A)p|LD$A+IB#>B-{dD2r?S65ibn38`QR`7mr ze;?&{-t(Lt6X&EeBPq9U%8acNA$XWBRW;cTlaj}4-l*tC|Kc)nEa03#o9Qr0^s-Qo zn1bS)^*HNit!t1dEoHf9e4aa$BR%IC$F+&+`B>S2DF7fI?(g&Pey>A-hbT(He9xjG z3xd0k?V+7xJzlm}jB}yF#eP6j->m=JFf8>I-<_6(_N6;kV$e#-SauHU4SN z4+HlEr<9xEuE?2Fu1!ovV@D7Y2dDYANwn`8Jy97j-NW_KYu?S8EaK?~P<>5lLBij= zR~DypV7}}5Er_3E3{B9#FxGWDT;Wt)kF`*3>#`Kob_lZ4gnT3AohN4S* z#i(n%T%J0Ulj^zPd$~Lr6H?om!v8eKetWBxtXRs#J)g$A^prZ4tp+(>fp6p-rm}9!8RKMvbm!}n$e?%`qsb`bTJEV ziUd>9qo<1-8#CZFlx-|5WS3-ELi5It@vZgd_d`1eKRp$;G)9_S0s^q!8J;&N(GuU5 zJ>#gBiJ&HmrR}OHO}}Kx>-q$(HCQk3<#Z=!W`w|!4v;$89}r~&-pX8KJAE`#q$UvR zd*)5|sruk^X3)dQKCzpjy(rQ6 zXq=!Y_?_6PC}6Q;8`eg9g~I|H9o2YZvhjq5quAi2xL!S!9a;~tMXuLmA^5&mbeG9C z3m*r^Q2LxQBdO8nuK`kzUAAWOo+Ty-{QkL6VbFz@6jeDXYFyMQBJRj_Az~A4sttZ? zy@rY$v}I*ui9=h-DX0^X3i1(WAAI8DM}1JxKMhI{uOWf%Y#7v?(KLZt~~>7CrNN zzMXFUtFcmTfS+RSv@{A+Z4GyHodQKlW@hHFYyfE3#NVAai%WFgSp0RAl4w3Td#pJjcX$mR64pt^>gp;e;Xh9LQ(f`o6P>5k>$unF zJI#)G9-2M<>B6b$=|7zQjsaiM*3r3y288uk>F)9I4Q%X(UU8pdI@Oz%j@P*K0w!)t z#QT5xv@Q4y>4yZ@@Z)$hRmHuxQKy+E)JXMq<*JT_S$n7pk zYxVdwm!G!KZK=0wu(8Yid@FXc@bdBeV{uH$;Wy6CN9zR<;kP*ot1|~eNv6(A2^Uo2 zpJC$Ra@h?rU&uPffxn4N?8gim8)037B=Jl$84mMx11 z+?Jf_+xCQKN~7!K7keQD#Zuf=EkYsDS(E&aSv~Bk1|O_6CoBIj7F+dJ&yq+|Z?*Gi zXW1{yC6V}54o(`y#hM6LmWcMRyWjnR`UD4e zumMdC7%ULL+1%I&BB1Pckw9KV75@|lDX&y#W_G$GT1hAV9BMuRqVz5GhE}A|=dA*0 zm%xo))q74ybZMWMy!@biY&D+tR(_+2Dd>;BfDtfoGvlF@@o~zsZzg`UHs63heqwqgZ zTr<4Xb9?zf)vrL3h_m#Y(_hTUvm*~qu1SAPjZiTMmhta2i-_l~=R`mdzQV@5I6n?&V_{@eE*S-l6%--C6t7`D^G6)X)4`!w z^oHGpK&sG&$vpU>9*VSnORk$ohcDr$0C7WTaO<-2X+=UWdIid}rm#^yH`Rk}-G;aO z`Ax-Ag}rZvAAoYJ1?YpaoxrBw|A^o7>;SY8=;PGb*q>myCOmqvo|~P%<$w(I9hN_ab&BXcBvPC?cOUUa{k(hZ{x7(4bh4?dWEo%kduk_7gbg~ZU+`+Hb8Ie>|XwP+}GDPPCwN0 zx4Gv$ar%p83q(vl00LJF%{>c-@J=i(bbbFv3Zo7x&}_g>-DT-KT~5?@<=>e+j6MAs zT58&_`hK2ST80f)9E8N}pto0#fV183GCAE_p7Yp`gpOEbWFHtc+6fn0?xoNOQ&v{y ze_*2R3>{%_o2efcwOwiHt)kX;UhvAxGxXv?SshKB*X+oG8#vJe>31yinGbk)4&m*u zS;Sjc*Y^dgHw}KNnK@dc_`GzRMDF_0c0tE9(aSTz^wqC#JhoHZ2pI*1!@a$q$0`rS z2vfmmNi&*lig&6ZWzj14$!5ffl2Vb({3I~}35Q#^Zb?gbpkWa}#9?I|DC0zmp`qcF z#OlgQ4ea*4kvr*3gua=dle077H|27M-kV#0GonqDZA|(f2gFseXLXHX8XA7=bGsm0 zAy(iCS*Od3b2KzGFuFV7+M=Ty*?a%=r9xl)l@EETjm>yQ!Ptn>cl&uskbwBg55xsj zimxr7>yFipA3SeXmBD5|PIQ0!HnhzMv9rIQD&S;=9csTZ@p-xS@zpz9zQSz}*<|@< z%1DZ?Lu>?)1drgIg#UOLJ7)m|rza=e7J~+G9f(FfzxCmQ0YLDwm-T1UuNv0QUw?(r z4j9TC0^X@{*xc0*%7h0-teuc+8b!oe5;VzN=89{4-4&(&{w?4W6>EM{96rEevsrEF z5jygscgM^48c>{mLjTBObt zFIOf{*t)me15IQC1hhTxau@}aSv{Q0XEf6l_E;bm7jk7yq`&F=d}bN!njK19@eFCc zab(neQLpV;VF-NjSdSe5-@5K3es}&=@6Ol44jr$Ce|zDIRu7Z?{C;uBOMH?O)FcpZ zw}b!iFIVYpp_v3nX|^qK_2E;B9u7HE6oun%senIRw(}a%v;bnG<~!yT6ns4 zMUP7TvwCdrYqq#P4Y3ZZ>*|q7uay3eiqqk*qnQ<21tgltJa--yr!aMqRc*2*eLj>S z-lq}|ssiT06k8 z6t(N;&SuV?Ha6y(YdcK@x8vgVv%V{2|S+wB%_BLCL%n1QbRz)FRsQo#v*#|5$ zu>iaR@fSy&{HgNs;bSP_O_rMpTDS&#{Ro0=X@mEpKYk40oOMe$uYLh2& zd9j8`6CS+u89d?Ny%t)ruNR_@$v(Jg9F7KLNfk=Nv zB$IDg3Q}#!o-J}0snYaqtm4U=JI<>k?-&QB- z)6=)^?%VK2X9sJLfP7JIfPy&Ko&&Lb3O?SLtR#rEtmChSPBLsSK0!-*?^PwSy}dm| z&^)A7@*oNYRTS#rC!719qC!W6HJhXDhKv$HStCmWKs46=(4My#Bmvv}U4pr9bTM=~jH9vW?!bj*DI zyLnHLu&2~Ia79vuEH)vb+Hs{H@Y%GRG+7=Kf!ETTE_@DG0%;6*N41jC85~f+RD|!N zF@x}K??EH#wJDg%j*j~`%Bn#O0X8akc(n+JNd1)JR`9k_AV|cBT>3x_5>YywN#L`V z?jLU>N?P#*Qpy@0PB(a8L0w1}bWx-w zPEJi_)BDdJU;Y&0$oOtkGn|U!Ms{%k6%tYb#2SwsmSi9?Ix{o#b9h(@!gTQ8p%Tq$ z1PZdU*m!s><_JYQzP z9X>mGIk~Ay^Ba5~UxYNubl*Y>`048ngi8fjdsDd4$mgg6t3;5DFo8q>qLE!*LZTHc z{l&$_zTV!iE~cE5AU>a+oWLHtW{2_CBZVB`$PsJ2>k58@!CeWhqhSIU&t`JF7!*XL zl-k*J#6fo2-0a8Jb(C zM%16mrzG1_Sf3)Q6LMFc;=Q9a@^<`B-b--&k8cs98=9KtdDcR;0{Vl6fpGz?fhjy% zQXDcIN$dw>@8wLLoq08@D=M-rzXh-pbi~z{p!KlReFT6yxH=yR3-XDn8|=g@s1V~J_QHh(i~bpN2m^89^F9pKwE}PL2&>i!T9-PA-zV19=J<6{_o(lEcxL2j8#>XM3Sh7(n9$7zuFdG^XnQ**qCMia8uV!q zj7&^8*x2)eU%F+n(9j_k0llky&)bBA2n6g|A*c0`O7j6YIFR1ezJBpmYqM{B`C1Xq-=ZHbRBr;H%5JB1p3cF}C;&$B&cZ5)<_` ztL8sIh-7VhJZN5j;)P_9?kAMRn92_hqz{xEHE(9_G zn$g#OE!eW6k3XR@72t4fm{OIhI@yAe;RFAlq5EiuG)jP)adWKybv&w#6 zJSD&2o|%#p%^j_(-+~{yBl0ab{}H!s#|*N98L(pa%pL9+xVf&_6&1zg!p*o{rBV2X zFJ`+q2*ootE$_)9l9pOKT3dH0pXD8HH(dR2-E9VE4f;&bRGV79abKjjv(psU&V4|f|;bEcGOXa>)qNqCy~!7pkZ4AeEUJoR=pJ=#7$tM!mgPYEbK0V z;shCSFYj8Y5>U>o*qof0Dcq#nWC9GaZ7dk%1CX~zK<7MX8m-{!#qkz%T=@emH~NRV z&&D2ODjf=i*1XK^{}4tY4T`vdu?Jj`7;15Oj9hz$Arp{=ZR2v6EiBs5_Sd&V2y8LT zu>&y~$qgtkX9pscp{%Z+>**s#r=g{F(iHsUy)^ueA+|gp@Vx^YVC4`by%r)69=i6% zpuR(8xd~b9KYyJ6jwzgfIff9o*CKZkSa?nw6C4Nl(|wT^KXbcsW#dlX6>{R@n9rD^^*((#^at?a z;_&Y{jBJ4H7!?V2hI9R|-I5*Ew^73E?ZZGz;~e>C&lUwg$X(m|kje+=GlkpoC+Ns} zkJCNS86^N1;Z}H!o&=sXa7|kckkXkaA#Kpv%}`Wh?!a1 z#>U2;>D$)U)`q={$W&&8VUF5SV3l8+J6P9%Oac_WyDuuj#pMVD5!Pr1ijJ}}VNV5M z7X)MZ;V#?Tv%8FO2nl&WuHgE{yKLP(-i3O79WjP|uS*6PD&o5H8xjuJSO!K4bd(`% zX9wmE1_lOfKdZ?KIzC}TGqXQ1kaCTU*W(Ln(2x71R?rWVE?~S|x5=(9FI*Ptysx~# zM1XiM1_Ef+I&5GF))%Q2&_e-GS3g)CdSJIc0=qgtHy6KnMVb*@wz| z3=9kzf-bDW!u1ffoFRF#{p}$+P?Evmx*6UEYqAYf5>ZhDP`H;IPBX>)Vf+h}^u@4}C`u!!diZyI1Q8KYnt;VXcoXZ(W|l3`>0a|J3|jP0)P*& zmQtXd2p9?I0g4t!TWiqdH!(56s`u@BL%`N_b7TCX!q6P74Vu~N#cvF8btWU?ts_<_ zP21a|z)E4VWOJ3g5LtN-66~-PR?Pv}V7Tgg?6v_B16YAM0N9XtcJs>~}M;;|Y5Pdkm%H0%i#7z+qQvb;u}O5OhJG6OIdt*vpdCl$6oc|eQ`gwM{7 zJ;L6}sjsI84dF)xho?l(e56oC&}Hibz`q#bjAI|og9q0GE7sN}>P(g?H*pSQD;jlQ_F^*HZrrZg&$s%^2|gV3_WZ`ZFg7(8p~syMJn8!HpSjqGI*(Ha zz*`W>5W(W(;mNbq$x?oQ_vs}0p>{`QsC=-M6eV@BJT%7zKPEm^A?1tyH9MRBeKsyB z$zt`l+z4POY@*l8Gc&rI$=6A>z@VU)4PNzl(?W5b&dzjY%||>7`pc5w^nZ1!_4LUuyirX6^IoH6dT}p0O%))x6%m0K zWC(}-{9p|ZgX7juOL_TTaM4hn6sg}F_yidxI8(FVUD(DI>RflTYiep>3;#>RvLrZ zpx-K$LxdAJ<=?j#5)}o0$R2n=J(V$=(=|QI4nc+s3k%R;%S^kxf|UR)4-E~46Y1dO z6d4u=|K2G;c zm6fC6BNsqY#MKpKJOtcoYvA*EczC!SZ+U|21mQk-rc?h1epXgim_LEZ6KL*_7C#qs zT7Un9iV!DoCcPa1r62YH5Q#h77Ujkr1kaxJfXWLLCn6(@PfT1{Utize_6ApHX=#ao zQb0>v+X;-E4DZXX-rn~;%K&Y&Rn{;q5f&a^@}j0JG!(l=LGwQkJB?|u&4od_11sj> z_zay9D2Fy`8@ly9o+IvcFvE1qoP(n{YX&5brabYx6 zRM^ui3kx(cGC;?w620Xtvl*&I)Ss`~=zH%iVY?h9xF@Hj^+i%ZQfY>gHsh}i@idWi zHz zwR=EjO|^M>c{M(2+-j)Maa)O-bY`UueM^|73Un8`ORnaI;Wr;$&si!n&NlGiQ>%X< zP>h7gO3~;mhGG~PTWn86^9c_RS7*8?#U80@j+@A|O2Q;!y5&dpWjN5g3`(zq{~A&8 z!JEMcrJA2+;d&81hDk8Tc=2z40LBDcL52zo3&YGW2e3a#=|*|}16Y_!J_amkOdI9T zpiR+MQcLgBsU@<$U$@(~Ly8j^)2nDC6w|rvkdp6293hvt;IhMou;5O)p~%V-FX!YC z9x{J%Fsvld@B){5QN36~Mr8yts-h<^z#Qb-#D5O)^>a`;wX*(z8tuP6P!sA3ZT=!&ex)`0a7)IO%`004fd+oP6J3Ouny{mLbgey0`<_ z0#FMKN~L7f=(8&{jKEyyg+5zNO@sTf9rn{2SR0rw)AH!bQ_rt|_SRbr-WDt($ZK?0 z5~YEc0va*nLsbtu#PxvRnB%22LktR17HPeIfFB#XxEzCGPS{*r&T#jo?E+r-eCEUP z4`>-Nk%p<{NX45!fN{n4VNBruhvS$ly~%24ShpM-qTRxlDj_c%o6E6L^UccNG|Fc* z5D0ltK!E)+NG1+E5<o>qIbaA!6Sg3%m2AYr-7x$E~6$B(8hFsnhY^?j<|8}+Skj5C8 zlbUZ67H3P)#H6KB&Z$tDPzjK;0DFN0W4ZZk>)tF#N2md702a8nSl>80GJ8KlQU^qf ztWps{Aklm0_I(r-l;yd{;z;UTO(CUkls9i4fq%4Sd^bn6Xapid>`mMX!$9TW?08_I zBBJBtO=0#LFR^2Pe;=Bw96~rxKYi~89T8#jVWh+04mbOMKgiw(BnlNP>$_ue!k;7n zst<0{XrbEp7ZNio+*uMOqJi%`5L#$(e+OUju!|(27Dr|=;r|~$KH|`c2a3leC3U)x zM!$bgNYU4kF`K;0SPN{gP zJy`$H2U+3^0?{Y}&IZ673_d~6EvQrgLZ*Jl$$82*DA1#{)#Rdbs6 z-@W%>=o_$$>nz;6?7bG=5ZfK&6I*{sWk^X$?Rj-xEp^2~X@II?2x%G!F9X={nrs4& zQ1ID(01F4IJDA+~+J}N;pY*pkH-CXPx|ZWmgKC}r-i&I3Z#!TW^*rOAmCk)Nf%zHV z3vb)DT#i;eUXJY6t&Vf_`=BpCVd?@b@>5ecMOQ(pcx-f(41ta+=H~Xz*P5-rt&L7V zASJxDsmZrP<}q-p?VowL4IJuj6PPH|zdLqkF-s(9gdlbS-uvl0R2T#Rw%r)t(%d}z z7ANPW3MmYn+L@WP*TdIIC@CuPYHe;Of%=DG9TY04ymWMQ?CjB>WDvgrzJLvI1YvP& z*!|~PBRIZ@QObYo`wnab3jJ|nc-sP?jELvi7myP$fXLn?1PaFg);#DVYeU1Pco12e zi;G+yC-yA?=s?&%v-Ggz)zlE7XK9q+Yn)?!W)Xy)4<~l-FG$s|5w;E&DAm{19gxDe zd4S=_t z2YI#KFNkubfk1$TEy&B$`{7XwyjiGS?X}^Zf3!CK_f2SR%Q7}RP^-i*HfMK$auVI151ND9vLxTn!2Nr&op4v z9N`M;M@MH0JOnZv4E-1rGc(A34t9520oo=fCqs?oGufQ;JpKh!tiaPf!9#`lQQW|~ zi{XsNudnmy{r&w+@l5>uwI^z6_`O3zL$ z?GBGs>$cNgYpcYT7w4AW&8Qf4*jTV>rFUS{%F##wQi&AG&4~B{NO~tn;Y2mdn zfH-5d*dtH6X7GqEp<6JJAquQs3_5om8s~pMHGp$n7=^G%(I^=`n{c=dwg1Oc*B%&E zVTza9WA0+VW&wgYv)53&7>YifnV_sR4FI=wR7G>Q4aVJM6f~l{a_ArDWKIn}O*jB`3c&O_o0b zZbhu|1G)#$Y5ZCNz6|sooQTmRRj##m7JRoyoLnApc$%-hNjF=IcYSCv+@+*Mq2lMCgs*d1qSxiDzXr*?NLebvLUpF=l&vFGyxhmq z(ytXT{4BB{Mn)m20?Mj3FifbHiHm`anQx%Cy(zD{z>i*?n^a}vG{^m=NJQsA51PUL zzQgSlm6G>XdCdjmUsLqjmGKZy(t{q0<^|F)Cudm(<@=aw2btZ*)M~fwHmf)_K%gQD z4^9HM>wL0&FXpTn(6L2zgHyb2(hx>x)_P~zD_mXa`6iwJ%@0(QqSzQvtpEW6*=#j6 zHAJC)!{tn^BL`8o!_RvJQqt0(eO^!zN=Qq;fC1U*Y4FzJh66d=*w}z`e7rh@h5)0e z&FQ)({onucnR9SY172(Q z_xrwI_j6p&>$;u^`uh3^h_4uSE^W7-H>PM-6Z4j9s;t~!cMl!t&0Du#k2?W$z>PrZ z1-~GK>z|`dehv-}Mn)eoEkLaU>53cB4#Eo*gFu{j5Zzx!svJ3knxQ#)`uw@Rfk8=O zVd0Y}Wfc`@9-^Y6hA0O5`kMEHkwtroBZjLT)8Y8|czV)x_`Ce^mJ8%6C@8?hifJ_k zA}oesIU%8Nl`L*bm$4T+ey~0|dYQa3y|mO1Lm`K9T+Zb^IIb<1uCMTJzm2!b<5R>i zRxT~(2!vEPzV5Ftam_UE3oK=;xJ#{JpG`1dl-wv;^t1<1rqs(*vDWq}` zybZiIOG0e|2-vC$l2=e-);Mtc#Kk*;Qf#Q{QzMMSs4&>{Qdi& zT;!ro1O?vJ)z#KE7pO5PZ!jG?_gH!dAv`LI{e%l&!0FR;q;s8isGKass!ea8w1~QU zx7xHMIr$SRr5aNa!Wh!*(93A{lZL_ux*sp*b-cUnao*r~V*2HSSaK*uaku9c%Nw&Svv$<%))^&2 z1>2PlKYJK=PmX?~jQ8)ZD!$R%gTWf3Y$YBFXrP0y&6t6#`1$#jwfXZaqU68}4^?8M z$bJ4o-usM9Oc#V0si~uYI{`joK4vQ+At7vBF^zcv!5{U%d5sGh!M|mkl9VLnET^2e zbWn5pyBPq6Fwzo^Qh7Hx*e3FG_O`{x+=iIZDM8tDMMVX*zdtvXlr)w#Mvt_9d-{lx zntZcy)!{Pt7^QB>z2qNKk?|WW%E1H$e8CJm%^gHP8X15j78SdRdAk+p)R)n855_=SMx@7P0hHY73ZCe zF1%E{^XT+(CNeLBJsRRdRAiyC$+^+KY4Sg!sF#myB;95EYtq!1N!3OW5#?Xs!9^et z(8I9NP?+3saA*}_Ouvy)QgT*Tm!9B+st_Q~4uZFj4_b0~x4rA{pLVT9*hK^i2-rb5 zefI1Eh+zamh7MS?8e|j$PP*)2Mgbi^7fcN=7#VyrG;gt#{57NwGv8;5Xv zJ<0omJ8v^YaQ>?nw?3Jw#2x#D?@*Lh)`gj8&s=sFP|fNbe|(C$m(t_V4hd6f^T%-| z8M#HdXz+V<9%Q`f#1;JAt}Leu*iilPm0@lkb#6gU4jG}nv-91%GirS26%=-u-D#Zb zUSC>T`t->dq63<+xO?};rl;E{78TE&5&X!yg8=R*KtY-J?v~cp*_D;6D63J~;(G

REjR z!i1T?l)hB)NssUU!j?r}DL>cWr|tIsltGq%C%FfS?+yQ@{1a~&AB*q1*O-30rpDtz zX?~`ywJn_QN-{D^MpRpLRPF!tHX3X25eSdEa<$g@O;_#h>@wfYc${jiu6}qnudtxN zcYP@*^J{Tb6!o1~wP2nBLw@n|ADY)Dj1LJffZKtEgpYCn9BpfB2_X0Ai{$;bypb2+ zmh|Pz7sI&Un5LC=ivkW@{Qdh~Z!g5~;t~=Ay}kMQ`RyGYpwT1813-3hbtNT)dnSZ` znt;#XJPbFUAYq0AY=RXfJ5`hNJ9|&5$9NTt1D<=*=X|GMiVT>leXWz zS5%CYOEj(D5L&qQb#T@==sEXC2V-8^;}=^Get)(9%fGgP&sKkQhWjdazKx=;O+Q2Z zsB0_2B+G1kdzkWCERE7zJ3BiJ9rX0|fgv%8nEuH!$$Cn2i=R{h@ftW6rY%ST-rk!? zd_F!tHKsRQT$YeO@n*>?s~6(O1y^x8$~~uTNJLHAOP>u;GEjhfzdSn~Ef9y`ESEtc|wF#H21gol8*g3Q#QJ z{q~?g7|NkH(sOVS4mfQHBp;aa^2&Uul$Htt%oP{EUU7Pxs$)<=FB5nx++PNm`G2o_(}cZa=C3+)`Er)hYUg*S z7>W_*qiz|ojVDsI$UI1*Pe0HdoRhM8*+u(aSqs-^X(*glgE?Z5R-4Ow9 zPs^26RriD_GAhHJas`dT>({Sw6uy7^W>LBna`mjNOP)@(l=hH&@ZkV~Jx(6}f&{iy zzQXi@!MEP6fAm#L1bmb=jKY6y7doIWM&OCeBncP;N;vcexW3gMbdhLTg(8WQ896x) zWzc564hRZ3R9BeyuIAj zZl2PR8{OoelI4GR?NZPQ%h>k2UI8^XGRRH>k(lT4%R*WC6>+rj@y9*Ya~!RF*?Bnz zj6NhTb{6mU95IRw;px`n&`PbUn4x{7+cZ2xJ)4)8Wp?30&~q<5yOnX@nXFtdKoo#! zaMcM334nd=QNbCeopTrj6{>ZAcj~`$@^%2g92EW=Eie&KFI_4rE zYdc34K))FN@!?3{%LjIhF*{$cJXIDiKYq0}!QjTd&w=TL`d)WAr^Nj6nA6WEDyCm+o zYP9ch0u>Fx{^{mfL&Jyna~VPuBktUpSoFX(4EFUgOFA4NoK{gG?+i>=<_WH)(h43q zb4~eu!rsIaXP8t54?jFOCNMW17evF-@gqK#zJfcyZ3E#C_=&~v$t=bAfzLy0$wBd+n%oG*oZ7CIKo?6{d4BZO)}HZ9xwA|TJu!B zMSgsA`kU4=#9Fpkw1a^A)OuDw$NH*NKD8LV1`{cnE`NT>yZVEdRgddkydmou-D&l` z@$N0o!>PM=a8cT^N=7&Xrv*m7gJ8#!T*kTEX27j0ipuH+*$SulDH)m7k8JFmn)Qk| z-vhKD5cCe|54aI2whUApcAO6>bI3FaXKi?FtH>T*2rjJp(b?WrdgRTgMi=h-?Oh6H z0{r}BkE)4Yu)ePuiUV$;J(4=hF3xPdv8;cnu4L(b5UZ!%5^!E}LjIE{Dd?HJCO{_P zW0o2?ip&2J^oQxwPT#>V9Gn+JFR`6fbLQVVULRA#?V#zOE^;1ja^Kof?>&FM5-qz? z>Z=$kPggA$ceL`R94ABfm1Co_jm4+QfeVon4zRL<_d?^970d9LJx9<#LMQgdXTPE% z$M0VyB1a6}eHkaa54;qk9wMbwFw4ea$84EE0RIV6Sg*caEwWO~VM~ttx%^hzn7ruH zQa@Q;AE2mHVH30PpWa1rFzFK=1JgrF{)z8F2gt5x8Qq9;@_0ECcJ#UO#|!K~_}Cdq zS44y&aTIrlyp+u=xtluuUGj!jcX@l=nK|AzhNwr+Xck8I5Dti05F~j{NOq)Up9z&A z(+V~_!rNrn;(D)*k%=dQk$&HWlx^q8y+0^-Dwt6y>mp@=+TSwn@oZK?`F(L?XT{yO zmm1Z5Slzcf>#h(8d-ukH(zuiGl=%ifyZAtj18OrP8eLcP_AoWuNi9Wz2EO%!|&o_s>>^` zEG#TQP{k}>5LlwfuQ2`$8QW)>chdcMYpdm_&7s^bgFn@aN}9?;?Bs#911R4>GuS~$ zm1A|^cP%=aaWfWC<3y3*u#I4_@V+v|B!<2S+f(DX}WFfrk9w6}k>`8TR;zTM^DQT(3r zy!#xZb?Wo-JL&QP9mGb{$JTe(@(7EFpsELT99AbYAx02fZ#8WAU)xrw`FY{Hr%0FA z_NTXM=D3Fqdm=M4ZTdbJH)tp*G=rC#q?OxK`H!P~YQIgr$w|HSp`WB=;`v2Im^={( z$R6rpyVNfUR??@+`LH~CQd8B@)8l?C%}mEIJ-c9?YxytTaI=$)Y5trU-0ncfLrs8! z21G@)-uGRV%%%Wop#Uc*C(lX|yt;h+`%|OJtk=qyPS#0W?f-Rdt`3krv!+3AUQX<;KNs>yUPrT)6cyD&%HXk!(x?Dy|#CnzzI3>kP?LFj$O0OM*DuHJIy!a4u42akL7@a<#m zxvqI#!qn{bZ;w+B90{lbN~p4^`}gh&YZ+hk=>HY`Oh@S3q6bMp`mx=4?`9Y_?*a>7 z(08}d(A4a(epXO0n{>E2bf?2yqj$2jjj_&)49u%e-Cj~MH#dic2<}5a<=`{!@x;&Y z%goQ8-dI6k!VI+F1?l^*N~lxt*y2|{04p&*-UrDClPPozA7*Vq z`;5MGXKz=|E#*ndUoR+3MFDaz*Sv9Asb6!J z>%Q9N;JfBkbNcRwM<-N8GlU8aZjGs| z{k?6$g#J(dnbSB_ptf7d+bu89N;e5|SqRBodc?&T4a`hn|LH-crDNGg7_Lc+yMH`4 zsVX!1L5L=70~!Q3S62c7{!-Vg=x^CONDg0ays~_U@pk-+uLDzhRy~^!@442Z>tS`b z0X$cxZnd(!%pvdl`E|!MruHW%N~r!rE6sYfz{~>a?o10Wr~a|!5ubOTt)tf+Sr+CB zGc%uoa_zv%5hKV?z1qJwLPzL9>S0AT8o}n$obpeU=$$wj8KHy#t)OgO+cPybmZ2@8 zWL=9bYLxNT>Y=xRI!yK(C2zr_#u(;#MFk4+*x2(ir|a}M%l*CuAA51-3=X(oXZkbTtPve(rCkNFH?uN{#&SR6JlwoXU~oq zX3a0+4$R#A{`-BG@V8)R+O(C*#u=R|F^s@)Y+#>G&-?ZLi{}5&1yGTek*-WN&qV2QH;u^T&v>IW~r4(=_ zQ>&>0Cf`@y@-N_Zq7T&}kR?K*MKikQF>{$=Fx$m0AlT9`NFBr6|hmpH^{Yu~Z ziNUn{KT8sxyc3fn^tBsj)x&lwH1whLc>MaDh0$4#5_);Z4;^dW-M~L$X?D+YTdrtWR_Mg& zD<=Uv35uFAJ&f+ZJQXrhqb(4}&rVOeA_FT$XSTVw&r&qK={*FOL+k6k>+*K{ryDez zWXZcuyr$j`Qe9vDqERZU+~6(3%gZ|x&XwQTlTtDb#)>+`-qKjwPNl_dmmbs4t&w!rJb#aG7y7GpxT z;>hbg7agFvi>09m^kfP+JqA0|?~nDMUvQgroHN*^ z(4Ygx8jyRhn&ueI-llwY>gTZkR40aYX9Zbit>A6w1Il>wM9Brct0l6ND4PyY&bpvU z16!QU>Ppj$?xkwpRH1L6p=!toXarq%4*jc1`b+BrZXe@fRj46B_$-RZESXhD5}5u z9mm4T3ieF5z#D_vsQkzpr=@S|yF*Cqc?N7+e_t>uvkHMN+H@1Z8H6%uSV)M@hFznh z=E7bxR2a$R^*k`QL?5cF#!1qFGBdw|hj{X`sFmea z@I`w}MwqsxquCbVWlVfPDZtFuyT)iSqcz1&IUWw<8bIH8YFKMAou*n*di(Sw+RF6LkPk`+b^JB*ZMV zIrn7p)hWfD>N7pN&ZMxli6v>3{rKzE|6BU$veL_OCy!I9)4}S2!HV25ZF?z6Q*@$2&@`DX3ZaX2R3IOy@CJ>g&G; z^Um1k8W3VC@`+)^2qj(euZteiYzwMJbObT5+8iAn9g>a_EFvJStt`{#-aC@Q z?I9CzHAC85_~iAVYPpMT>gs$jTvNQG2h~OtLqN!}44mcft(0Zi(N_{Hjw=fDxDXE(^re#0=v6_lp)bJct)8`Y4)HuiP zz^m;_nxwTgxcWU!BN*=mkwUeSV-Qn-x|28$6OZfAyLizt%LxGqtktw#XYLjh7OGn> zy}Ys-dwe_!={S^ zCJWGCRnxp1Vs>K9*+m-Ix7XBfkLp^I!B<1+f@I_2D<$j-JomHt{;8Sp_1~jh^f(on znX>$RQ_J_?6Q%IB3BPOMyqKVgC*9Sr2eytXb9>XP?ltQ0MU^wRc{ja0*(qQg5H`Sn1x$*os9-aAmBlP*H@1UMOHdDJ0rY-j9-|tma?+Jyr(fZ zBqwr#%8i^PfQwqyo6G0En!Z@Ml+eC1va_3^B?b8#LM@?#2LX030D*(57I1D{ zW$0g+gnX_KoK+9l*?q2*^Ogc4T0gx)Zxh%ufBMXl42%p7;bnDX5O9n4i;EgSUW|sl ze?2np=ZyJITjQ&}n>W~#c9HKTn8;wb=1;O9pd%9CnjvLSfB$&_ZyAb)qn9{83C=UB zzX&j*n`L%!aRK)?J|*QJ;t8<%L-Dg;zb=7W0vg?C)ED43#Jc)#zyJJs6a1^;s-I(3 zP4OqEUVqR_z?3WoTZpPP=`@itX#Hb8f79AsFcD>a5AAPUA54{3ofY|0bY6KIx zJJ-Wr72ef9VX*~_<@Sv!ACjLWUr+Ubx2b7#X3cp#JTWoH4xIW9i9FTr# z3Pk|O!Y7Bsrq5~e%*v<(@+3ZfzvToaq0|1OhsVdcNzcrrve5^69>SCT#`1FR)bbJ~ z&wLd-F~x2bPQ?@kakiH`|<0BbGATSyO^9 zA&Pn5Xn(6=F9-9XeY`+y9ph3`uEP}GYaxJ#Oslu~xP*inv2^8>Jw^Y_un6Zh`N#?B zhZ~EFAv^y|>RD)J4jf<^iNE>b?ytkb!cb=iCxRXwTc9w1keG<@PPbf~`q?PAw)QVHm=ouk4qRB9)ju10vDNIq z+a3)P4Z(g`%)*BLR9{fysrI0#sOX_X z;v*iNWwS;a&dncGX|t^l$x+%pY5T%*L_~xJa@v^TX&#T28ZJlB-v2ukiFj8~$JF$g zv32!rd3big^UPv@aQ1~5?O@1*-c5<#{+w<~SK$MG-7goqz}g|yClF@asg$O)-<@(k zfvLtle4K$=jz^S$>EBB@Z}4pu#(sT?x)2pP<*Xan1su8i3>x*J$xFrKLPqeDA`jxPj_4u(dnQ|(31PwQk@z&nWe?Y{YKw&DLcbwu~N z&fas5r6CRUNbG_A?nB=|Rg3_dCsL z|7u;T=t%!1*425PPI`T4ukXzZ4a5`Z$DOKlfSw-esGkc9BE~o8KOH}F<_xxEK=D`a zwRTq6Yn>VJDENYM_XG^DQaw0$HF2IPkdH>W52}woKV%VnVW}B4HIA;?eibCt_;`J% z^a`3A4#>JMb~n_3x7gRcSCq*= zEU1alo7cIzMxiPFACR&D!>JXHh1o6k6NwGjr8=Yv zGy9>?5E(vB`aH&ifg5ZQ-%0G5aHNz@Q*Z z(!fT8-oGR#XCc5c9;!FuoQzcHIfsd zd-G=bgH!E7`rvX=l6IWi#lyH8Jv+zfLD~EcYZs<8no1=;1`rX27%w zotceo6F%Dxw^1#@mwW!pHP9a4lHDmyyvKH#278QLkOtv`F)yNh7dGGs}LmT}Z@3pFa}?dHO03 zO^i71gomH*Gx+Nl#dhMnzW&O_3erQ%>^4#iQa5p4!ntn?@^kn)5ljrupSN?*f(>#r zO3JB{#&54nXm-`bJ$i(0RWs-O@Wo`0453f)Df>Qg$=)AZ^q6@P1Y(a$mii%o6g*g+ z(b+jSwMSXk73Nh??!6xvKY>gl`;kNKE$4@6S-`3Cqi^O0SR$8&2QMkAe(?pOZ27?_PQ2?nqwvVb{dO< z@pJ^LO7)H-9+b1vr;*skGuYHGCo zPah2;6bm9lyrdB!9wn6QY;282XL#;oEH04n(y#x%t98x^-{ssUOZV&fOe}!f0chdz z?jRWM9e8l?0+It!yNXLK%+Gh7{Q`eWP`u=Qf8VBNF+9=vN>90f)s%ahRo1$w&rXE$ zu!@|v;K2(HeC&Hj)lV1pEE(e#^;5sUb?oTT6StSLFAIOQj~MLlH?DTFWctAv*qSrk zBd3cQa9Wz*QHrF@u#2I28FT}aPgmO=B%cS9G49@rd6W?H(e$9}#^IN6Y1Dl(U=noe z%Yh^7XJ60_?1AA{a9G$$)2=nimO=vm&^hYaTeWNZxF2=V?A|TzY?P_CBjnr-%?@dC zx`9IYCX|fS|GRWl(aSBq%NFH}4zW5Yo%s^+3-cLN4t1RaDTU2VrB7Kh+vkE`#{67Hc(POGYM-JC>2F*is#Skavf<%BQ8qa;Pqaa)xUgsLa?~+Sao~X${?o0d_;gk1aCb8Aa6ru<6*LM zrKW~;w~9f+%MP05fm(oEeTyI8a9+Hd5Z1yS*KuCFl_jlb%B5yrWJgS!`?spWsJ-LJ zjrH{jsogE>ltf3{l{pH*?67nZ#Cr)W!Z}Ipllvy|$K9BpZ;WX9MWm&NK#lIqUGe|) zdeVhx2>|Q%)rN*TftPrtP+O9<+FbMAvB)HH+tZT>yUm*EFq*+aB|>x1yYs$p8VNd4 zulHk_z{w$304qD*6+B(SpN)-j!=hVF_=^A@i{^Y^iO&}Ml==4DFx%jxASCiys+fedK?gxeZshN_C z!36)qY4_G&djKm!R#k5!K@|!Z8~C&66i=#v$Dw!`haoXynKVE6j{!1&e3y2q@?3JW znzr-^1F73_xikcb$%5SP6lwJKmp*{COhlI4)Yj^bvv6?u&h+!|F!D z^~)i}vmHJ`WBg}_ii|5kH7=mGYQbIE1N`!u-`~R+xg@(du_1?#od)<6m@2R@tTFJv znaj9do?2G^$oxJe`OYMd(7(ZI;ErxhBsE%i(>gg zKfT*U;r(-H1c*cZFJFkB6!07YQ5sQ`h9D#)L={_nmf$HSq#lhj6 zff<>4Of|LAJDY!6la=tTLn%Qe(3o||l=9v>7rTj_?ygvuO@Tigp|7@%z zK9FRIjt?L96F_GVL2hlru>z!p{~;kJ_F`fYP}(q@t!qzy@5tg8)=4Ny;(RkipW)=@&dP9QS|Cy(+OTd%OS9$?N7#zCInX{7$FX z4>nJ_fOG>DqJTOX9Se)Go?dX<1U9(=+eMB8gZCz~A-Wh60wNueZ-rLL-CLEYsw^vW z_Q?4mymMK}#l^zHqOJ4Cf(uKp<(l98x22ptWoa+{j#*-_9g4%#r-^jnr^hjrhFt^~ z!2*`+K->}%6W5lO7%soou%2N$4Zys@Ah+y%!cBGj7+2*b08PN zAFQ9q5*>ucOZ6Vy3G69yKDD)#U%nr7Q5sK5Nceyz^3{2VD%%qO!j$2RBpn?0kN=RdttDt8`G);+&S*fFKVM zTSZycN%LJ<1NWa!NWo(0p_VPHF5IU9gYh>Mx5`mIwx+G?UGq3~sZd@G?-5y*A0j@! zb`T(C0-;(SJ~McSm87VEEP^aOoF7C)b$$QIRGuOhm+zZF4c;9R`!ZoLi+B8DkFkkK z)xVbgn@y{JVvCk#hF=@*9A3Ru2Wf8AzY87h?N7mH`+3%jff#XVW6$PIEN=H7S$4!Z6}El2;(_{ z{QrL0QqX{X(edU@Wf{N0PUUsOuTLG-O|z$NDjfK6W%Tyswide*oi+Maz@Z;{d$H*4 z+v3(cmsb;ul%$HRZbSwrppWK+YZ9Mg>^`7kMcu#m$%0%Psg!RI2QD;fmMb#zy6Q$D zD7S`h(C!FDk`NlNUEZH^8ZeMP)OwV$iHQ_@Em_;GUFx=5<>M*g0EG`U`kEY35xD`0 z3QMDKJ3;+svauD2yJlh0+c?&!`--fzu-7ZvLxz5fCH;%Cws!G*10GDgG+tl~499^- zSda-nFeJz6-SHjPl_O2Y3JT*_4n2X4_#$Nj+3&UCSc$8lOY76kdpKz zCV`-GKecZLp%}HsEhi?F?UY}61+Sp7& zqKfiOHACm=a0aRhFr?V2xI{%uQBfGTElAmXuFH!-P#X6aXXoSPB_=A6U;MWF0ir42 z2GYNX8Mqm%aVc+T(Gqr;qF3yzc#NB{lYnUfDLFX=iBd>5%cr&KAeV**pv0%o9_mdu zBuNOQZUNKNSAiP*pkXD9;Tw_d0o=<${|@*8^p~=QEcP<2<3T{`Uy{#d9$TSn@Vo9B zkbXmnifDI^uf<34xNG<#_`_<#^{IEjX*y`+es8Xn@GYnN`Whu6N%_ba3XrbwVdw`i zZnVhN-Tg|LA>eZOpUip$8326=;N6zK!nLRUr@VJZP*2|(S#>GtVvH=0Nl0iZxUoI| zwQ=i)*RQbvBG9s+prT@n0L0!Ch|okzlA(9a|M~8PZ(AwI%=G*Vrvd9fD=0y-^E7f3tE&J;=20# z`i$U%fI^LEdV{+Tt5y;MF+W-~lk_EIc8BSdPR`Apg!phcB)^UWW0T1>cmMQTn%dfK zH*b3V{!Sg%QUkI95?9%WoHW=Jmf;ry00a^UkKnpVO+dxS5lo#P0`ku=9H^wY#PoAV zHmDjSBR+WXtc3Tf`Jp=w!7S~vV2y0y>B-#;M}&oO82H)s&T;YbPQd&nH&^mnKOfEu zCJby?@eKPk&@g6TAs$YB-Pa1?o33pMqK?=7P@%^NfsxUvmr8b z8bbKi>XX{qT?TbbZucq#Uj_%0!br#Uakmro=LrhSUoKDMy#M`;@0o%rVP#Pf1`hfX zWN(p>z;h8@`@3SlJLLQ!M>!p4_ld6uPo6xft~zjY6Z=PTYc7JaS#0#odRLUnMvmfB z)&g+WU_ks>+~$o`>8(|Zr*JvCmi11>!=rjS$>o~Bd=C%6TSs`|SiYLJmi^pKLHl3L ze&2PRJ3TmLtKX!;#_8M7Za?^_AU|MB?+(IHNSP_n zVz*e5A(qCrNh>Sxb&kPl8TSOA8OBFNml{|U?t!O_nG}jFGq|zT*Kb0v1q)};GOV?= zsU_8Vk!YbuLtUV>(~yLK2z5HEPnwMeTwN>(#ast=V?@d1v(p?gRk%FJ%^8*T_4L5y zLxxfq=j&>3}w69*Oeb<|dJAkS7NxsnwsO=BB1>bad(2*-Clu z%uG$OwE-awGYvd#5J-U5;Ry7hI$9Vi+sVql27D9m94cAtaWHA}mBqATZ7t^R-3RgU z7%Lt>c5JR*@L-NCeRt9^)vTD+jJ$S|B`@+IODij$nec2)?!D}s{F^ydN+%R{mV!Lm+UZsVgbyUB=dFI?6`Qd+oiw_gcnzZ4M-IXmD4Lk4VZ9UvT(IN9FC_Tm z{e=>FY|Mi20tRRxB#+t!;0ECsfTm)Lh1xK7Lj45~hWoLxid=L6=oUtsL|K^G2t+b4 zL`*-D9Z~CJJp>G&Fem~QA74ojL<0)eWqkFrv$OLg^vjBhBm@i8x!BOy*w_e(FiynM z>gtESzK^gXKz@P&CxL*D1O6p^3>`p;v$L&XiiB1J6VtiuJ*L;r1*43NIHb_EFQg~y z<+xAZ&!jgKOIKW{W5TKb6lP_;bIBpS!Wftxb^m@}oQvK?$d^B@!&3gm?Zie)JGo%gP+;E<^3_P?RZ z%nS`Uj6=mA^)E*$ZZ8^spb;QL?I&1Vziwt^^rWPOIJiIbw5$wj5NLLBq6mb?*bbSJ z0@XOGCI~_Z1Vp4!_tX2>4~Hvd!gb|C*cm-zP2{{x0+Or}?It#5h5PEEXF`B3Kyr`F zBy#cmfk}mgg}s|`r3vF@+a|EOZdPA)+ zhnfaXM$3l}k{Zp|G|cN$>by2=n}dy&`53NrC5LFA5i*X*yb8m?&`_1mT{|)4G>`?1 z=8$~ffS!`HZ+bV+e~coFF3zE*;{9qo_l=`x&gH#(qBkq`)9!MU{BF@BN21Q=;lcV{ zU<^^bcH54)4>o_x;s?us8qQy1-y1LboKsfx4002RxT^f)ftqm z863;(mYHj)M6bRWO^{dB<=4~I#qj_UO#Wu~o|RL``>+{b0KOp*L5a&v$y4x9aomqpNvT0AM3KAF7kLCfWOxX> zj$iz(>q5*u%L_^bY4 z3|;-YpbclO5}EYXtDF)3r-nsn51ngF1>#`Z_G!tEX!(ea&_U7*=K7E*0qsNS30?Ym z@hE>r4Z!oeLyEWV=9sg(J;i_q(SZN-qV0~^C!007xhm@FGgF$OZokZRpXgyC1qcwa z62S_yG^8(d0VDxLgJIN*XoS!K+M)lFU^&9YMNrlSEsl#0TV8=ggBk|!JaC$!RL#h+ zw6NfcX^TrtgaivtS~20_`6VSc?Cdnu)$hl}A$sE>eMB@zGZCR=g?As0I&`G#xi3(3 z0VIQm6ytx4Y~0hckEItuF+iS8llS95Nzi%hJ1P&r5`zN+R0qU}0tlfHMaMbckix6KxO8LOvIyCeL z=K;k!ay4o@BnA*p0jT00AcGuXb(Jxsa2_>RaL)mDbJ;AY(QgKrDaItTcg{)t>TNmW zdxBo$J);;zy7f1gVLP~hPv~TI=^&;tqn|(5!Z-u_toI^OF4~=UxtS_wyszgk(V|$ zZY4hBB6Z_EpG*m=Ji>(5^&h`kbZ>bIJ0tA zw7_4#7N}i*zStMT0XwSm{q!q_LXoiJ73Al)j2ap}XZ5QT6iIxU_+SHpCH0o93dAO6 zW+~ZO)4SoUGLKDH6Q$W@pZ;u&TROdAdZK^f0`ICqehzgY7oJgy)SLJR0xyZ!0iIhQ zlWRSV2%K-DuMeq1?hhXwG)Z~s=}t%(fKYx8l}#jS3^n=rfx{LSrmCaU3sA?naR6L; z)B`g!2S41~Zi63Rnc!nzrTUhWs0oOo#z!$9FJQ3wHceiB>6nZHMPhPtEh<~oTJP|7 zAoD`wCW`Mv83M`&{JqqztiFowTY{OHmzUQMGuVWhRyr85*6lq4T+tKI!+i86?#^}E zZ(S&-hd20SWIO>aVUw_xa5?`qt>9A}qNs;DZhM}jBcS!%zOimj8{{2KpUEoi|joE(N4U?J}?Tk`!^54Qh09}^6=RRd!hw6jQ zeZ9Sy@$A1GaQ;r`q#TN-*x`Y^$#8?~GLO0{;q&gQBMYG{O#bXqK><1cj>hVL6kWV7 z6*3mXF;uLQS`y6RB_#{6cl52W>FoQB;3L!DVRYt9);zNsg67DQWotAf(_7zv@m%VE zNQt)zFXZ1{`p08i8lwG%Znw-XyTFhI9!JoWWgg$NU-*T}%|Evd?kId3`nW5<3ce&LrnWoK{S z^YCca`|FTSOJ1T<(bO#DJu`_n&+>B3^WDmQVD4Ae{o)xkC!C zuT5~6d!GWa3$qZ7o6*hRw3v+{8sClz~brDhnUWRUx9E^lw59@)w5-Z}ScUf+k!kK)6C5 zm-!O1<&5m?R@-til7OnJHT06``f34HCy;@wiBcDlcXx;)z8B8!c*GSI#EKaC2Uc-o zi`%kpXq#5B69(uuiXfDo*y2(H`{Gj^ARqk00sS^L1z#EiV)roSh}k7pj0XlPBCntr zMrW}fm?ZK8s7vr)Qier&TN~`8+%clV{d#@F1T=(;+S*5Wd0TpWbtLxNI}({hApf+B zwv(E9Tb+JW+I?@Z&`J6+5!@v@X2OIhggPFUWY_KNB3s5Az1K6dv(bPD$zCkU%?*u; z>YruYi$dc#hKS|mUf{RDmImGAgXH7_G z>_hRU>>z+5LYyMOiH?MTnK~Q|NAUYG)J0cCPD<*H%f}mu^MO}i9{dl4$Ff}sS(%xr z;Bj05cLLQ_$z07C>9E zf~5i-EqV+bjZ3wz>cEWv-yk2twim4p4jYmonk`^-cx~~F!302QF$^vTUMpl^q@hGi zw!u7x6a?lf{QNn1s6Y(tF;)I+PQQ-=KrjZLDC?lcwy=lqFf3gMhlU=;#l_vfznCTU zcVh$V2cEo7B>y}vnBasLRM}>WK`yOl-09Qe>+as6A)zlMxSYAk$&k!xpl3xU$4~+# z467r-2<;}l{~YDvu|oJC^2CBbgk(@hEORo1WY9e!abyS`V%&{7 zPM7f?sFyg|WswtOV>*V0;hFhpVXo#;1}@JH(kl1CiM}W_Qz$tuZsPm*(9ElFKqxCa zn%Ie%GP-vZOdDGhmmO7g8SOu{uHLwMHM|&ir7d}vE+d>}NAPy&Fitv2WBhkUR|4+F zcwm%)t*>5X2qAjSjgJRh_r{EeQ5k8g9t%ksl@ZbAr>4RhFx?OEi7TG9U}tRE2cjZz zQAA$e6x`uqmgKy)V1|%@ri_R`;A4Q-cM{~y(gdT|DxE`DIWMgPX7{}&cIV6M_I4qD z{z-GGo6gSbwv<#Tr69feNwEp{){V{eWAPfVx_3O$vjdo~qVgIo20EIuk{~b|H_F_V zTc|Lly(n=u3R8tjI%a0|ob#x!b`ann@%;I7+|?s{T2<(1&=pWqlM$p5Do&h0FvCtj ztTNk4K)KnN_EfLp0yH1l*}>X$hXe)J{&}<>J$e+qrXjIw8hV(7i8i#X6^0=RXRN*@s10pWODFA#K z3F(;CjJUgHY4#p9_Bive5l7^wM{w^NkY0r3mX}*sN-mQs;5Gp&SeAMiV-foW6#mnXGxEdYeGP&{E^FXkZsP-J3B+p-&su(v>E^jfP@q`5eI?nHc^8 zpc1ViG9bkll$JVvjL$e0o>Bf5p$J2dt+cPdJo*NW@vNAbE?&CyC0d*_G4rhP(;%%W ztQmoReE+S65P&Jp(7!)_j)B4uUnUqJ)8K_JdH8@8^pJjtRrr09~1W@PzDpkoxQD&7N1_r$SZgfulE<+C* zUEHdHdK@*sA{{9!tMJK_tS77RmbQpwl|7Qr><4fKM<Ke}j5rL9$yx_x-vtdJtZn%3jE!t`qD zx3YO7TxYeks8%nWQ&U5KZHNoRe-a4Feef?XNrhdCAA{#OAu4(MJ=8k<|6vS<5}RL=-}3r1%)FVdE%Tt|9Pf z>3$p;_5Mb-XTNEqx6yK414=r~4KY*u5E{syWxSyFqg8BLGN}2**`o*!D2`F+*n95LhnSgaU*!(aE+_3!KT&5q?`Iy#r?&`yga-uc(S zPjhlBvg`%Y4=##^=8tBEj&1kj+sOT`t*vNq(G^orP~d8zyq}$g+g=yS8!%E(oh2$} zJeF~TDsYFFjNFkU#ZR87pE<)TEPPib3(kb-OU>(UMX6-nSK~v@#a@;0@NiT|OY57+ zeB?n#m$yIFtv%%1mZcdp7gOq3ZxR0;7|KyD{=DK&qG@ioxu?tI%tJpt1Cv{i?B}ZA zNIvR(oG1TjHpcbD^CwBR9}ZO|%n8^OFC`@CXctrLJ;+Hfu18AW&DrzX**}-(6lvO^ zxlsG>kKg9ordHQOL$HtA5+8;uF8uQRst+ zj6GedHa7gx6N(u+_G6A0F1$nhcx6@A)FgkO2aDd~I#3$IaHiF&3@y(eh)gqd4#17k zS{#P8;Q>v6L%{n(Q6#}a-SQM90{aqE%q4{P-;wrsZuKe>H*kA%9UkR)`5AWlMKo{1 zVKsTD24zNQD0x?h-Vrn-y|ukB4!~3exR@ek+ajcofY30-7ZbC%a-|A`U{o)tGOVO7 z4gvy2rs_#KnTYZ%E{h`Q4zh`=NjYYYwEweGy9N-mpXOZg@CY=CGqGFMkIz zcX9D`G#E%kRODEFS7KByASc&YRCKp`&=+}^k{dLCny`t<$yc`JNIA?@p^6!Xn0kZb zhhop&y~!AdZINtCAS{$5KEM#xOk>PAubb6{Fv-PQ}|fj5k=+lF^{1woppCn)7QP|^#!F~e7`_eepVQO%whG)m zV+bM3?92@IGT~==wkR=Pn{Mb1Z?{17S#Z>_zsBBRE8%J)0)B7eoX?@7M|F?6e1$zH ziYw!tyTJBC96Bi@v%KYT=@Lh%WwKX!`0;5IrP-aL0@d4+(@&33ZITIG4$K&CU8#s; z8ts$5_2WQaJa&+vLxY73JnICT#!!SJc4E}-IcoH?_j@x z7XA)VHPEGH8%jhhh>f)}oF_c2OX8?-P9`Z69PaHMo10_ZyO%eyGvZ`NtJTr6ormct zVfrFpNV>3;`@N(bYug80jR-H6OU#=%7#Wpv&Lbmp%6ckg=ztK6(FC0$z+rTqx1 zSov2tr%^;<(x`HnnQ-$c?Zv2DuCD0%>_yH4vIJQYxTL*QLVRabk^2=1Dlo4agfW$}!KaE!R=75H1RBT9al(!a5O_N3ch3wRUQrMKh? zcE0$`yyz);?p#lAZ*NV?q0-XQG+R$$f`tJ1twGHXhkF=y&VH}rAiRw`K{KL4iD-Q$ zG;1e;3erm_wlrI9PJ(d6iRim%yN{)uP4m7;-FsY2toPF=ja&N&+`%~XM@PD>P=lcD zZO+{LE{H0bcIUCH(NP$D9Pl+=#ylNd0Jf6btel)ol%%a%f7hQusf@FaO#XJ2V&Co% zZuihG-kkF&LbaXB*B$>`Luu2_1}n${=7ZOmF}ra5(YAe4+n6GUa~V254wUR-bJ*{B&}1l#lHUP{Lex?A?Nq;pox2n;DVQaWwx^5qA5(MTRid;!pyEkQII!!^EAco@(!3qar28EN znE9IMN$`95c>LDb@Fy3&gd+Yhpl0F5ABE&mpwId1muIR7#dSa!`cXd}?aji%jd%^4 zzuGZTQ2=+IsJG} zv%YD~7NA`Kse*YCe3u{eOu=XlV_Kl-y!x5vpNLUCdib!{hv?x=IgS?u@QR3eqeC{! zGH~RZPEH3nIBvPPC{p5u?7!&!wx1K~vK?*K z1yXJ8x7HCo){)LHp(yzXm~5V>Bt}ls6om$bh@JSldvRC1#t8Tf5~44}iptBu%iGC~ z@e$MocJ)>j7}06+G4KRaa?8uh$6YtRh=&6{3y5C`Ttv+A`ml-yaU_JX9c2dOTOi7u z{e6%G5vQBfS|>5dx|e}ve7U}%0h_kvyr-e3n3=iXVQsC?!w_-){Q3Vk;2J961cHYO z@ahS(r`LdEQouC_#lR*%aMS>pz;z=WZmr+9ZJU|4_G#eR{Cd;5TMq%Z;cwCcUJ|Y8 z1Uyw1xL{G{_|c?~z_W}`XdP>EWaFIR7@Ic9F*~xh>G)&dN(#jm;Nhi@fTKuWyukS* z;5xKf`k)n$I+01In1Dwo#>dO|x^;JU8a^%CnjoPj)@>ROT=qEa&gUn<^QLyp@3C0u znxs07*Qi6GIatBA`PwDmm7EL=d)GYr20VL{AwlZ`h#z5^0V<5vNMS1sn%y`U7!Is* ngp^YYK*h&~EEH|G{xgKm+5P|h!R19jQy4s5{an^LB{Ts5L^5Ao literal 0 HcmV?d00001 diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 171a4f8..1b3233a 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -391,6 +391,33 @@ julia> fit_WLS = curve_fit(m, tdata, ydata, wt, p0) julia> cov = estimate_covar(fit_WLS) ``` +## Visualization +A `Plots.jl` plot recipe is provided for visualizing a fit. The simplest case +is plotting just the fit curve by itself, which can be done by supplying the +model and the fit object to a `plot` call: +```julia +plot(m, fit) +``` + +The plot recipe can also show the confidence and prediction intervals by +supplying the `purpose` keyword (default `:neither`): +```julia +plot( + plot(m, fit; purpose=:neither), + plot(m, fit; purpose=:prediction), + plot(m, fit; purpose=:confidence), + plot(m, fit; purpose=:both), + ; + layout=(2,2), +) +``` + +![Plots with different intervals marked.](./img/plots.png) + +Changing the keyword `significance` (default `0.05`) changes which confidence +is being illustrated. + + ## References Hansen, P. C., Pereyra, V. and Scherer, G. (2013) Least squares data fitting with applications. Baltimore, Md: Johns Hopkins University Press, p. 147-155.