From 02a1747fe1518d831b54b99a294b76a2d46e0a66 Mon Sep 17 00:00:00 2001 From: patternizer Date: Sat, 18 Jul 2020 22:51:02 +0100 Subject: [PATCH] Add block to loop over 50x50 matrices and perform sensitivity analysis --- E2.py | 359 ++++++++++++++++++++++----------------- sensitivity_analysis.png | Bin 0 -> 21764 bytes 2 files changed, 200 insertions(+), 159 deletions(-) create mode 100755 sensitivity_analysis.png diff --git a/E2.py b/E2.py index 9e13d40..22ce2d5 100644 --- a/E2.py +++ b/E2.py @@ -1,159 +1,200 @@ -import numpy as np -from functools import reduce - - -# Calculate the norm square of the jth element of the ith eigenvector using the alternate algorithm -# H is a square np.matrix -def vij_sq(H, i, j): - # Get the size of the matrix - n = H.shape[0] - - # Get the eigenvalues of the matrix using numpy - eigenvalues = np.linalg.eigvalsh(H) - - # Determine the jth minor - Mj = np.delete(np.delete(H, j, 0), j, 1) - - # Get the eigenvalues of the minor - minor_eigenvalues = np.linalg.eigvalsh(Mj) - - # Calculate the numerator - numerator = 1. - for k in range(n - 1): - numerator *= eigenvalues[i] - minor_eigenvalues[k] - - # Calculate the denominator - denominator = 1. - for k in range(n): - if k == i: - continue - denominator *= eigenvalues[i] - eigenvalues[k] - - return numerator / denominator - - -def eig_sq(H): - # Get the size of the matrix - n = H.shape[0] - res = np.zeros((n, n)) - - for i in range(n): - for j in range(n): - res[i, j] = vij_sq(H, j, i) - - return res - - -# Calculate the norm square with performance improvements -def vij_sq_perf(H, i, j): - # Get the size of the matrix - n = H.shape[0] - - # Get the eigenvalues of the matrix using numpy - eigenvalues = np.linalg.eigvalsh(H) - - # Determine the jth minor - Mj = np.delete(np.delete(H, j, 0), j, 1) - - # Get the eigenvalues of the minor - minor_eigenvalues = np.linalg.eigvalsh(Mj) - - # Calculate the numerator - eigenvaluei = eigenvalues[i] - numerator = reduce(lambda x, y: x * y, [eigenvaluei - minor_eigenvalues[k] for k in range(n - 1)]) - - # Calculate the denominator - denominator = reduce(lambda x, y: x * y, [eigenvaluei - eigenvalues[k] for k in range(n) if k != i]) - - return numerator / denominator - - -def eig_sq_perf(H): - # Get the size of the matrix - n = H.shape[0] - res = np.zeros((n, n)) - - for i in range(n): - for j in range(n): - res[i, j] = vij_sq_perf(H, j, i) - - return res - - -# Calculate the norm square of all the elements of all the eigenvectors using the alternate algorithm -# H is a square np.matrix -def eig_sq_vect(H): - # Get size of the matrix - n = H.shape[0] - - # Get eigenvalues of the matrix - eigvals = np.linalg.eigvalsh(H) - - # Create a matrix of all the minors minors and a matrix of eigenvalues as required in the algorithm - minors = np.empty((n, n - 1, n - 1), dtype=H.dtype) - eig_matrix = np.empty((n, n - 1)) - for i in range(n): - minors[i] = np.delete(np.delete(H, i, 0), i, 1) - eig_matrix[i] = np.delete(eigvals, i) - - # Get eigenvalues of all the minors - minor_eigvals = np.linalg.eigvalsh(minors) - - # Calculate the numerators and the denominators - numerator = np.empty((n, n)) - denominator = np.empty((n,)) - for i in range(n): - numerator[i] = np.prod(eigvals[i] - minor_eigvals, axis=1) - denominator[i] = np.prod(eigvals[i] - eig_matrix[i]) - - return np.divide(numerator.T, denominator.reshape((1, n))) - - -# Calculate the norm square of all the elements of all the eigenvectors using numpy's built in eigenvector calculator -# H is a square np.matrix -def eig_sq_np(H): - # Get the eigenvalues and eigenvalues of the matrix - eigenvalues, eigenvectors = np.linalg.eigh(H) - return abs(eigenvectors) ** 2 - - -# Prints true for each element that agrees between the E2 method and the method built into numpy -# H is a square np.matrix -# mode determines comparison with normal or improved function -def compare(H, mode=0): - # The tolerance for comparison purposes - tol = 1e-12 - - # Select function depending upon the value of mode - if mode == 0: - res = eig_sq(H) - elif mode == 1: - res = eig_sq_perf(H) - elif mode == 2: - res = eig_sq_vect(H) - else: - raise IndexError - - # Return the comparison matrix - return abs(res - eig_sq_np(H)) < tol - - -if __name__ == "__main__": - # Generate a random Hermitian matrix - # Dimension of the matrix - n = 8 - H = np.random.rand(n, n) + 1j * np.random.rand(n, n) - H = 0.5 * (H + H.T) - - # Print the second element of the first eigenvector as calculated in both methods - print("\n--- using the algorithm ---") - basic = eig_sq(H) - print("basic implementation: ", basic[0, 1]) - perf = eig_sq_perf(H) - print("performant implementation: ", perf[0, 1]) - vect = eig_sq_vect(H) - print("vectorized implementation implementation: ", vect[0, 1]) - - print("\n--- using numpy ---") - nump = eig_sq_np(H) - print(nump[0, 1]) +import numpy as np +from functools import reduce +from numpy import linalg as LA +import matplotlib.pyplot as plt + + +# Calculate the norm square of the jth element of the ith eigenvector using the alternate algorithm +# H is a square np.matrix +def vij_sq(H, i, j): + # Get the size of the matrix + n = H.shape[0] + + # Get the eigenvalues of the matrix using numpy + eigenvalues = np.linalg.eigvalsh(H) + + # Determine the jth minor + Mj = np.delete(np.delete(H, j, 0), j, 1) + + # Get the eigenvalues of the minor + minor_eigenvalues = np.linalg.eigvalsh(Mj) + + # Calculate the numerator + numerator = 1. + for k in range(n - 1): + numerator *= eigenvalues[i] - minor_eigenvalues[k] + + # Calculate the denominator + denominator = 1. + for k in range(n): + if k == i: + continue + denominator *= eigenvalues[i] - eigenvalues[k] + + return numerator / denominator + + +def eig_sq(H): + # Get the size of the matrix + n = H.shape[0] + res = np.zeros((n, n)) + + for i in range(n): + for j in range(n): + res[i, j] = vij_sq(H, j, i) + + return res + + +# Calculate the norm square with performance improvements +def vij_sq_perf(H, i, j): + # Get the size of the matrix + n = H.shape[0] + + # Get the eigenvalues of the matrix using numpy + eigenvalues = np.linalg.eigvalsh(H) + + # Determine the jth minor + Mj = np.delete(np.delete(H, j, 0), j, 1) + + # Get the eigenvalues of the minor + minor_eigenvalues = np.linalg.eigvalsh(Mj) + + # Calculate the numerator + eigenvaluei = eigenvalues[i] + numerator = reduce(lambda x, y: x * y, [eigenvaluei - minor_eigenvalues[k] for k in range(n - 1)]) + + # Calculate the denominator + denominator = reduce(lambda x, y: x * y, [eigenvaluei - eigenvalues[k] for k in range(n) if k != i]) + + return numerator / denominator + + +def eig_sq_perf(H): + # Get the size of the matrix + n = H.shape[0] + res = np.zeros((n, n)) + + for i in range(n): + for j in range(n): + res[i, j] = vij_sq_perf(H, j, i) + + return res + + +# Calculate the norm square of all the elements of all the eigenvectors using the alternate algorithm +# H is a square np.matrix +def eig_sq_vect(H): + # Get size of the matrix + n = H.shape[0] + + # Get eigenvalues of the matrix + eigvals = np.linalg.eigvalsh(H) + + # Create a matrix of all the minors minors and a matrix of eigenvalues as required in the algorithm + minors = np.empty((n, n - 1, n - 1), dtype=H.dtype) + eig_matrix = np.empty((n, n - 1)) + for i in range(n): + minors[i] = np.delete(np.delete(H, i, 0), i, 1) + eig_matrix[i] = np.delete(eigvals, i) + + # Get eigenvalues of all the minors + minor_eigvals = np.linalg.eigvalsh(minors) + + # Calculate the numerators and the denominators + numerator = np.empty((n, n)) + denominator = np.empty((n,)) + for i in range(n): + numerator[i] = np.prod(eigvals[i] - minor_eigvals, axis=1) + denominator[i] = np.prod(eigvals[i] - eig_matrix[i]) + + return np.divide(numerator.T, denominator.reshape((1, n))) + + +# Calculate the norm square of all the elements of all the eigenvectors using numpy's built in eigenvector calculator +# H is a square np.matrix +def eig_sq_np(H): + # Get the eigenvalues and eigenvalues of the matrix + eigenvalues, eigenvectors = np.linalg.eigh(H) + return abs(eigenvectors) ** 2 + + +# Prints true for each element that agrees between the E2 method and the method built into numpy +# H is a square np.matrix +# mode determines comparison with normal or improved function +def compare(H, mode=0): + # The tolerance for comparison purposes + tol = 1e-12 + + # Select function depending upon the value of mode + if mode == 0: + res = eig_sq(H) + elif mode == 1: + res = eig_sq_perf(H) + elif mode == 2: + res = eig_sq_vect(H) + else: + raise IndexError + + # Return the comparison matrix + return abs(res - eig_sq_np(H)) < tol + + +if __name__ == "__main__": + + # FPE as a function of matrix size or not + sensitivity_analysis = True + + if not sensitivity_analysis: + + # Generate a random Hermitian matrix + # Dimension of the matrix + + n = 8 + H = np.random.rand(n, n) + 1j * np.random.rand(n, n) + H = 0.5 * (H + H.T) + + # Print the second element of the first eigenvector as calculated in both methods + print("\n--- using the algorithm ---") + basic = eig_sq(H) + print("basic implementation: ", basic[0, 1]) + perf = eig_sq_perf(H) + print("performant implementation: ", perf[0, 1]) + vect = eig_sq_vect(H) + print("vectorized implementation implementation: ", vect[0, 1]) + + print("\n--- using numpy ---") + nump = eig_sq_np(H) + print(nump[0, 1]) + + else: + nmin = 2 + nmax = 50 + fpe_basic_array = np.zeros([nmax-nmin+1]) + fpe_perf_array = np.zeros([nmax-nmin+1]) + fpe_vect_array = np.zeros([nmax-nmin+1]) + for i in range(nmin,nmax+1): + n = i + H = np.random.rand(n, n) + 1j * np.random.rand(n, n) + H = 0.5 * (H + H.T) + basic = eig_sq(H) + perf = eig_sq_perf(H) + vect = eig_sq_vect(H) + nump = eig_sq_np(H) + fpe_basic = 100.*(LA.norm(basic)-LA.norm(nump))/LA.norm(basic) + fpe_perf = 100.*(LA.norm(perf)-LA.norm(nump))/LA.norm(perf) + fpe_vect = 100.*(LA.norm(vect)-LA.norm(nump))/LA.norm(vect) + fpe_basic_array[i-nmin] = fpe_basic + fpe_perf_array[i-nmin] = fpe_perf + fpe_vect_array[i-nmin] = fpe_vect + + fig,ax = plt.subplots(figsize=(15,10)) + plt.scatter(np.arange(nmin, nmax+1), fpe_basic_array, label='f(basic-numpy)', marker='o', color='black') + plt.scatter(np.arange(nmin, nmax+1), fpe_perf_array, label='f(perf-numpy)', marker='.', color='red') + plt.scatter(np.arange(nmin, nmax+1), fpe_vect_array, label='f(vect-numpy)', marker='o', color='grey') + plt.legend() + plt.xlabel('Matrix size, n') + plt.ylabel('FPE, %') + plt.title('Sensitivity analysis') + plt.savefig('sensitivity_analysis.png') + diff --git a/sensitivity_analysis.png b/sensitivity_analysis.png new file mode 100755 index 0000000000000000000000000000000000000000..e220b4a52305d051db1a2728d144fa3f893e7d91 GIT binary patch literal 21764 zcmeHvbzD{Jw)VurK$LBON-3p+goHF$;1U7plI}*j!2*@ir6ApiG)PJ-2qN7j4I(8V z-QRe@ea^Z2o^$Va@Av(A{QdUYi?!yO^PTTI-Z7r>jAzVecO^s*A3S{!!?44+TQ{UJ zY!5ev?P}P+4?eMJFI9&B?XkLslid&hocHTIg})D2+)}c_utR^K{|K6TT-V`4ZtI(h z)-vXL)^?hfx|pe^wS|efwTYn?wXLqDm7%#AJM(2`HYREVYikQ$7M8z%&ung~&l2;C zw+q9lG2D%-vi6~qJ#Nl23q_T)b_WT%JWf^|`Mm4W{fF6EtmhBP9AqZ8*rl-Ns&`Fv zhg1$@zIa+wZGy_Ej$`Wwt2KArmYV8QQU78rn2Imds9`Oyx`6D_{igKSewa=8=RWf96Kzk#zvrlb zdmgZVWk&{_Z>s6aSbbOQ<>dn|GBQfqhc-1yiwARDcQt${0czSrd)@@8)-urqhrTI*|X#D%Hkm^3~Y9gbY*CJ!o?QG;v+oN z*EZLtlg#^y5AJ*9INvC0)9Wyb!~gpAixx*C;{Bp}q$yD@iNkr>!r$NjiTW9$eXY6X z>T(JSSDaWDfA%CL$iyeAm$;7i6y&~@2uqym&h?~<;JzR2y4oBu{>`dhu0C7<-PFbq z)z*F?*Ri-a<1I3haz;jpaEbW=PpX{daH_PoFs(nHGp62+5OirtQeg2p#V%4{J9RK( z%dTXt=RltdhBXKxzo6LA5OuvZPK;F5VZ`5Yp;aNxN?NSy1>#|j zjn|iFgZb>t%B{XuS1^s{wws(T?~YbBz# zPW|p(>a_MGwEV6mLW33RcnQ8K4xckuB6b5(|en5JC#T$7Vs`VpI7T-rTCSeShXTY!1#>r0Y?=7vMG^uf^A z`ucXN%^A;k6w@AO%o?=)1lmm2!<~Lg!3wsesi(pAoX9QNFm0(jfX!ZnZRaPr`Yy-3 zFZGDvvM`_RWP?q+Dg!PtEv%w@enEF~*`|*QMS{fRTx^5)Iho9c$h3k?t zv#93FM%}%oAK{Yo{{8!T*vJvGE}43@O)%&<3O?H(K395Qlnpd2&km_9bTv!xiI0Xm zrG~pKD#7lNOe0>RStMYHs$&>DzlN*wC*N0FcNxZyRh%)GWhFDQ^)@KV1^31SW}={` z7GAuVop#gPv?pIq=Ix!vAB{1buIu)u-|{kAE9dO_VUG1PXJbl64nDcII+Rg(1nGw-fD?6@`ACS5z4~b1*hfOD7H~PA+yvJq*@FR2&bL$mhi zh~T@MUMCH!11{Bu+eN!2IW3IG;a<^?!E7~s^=IU?`}vJ)W6c$|%=m02ReQ6fP|E$G ziVK%6rFu1%MyyRB5KTeqv) z88$}Wl;Zgb`_pDTA>L@H!Y3KG&;3|O-u;-L#cMw;*lQ&PS8ho1&pxj3rC`J<=UGJk z%p3BxD66Tc5E0mzZ%*M|gdtN1noai>jaYDl^}W%n=Now~qeE74-_K@c!6?nxg>d;VO_vj z3zMBP<|P|Bbv&cwbPNo}lN}kS1)cL?BFu12-ly2-)wdSZkzu_CbA*#kv?{2Zwr0Sc zsD@e9@mPMi7cA(K53W}5@TtqOopY(ht43%siX0c3K9j0{8_WRTl3E*bW#YArL?bvO zUxzxhixt!y;vB&&;=>(p<6j2`z5`z{KIA9VM)&5PWlg)9U5^bi!WM`+J`{ZNV6cX* zHn&Ajd%lg2XS{StTtCc4o|~FCK)ZDwcCa|?GwIv6pI>HY=hXbPNA>>a6SZo6UC0q+ zh0(MGuY{y$gC|3gnaX+Q!QymZnr4~Dlj}aOq?D8r3V!B`jf{+35-hfzQr>{TF8F;2 z!`AXP7Q3?zo8La|5vw`r#R(8r&mz(XoF|9Qmlrb}qu}O(gM)1;*ufjLQ+ZxH3!!qX zhP~dzE(ODk*XElg7aFb~Nr;Rros|Bhl|heT7JCQ?eT^aPpT1%ix(OFD>MwoZ_#?^_ z1#C}B$6FA3?*fKV1A~#LGKXlEKi?!KCnXj2U^n4-cDtI)`QuICOA{t1O?373jYsRk z@<#$xcxMtZEW>Q&*T`uOgBJ+pz}*C?`o?x7AHcBnGh(va_{n{Q?a^Q1cVVg?se1>C z?YcTTfvPBj!2+XGQzf>RD}`#9ccd9G zm1NI)Cc&fI35YQHO0BtT3_v-DBDJEKPlxhhJ?U^OD=UuQgY>&x)~B6E`MTiF3;dcp z5njlQ*r=Jj-0`GLvW?%werLaZkNFL4!u$8q$LQ}gKoI4${&Ba*x>@|< zrAw`y`VlE$JrD#q*H@n$BIPOwhan`v=5PksCZI6Up7sWURQa=GW%bwBckT##^*SKn zB7~a*2M!=xFf=qQ@Xa;pyc}lRWvHB~=d-#rod&q&CchfoQMPgW_{baewogwE+4PqZ zEKK*QyR8pVE%Z7~8@49OwaSjd+%Eitr2hxUq%o-URFLY@8!R8cq91XL}()hTm$8Bq)vM7A- zp_4xvMXCI0O}nyi5pJ6wEG;c*aUVW>7>L|A>$zibckaS=d-!lyYkNE6@q-D)#cF%^ z?wum$@$LFx6~DSLLCVJEv}i&{OZ&r{z3xT#R9l75nG;0&Y-%~%&-3xAyis#Z;4upej1A_HG?mxspPtVo8I-gYVgqxhW2Ue07;*qS3OiW_pMGAiVyZF~& zR9AZKufyU6Uoj7tPgW$E?YeuuWzrtdMbnp;G@O>hw|XqAFX3t9n@_o_y?!apk+^khzgxtM7XW}VQW*fj^=&TZE)06HKs*g^g2ktF zT`zsG8!#OI&$WFM=lauMx~wf4wk0d^&egE{L;T6@b*!GmqwUh&m1C+BQqhgz7o;@L z$edQk;%$~^G!ahm&@|{s*95rN;!d#p^@|sgm6T3*V7z!*Si2;*t1JzdYhr3D@pP{x)2UOZR;$%F&*2e6@)0wI#Jn)6gwnDGT4bLZ=jWF((tQ6s;zo=(e;p9#Y2P1;-fadVsl{rd=GW79aOY@u)FX z0e)^tkl_RXqUSf&k#YSPgG38#2LBfp!%M4>siCz_eg6FUNx$|c{YrtC*{qMs=!wei zmbh-f*qmL4&h#}r+F90>ZLQSP zgeyo(N5gq>O<kvD4vfn9kE;fsP;lt$IcSl35d})y})&}pv zNngGC0M0BcC)a8sb^(5bVOl2$JjFvD#}iOAMN#@T{^G@p)@@4qLHBk-UgOnLO&VD9 zl0vmK@Gr99#UNJ-R7GT1bZl0ZEEwADTet3th|MjDLQ!YAnk?3!;A#f0Xrc&Juf=?+z1C}!nLM)I=Ne44UMEBL8|iq zNK^tr4Q%iN%p(}tRWJ!BPjWtaRohP8u0n^*PTkTyHT`e_jLED!$26ERA;fq6i0iGU z8#p_&ZaqCE)~ek{_#+k6UF0CLHeM-N5ibIVGbQ*Pf)WCcuPhasEQd&$Uu%+W!@t5g37hF&RoLYIv?HUV#si<8 zoECB|^dP2yeBTl$`pN5zS>9Fc3U4V1313R*nX5QyUER1kp)F@s*Oig4b#=bnivmmE zj@|DSrE-x41o-UjouFv;I_~S>vb2}k-c+C1T!lPD%GK31FC4}f=C-+FSF%rRL=+PQ1bozqIzE^^0A=;lKvDK3ulW(X9>`%W9Tr1`431#BJ;8QDOY} z_;?lo#iI3J)hIVWhy-D`goFggr5-y(I~V~JTU{RVqr*)X&iSI4259rCwz}KKLT#NY zSPQ`Kgq$4tB>9vk$doNVIRZvVf&_B{lGDa4!)EJ^h4zWEBSNwO`yth#!#xldjuG)` z0!}Laz%lwCIo)Q>z~ZE~w6x4vK?r(jaUKue5&ZQ;VnTvpy~=?@TF($ct#My1O#!l~ zH=1PQ;i1WMy9iv`r+mi*4=m$zd+`?;uU&gs$2~;! zp~noeL#eH`lC7kLiS`LtfS0tKTwI8&7;j4nse}Aao<(%lp8vMExIb-}O`8(Rl_Re{ z@?8st*jxo$9O7kssO9tbnVEDrB_$=enP%{O{xDb`tM9k)sdiA3$`35+EAL+t_3aZZ2yhHVFSZG!)Vl zzzPlshiJuDpcqFr-4rKg|Ln_~`ju8k-(vc+XU`5U0RE}Gr)4>B&#qtBl9X1Aj2lID zm(}@a+>%;aT3%_Q>*s!3RCM&Uj~c96Ok_lf(b1n`pEnr!`!{`j)D(4HbW=s?RQ7d! zJv}2xkx8^Ag*LK?sRWZ#RC2$v5*|g>WI+Asp)DeWVjc+h=T(NYX^U_lt*WgR z_-Q4{GjiMoI6GeF<@9ja2wduwkeQWbipi(mOPH(Ulf<{gh#nCS=U0022tlZVq3v%lD&OQo3JS#Hi%b^ointCUBjPp}GJ$jlKX8 z38`43n268mILLooH&(}5WR=38V`4DL$fs zi`Fpj%D8al$`$?|B3SH<2|$geY&Eyw*TBW-U%R=u0MY9LJA0-n%XlOuIX=sNR!cmD zOPutc#e-jSbDS^)guW6-e;6>onFx&fqH5O@2YLI(`l!$mKLm_!<8Y0T#P+$(OAF{e|nCNK)t$q1Mr zk|C7j7COw`Nj#{hwU5yESr!C3wQGz7XNv?kO2wXrDr8{FjvwYTxzmS4yOSl!I`i29#$Hj75rX5h$4X(%naHOxaV6 z^nZK=DNJ+8<^l_3A`Ck==+6BMzDMMnl%Zil0JCZ$urfENI^+*<=i5xs`&Gt(9|}@E z$V;+QjVXTy2?1E3c!Z#;U9ZD7;|3Q=%*c>}fQGD%7MvXJ(Z0jdtYj^~nRDuXJ&x>u zsM4<;U;!NtK>{ERIJC<5Blv-63A6k$dIaTYM15NhU-{NBxekyPB>@YNoynNxs?_t( zd2eSc0Fn<2O%n<&KcIA(m%>Ry^bBB9ibl-Pl`9IspM=v(h9@Jy47=t~i69G$B)eWs zgOYw&9E^K>W}pnE(#Ums`V}~@SZ>UC3caR9rpfH1i6(snGXS3uDyLkS5!+YIo z>xxtQKc;f~+M;iIWdZP4b(+c!X4QNMyV8LI=sKWg(8ibt#tfL$)j-b4Zy#hP+ESWf zwm4l@t-(|$9@`v&J!4+95Rh83=?Kwt0=S_}6ji_)FWt+#i=T%q(xA-a@NG%ScHnpC zAyFH9U(GB7u@n~dVtjLgD5b*<;4C9)4IyVlQAl#N=<&i%9_)Sq8Kk`YTMgFpIDpc! zN=iYr*^I;7MnE3``K!LX2Cmoub|<2C7l3L;neMN(4ad3KD>lH~dqga)j7Esz<v+Xr`HkdKn&TzkSP$oKxRy9BP9}s4I3N zDAKuhiiE6JZ{EDg{f+zLwo5&Ckw8>SD;nHbS)U|W?3F&}OoaJ2Mp)b$bBh-)?cGHH zyL2Lt{jzWTyLY!(O4i@gKJjaW(#p&b5Va_=X`lH;h68p&xcvS<`1k>lC}K~VK~B*U z^>mE~u0R`(1Fh2ngdt7%985(L+z*3|X*#9jHZpQQP)shErB|2&>{B-2g_$MhH)4ne z23Vhr3MFE1n32E|p6lLDoUt`fKH$#UDAix|*8NlFR{ z`RZc=oTFntq!f|$LRG|%R!Q4V2cz6UZP>z(MoP8VOWH{|us8wAfK#6F7Mg=1v zm>gMFyFl3#5tXpXHJHAQ&vOB=YXYx64sZA zh4C|e%OW9Z`rq>I=a0P&K}96s25-P#-7ka^>oBLEAE0n#XKw;$U9Gb#ZvH-3nol0Z zvX!-49Ul^s1SxA~HQK05Lot4W(t3J&N#Sd);OvL3c_BT?ddDR7G8Uy}P~6Jze|RX` zl0Q>cULF-wq=6;y@{K-QSc`%TyxkjRd+ z#|%PpYWjsuH1~Nf{FPk0LGo%ihbPR zEx`Y*vaTJXh|Z9HA)=n4j1`bjvM@DuZg&bm?9A$25(j_JlLr{;>C>lyZ32c4UZMCw z9Z|1uYgmf=O%a{BOmtl$o?WYE%o-xdtBbq3ev>UG=XrQOhz|)T|5gK>Ik;z+>1jeO zik+C#0sM+QE15Uhp?!}4(9_H9$vX1;o4cBtQB6%v!%J%?4nBcbFwP3*i|I@`&$e$n zAX=l_?*6snY-C^n)lPib z0K7_Hzy1i0?u;WapxJRbIl5#XsPY#Gm_U;ZM%S;iqe@6m86F$KoY~zr!}65`Or(d0$5Y*lA2sfp<#v}RVxdGk;h)=3-n>akQSn_?7RSnl z>yv9|;WaKUPKiE}V8(n0kytblmBFB%lVvKulK4bKEy%jfs+Qtna18E>SGP44{jrU) z$Apu2vKWh|Z)UjU1B>SPtBPR(LdQhjz&2@( z$zS_ZWReEkDRNFclv@hu3XTGVBOxvKIQ4IR37UxiXAMq!>*Ovpp;e=ERjfo%GXl&1 zGG&{Tp58F|9!!c*8*Q$qXFvjR`0!yU5Lt4Q!zr2g}k z+6LHMr3m^_rykV<_$Rruzp#HPD1qDi*Kquw3X^_M3WhOM{?loZ+1a}o)0fe?gs7;i zfDnj1ufr=$5jY{eH~Rt~qFaac`NB#BCBCO@@y!>?Snnk6XC>Y7hl=3BEJG?D*#A~j zodi~`i-WkFZ^u)mBe$)=O2&4QXy1r!f%rfE@Si)q~(u`7OxOx8gxRvDN5t_{>~6dUKi<8Mr~$Yb+U;sTTmZt!om7@rU1~wYxpJcKWV* z8gTR7--f$4Rv#ExK89G?c{9=qXxs~c_)+-|m`~r8!f)U1$;rvtXz8NW3V|!a2^5~h zn&F9W|GL!`6rKt-4H@GJBHNNp#oyGV9>_NI>OvEMM#=qb5fdAG1N_t1+S*pg_{HBH zMT2>7+RY{2@-k#C6(5SHcV#mpv@LO^4&Kjgf@)cy`S|5nls}9YQ+Ag48*>;3X|Rpu zlkIT;Oh}Rz)nndMvMjUmIgg?P-30j3XOJxeW1t)q6tp^U#Ld8#2M%|qP53+Sa_k%27sw?#oht0{dmtompBksHT$?6H}Hiz~;Fzjyl%x7jA)<15aBQ1}`i)&cAd z?9qvnCxcZ_9?g`Nkx?lOk`x!`j0&#Xxx}_NmKvN}K~MtZTb^ZS-&DtWc=Qc>T(w9pY^0sJv|u!7z?|BL%ztKKT7HxR%)^x`=thL3sgcs_oGesZjL^Yx52!6XV>e~ko9fGYuU(%O$_dG4ULcmv*A0kiz9&}7QXip1Dm`6iYi4Oaq-~uhGDUUcs$hQ&!x2+TiRe4 zp9bRmo?a|tRlb*~;h7>ml?-_g9Gm!J-U@dD4}3Eei*!J1iWO{=X-_Zmnfn{*15w<3 z57%Qo?8bT>Pk%e#s+)_w`u%Y}mZ@Ft2LE^M@#rle&^ zibdM&47U@8%g4yf+_oA^>pDkIF2>~mLC6C!T@cT0w%@dJ6j>LxwWk= z5l92TdH9x&j(6|hU&fsqF2!kXldxE1Gq4Sh%Y0%UevE`t0yuT~-%k98HRN>whb12G zc>jhm*el%}#G(8DJTp+)e})6UN^R+z>Ab<-yD;~_p}7?uG7tQ#S10!G8MZBuQ&cpt zReEYQ`5m)Pg#s^dRH>aMOY22tRpl~ zGxyQM+lPV*4o+;w`#!7mg5ecm{u=%v{$dFkp~OCz66+TR^CXc26q5*XL`9v|3cAfo_PTqDE92q>wW0tI zr5{5yp=UYVlPDOWJh*6n_4)pHud@U%;OIVo{(f>_HF8zfFBTU`xVvV4{xCo8$5oh= zlH%i-BFuDwi3y7R!4M@N=)fyvA9*HX{yzM5pU0ZR>$*)|=u1gxIz6C=}P+qkshs&6Yx zWVXv51AD^A6oUrGUvL1HjMUP+PrIkg?(Y*ZUfWNeltFQL8Y35TdL+lEeT5K97lhFN zgL-O;rhnF6d2--t86Q^J)T_~2x9+Dxj0Yq1!4$rf5B58 z!dr*2N@+#Kj-RTsllg}Z74Ts03*ZBVy`u*x++Q_@vz|NMt)kM}=%(QH{JHgLj_GNB z`}Et-y;eX~!&FWJDlt9M2OY2fmam&IfT_=H>(jq#@={}_XJAB`}}8}QI8NWX>e zn{J0Hh4HbAyD%$6zRve!?qoM%;9)+hkc!HS0E%aO?qP@Vo+(Jjuw-vaa*Q`r^-Nnf zDcrj?*YA{yBDCYlUQFWzAi0L9;&1770~R+plFMd+U53oVK;JR{r!{vLyMQVk9UVvE z?gq`Sf+9{3@Id0{POM8NlIdadJ50RHwl7lKbr=hK4lYAOZu^7=iYQ%)WCu+7<1R-s z!etGoU>8)g7D&Z)j^vkn^eO%m9pHtB5$(#z$ zH?{&-IJ94?DabG-OpLQ|jNVZmC#lCUz0Ql>pl zFzB{2=8IhSe>9v$EC|wYZ~{|?G`dI|SI0Y%tX=Hk`!S+Z{I z%w~UhDKH9Pimt9M3HvB89&g1%M7|JU^!?#qlso9#=m9)(7&JuAEFoQXw|R-F)=hWL zz@nRMHP-_P*dfSeevBH>R{;GuHk#sd^_)1SR4K=U=$xHI)8B?az3z5CAD)&4k3n+d!#J)co-v@J1#M3%kI{ znJqp9tJ??oPzzKYs054ZN=R`5BEneEzM0u73D;OKL_!5U`*@hJKukPyB-T+ zh>C#F z7_0z6Ijrpu$)y{%!k}<_o?BA0XPYSKI@kO2*&s!C`6sUZ?C2nB*T(;(n)`dH2He-t zB?Hz&aA3I*`245FqPhER|4GgZDb$Tp0_xOL9%n&b!}i^p7nmj7iETHsA7JF5@;2~9 ze=7pOcF#Mxo!!BorI37J<(z-b?%*3?QQ&=!w*O@tBZO%yg-$d`)Lp*osng?A=-Bn` zv>e|3LuSh)MPO3;AM@2r6XflL?*XxPbo5`besBbKloMR#V5lY;$krvIrXAu)LasJ?+xI{5}_tZ(m^AJLnJZ%ZZib8EGLcAgK_C_%dN66!o!wmj+S z!qg4us4`YydTA;Gi(oA7_BTsOXaK#RD~9rUdp442j!C*TRkM_sNXtDXgbLE{GV9kV z+Ynh2-a|VhVAw;>77enoty!wADA1VUq@i2_fCg2Z!LJ0}qexCpPe)odD0M(EgI7S; zF-IO9W_3c|-}c(8MzYyg5$iYH*1G$wt1so=!kd+B*aA-n3iJyM43)#<8a9tlEduyz z1ub~SohO}$411~ow|=&&stTnq3n2Va2XUF>;4=meP~icmsv%X=3yb3(#`=gxWE~n) z8lm$5(yVt-g!vKoMr~VA=*66W3(90jsMSfgE1miivdxaY?A%&692`FMLKi~O%J&fE z46P@$hM*R9oXT#W2R%UWJqkDuay-c9*|eY(gie4&+Mr1~%ijw^u1ji#9pxx%bnp27 z9|2_0=lvZoepg^cTiV({q=~8?90kX5(q?9ZqpOfZ0+pSwC5OlOZDG=xeSn9F#}ToW6ayMLn`6^MXIA@Y3s!=EezG{-aTRBr3{%eA4Z;vDaJ z2Gg+d(7f^Z%pa-?ur7nA2t>V?XZ_sBHyt6{C-K`DBrW+dX^EBo`jL19ESm-FYOtE| zNpZf>Q@I;hgRONdm`gh8E?=?%3vT80{J3cXj6p=f1*XGElhmI~lgtB2paCgtD)U}; z8cGs?)E^yC>bT5z2s1chEA@;7BwhkPhoIK}Blh;DkMBoJL;agUbY{jJiMwL?Hv35< zatDNE(6tn8>&OMxCNuMYM0;gHjWknim%HuCFY)(f15fjCA|~2is(}x}{4i|l zGm$S@hEeN9==os;g&*p)fZC%qXdnY2l!tP-o-tPaX7GvVfpBDJ$th9D$_x&JEj)#q z#(6j%-x4oL1+u}8;pbSSOCBzVIj7^vsg9NXSc4jD5E78)_;m1$RQHJ?gV&Ka#4MyI zH@g)Xm#FR+^jzX{8WuT5@oZ-#%oo`ns>)*JXZWU{9zD`t=R=QSP{cV5+o2qaY8)U; zf_jnryNo$=?V+fxrxbz=;GfjN$&`ZxJsI|)B%Z_S`yF)S^a^{hDgE{komIF=%vvzm zyW%NJ+VS^I=F1pX3uHRrzBv+$^8KD=jf3p~oi(rQhuu^>)*qL_iO5axYdVXy!PT_76 zU`AngljvD=)PV9R`4CHC!{lZ3888a4u+I)L=k;%O7c6z;PptY4rR|>ZFy4sTe_g(O z1bgsAjw_*#@uBwluN1>M=XH+8Hd3G7^^W9)&U4-qH>73u(h(5Qjqj&E_b5v1JCpRO zkJaAl;Vuy>r?19>vYgipsr+=RR(awXrnw?Qo9AkloQfSl=8~(Zo@F5A z%RTsXU`rC3`aptW|FLpxYz!wUDT&va>dKB`XM$mRrz>Z6n#c?6iJ4AJO8p8Y;1nj+ zd=8KmxZwJlJ3E7VdwZqAiTa<|*xAK58pLL2b3S_XD4~5qLTMMK#@2hx|MP{5$4XOj zdxKeA#>76%R-93{{D!?#wzIQaU0Z8zYYR|L3CPLGnO|CpFDNJ&Z^SUNM|*q^^aCYh z4|0e%(7U7$bv+*F)Jvp1&|hN2aT?-Fr%6g&9G#CyO!PDOLACR~$~`l0O4j4XNw6Gc zW8(?=l-(F>2G7IY?req@78Xg6TLgCOA;v0V9gD)FqM~Zspok=$*3KBb9}Dl}Y`^#F z6uX{Ih-pq0NG&X3L&^&Lap%@9LQGp)K|!IbyPKVr^@Uz$(1QmLvT}0^o*%@tzoK7( zazdbhe_5l>lq<9P(UI>^izelX+{JvX`j91#7Q^iCLy7ha^AW+oI2=7g=L;P*ylC%o zMY38^ep9FYN+H8VFB+g^uuV}nXhF0R-lVufPiKg;EVY^_Wfbd&ZwWphrSCa)14Utm&9j~s8*_(TGiCL_%+C}=j|XzJ<;{`m1@%Y==BnpR87 zy$_cGZ0mG&b>V4EO-&u3Gz5)yU#+dJW%()YFi8LyinQ3#Q3fqnuk~5ISxDfMAvQn`gJYjJHbf(YSTU$piUX*8Vx<}s;*PK1$;Fp5$uICZWO#b#MA z!8Tuak%;d+==_BwT(?&QIWvz z-MbCBLQav9-Ewj&VAPyZM^{eb45`&?z0S&IV{7|1H#c{j4Q}Y-fPj<|8u)*@;s1V~ zdSeq|+HzWYdcs6vUbeNh5fKx=RZELcN=h2faKHc6ikOv>`{8b^BJf|%3Bk~dQF4vK zquiK=x8Yy6z58qpoI(wZgS8?nC&#Rnx3qeIhmec~Zj}shIU&qb|Hw?}Zf{bUId8TF zM}h0tuZN^SaVQSReGU@VaECgAibJrf6@jo$u&wnqG^lI!*wOr$o7>&ots!d+3r@IVDlT=aGLPH%Ib6G@1W$2DNu7724m2O-1Q6nUN?8 z3*b&dXL0A~5rO%;uApbp&+V~Iay$;}LUw;w<(}CpN!QEjo4?-3B()||?uCylC@nQe_-z=DAcPD8yr?_gl7PcQ8wEFLlIo#%<`<~BH*c2YW@fS%6&0nG z2f~H4eV1z(-TpZzcP>Ob{q|K&_-Kg@&v#rM$Rp&Q|{SOvWUG7HFcvPL5d zioR&f$p1%XOD<^&qKBJ$H@t=lb61N=x+|> zb$XxXM2$gU)a*rPG8fiSl9N0z9)o z2xMMr#N7Svv{fK54)1FwF*$k1+6IV+fpJX!KYJphuc6T(D<7?ore<*f4{4XOn{yFhML z!?yDAi0cdJ$>M}gm9;rR34rX5j*fU*NO&n)gJ73w#I^+2>Sd0fB~wKZC>;83T8y$e zLTXJ)HWtg`lM4gu>gvE77wn zB#P~(rnUY38siI&`T6;YlwJqEZpTxT?0*Rk*sgAaG?Q#V7XbkQwQb7E%7Ah&0s@ri zhcjy=WpASbZ!Uz}+*7!?x%1>ESi80{U5;{*T75%9gPt*90?1ZQh?n2khH~kxaQ&v< z-f#f=FkeJOM6pFhJ#w<|ysOJ$s%C^p*vZ8Hmx$f zc0x0P6+Apd9|)@xl3(pnuX#-nI|DQm8T{eo*LVeB2;?w2kWr3-lON0uczM=h1`h(* z3sh?zmA{}hnBxfljBqNHQBbe4oZn||o39t*huc)ldrYRF9kULa6W$xOCPE@2siR$v zo`|EqF|IJ@DYSkTy1VK@+ZRqqX$na=LpeNk6P?c%|YK!c`n0E>D!jPXhY z(3F{QnnHMfMC}QMdr}MUyk(#Z))2HbJnPv-m5c)rL+_}o_szs<7PVdeV-K;~dr*q& zz|%f%f$>74uhpjm3_4YKIyA=Yt<5!kl=YaH!Y4c{!>?Y^qAn=cU#}u$*-j62Fc(bX zp);)(hgX17RS=Jr9;Dw>0wQuL@?wy`;dKlQ&ZE81YJ~Q_9l z4~4E|99}Qea9Zj>TDJF|2anz)my`&Z+n1g3DCvbZlIFfX^?@6AI3RipWmYp6(5@_< zX$0ga%TF@Y2~T9Of?fkXSfEA97&GYD3=(j%g#uIn%I-lWS#PX17kme%@>%;mk$UJp z@MnQ1OLXdmxRqS8hvkHUs)92G$$uyiZY9Lee3O>=v=TwGinppn+HMM`w!h}&}6$)(xg>RSaMrM?VnKM9Uh zDaSWmFe#c^}(F$kGXF=wM_@%M-r=aGJ$A`rt6T zULHUGU$>xN@C`pai6#it&#g3U#E@!?pwWQT(-JrrhKwcqk;3nsX;lP-S>&;d ztbJAQuYg<~=rUO8?%NCiNHUy1ky;I6a+p~-y?o+U(w+k{%?$YQs7J!~2)*nZ!ED!ZlC$HrT{{F&uD zoVgCIWsgn!NYS*5dFqGnr(tlMknGw_v>o2s0;ZxG=B=>p@;Rg0Q9k0sYNeNQzu}>e zfiZ;3^eeQl&vn+j=8V52iv%b!Ru+x@%;m$w3tIRY544h-rS1k>;EV3BHQ zDTsi`a7=>xEL2u+fERfMQZZgYoC%KiDiWZjKzDZV^9L}1N!SE2jbEfhDg0?W^X_X) z@Q%egFJ?y6f^7C@cEB5P6iu+_qsLLa`ftgA`Xu?J5omdtnQ2FH5*~Z>`Z{E7H%rox z1NdfHeWLt^kcU6L|AkNL)dFTv_aSl8SVHB3fUdfLD^Jal0cjh56&+vR$`$rm(56KN%a(4APf(5@Z<>?CV^gsUU5;Bvq7 zFb!1|G~;;q?c#yIcQxa@y;WR>QyRWu7zXb=3n=c-kVW(YXiz{_~9YuFCl&F-@z#l z!)y#yZ|&WMrIA6*w!g{{YHIIpg*T7u7D!2N`n!ld+`wr^{5>59`jJnDuUIv8AjLbpdKkYX!8?Mi_0 z>s+sY1&t4*zp7az)z#HeYsz=IG{;uGT(?oaJ{)`ml{9>N(N7y5X#cka@tiK~IM{_Z z(S_;-mWD}7mOr~~IDo6KT!(Xt0b>YlD(fYiU6}j9JBf|}MkXDbLwa9BmOI*IoK{OH zIP;*^3^nCj0k5VHrvtxx>hx*HkCi;@GiAB^kKr!cjJG&Lk*p3LDHefxP&?rnKDl}I ze0?i*&f_u>Uw=CleR%jz&!8ZRik^b=QO@w>pc?2EskjHu=mRi^mtkpWt4d3Sd%XN4`F;n6qYD@X;BUJEYleu13*u^7{sMa-*|0mM7PD2$QY8 za1kr@f>vNn=n-27<1oqD1a9?++m9arnTrM%P!t1;>k)IDB&2d?MIji4nD1cRhORIP zCif!%EvJqB`0)j-AO!I4CBZX8u*UTu)v<(;9!>)fHeC&ke2CPg#ksk;LE`p`HVBwS z2#xAZ&X_RI_N+{SQ(Y-%5X|_z3UL1wJn^jA0W2d3S|ksrQQ7xDaDWN&5Ym}{0*Ddy zxFhZc%2VNtsd?uBiYdcmu}V5Om_H2TYoLu1!i%}uLZZURi|fRllZ#+_`tUfkGAfQ$ zklr)G%~&80R0vtGq^PKM7m-JOAOiiSo~XPFnu|Ns<0