From 093faf3fe849ab8eabdbd3a6714aec26e8f8f8c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gi=C3=B9=20Marcer?= Date: Mon, 1 Jun 2020 15:46:25 +0200 Subject: [PATCH] =?UTF-8?q?ex-5:=20implement=20the=20code=20for=20plotting?= =?UTF-8?q?=20=CF=83=5FMC=20vs=20calls?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit move all the functions to lib.c and lib.h; add this two librearies to makefile; make plot.py more easy to use for passing to one mode to the other (show or save figure); create fit_plot and create the figure 5-fit.pdf. --- ex-5/lib.c | 86 ++++++++++++++++++++++++++++++++++++++ ex-5/lib.h | 35 ++++++++++++++++ ex-5/main.c | 92 ++++++++++++++++++++++------------------- ex-5/plot.py | 25 ++++++----- ex-5/plots/fit_plot.py | 31 ++++++++++++++ makefile | 2 +- notes/images/5-fit.pdf | Bin 0 -> 14478 bytes 7 files changed, 216 insertions(+), 55 deletions(-) create mode 100644 ex-5/lib.c create mode 100644 ex-5/lib.h create mode 100755 ex-5/plots/fit_plot.py create mode 100644 notes/images/5-fit.pdf diff --git a/ex-5/lib.c b/ex-5/lib.c new file mode 100644 index 0000000..92528c4 --- /dev/null +++ b/ex-5/lib.c @@ -0,0 +1,86 @@ +#include +#include +#include +#include "lib.h" + + +// Wrapper for the gsl_function structure +// +double function (double * x, size_t dim, void * params) + { + return exp(x[0]); + } + +/////////////////////////////////////////////////////////////////////////////// + +// Results printer. +// +void results (size_t calls, double result, double error, double chi) + { + if (calls != 0) + { + printf ("%6.0e | ", (double)calls); + } + printf ("%5f | ", result); + printf ("%5f | ", error); + if (chi != 0) + { + printf ("%5f", chi); + } + } + +/////////////////////////////////////////////////////////////////////////////// + +// Perform a fit in order to compare the data with the expected function: +// +// y = a⋅x^b → ln(y) = ln(a) + b⋅ln(x) → Y = A + B⋅X +// +// with: +// +// A = ln(a) → a = e^A +// B = b → b = B +// +// using a linear regression. For b, the results is hence compared with the +// expected one: +// +// b_exp = - 0.5 +// +void fit (struct bag full_bag, double* p, + double* a, double* a_err, + double* b, double* b_err) + { + // Expected value. + // + double b_exp = -0.5; + + // Parse arguments. + // + double size = full_bag.size; + double* x = full_bag.pokets.x; + double* y = full_bag.pokets.y; + for (size_t i = 0; i < size; i++) + { + x[i] = log(x[i]); + y[i] = log(y[i]); + } + + // Do fit. + // + double A, B, A_err, B_err, AB_cov, sum2; + gsl_fit_linear(x, 1, y, 1, size, &A, &B, &A_err, &AB_cov, &B_err, &sum2); + + // Parse results. + // + A_err = sqrt(A_err); + B_err = sqrt(B_err); + *a = exp(A); + *b = B; + *a_err = *a * A_err; + *b_err = B_err; + + // Check compatibility with expected values. + // + double t = fabs(*b - b_exp)/ *b_err; + *p = 1 - erf(t/sqrt(2)); + } + diff --git a/ex-5/lib.h b/ex-5/lib.h new file mode 100644 index 0000000..6c7bf8f --- /dev/null +++ b/ex-5/lib.h @@ -0,0 +1,35 @@ +#pragma once + + +// Subtruct for the results storage. +// +struct poket + { + double * x; // Vector with the numbers of function calls. + double * y; // Vector with Plain MC error estimates. + }; + +// Struct for the results storage. +// +struct bag + { + struct poket pokets; // Values. + size_t size; // Quantity of values. + }; + +/////////////////////////////////////////////////////////////////////////////// + +// Wrapper for the gsl_function structure. +// +double function (double * x, size_t dim, void * params); + +// Results printer. +// +void results (size_t calls, double result, double error, double chi); + +// Do a linear minimization fitting the model y = a⋅x^b. The p-value of the +// compatibility test of b with the expected value -0.5 is returned in p. +// +void fit (struct bag full_bag, double* p, + double* a, double* a_err, + double* b, double* b_err); diff --git a/ex-5/main.c b/ex-5/main.c index 115cb1c..4b58526 100644 --- a/ex-5/main.c +++ b/ex-5/main.c @@ -2,43 +2,34 @@ #include #include #include +#include #include - - -// exp() wrapper for the gsl_function structure -// -double function (double * x, size_t dim, void * params) - { - return exp(x[0]); - } - - -// Function which prints out the results when called. -// -void results (size_t calls, double result, double error, double chi) - { - if (calls != 0) - { - printf ("%6.0e | ", (double)calls); - } - printf ("%5f | ", result); - printf ("%5f | ", error); - if (chi != 0) - { - printf ("%5f", chi); - } - } +#include "lib.h" int main(int argc, char** argv) { // Some useful variables. // - size_t dims = 1; // Integral dimension - double lower[1] = {0}; // Integration lower limit - double upper[1] = {1}; // Integration upper limit - size_t calls = 50; // Initial number of function calls - double integral, error; + size_t dims = 1; // Integral dimension + double lower[1] = {0}; // Integration range lower limit + double upper[1] = {1}; // Integration range upper limit + double integral, error; // Result and error + + // An integral is estimated with a number c_i of function calls. + // Different number of function calls are tested. + // + size_t c_0 = 50; // Initial number of function calls + size_t c_f = 50000000; // Final number of function calls + size_t factor = 10; // c_(i+1) = c_i * factor + size_t size = round(1 + log(c_f/c_0)/log(factor)); // Number of integrations + + // A fit will be performed. This struct is needed to accomplish it. + // + struct bag full_bag; + full_bag.pokets.x = calloc(size, sizeof(double)); + full_bag.pokets.y = calloc(size, sizeof(double)); + full_bag.size = size; // Initialize an RNG. // @@ -52,23 +43,18 @@ int main(int argc, char** argv) expo.dim = 1; expo.params = NULL; - // Print the results table header + // Print the results table header. + // printf(" Calls | Plain | Error | MISER |" " Error | VEGAS | Error | χ²\n"); printf(" ------|----------|----------|----------|" "----------|----------|----------|---------\n"); - // Compute the integral for the following number of function calls: - // 50 - // 500 - // 5'000 - // 50'000 - // 500'000 - // 5'000'000 - // 50'000'000 - // with the three MC methods: plain MC, MISER and VEGAS. + // Compute the integral with the three MC methods: plain MC, MISER and VEGAS. // - while (calls <= 5000000) + size_t calls = c_0; + size_t i = 0; + while (calls <= c_f) { // Plain MC gsl_monte_plain_state *sMC = gsl_monte_plain_alloc (dims); @@ -77,6 +63,11 @@ int main(int argc, char** argv) gsl_monte_plain_free(sMC); results(calls, integral, error, 0); + // Update the struct for the fit. + // + full_bag.pokets.x[i] = calls; + full_bag.pokets.y[i++] = error; + // MISER gsl_monte_miser_state *sMI = gsl_monte_miser_alloc (dims); gsl_monte_miser_integrate (&expo, lower, upper, dims, calls, @@ -106,13 +97,28 @@ int main(int argc, char** argv) double chi = sVE->chisq; results(0, integral, error, chi); - // Update function calls + // printf ("\n"); - calls = calls*10; + calls = calls*factor; } + // Do a fit of the Plain MC errors. + // + double p, a, a_err, b, b_err; + fit(full_bag, &p, &a, &a_err, &b, &b_err); + + // Print the fit results. + // + fprintf (stderr, "\n## Fit results:\n\n"); + fprintf (stderr, "a = %.5f\t", a); + fprintf (stderr, "δa = %.5f\n", a_err); + fprintf (stderr, "b = %.5f\t", b); + fprintf (stderr, "δb = %.5f\n", b_err); + fprintf (stderr, "p-value = %.5f\n", p); + // Free memory. + // gsl_rng_free(r); return EXIT_SUCCESS; diff --git a/ex-5/plot.py b/ex-5/plot.py index d581431..7ca0ba4 100755 --- a/ex-5/plot.py +++ b/ex-5/plot.py @@ -4,26 +4,29 @@ import matplotlib.pyplot as plt import numpy as np import sys + table = np.loadtxt(sys.stdin, unpack=True, skiprows=2, delimiter='|') calls, I_MC, σ_MC, I_MI, σ_MI, I_VE, σ_VE, chi = table exact = 1.7182818285 -plt.rcParams['font.size'] = 17 - plt.figure() +# plt.figure(figsize=(5, 3)) +# plt.rcParams['font.size'] = 8 + plt.title('Plain MC', loc='right') plt.ylabel('$I^{oss}$') plt.xlabel('calls') plt.xscale('log') -plt.axhline(y=exact, color='#c556ea', linestyle='-', - label='Exact value') -plt.errorbar(calls, I_MC, linestyle='', marker='o', - yerr=σ_MC, color='#92182b', label='Plain MC') -plt.errorbar(calls, I_MI, linestyle='', marker='o', - yerr=σ_MI, color='black', label='MISER') -plt.errorbar(calls, I_VE, linestyle='', marker='o', - yerr=σ_VE, color='gray', label='VEGAS') - +plt.axhline(y=exact, color='#c556ea', linestyle='-', label='Exact value') +plt.errorbar(calls, I_MC, linestyle='', marker='o', yerr=σ_MC, + color='#92182b', label='Plain MC') +plt.errorbar(calls, I_MI, linestyle='', marker='o', yerr=σ_MI, + color='black', label='MISER') +plt.errorbar(calls, I_VE, linestyle='', marker='o', yerr=σ_VE, + color='gray', label='VEGAS') plt.legend() + +# plt.tight_layout() +# plt.savefig('notes/images/5-test.pdf') plt.show() diff --git a/ex-5/plots/fit_plot.py b/ex-5/plots/fit_plot.py new file mode 100755 index 0000000..a2d40b0 --- /dev/null +++ b/ex-5/plots/fit_plot.py @@ -0,0 +1,31 @@ +#!/usr/bin/env python + +import matplotlib.pyplot as plt +import numpy as np +import sys + + +table = np.loadtxt(sys.stdin, unpack=True, skiprows=2, delimiter='|') +calls, I_MC, σ_MC, I_MI, σ_MI, I_VE, σ_VE, chi = table +exact = 1.7182818285 + +plt.figure() +# plt.figure(figsize=(5, 3)) +# plt.rcParams['font.size'] = 8 + +plt.title('Integral uncertainty', loc='right') +plt.ylabel('$\sigma_I$') +plt.xlabel('calls') +plt.xscale('log') +plt.yscale('log') + +plt.plot(calls, σ_MC, linestyle='', marker='o', + color='#92182b', label='$\sigma_I$') +x = np.logspace(1, 8, 5) +a = 0.48734 +b = -0.49938 +plt.plot(x, a*(x**b), color='gray') + +plt.tight_layout() +# plt.savefig('notes/images/5-fit.pdf') +plt.show() diff --git a/makefile b/makefile index 783009b..b0ff4d0 100644 --- a/makefile +++ b/makefile @@ -27,7 +27,7 @@ ex-4/bin/main: ex-4/main.c ex-4/lib.c $(CCOMPILE) ex-5: ex-5/bin/main -ex-5/bin/%: ex-5/%.c +ex-5/bin/%: ex-5/main.c ex-5/lib.c $(CCOMPILE) ex-6: ex-6/bin/main ex-6/bin/test diff --git a/notes/images/5-fit.pdf b/notes/images/5-fit.pdf new file mode 100644 index 0000000000000000000000000000000000000000..d8bac2e762004af59e827c17b6b798260183d1db GIT binary patch literal 14478 zcmeHOc|28L*H`9CGE1qhc?fsD+{`l}^H7FN*F3uqPW}qd|!w9ASLOn+Qj!*ph7By`13k@-U*O15{CZMFX`!l7%Z)+1A0u(+Q4WYHqG)=Ss9E!KKKZ|6&J)1O>|jQ1bNjBKZLlpuGxUh3E@gNAe}wdcXoR@1*M*@klX8PR#5gk(lVN%J1(&Nj^+A zO4-v}kf3Xy54kxA)Ogp0i(uv3IgR-=htp9I!drJCwmD$ksE|>fz0| zEpT=|WH#Nk3^m}lvqH14Ys!^&>l(Em5%j=_)#Br&I8q@M0Gpo^!!7m$sJH z8;Z0Y@x0q1Zy)BLx-A7393-|l?6Z;A;cQw(Qu%)8Sv@#alm)i=Yw*CZnbbYa8Sc+@ zu9X^er(Xu~)_YVOooec1brVY94;Q%aU|-a(Ri@3V+R#iNY;wdZ5zbg_EO1?9V~b`m zPa>C$mgYurx>89Wj$J|m^g`P<&W}*{2-2Hs)B7Ok%Ct9L+R?|W`c92k3;BGHyfy-+ zb*$G=!X(H*lT}qeH}8fT%N7>Xs}2(p2lEQ_se%q1e<5p_RQQnAp;u&x9~Z%SHe+1y zIWC%};?ithJo}b%ckwYx_K0%g1c68&pqoS8dGU*Q2~&<1`>;TKiS;XMXIOn=-kU-7 z#0*u!{^p1jjfjE{cel%VvSQe1$BOVoWey`{Bj>X8b0x5`4CUVH0OtGypi0+9v$o>% zBN?i@SiMhN`X+{rup^>ND8+C`gSf(7xGTf~NO8BP}oVW(7!VC1(_xI;D zsu|i`O*tiYE?>xiM~dCGe2=j+EcTkQz*B+FE27CP30958AKz!FVobL=ik8B3b3}H8 z2wW`<7wYVCjm$_fivZnZ_*f^L%Byju95_}i93f_96dI@C&zkAP_R%@8G4Gq4{H?+RD17yJrea5u~(jF!T0Ouk!iYCV3Cd_x4C!r>&cDUn@6Ty6;(G48J(;>qGr7*%}t)^?Gd#dv~L(q zRS#_ydPz(u+#mPct@mW2j&Zn0Y=K8oPW-WHA*~O*Z0;ow7tVw=;YvwYM74|4U+=TN zcqV&MTH#i;T<+cPO-=Jj;n1;0Zx0`=6%y^=W- zb#-mu^{(jpL`88XQFPa^w`FNOuhaOTH8`))zrjn6f8p!Herl$Fp6-`#-}13#gtP3) zPX8&(_$0aSa^nLrHwPYDGE8`1x--AiZKCZo&3&f6EpAU!TX?Cc{A|iPzst+WAM%@g z!Sh;Yp;t7gmdUudWTy1^<(H3kq*l?T7+fge?CKD;rL{@AZWOPtTGcUrbK%Xq5O~tm z)mj>tOMzv~28MP`qdf3e?QE`}Z$80Td%W1LmpT#X8^iE^xQ(l1xJ_}#RPug3nE$O#)# zPlbs_T1Gmi+Biyb=^&$8IvN=ZMWPrk;c2f3MD5x{Q9ZOtWF3xH;x*+Z6eEE`0{kL@ z!J*(NJU}xjBwiYVGk=03Xp}SzVdQT|BEv#;FJF&kptOw4Rw5OgzDyi7FXJd70#U{~ zzj2t`r6hFL>lwWljGZmB$1x2twlAk5j%~%T;X{!h!abt$?u)~DJBXj`2HU88 zh*D>HD}o!!gpZzmIm|3{oneXM)GXKBCj~V^-6a7_MNgoL9cQh z`9kn^S-Hx#enc_^2Kds|*3{n!z)A@{J9if!e+WxE+xkK%*bhh@!ger(nkRsK5XgmK zsj7>kBN2k206<#6r7!>tLrPh~;XiB2%R_qd^&>&O!I3ZoMDJhe6;x3r`q}%sc$2(* z;V3d#UuvmqOY(IgBU?!%5<;Bd@8=buu4DUy5*mktBlg;QfN}gjwiPo2Vi&N!#PyOa z%v>Bu&S20>`+^nF9pYnIbV`2q5R@jM3DBm)Ta)PG{5Sq~TIh;2wj9;{b6fDB;j>Gyw+#3696% z;W!-V1doABVL@GKzyK0lLpk&u0}@zSr9gYAPl7Z7R7%UKkG&$j<2;`g~1^5N!coZ58aEUG!(8ZyFUPu&l1-9H|>f zIRVXwl?mdDd=L77=Eib@Mvg~H104chp*gTlA^;lXd4Yl?^(!sUiRB#31;8H`M$Vz@ z6^X3trRPwcp9z{5kXE7G>L&<*mi*<#P76-P1mxZQ59EMeUOXUE@FzozC9DGy%Bm<> zBU6W#QG%lUiIf&LlJe>Y(0m8?B_ z6~Wdb`Y2~t-3zt(iuBm(irw|awa9Z89}6>2Md~JH42x-VUuVdwy&6578d|qm3R_6) zAFKUf)5|B$`(SO}{);kqGWVy~7zwhQc&1JwbG|#JdA?n?pfe={{Udv#lIw$L( zQKw>DrRC*mYvs>O*SA(YP;%88t`wH-@ek6I^yF$|-t2JGcy6!3hp`UK-4A={F*o6Ai=9nwaykU)I!lQyVFN}oLSRM9WGxt38Bx^iJq`)uIj~o8J z_hh8ukeyk{VMZA*aYV0pV{$74TK%GKmULiaSZv0}Ln*_73me>L{AAV<{7U5TC&8n^ zJYEH^475Cf;BSv|%p9c{ki7R$NpEbibbA-Qjg7N*W*f`4Z7hM&_a7+@Hds83$(TFa zWHoh5mSYF~^<%;YRMwHfz4Ve^Th1&kYfeu&3(s`$uDyl2v1FdeXgPMZzh1@GOgpjb~!b zbU5v~KdSyCS^rU~xH2LueELENQ^kBCFRRss%%j3Iw7H5hWy9C9qs=1D`}=2Ki_O_8 zx?OJ(kHu;64YCPl_IBN{nERd;yzq^(z}P5z9dWIEf&M5iEN+#!)bww~s-PH0<*2#n zY!<|Wvkq#e^xc!&qwXESyI&!3KWo~F8g0h%WUm1YbtO55dt4C%n7T(?5lZdhd5W2z zxS7)fBcH2#js$qO->${P5H@nv@Aocvoc`IeXE2jin7Q2}c`8P^TFaQj$2i+ZXR5g3 zQS7s$`B&X5A18d^jSbh<5)?WUW+b*)C2Aoy@?4FrOmpyPk}kGiTtOw#uJwGf!7TD> z+((~^1P%X=rru0~M0)S|#W9_b&E`Vs2LFhj8%E0u7tL7i4IQuC-Z^+7-{foN!mEVT zX}+s{IZxi12^yrwKDiiwLv}GHd+Ix7-1&Xvv)a`g6ZGHiK_El7aEHTNm8*nUmm7|GKwb?=`Lftz(-ipU|^o!FF0A$qrDQ*>Hk@A;_kl*Hog zT~jf6lNR@_qHXmpEc&Pd9XptNUnRt{V8%MRefF`qDagIcH9%iA;NGn36f*S^Z8*~j zy6tLw*{tI9NLkoMU`Q)%*qgsUt51mm_0gv^RnL7kP1H1Dl3J#D1Ld8CSBbHYg5QY8 zSoVF;5uWUHxyZwI`R)$8h2skD+Gi~v?ut@lFT5>67xqkv%_qdXe%2vZ^YS^(qdSg8PNbs#pG=Oji{CD_mGfk677F&NI-T6w zS{TmJ9S~G_?UDWV?t@uhzkJ&uNs+_7&g-ypU-#$NL28vkbq#KLsDjtwHM`{-icz;i zSqz3&pS@I?=;9y1a}c5z1vx)AXfyWHT5P^Mm|&uEc!vV&14dk~X|P*i=w_2K`(2LW zvuXNVcKr*^RJS=v2)e#O2 zr}U%?PwT2%-;@l04vH5k{>J{az%M(>qi}#gIfG|CJWeD95+)N zrk_hl*)3-mdR=ibo49N4Yty5d1qKcgOR;t2_qXTp0CBB+A9Q1tSZm;PloPl4^z>Cb z7D+)}iw8=&neNo=-5AZOcI+e0RU0{e1I^*4FPfHorp5$m$?$yqLi2aDmd?4z$HCVr zk5^3j9PiTl`q?LlqGTb6=H=!t9zOHmIRvZzFU0UC4PyWfTFIiCp|)>5lCE!<-4YJY zz2DB<-fQ-5;d+Yq?w7J38R_WJazxF77x+|%@IjvGd^3j=PvW!H<5cUIUyJ6{T$#Rk zJGIWT(KC9_RL*vVGx<{JwksXQ*EKHUPpRW}@ejjQ!h=jlTE`9g9BM}cc;D?d;*9bZ zUsR#4NcL|QO{Vftdfmf6;>9V;@`9#v$hY#Qj_0R);}#5A&vsjsd(bM#Ob_>C#Ht{mkGiNV~TKUZ%FR0gl z<~Uy{4*%ZX>v;bas!6j+k~k>z0A# zpLkhg7T4t^mDv;=E{W5&ZtP6;hVf<5P-KTPUR{UiO0@dx^a1QcB(Pu?hms;-;V3i) z^RG^oz&_;f50&1cu}C8%9%s_{4_XSpvC_9joZpWh&LD8gM#OB+pF7xVE9=TYqNOMq_s^}^;TPHQbL%W(m z-h#fpA#eomiDMn4@;_pl_trVtF!W}yHxwWYVNv9fIU_%uq|_s4bk z`di?J0%riL7>mEjqJpAR^)peTXDy?9U;eQ3;9?Qi{GKW=3W=yvL8IcG!h8iUo;iO+ z9%9PY6*Lv^5xSPN)pr}G(jHV=;9gWo?)9AXsP54lg~d1g9aA2=v^3Wro^dwAB*{I% zy%(O;9hu9PttuL^IA6+;9IAL#`ka^xx~X1OUG z@0#6ccuW0qgVdpQM)}tO3+z-^j#9#q5GR0QY*~PX*-(s5Pi5`U>HfM%(lVi@E_QE` z9rUeL6$}%P?=x^4c&6@n{=sBvNo&M}W?$p(6K;)amQAS-oK(dJQFXF)$0qvbxOdE7 zYHg~dZ9w$Sa+`iAqHj@dIEAUI6^pACYb>`JxwN+oYLAMWJ$yS42lz!s+12dG1?0hoLy~c5 z52C7@t<&+mhWgJ#J*V@jtVSA$dD0u_)5RW8?#XC%zX1Q<_>T7?nz`W0yFJ3!4LAa+4A?Q&2a}`quO-t?fYK0^a%Cqg!{V$ zL+}AXzSDOy_HAw|($wU5Tn*cs+sLlnlRB4iGNY@zW%BizS?of*$1lDQfM%IFY$lo{Bt>Udm#rYZHQ+Xs+80?WAcJ6sLj&?K-``{>ebr*&29BaZl zy!yQIb3*GnS`RUBw}Bec38JyUw-YN z+0I$Pa%K6`V;QZ6^|PFn?2jZ$H6F2Cmw-_ju*uD#g*o`TA9auH4zxWZPg%;mSCli5 zHi?3!wEsY0^4{Y$GTYFvSdtR9G7H?#&_1QBnJlOEZP!1a@=~y*Fri(3^TEDF?5?>@ zbR2N=n?1C7yF_;gu*pShR`}kLGClDs-YoHSYOlQjtgGw@%TdCevXgg)7zms-Cy3`V#2G!D1#bn# zMbva;R_5C_Jj;4E^PN`BA~^SVEnxo^RIO1024AI*JW|W9B1y3|dEnL~|R;f`skCt+5m(VSsQQ|C1j7AXC#{sTercLjGl_Pgme_2vlq$WCPp0O$za4j zcfvfty0@|J|z z#domg4aXx@;JYP!^c(b;Q+K_ZZZxUtpGe4LwS0} z*%Iuk0w(7Nj(ybTRk>J8xm{yGT_du6B(SUr#dMniEo2=#lN~cE5J&%>!{dH+mCSgX zI7(QYjo-#C`Ndfc_PPEv2Bl*krfa#jsOGt)dEv7ClfuGoJ+|u3_n&>t56f(f(`Ve~ zP2KZYj_Cz!`$UG&J!4Vsv(m|%jkyvkdg&ty^X<+UnK;J`n95C@vL2h?aR0N+*LBAH zx5yB+W+y3kT8W7g1va~fA{3ZK1C-5QRhsDaEgsO$cxl~17lYzdR`DuyWjcFo)5BBw z+J#RiMNdZWVQNgbVam_z+U}GPWGB7l#g3|N?E~fSzOdX%G!5+-Qm+cD~5k}>NH>2EA8L;u*Dr9|xQ_5zA@#h)^10Jy{Q|n>{WZx0U03R4k90c+tm|fSU3W7>IaCx5&mF537Q-N9+lr+;=mjH%OwunC#ziI z>u*X|yTs8DEC3F9sOA68CB9-m`v+GzIVKgFm&@Ltm9XT0Ru z;L08j;ugq}^F%lrGVq{m$$*gESKy9pU_g)qd9fzj91t%=uIB_NpO=F=E+Fh-smvV? zY!~^dC!Bl^2%bWg3;5>ArM^JX&=m~W7!brlE?cr&U=A$Z0aJ$gLE(oDBwGr|Bw>G0 zS+am&8k2j1j3O2S>426$91EdMaCHa{PG8AaWc!I-8r-k-e-)}nfzJz;DyyU6Kp{}O zUsYG)S0&Wa2MtP~&azef`AG)C-WHrkfQASq^l~*|L6_p0eto&ITx|sjh~y7F{i3@9 zaQ!M=dFb_fSeS&Li<1Xn;)gI-RPb=|^anEkAtzlB{j#K}Upd2{T3#2bX~U zEvJ7GvGpfzXsUaq{wAF~CU2Zqhbp{r$Q~D_Cyu^0gl?^la zo;bNTzh+K=UVCrlV9?|Eu4`_uka*heFJ1Vk!`q&8s+_%``{_tv?%8+v zS4j$A%J%W9EATLth6i)12w94vXrf*-Ao^+-S+cCVYj$3U4YS{`vak-8`=tr+L5n{K zkOLajc5#3r1j%2%82$5(FmjtmZm7@YaAtAYVd8$Gm>xiiCMiHbyydU_3eyIJT zFa)^E z;dZw6ZUEUVwHIA6vjTmIjxZzuk}xu1@;?y&g#&M-Bm8F>6dnov;GRFqP#7o%Voezu z51b5Z%cP+982bHci4bK!DcTGBkj$Ys%0V z5Z|=A4233ulg8C$NDLMQRv75_iv|V_Vwl#Hp`{6H$AQFRrNE1`x-EL0{?IsRC0bPv z1)|T^@`6I4elJ5xK@0clwt&WO{9$n*KxlP6Bo+^?=&RZyLCDzJ{*X8f*a57qhlec0 z>b4-PYc0=6Job0l<8dIEZf#qF^lv;%L7R#-^`y|y>c5(<)H*bz@X(sSrXJ+@T2m&C zgzVeudeT_PZmlW90l=`fjDY*CjDQ6Dur>7vXb@w#whRwEacj#6AP#SJ848Jp3@Y^d zIsZ^dEM)aomtoMr75=k~l<52#g-I+yGN&=;Z~iQsiqRmmu=45iCl6 WBwJq+`SnL(upl;GP*6=z9rj;DusBcv literal 0 HcmV?d00001