From fd454916fe058251e34a8ac01e290ad4b039393a Mon Sep 17 00:00:00 2001 From: zaplotnik Date: Mon, 6 Apr 2020 02:36:26 +0200 Subject: [PATCH] Start full network --- full_network/generate_connections2.f90 | 518 ++++++++++++++++++ full_network/generate_connections2.so | Bin 0 -> 167760 bytes full_network/korona_full_network.py | 702 +++++++++++++++++++++++++ 3 files changed, 1220 insertions(+) create mode 100644 full_network/generate_connections2.f90 create mode 100644 full_network/generate_connections2.so create mode 100644 full_network/korona_full_network.py diff --git a/full_network/generate_connections2.f90 b/full_network/generate_connections2.f90 new file mode 100644 index 0000000..de88ec0 --- /dev/null +++ b/full_network/generate_connections2.f90 @@ -0,0 +1,518 @@ +! File generate_connections.f90 +SUBROUTINE others(conmat,n,maxcother,randsinput,randsinputsorted,kk,th,connectionothermax) + +IMPLICIT NONE +!Pythonic array indices, from 0 to n-1. + +integer (kind=4), dimension(0:n-1), intent(in) :: randsinput, randsinputsorted +integer (kind=4), dimension(0:n-1, 0:maxcother-1), intent(inout) :: conmat +integer (kind=4), dimension(0:n-1), intent(out) :: connectionothermax +integer (kind=4), dimension(0:maxcother) :: tpind +integer (kind=4), dimension(0:maxcother-1) :: cons +integer, intent(in) :: n,maxcother +integer :: i,j,pn,iad,pnj,diff,indmax,indmin,ind_cind,suma +real (kind=4) :: pnj_real +real (kind=8) :: rnd +real, intent(in) :: kk,th +logical :: first +REAL :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0, & + vsmall = TINY(1.0), vlarge = HUGE(1.0) + +INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) +real (kind=4) :: random_gamma + +!f2py intent(in) :: n,maxcother,randsinput,randsinputsorted,k,th +!f2py intent(in,out) :: conmat +!f2py intent(out) :: connectionothermax +!f2py intent(hide) :: n,maxcother + +!write(6,*) kk,th + +cons(:) = 0 + +first = .true. + +! compute number of nodes with j contacts +do i = 0,n-1 + j = randsinput(i) + cons(j) = cons(j) + 1 +enddo +!write(6,*) cons +! compute max,min indices of nodes with j contacts in randsinputsorted +tpind(:) = 0 +suma = 0 +do j = 0,maxcother-1 + suma = suma + cons(j) + tpind(j+1) = suma +enddo + +!OPEN(UNIT=12, FILE="numbers.txt", ACTION="write", STATUS="replace") +do i = 0,n-1 + pn = randsinput(i) ! number of connections of node i + + ! randomly draw the connected nodes based on discrete probability distribution + do j = 0,pn-1 + ! how many connections has the neigbouring node? + pnj_real = (random_gamma(kk, first)-1.)*th + first = .false. + pnj = NINT(pnj_real) + + !WRITE(12,FMT="(F7.2)") pnj_real + + + ! dont allow pnj larger than maxcother + if ( pnj > maxcother - 1 ) then + pnj = maxcother - 1 + endif + + ! the indices that have "pnj" connections + indmax = tpind(pnj+1) + indmin = tpind(pnj) + + diff = indmax - indmin + if ( diff /= 0 .and. indmax < n ) then + ! uniform sampling over cind, which includes indices "pnj" connections + call random_number(rnd) + + ind_cind = FLOOR(diff*rnd) + iad = randsinputsorted(indmin+ind_cind) + + ! wire the connections + if ( connectionothermax(i) < maxcother .and. connectionothermax(iad) < maxcother ) then + conmat(i,connectionothermax(i)) = iad + conmat(iad,connectionothermax(iad)) = i + connectionothermax(i) = connectionothermax(i) + 1 + connectionothermax(iad) = connectionothermax(iad) + 1 + endif + endif + end do +end do +! CLOSE(UNIT=12) + +END SUBROUTINE + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +subroutine household(conmat,n,maxfother,connectionfamilymax) +implicit none + +integer (kind=4), dimension(0:n-1, 0:maxfother-1), intent(inout) :: conmat +integer (kind=4), dimension(0:n-1), intent(out) :: connectionfamilymax +integer, intent(in) :: n,maxfother +integer :: h1,h2,h3,h4,h5,h6,h7,h8,ec_size,ec_centers,pp_center,group_size,gp_center +integer :: i,j,k,l,ps,en,a,b + +!f2py intent(in) :: n,maxfother +!f2py intent(in,out) :: conmat +!f2py intent(out) :: connectionfamilymax +!f2py intent(hide) :: n,maxfother + + +connectionfamilymax(:) = 0 + + +! households in Slovenia: https://www.stat.si/StatWeb/News/Index/7725 +h1 = 269898 ! 1 person +h2 = 209573 ! 2 person +h3 = 152959 ! 3 person +h4 = 122195 ! 4 person +h5 = 43327 ! 5 person +h6 = 17398 ! 6 person +h7 = 6073 ! 7 person +h8 = 3195 ! 8 person + +! elderly care +ec_size = 20000 ! persons in elderly care centers +ec_centers = 100 ! elderly care centers (100 for simplicity, actually 102) +pp_center = 200 ! people per center +group_size = 25 ! number of poeple in one group +gp_center = 8 ! #groups per center + +! generate h1 +i = h1 +write(6,*) i + +! generate h2 +ps = 2 +en = i + ps*h2 + +do while (i < en) + do j = 0,ps-1 + l = 0 + do k = 0,ps-1 + if (j /= k) then + conmat(i+j,l) = i+k + connectionfamilymax(i+j) = connectionfamilymax(i+j) + 1 + l = l + 1 + endif + enddo + enddo + i = i + ps +enddo +write(6,*) i + +! generate h3 +ps = 3 +en = i + ps*h3 + +do while (i < en) + do j = 0,ps-1 + l = 0 + do k = 0,ps-1 + if (j /= k) then + conmat(i+j,l) = i+k + connectionfamilymax(i+j) = connectionfamilymax(i+j) + 1 + l = l + 1 + endif + enddo + enddo + i = i + ps +enddo +write(6,*) i + +! generate h4 +ps = 4 +en = i + ps*h4 + +do while (i < en) + do j = 0,ps-1 + l = 0 + do k = 0,ps-1 + if (j /= k) then + conmat(i+j,l) = i+k + connectionfamilymax(i+j) = connectionfamilymax(i+j) + 1 + l = l + 1 + endif + enddo + enddo + i = i + ps +enddo +write(6,*) i + +! generate h5 +ps = 5 +en = i + ps*h5 + +do while (i < en) + do j = 0,ps-1 + l = 0 + do k = 0,ps-1 + if (j /= k) then + conmat(i+j,l) = i+k + connectionfamilymax(i+j) = connectionfamilymax(i+j) + 1 + l = l + 1 + endif + enddo + enddo + i = i + ps +enddo +write(6,*) i + +! generate h6 +ps = 6 +en = i + ps*h6 + +do while (i < en) + do j = 0,ps-1 + l = 0 + do k = 0,ps-1 + if (j /= k) then + conmat(i+j,l) = i+k + connectionfamilymax(i+j) = connectionfamilymax(i+j) + 1 + l = l + 1 + endif + enddo + enddo + i = i + ps +enddo +write(6,*) i + +! generate h7 +ps = 7 +en = i + ps*h7 + +do while (i < en) + do j = 0,ps-1 + l = 0 + do k = 0,ps-1 + if (j /= k) then + conmat(i+j,l) = i+k + connectionfamilymax(i+j) = connectionfamilymax(i+j) + 1 + l = l + 1 + endif + enddo + enddo + i = i + ps +enddo +write(6,*) i + +! generate h8 +ps = 8 +en = i + ps*h8 + +do while (i < en) + do j = 0,ps-1 + l = 0 + do k = 0,ps-1 + if (j /= k) then + conmat(i+j,l) = i+k + connectionfamilymax(i+j) = connectionfamilymax(i+j) + 1 + l = l + 1 + endif + enddo + enddo + i = i + ps +enddo +write(6,*) i + + +! elderly centers - groups of 25 in a center of 200 +en = i + ec_size +do while (i < en) + do a = 0,ec_centers-1 + do b = 0,gp_center-1 + do j = 0,group_size-1 + l = 0 + do k = 0,group_size-1 + if (j /= k) then + conmat(i+j,l) = i+k + connectionfamilymax(i+j) = connectionfamilymax(i+j) + 1 + l = l +1 + endif + enddo + enddo + i = i + group_size + enddo + enddo +enddo +write(6,*) i +end subroutine + + +FUNCTION random_normal() RESULT(fn_val) + +! Adapted from the following Fortran 77 code +! ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM. +! THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, +! VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435. + +! The function random_normal() returns a normally distributed pseudo-random +! number with zero mean and unit variance. + +! The algorithm uses the ratio of uniforms method of A.J. Kinderman +! and J.F. Monahan augmented with quadratic bounding curves. + +REAL :: fn_val + +! Local variables +REAL :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472, & + r1 = 0.27597, r2 = 0.27846, u, v, x, y, q + +! Generate P = (u,v) uniform in rectangle enclosing acceptance region + +DO + CALL RANDOM_NUMBER(u) + CALL RANDOM_NUMBER(v) + v = 1.7156 * (v - half) + +! Evaluate the quadratic form + x = u - s + y = ABS(v) - t + q = x**2 + y*(a*y - b*x) + +! Accept P if inside inner ellipse + IF (q < r1) EXIT +! Reject P if outside outer ellipse + IF (q > r2) CYCLE +! Reject P if outside acceptance region + IF (v**2 < -4.0*LOG(u)*u**2) EXIT +END DO + +! Return ratio of P's coordinates as the normal deviate +fn_val = v/u +RETURN + +END FUNCTION random_normal + + +FUNCTION random_gamma(s, first) RESULT(fn_val) + +! Adapted from Fortran 77 code from the book: +! Dagpunar, J. 'Principles of random variate generation' +! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 + +! FUNCTION GENERATES A RANDOM GAMMA VARIATE. +! CALLS EITHER random_gamma1 (S > 1.0) +! OR random_exponential (S = 1.0) +! OR random_gamma2 (S < 1.0). + +! S = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL). + +REAL, INTENT(IN) :: s +LOGICAL, INTENT(IN) :: first +REAL :: fn_val +REAL :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0 + +IF (s <= zero) THEN + WRITE(*, *) 'SHAPE PARAMETER VALUE MUST BE POSITIVE' + STOP +END IF + +IF (s > one) THEN + fn_val = random_gamma1(s, first) +ELSE IF (s < one) THEN + fn_val = random_gamma2(s, first) +ELSE + fn_val = random_exponential() +END IF + +RETURN +END FUNCTION random_gamma + + + +FUNCTION random_gamma1(s, first) RESULT(fn_val) + +! Uses the algorithm in +! Marsaglia, G. and Tsang, W.W. (2000) `A simple method for generating +! gamma variables', Trans. om Math. Software (TOMS), vol.26(3), pp.363-372. + +! Generates a random gamma deviate for shape parameter s >= 1. + +REAL, INTENT(IN) :: s +LOGICAL, INTENT(IN) :: first +REAL :: fn_val +REAL :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0 +! Local variables +REAL, SAVE :: c, d +REAL :: u, v, x + +IF (first) THEN + d = s - one/3. + c = one/SQRT(9.0*d) +END IF + +! Start of main loop +DO + +! Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0. + + DO + x = random_normal() + v = (one + c*x)**3 + IF (v > zero) EXIT + END DO + +! Generate uniform variable U + + CALL RANDOM_NUMBER(u) + IF (u < one - 0.0331*x**4) THEN + fn_val = d*v + EXIT + ELSE IF (LOG(u) < half*x**2 + d*(one - v + LOG(v))) THEN + fn_val = d*v + EXIT + END IF +END DO + +RETURN +END FUNCTION random_gamma1 + + + + + +FUNCTION random_gamma2(s, first) RESULT(fn_val) + +! Adapted from Fortran 77 code from the book: +! Dagpunar, J. 'Principles of random variate generation' +! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 + +! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM +! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO +! GAMMA2**(S-1) * EXP(-GAMMA2), +! USING A SWITCHING METHOD. + +! S = SHAPE PARAMETER OF DISTRIBUTION +! (REAL < 1.0) + +REAL, INTENT(IN) :: s +LOGICAL, INTENT(IN) :: first +REAL :: fn_val + +! Local variables +REAL :: r, x, w +REAL, SAVE :: a, p, c, uf, vr, d +REAL :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0 + +IF (s <= zero .OR. s >= one) THEN + WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE' + STOP +END IF + +IF (first) THEN ! Initialization, if necessary + a = one - s + p = a/(a + s*EXP(-a)) + IF (s < vsmall) THEN + WRITE(*, *) 'SHAPE PARAMETER VALUE TOO SMALL' + STOP + END IF + c = one/s + uf = p*(vsmall/a)**s + vr = one - vsmall + d = a*LOG(a) +END IF + +DO + CALL RANDOM_NUMBER(r) + IF (r >= vr) THEN + CYCLE + ELSE IF (r > p) THEN + x = a - LOG((one - r)/(one - p)) + w = a*LOG(x)-d + ELSE IF (r > uf) THEN + x = a*(r/p)**c + w = x + ELSE + fn_val = zero + RETURN + END IF + + CALL RANDOM_NUMBER(r) + IF (one-r <= w .AND. r > zero) THEN + IF (r*(w + one) >= one) CYCLE + IF (-LOG(r) <= w) CYCLE + END IF + EXIT +END DO + +fn_val = x +RETURN + +END FUNCTION random_gamma2 + +FUNCTION random_exponential() RESULT(fn_val) + +! Adapted from Fortran 77 code from the book: +! Dagpunar, J. 'Principles of random variate generation' +! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 + +! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM +! A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL +! TO EXP(-random_exponential), USING INVERSION. + +REAL :: fn_val + +! Local variable +REAL :: r +REAL :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0 + +DO + CALL RANDOM_NUMBER(r) + IF (r > zero) EXIT +END DO + +fn_val = -LOG(r) +RETURN + +END FUNCTION random_exponential diff --git a/full_network/generate_connections2.so b/full_network/generate_connections2.so new file mode 100644 index 0000000000000000000000000000000000000000..125f98c259e46ddf244444e8e1cf5fd3883319f8 GIT binary patch literal 167760 zcmeFadwdi{);Hdh3ycucgMvmyZQVf=4KYE`M2%*`4D`SRA_)*Qn#*KD0wIaX1cIUl zCxf&@Lv+<;cfGK?t`EDqE4%V2culwkFRKAp#mhoe)E)y80U;nV@Ap*KOgb^q_1WkB z{r-6Ad}gZ8sZ*y;ojP^uRCRUcZij19l*uH7uV~?VfuL@yj>+N~Ez6?#ogich>B3Na z|3bK!(Z)t{et6UUn2rmlf*>&j)iDI0+`ljosl;x=X?lW=Pd`Ub!{tZ{Dp$l;ij>v%VWQ(ir5+eq+;#%J=hnL=Ljb6?)Nct`$+!kNE%e2d%s(ACS&rYc6@GX)=# zlXTPkiD0rKo!1n5R<0m4HC4q3X1VD!QL&n$dy!pC?O z0=N~QB91IU>c(fM@JQhYYb$U6>YW4exz;QH`q6qJbEEoA{7tXQn?AU$VDwwHmYx?s zdtk}3qtCp1nrGqao=^W8C|^47iy04;=RdPDscjeFiH>*oUA*ZR1=B}A{;$&~EA8=P zuDLW&ZNKlnGb)nLiNEb%hf&9$zy4s-oJo!M^XdmohC+X&gZ|eSP}~{)k$`=J;QtPz zABa8^iX905>LB=s(SZZeZ$`ZX$r&{WpFa&EC)tNIK)dc4giq=q?K(CHea;~AKR*bc ze-6^FwS&-G2B~)!+BA@UUKoUa`ylw$gW%tXoeb13qX*HaeGoo(f_@-AZw%7jZG+&S z9E5-RAo7$9Qtwm152Uy9LF)a>Aov;Jb1^>suh)PWsJ$}=(dR>h*sX04`pXBgw{r%e z*FZl|d-n{|t__38zjhEl=MEyzib42H9Hic(gU~-X2>ovcp>G{T57!SOfB7KvHx6QN zuM8s3F9+fC-5~ne3;Kcj<+DNREgA%W*C6t29;Dv*LGWJ=qR;9<=)W7pkF6iXo<9Je zf!h0*LE2R?2>yaW@ENS$I85kh2(D3jA{=A=5+~$(bi^G)_4F(xX9%YY`JA4U3l}hb z#t9Oi#KU8GFXwNWtW(nXr;Fm2GlbKG?Giz~*DXwMaYC(A$C-!eX%>Pai*UBEe5=$; zj5{;0*TU-s$Szlyew)NGdUAPK9^S6&d2SrXxB7HMSmy;C-^BTgfa#jc@mX9CVfwMm zKF=1Kc)j5`a&1_jv-R>}Ic=QJ8n=$0kM`4*!Q?;d2krd@r|TB&YxmOb`=%PT~Je9RPQaR^A;5eMU~Z+UZJQQ ziBROs&nYVN)OqGs)_Xm5`8gA-YN|c?C8bp!J#XNgB6nj6F(|33yxk+@F1FRpEy^va ztM}wL)K++>2)y*grA!-7&zT@fQf{=}LEsIDnU>A zR98cw%9`r>abQ&BL}pG+Swod)e9=|ofnn|CZJ$(A=dCNL4!6y@pcaTCJ<$uHz(n3e zUdGYrE}GF$?X6tkK_hFZfz;Ssl85e*tuqMchRfUPXVp}e<$CLcqUu^SXB@t4GiKQ4 z6xni}D5+OSnzUCkCFsg5nuK=u(d#tNBGlrkrADG0Zv>YTnYMcB3gi$U;zv5>943uJ zgGy@ap%2D5EK;N+(YHV7NYt>PB-}X4tt@Tu^48aQD?D}eMM6bQL%pY>rmCz+pq4GF zE2)KM%g|-jRJ!U1=(!~e7L*K(O&b_J?gwa3V=b%(eL?(TXv`>!Z1oX!aSW+2(uf?- z0yN#jMn<*4jP$-^gC^W{)S1-nI!wA_Dus>=Dy{coY@<_)$_>n5(Y+lK?h5qHK>Z~w z@GPkJz*}IX*10h<$*AdeK`mvuYZr$Z*gaK-irqD^plS9ZHatKQHZ*E&Z@(5zs;n-{ z@pvn0%7|%|hiry)kW*4sRpSOQJ+G*qD*})~W%?~VO`h(7c-F+4YDft*#6bZl2ZI+J zYw8vYM*HeDcq3MK0!4Gn`Pi%X7A>l)^rAZ=C>oRXq&V5jn!e(haO4 zvSGM;=0+r|sw}Nt3>{aGOSw9wz9uCtB~<_}%!*K`ys~rw-Cad@ZvC#vQ&!?FLEo3w z*Xy%~P~@pD6DCiZJ|lmIZCcUQ!eo~-dt%YJl(g}Y+upk?!8c4WE@eXR{ng=_Nr=L( z38pz*e9{y2S)Fj0&iY(L5gVbWdn(=c3daPLOP}yQP&9VUdap2549e^EQjB5x5B^57 zSVLd$&pq}4DvK3v^-OQG}6ay9Bw$5 z;Ar7az;v0=?lC-%U_E~*%cq#|UdBI4P&j>fUoRX-NtaM5K`8i(#L_Una|{Co;qQI% zVYu`8zW7g9>hxMae30Y2`{56Bd`~}o>?&UGhQ9oVae7NXe7IhzAO1|vC$S&i%<)L|Q$MbqSH}-AUx^Y}@{qVwbI$d`^{3y;R_r<<^Ixf@s z6!gP)uGi^``{DgupRF(Th7tO?~;K^6|2~ zA70{o61McEAJQVQE43efU8auD>W44p>j$|Xz9~}CkQ z`O<*F1qOV0%u`0O0Y8gT>R%NG{BQ%l)_@PM*(kfwfG;r6HyQA>PaeLO8}R3Zkw{w% z_zMmA)du`c2K*WW{vrc@tpRT};MW=OWd{5P1D^KT!&j>Te|Z>*bh`ncV!(G8@V6T9 zyA1eZ1HRLMzsi8u4EXB|_-+H$ApoT{!L;K+0G^GPX96+@O+A885RRR(LkSI!1Jk;@@~FpL$uTwFdn420rTy_!|uP z4F-Ic0pDuC+YI>a27I;w-(kQ{G~jm`@D2mM(}16Bz-tD)Y`}LL@J<81$AF(>z{|dU zW;qc1&`L>={jJ`pP>1Z>V%{u-CS1Hq5<+9F_!ckB!X3dCR5I8Z!Zqexf@vBEb};x1 zg30v7vbqtOom|T5u4TF!P$_jGv!4?J|B$!-#u!+H65lq7}Sj*r~38v5| zSj^xL2`1Mb%w_OSg2|-^vl#q3!4&!gQyIL8U~=t2iNP-rOfEfWVesDxjwe`P@Sh15 z3GO}##F*a`Y$3Rl!M`DR1i>8)et=+d;lWl0-$yXH?%+BG-$O8kP{B0}zMWul;lUOL z*AYyvJJ`hFc?6To4%RZblwfkz!D0sAL@W?48DwD3gLp?-&6Z9CRif4lfmZ_d=bGN3_gQka?!z71`i{c zTyt<8gQEzx61;}N$5p`Oii0f-K1eXR;9wJjzap4iZ?KlZpAt+iH(1Qz4+$n$8_Z?! zPJ+qB2D2FaI>8k31ydQkiQr^{B?iAhFuB;Eg~5L(IE7$=!G9+B3WB>&u>L1FmEcYW z|At_4t-%fkKR__K)L<)v?<1I8W^f&Y?;)66WpE9HZzq^sWUz(7bp%rg8Ej(kJc7wp z25T8yN-(*|U@?PlBA8rbFqgqI2qu>p%wq5qg2@#IQyH91FomK)iNU|QPZG*oLStU% z=}w-0$kVU-qja82hz z4suPW-xBI9PZ2L{1~bIGo5?Hn!YKdE;~-Qge1&ADoOBZ9<+hX2veGTLX*bBGZSszj z-UKk*q38_bh~Z%Tnld+nkI-<1?8|%}1;7&W#L0oos1=f+{R4t)GR$II?ERoKkvL5d zOE|PyZx4~*5-bOW4(G|rQCWFIR^F4<)9+=Hj!L*y5a55MGC2@4hAf*DnUK~B9mWwa znUqO4NX^ixpdI-x6oM+|$?6hoT<9I`Nx-BVGf@y_0I9#V;UXC~7XDReO(>JE9L;KI zxDxJNx_NU;xp<%zI(w#E5H=9&L*h-SRr64P;P4QG~mPwVZXg(3wz0k zVlIJDgNaElHZv~Dp0p#{GbF9DIk=edXrhMnv8PutRluHdZ-&6WeOZ0RM|v9Mi@jwT z8mYRhi7wS{m1v-_MGR6la>4@!DH}j!0up^Cve8Aeu%P{^hZ-w3e*-`cEU_kP?MFkQ za8VGCK!&FKWW@Ohp5pzX`cuNVA5Z{Tf~m zJ?H^vv(DN3q^#y!N6M-fgOv~IcOb(6Lt1QJ1ZDA2d#tptl|O@)Q0o$cGnhtUf1rGg z4TOW@O)hKdA#sEYgSr{UrTfpcRu&d%TTqldDcL_<6xsiKPUTGqKZArLO z7`q%8`+2h@7`mg`Dg&R^%6p*|wTB{|K@z3XS9-7J<7mvtN|jX#cMVd#JI+PJsXIOd zg$w$D7QNQgjT9liKw;&pbf;N$^}u?><_(O}WleVlDlATSn_Y2PM>;dwME_kpCzWF@ zgl(oVBddkhESDN<*IT!vL=axlW!m!sZMdbiYD-C*fy{G%gN7@+wNeK58o&zHX4pU+ z^QfY~o7gVhLK=#AfJi-$AlQZ$r=gy%hcRB@PhTQ^YhB;6er|XIRXy`-RHaTZgH}6) z26nB&c=cKnor<=R?D++fkhR~1CMIeZ9${r=xT0J+Aj5XTR^(C+*h2r*R_e^zyubzo z^chFWy^=ae&D>uMso1DTdv7Cd%z8#4SDT@i75zhrHjsISF5E;CZWe?K9-z5z1qzk7 z^v(yie9U-9d#z`Q-tJKWv~1uZ{F{6J^CA4jVclj@pddL=2%8B^hyvfppkrlY%gV->qF6sj50OQ=N@G(BBsP7nm15MR?sSN$nLXI_&P?-HE>hp^@g2tOK(IonJ zqlvw946>0EJRRkXR%t;}#SO)7h69-2-d>R$|*X$;(;!;%l>v#X>_I%dCAKLRc@EyROf6GY!FWd8T-T$^dpZAS!jAoiZVIn*fjo7oJ zu&+ITNA%3=MC`eMv-n@M=XdXf@|Zo(2XoB6gV^&eXedQHu;+^(L=km@d#1qb*$qAS zwdXr2EP*}EVJ+bHoD8zCJtwlFyS25<7pM^Ub=4`^P=aQptC-6{2Jg9b+t#P_jKxC^m8JmveIl_ z&eIl_+O4(JvyGI~L>Rh#t0tx`oLYZY3=d$Abn?enj zwV2|K0!{lG&Bbztl(3J9R*=`hGV*X zYcXc5d}|^U&4m@4$wxdIXT-1b2@#HGSvLnJ ze7+c^(AyQ-9S{Z`SR@A=eb-B8Aeht^=(W*!AE#Gz=thW$`KzAlI)+dktPw{Wv@;@? zz|+@bCSY5DJb%IhHhmq?oYZ_}tL)$JeOwz&UCeuQHDmp?wwk#h2qcKjRiu&doaGQ2 zE@W~l2g0*laj=ABi0pT;`K~zV0+QG3wKh7HT?532jasOS>5>J>4pig{)LLXW524+Dw6u zl}bytSy*7wP#&L4EG;yS(`})5wNJkag>J%VnBzhKXJ#SmXvjwM2S%Xi?;zXEfYb0= zoW9Gkcy?k1>|#^9D`SaO7MlyOJ>m3Sfw*TH-32BiJ47@84Q$_{iMSbbE&9KQ7JKy_ z{F*@6ZvmP>zGSx!BS08JjZ(G-6Jcsv1_XtiqF=+B85VksRz~WC<1h;?6ZRDB0A4#Y zOtnd;dPk?a4mO~jO;mp-yqY&cjGl(E05r1Jm1k zfbqyYjrc|D{66T--JxTSG0XuiO6PYd!jJ3Y(cTh|>-^e*(VnF~Kqsjamgxj!&A~fC z5wc(o4a7bGM##N#8vLzpqS(UH-uxURw;TRRw^K%0p;OKwN_sGbF#*co6XiBgg4ZOS z;(1W$nSRRrlb(5_o_RYTe1C=r!IgSeyq+~ZQZ9$e-KJ-LQ^y3H&CAtMxmn2SoqX6n z-}S-<=14_*4@$|os6Ep36d0Hr zk};-&G-Mhf1WlOz$2GP@W!g;zW7M%nYxxXxYUWzvu55!!rt1Z}w)Toh`>3F~A2c`L z61)HeeV-5H$${CDtb8Z?cX~(2zROGv!{xx;l8ng$GX-k7IjEo-3_5*(%!YWGf@Pa> znD)fprJ_ss?Lq?ri?E#@*|2Z6ZU-){)qg~6eiPmWs+n^eR1<&$+6CnV>JyyGCR@f~ zaph|u1lh_&GYaf{6`4ZU<>C5V8ExKSPVz~st|B)6*Tm*p?0CuIGsk?JOc@97947~+ zO0s*CD-hF$d2}o4FoSII)*SL^D?sZMpXqTaR3$614=<$s@ z$>ajk-GLeUU7;&pZjT~?WZx#UuhnGB*rIF`SI(uPfnj>7*8w>Et&87vxm&y57n{d} ziz98n|HxhAigsd(_~RJt5AiLrRcmhsJ%_7 z+fp_QRWppE~ZP@D!47ag3nzX#!fjrzGu4t!N zoU702+ibRF{IlU+=(# zL7{SBN{>!+1T$5y95&^f9AyY~SGy_W76dw zy)%r6BE+=@L^OKhuFDLF7($!_gs(MLF9nOm1Dn+ZY!z}cc8JZt!PM?lzQD4kXFA@y zr@_3D#;R=w+fIFj?p?~!v{nZmL_YdyC}ew?P_&?QDcu-bzSbD^{?N&jCy(y9==G(a zHsSJW%8qFP*Z1^*r7dVVx~ZuNOw`!Ra}YP!cgU(7rM!UBfpB=TdIqE<)7iMv<*XoXDDNNqo;MjiEc7m2+qx*vCCcPJE?m;h5* z2`MpF_p%25Xk7?L(VRjs9K1pdcB;+ke$(dLqBqGEhS>hc-O3c^ig- zRXh3eKC5<85_a0yv0DtvpdZhzBL0MY200r4YmFn}t404d_zHNFU>m}D%u~3{QKG5X zt3Z^u=z}u206cX6qoUyQ~bq4@@ZW28Qi z-VajcAcz|y#8*X#w}2S#Mi5ukqoHE+o%kk~?cWQRy@_15sC*B7rR^lgz5F1lE)aV? zvHS)sVC{E&{5ULTHMlY_H=eYIWnn0-m061H+h*#rLNDa6?>dcXIGIL||23}* z>ezwNXp%Fwd95@DuqlwuXPE!~;UY~@KIXp;&@_DP%GinVkGe=_>V#CZFSr>X#S`ry zM&&}+U+l_T!Nv3}G!VP2l4HomwK3^2 zO~L=zly5OEaUVkJFY!go&SBc4kYCeC=@6fkU}jNe@BBqU z*!Mcq%v3Pdrt-3v&=Myec1tTbWuoLA)jrWe-6bpY6Nr(l>_RA+_EGOW&EvBFh#ue( z@ix9qLtTTY-qS+y+L|+ zx7o0O$L_}_RW`|$t+aHLeb+?`4SRjZLHr_!o!I0|5LZ5nTSwXvgiT#bT!8@fN)GJ3 z6if9{kR%YurLuu*e8=CgxxV0@x>^)A>=Am3`S`L-g)d-u-8PVBiPINkf@%Kc1GQ&*tg z>^m6}{qw0G5PbV|J%ee(U(r6yA~(B~PIgqmb}K#nYKA0;YQ8p6??4>c{T7!9+oo?= zHlBgp++NZ6c?6R1Ij1kh)UeJOn1l_RqX@%|oS9Yi_CLkrPUQ&A;xqGa!{PbCfMM-{;Vogy_q4ZR|zmIMUc0K&C5@`Lst8oPnXC zmI}Z=6gvYwC;G=y{23T3V>VZyuTV=5(>!j;(9dNHh>>!j4%Wpw`~A=?wo8ocRm2tsyY)ghnG7LlL)3Tx1hr5fAxOj?UB*$d`5#c6B_#`h|rbyq3s<1ezDn$ z3T*8TCaZ0m;~Y1br|tB8o+UoG2|*4#n_azU1k8&facqi&S`tYUG(NC-eHHFuT*_|R zEqYCVc;&wcNcf8tR_cz(Zn;a*ME`LJiK&-5uFuq)CHOzO`*w84`DGGw9ceqWj%>|R zc0&q?0IQ%S9ArUD72 zLoMAgs9bzz7mya(8Tn9~{VtZeGhpZx1ny@R=M0S8#^#!~uWrCU3=RMD$S-Tttpqwhe0_Jd8${H+rQ+fy%v!a+IA) zz}p9-(s=i6r1o_aqN%E|bm0Y^?Iyqb+2hV5jG7~nHV3E!j7 zSo^b@_lwHj6PP-$jqNFwV0P~V!1{KlxT=$8)IjE9plPg!&E$t2B`AgX} zoCQC!Ke!81YuJB4$=<2to6u)?asqpevmL?CbF|pJ4E=`wwq(eXxMBx}7PC0oprhwq z*yDfJ_j!yr*6E9}nsSuSgG<4-B|>whPP59UXf|JG7R*W#{nw+pIpCR2JJGfb8&H=6 zWeeCwtGOZH@zWN?A=j+lNcRh`RCB-e9gh;1xgd$J=WMZgB5u^&gFuD}>3-p`;db&n zoiM|rTj@0zYKukmj{`Xn_!j+t)-y4dH_{X$`X2+r2D&VSfN(v{BFL70YQm6PdYrz+ z=J~)t?kKUj1+>hx(1|-?a2u1rB1d@}{dnk!F+`ooiuhkb5nI#U=Uy*1{|&b}?&*@9 zWN~1dAzAjlPV+h}OxeNg=-^*v1&h|x-v>+iD7FNaM1w$Ugk?DKXevE0X&|-SrmJOP zEL74kLf04RX<>oy*xBB^0rYbN5uMiNz1k>&U3_KIIrwr}FFg2HJwjl@)b$mn>z|=e zAEj=CC4`mwU$|fTIsIcc?^SC>+Xm9?kF@i{8~U`f1*{R@EG>dv!fk_7jDTW+j;527u0fQ6T(YsW!!^OZkd@yY4N6 zN*gnL-CqX7{t`QAQ=;H4cPMS}m+I*=vrxP5!)O{7Sh+L*sWt9|6?_IOz+4yHix`|F zm`HMxckdOJxqX?0C4zBb5eb+;0G@-&+a}6v9vIIM6Vrp&!bA6t8FASgi1{&Y#AbSK z6nq?C8!>Q&;N7Ty>8BLnEj>=(-aOqE_<;|*f_B(*Usphsne^2#Tz3QhW8es{@g_(k zu7CF1sHzhgwx-2!75??Gw5~NA_%8rKKP!dS!njfT)^J74K5`JZ;A|gUkLbS}gk8%q zit-9Eiqr`ogAzSVgXy9jV6+b>$F$iUq0ri)V)H6kfo%f~tsp(STMH+@grQiA9k>OD z(bpLy{mWb6aTV$a&llnSSqk63MJ(8eq0%Q7q)hUW)46mTp>i_56`PB3;|yGtqa4OVA7|3X+N;nTog`>)*kudM zGW(8)#O8OSXi9mQXmiL_{1z~_iLLUrvFmP@FlS8&VqX}~`+@iGlvNyoylX5=*a}{5 z-$rwGWW2}`Or5fdwP1_hn7e@ZQg%-U?Z~Yt2jg#troPtJqiJ;90r6-f zqw&2;>kpLO;(qODD+Qo^AK=+u zr6zZ!vh^6R6R{kWxspkAofXyofo=ozZo${uqBmeS*bogoVx{hl?CIE*9?}0GFo-S5 z5cX!Ew@#PcGbBWX9oh?|N|JL8lk+#d4;s5xmbYk6 zQ8v|0y7A+-Q4;+HkUuy9Uo9kF&p)_Xh{_z`Fc@B$Bhg^EFq7&^*Uly@z8+cSEq?tZ zhi(UF<9s&uB`7Qwu#1EKwq?vR-`(X@v8Ty8Qyp~&d!+I;Eu@UygD&*9-euJ{QFI$( zaUAUD-5%;-?10eIfbDnhPx~ksADO=l#^h8NT9`p;?a(i`CjU|Rm*7q0S_8313s^H= zN{4vp3A7{!O+}~Jl=p+Qr&Hd7$ximYO2#GFF!$Mxe&A5v4$eTHR`Mq4ClG5RF4)1M z7a_3PhDMrbx}!aaEFpN_>!9i#f8p`1pUKCQ7q|v`VAE%e(?6#qk4-l34|ykQlTel3 z0PJ)4qA|`Nv=R|4Wh6VTi_HMhf5TAl!UiBR43?tL>*Uxh2*fb4iYsZC3f{trMOT7L zDZ}EAujlEcf>XiMi>G6mmr%a$I!iC_T)IW~8()ac-B5-zFy4xw6kgJqgq7%*nh!cg z%=*49FdGwTmn#6a%z?CnIwx?i?m$|_=1kh=-+5OG-&YOqXSzFJ9!|yezJ=Uxg449e zcI{;Ik}r8tmy)@30ri&4TH#a|S)A(KFriVIdoglC9nNSh_p!>YKrHz>bdSY;@5Y1J zyaN*kI%@SSFhGFMHz2U}K+AgSFj|1|%?3kgj}vFVL`O5bW9^~yOsF+;Ena8ZcpBPu z=t=Qct)qozG|SE0h@PJ?0r@ly@O(spJq~o_lrV!NA{f?70MQdn;w!BUFnY1Mi`hgz ztLo`0^08lNW!qr^5GaZXlm$g({t-;HxYwWn%(rlIIIV)`J}YURP3KW6?3f@d)V+z6 zmqoZ*y2J6UW-I+4Hv(h=YB_+jTFur@av;bl2%e8ST8_w6w^~s~=i5{1e86njG#^jR zPoX*zAKtBp*aTwJW=aCq@$69{xij$<`5biBFCi{1+38>i%totghXq1>Su(Aw1PVcRmzY;iWOHrgDO$a z2K*~Eq1_?WyGyFd6-SKIw||K2+d9PQ``#q3%q0G5gK3*X7)A?*q}|#6A`oN z7FL)=-tQ}FJq5Wwj35&Yqib!5&U~aM920`$MGU$zg(P5c10PD!*$J+|8frg#CuW-~ zkZupGrHSs4{71AAF(?Suk-HZC52IIKE+ZCLTc^@Z@Iltm0-i6cwWeSSycV-V3TY+T zCt&a*dI=t4p@UtoTPE{LXQH4E0wn1DqVc*jL5Hlc7Ngx-&sG|$;&MM#xLVEt$Gfmm z2lZhOJXUZfwWE>*Y*>rUd#IIYAXvMsxjrq5YydvvTbDX7K~`h%aNuhwyX)d`{DMbS z89jHV(9&EsZ3Qb9X|#ZiGY}7H9=(#JaVdC>?SNLp1X=tB9h^#aDWAZ2s`4as%&7es zVJx&|p%&I`>$kAp(OyEzUK9=c4HT`wiydq~P#<2(rnHYVE9lTh|4u}Ebb`JjxE*&i z-!Vp8gQIo*%I(e8?U0Tx>c0ZslTakY_2J;%5=&L`DKq{2;i z70F4PTxwM;o$taT<~`c)_mI{Kpf#7eIKicswMkv)hvkNzy6?OU56{d_($lt_j1%4o z&VYzo7o%3v(Ge_2w5hDB#c#8ja05x|g4>p4)4R^3t;{YkEbjWmiH%C9>M&uiOGtA1 zPKtNEB`0YVs{a|Plo7{>54MWWwc@od2af4%>e}eUyHqp*<`{-UBb$%LE)9zA^k$d@CJ99J8tKp@@<#radU|ooQ5?5YD zP!cwllB4nRTW}vnP}*xbN_X%sj4^Em#5cbG5$pmV!~6{s-H|KdWY}c0uF%cBb34uV zd6-}B+hBWofmnDU!cN|4v7WZ#FP&XiMMv5%$ORK74Ibn3#!;Zp)hh# zVnmrW4&Mk8@Qu>=j=h0s6E}lkoAx$NQXjJODlB-FlYrPtR(?*;A%|E<)!JFX%qSAj z4i4XXBtf&|p)yfJULCPH)=k8Hvu-W@D0iuy8_g>8Bt}8v5AS>sEyi7)CAE26f~-diCfXe zDyvz0f;X*9m#>ppj=35Pm_Wn@IokCRxxiyoE!8F#cHo{=opP|w24(LRQ?@Kj2k!h_ zmS3CCIs~##^;s9bYqw)@q3N& zW1$OXT}+2s=#^0Dg$s$BI^hn`cm1ga>tNr$L|RT@x+%1dCUqas-lH0Dw?}&ucW9a= zb|~oz4aMnNGz;Wn;=8yx|0wvN6GnBh)RM{5*e`W>!e>a)``bxBC@bclp22eFU=+>S ziGhJpG3PsQ=+V=H<6_FuMTKHY=3j`8l`@O>mO-syYl+hqq9IHQ43{~vq)##Na}pEk z2Ig}IRp=1s;AH-sE$Hay;JjeBZQ`=g-~+;2ae8b^BBFRPv30Z2ACWpF8vfq#s9-4M ziXF@39TEG-QVA`c6MGSZZv<~k?0zujV>fmw@|h0FPicopanR3vz0@x9J813sFQSB) z_!$_{=VLMDUkl&SZJhK&1HTZ4J`>%K(g=6O#FNCPnN@`LLHlz~Qd2RK?P6j(Zs)MC zSd~aH_+ab1#CrqiW-!y81l*Kav%O3_Q3-|~DN0Kv1UQYkluF_T1Cg8v(- z5!Zq(sEo#V>}7g2s1yo5M2bf4TVjd}qi9(y)C--GcBPw&y7l!ffZ#CeiQe{ zPz!nOp<>GR`b6A}`fV!old@^iSz=1ZqSJ!wp$N9l^RXKtNGUBcwBO?f=Oy@gyP0|H zjXZev)nwn;R7(%U)1=pI-Ng_~C^$cgXa0B#1i7hpwIqaeq`nFS8pq4mTRTAObbl=S zHk)Ko#^b3%ypG~N7msS3${rA9;TR6Z*c5}{S*TNVelNFu8AS_w(f+;m){77hF$?yM zFyyd-iX9JPblgbvCnBPd9%627#M3~$LOewsN!u_ud8-|N0ZZg1){zLDnt{O}7AS%^ z2R9TEq&t*Pa@5?I9A)ctb@Y8MYf8$L@$Go@WI|`46 zeaE7`Q}G>xbHcu3vEI@69wPdW#Zr4`GR}!;QH;$O(_+eT=Kw=nzralN#oT9d0>cC) zj4x&t5Dh@k7BNGt17aX9EO;jpr?NAf80BPid2b|cPWPQyGi;T>NOQU!`*PflZq$1i z0#J1`dvNmt)q4Urf5#EFn0rhRg?-&e-~;r<7&bhEv!ac7);;pk7^XvVVb5-)BO_28 zJp@I7nnlm*cDmHjjD|Yc_Zbv577CDk-N_ZEFE}#= z*rP|{DQbRf47OWhC?HM*2B-WuhPo)x5j_$svTk3{wh9MHw)h4-Pv-kgkkf z-U5*61?|!KrWndI)18^_;^;1p?ksd?p}UcEH|OL?s|%AyfeS-H#(oo-Aie2ipA<-5!ya)~ zgGJg$wj*0j-da*oDqALKU%|}DeZHQIb)s&|=p6DTbB1VJVdS=qcDyjMiuMzTYXXX9 zVYmbe{esKgJ}KbxLy4~^edrx-&%i(xSJ7|q;JLMlovIO=kANF3za`o`NDZz62>Djq zFa`p4KMZUKcJLnpk8MfYhzDL67gs$EVVQi0rC|e_ol8mYmHj}<3w0@N%BH*TW+K`# zD6(wW7bD@E?w$5PAyu#|X**Q_m5Qr=3o*dhmhmqy&JM$YQ_#K}j7S%+Ts}Z88)2!? z1x|vAt1?66K%fjeRyJAeReL?ol~;G&OMZ~rDmm%=Tt=By7X4Lht#Cj){5%zeCqd{! z9ueOciOZ@P1X=)5wIB6%;tp6h$x7z!SQx(n0aEFD3>?1`djI zo?td{P)yzeBo;j>pZj(Rwv6{NxF9cTts=88HubV_cyS!BV8JN4n<*aSTjQ2D#!~lwG+>= z3Q69b__ZkrT{`yIsA-?M|* ztw3&v*liHok@2>;atN2I*tCoO84R>M(kH9q@0Z;}n3li5GLLE5={`nWYiSTnQsVyL z8$v2RuB-V=+0>zo+-O&HcVnzz}J_m|1Q=oQhEn1(m7WfVO#S&h{x?IpWJ=T}%S zHpa_spGMI^inq1+^@f7?Nhg1Elyx4+y6SM1s zY+^@dCo=u`&}41a&SNasvm*54L^|~@GFa+?4xSc5Om zI(%Bi5w~!ayJv8llZb0b%01_>^kL)??wOH7^}{W zzw`|g9>B+i|LZaaJ@#S~`jpm9o91s}|<=s~ce z=dKoRIOHp}FaL&h5!njpWNf@^<s)Oi^Q4pymkNTL5q) z@yk(vRfv11nmh#h#&dS{rYQBUC|mn7D)Mu>`;h}T@`o@b~T2vha@Mml_jbf(xco7dFgi} znnz&&R>U;J`Wsd-? zE)bV}PJuSLM7SoqWs)#eedKFOo$9YB$x)xQf9@7w<29Kj5IjeHFrM&)7}Ep<-JE|f z-M}TS(ort~P~@1d{v z0iHZT{P8-)ye|lK5??1j0&3p37`b=_qx5@J+l3+9O^yP&rTg(U?{kopZo`VOW=1P+ zYxm=fS?#yS2{DEF0Tiv>N65Y8vrfK-ao&mn-$0qQ2Xooj$7r5MrE9*zkS|3AwVx7k zX&0nz=3*wSkr3a~+>`qX9dPuCfy_+g%1L-`VP976S+)Opqci;o} z;_e18gD_Lb<75MMd$%)VpXh&qyodVojAIE<`lESBoT@(`35F$g31v)dZ~iOY>$#84 z)KM!5#VKXn-UxJ#`p_&^Hs%!~c9{-10;`D5)*kZ;r*%bW|2Ug*ig}FHZt8-wAVw40 zV;*C)f#-=9(iJeC%W`>1Ti_YG#pE3gzRRezEwCydD3?itA1a_sWRM6pUn_`~Q6`t{ zF{;&N>H@`MoG}%kR{6F?qw*WJWsczI${b9RRGLED+8;lh&1r4T1wCmqS4teBN}RkuB4b{R1_+N8c2Pdh>8;t z)Ilnv+MqpM_l%LQ8Dh^>@-Wdrflvua4V5ry6(%B7Le)YS#F*601Xh@Q%z5For+NV# zN5pxN-y(m-=cZf0{YCncm8Xrt%_|_8s#c99zouRLR6h^)!qd73yBO&-wQwx#YMB=j zzyndD-eFSfqwMO$7|a&-z%V-|3?)5V8NU`TW26$dc9IgaS}BOhR=jJm`frO-Tx;R< z?W+N7ict#8aL8-TK$;0@Rtn+>5eMoxLc}Qr7JDGO$G16_@Rq>sC+vaJ!*+OmkXSJB zVfLBPtu8_ZuGKio(K-?{18`G%>`LijrS8Pi?^aJ94IiHu_VIKm_WV9GJ)gz7nuxrg zhIZ&qU%TMX{gm@4v_)ErMLJ#0hkI9oP)zX!w~?Ay()LM=ZNsDBIaOWpFk#TauRa18 z4^kciU*VPC0 z+E}VV7yFdk_9;seNvCGgCy`(dg-&n7ix;K)K+=tTQm>iROL=ibWLIuHj0`OBOS?&@ zL_q2$0<86K#1086^Cse6&&2Wk)*+zmQ3^&95w-zt#O}%FAs~X*Hw{4#nMwVz=!OK} zHtbrI_JIC`uYHJKVd)9P5kqWCPu#rce7wSwf*9Ved=?zex`wS3g0Ik6a31{tk93{t zy=`GQAtshf;?a?fzu_>^ge`f9`gXJHveH%=KUoEv_BZ5G#;i3ENsgI zH9w4+T@N9;vj2(%c{mQ>{rE2eSW4RiSv`vP{@X{==B9 zh+BGaCstMgGfJ?{B@4|~h!@VFIL1$RXX4J*eh=M_#I1T`v`uY{rfO1ky|nvePNrCt z#uz*u4>ZQvFm?z{MVP+Ym$?8kwoB>tlD25ENX3AP>!$7l2perQ(~LUkh9iw|wt6+dYhy@c5F zUqB$8j5Npn=2X=_R!{?@o)m;>>PQu(1hsy+ddDz$x6`nLnuzX%4W5-_I*z5|I%Lz% zRjx}3I+RbbdtiYbO3J_3P8NyM1d1H|&qDr>h+N7`*0n%6u)!O6gyqr(PNGe{-6Tij z4KjzCJ%sI#fRNS=*io|qj}t z39-2~mcmZ_7ZW?a3VoYmZ3=7v2T*(`v0ZrG^DZ@Rh}e7wsOX1~>3=Sz<6dA3djJK| z0PF$18M5yrqWkNfw9_`4L-~c+d#GNXPr4&iF?JHNM*1|V?h=ggzf7yU8YtL4Z*Np z=1rwD9_h5D#^_#Riln#tK!C6ayNkajE3yxIV ztui*Zv7bL0+hjxS*FMNBL0m!qi-n9E<~erA*_dDOpR6K#4G!g>E~taL-ww;5{f6-N z1KVuqGanyCuJ-I7`WWjTOwbGP|Ly)r7Q07#=m}=A4s4?iq4yPI>6wGQvpwV0s zxCyj0%GM*LQxwEQjz*PJ)zMe6Jsnqc6+naVH`)GJ-bF14$4smF-)DR#!0YU7&7?}jKqu7 zy4dMz+{KW*01cuMNi8Xc{ltjPU&9E=X0y;Kwb{z(H7JV%S=3XJ4wuwE$2y$5ft?^g1$Rcv|#6~~IrcA8Pajs^-aLuB{+Hl;zp)c&a8y$HPD zb0|52U7c-`O|OFhj~bRdgjW)<*NoVV_VV7dsduC37F2opR9l!hfzS!J> zTJaubHo_2nmw79;KcO}F8{hGBY|7|mICqUZyYha+^HAqGj*Jf){$S7ehyI>!!~MbA zX*(MVo2bORg0ukp9@x+=r70OxhfUc@PKtIh$k|!I*P+BM3uaJ0W=*)nhe0zDdnv&o zz>$}#z*{!2lKt6~>hzEA z3ZNAsM>)=RVQp${)FfreqiFAY!7o_n6f#h1y4|H>O}aJylCtd!(d^jFeJYa zT+c|!NXZBhrmE;dgbpy&x16TeV25Pb-vsMG$1NUf1=K_HUaDPb&!!@vuq*G8al^2K z$&4qoMyF=Uci^drHFuIS`u93J){5Y6Y&c`@+QuNZ{laD@HcqLTXkj6>Fqj9z-t#DW zf9(|g=P}D>eJ1+5P-So;FplPryhFKtu`}#0bxhI-4D*+pz>7AB$f0#0<*|kY8*#x= zGUHW_YWp75Z4xh%GNH8jkJF0bsP@id$XL3Cwy}GkrEG4=qi6NsX`^7Vc$OX=eu^FW ze{1wsuzpfmqu;C0X$BLDsY^Y!{4jpU(e^f(;P)$!{af^f}%KRlble>u|d z`=p?y>&NN!pGqaV>C=CllKpxr=zg?%MwFN&rPX`!y3X8&1^B5{DX~6Da+g$J>Xk}8 z5`Harp{LI4DU-Z4QaSxBFjbCgu<9c}^Qw~x<@iNml&-7}v%(K=_N}iNb>P=yFO_f# zYhY<3JV}U&SMcAn)n72;)d-386(zNvM4eZ1T2f}_f|AB^_T#fj5*Gr$C@NJ{mX&#` zrAX$$(|D%rn zx9Z0MoS&|L-2brt>4S=(|Fp2)BLmm3Q=sAr-U@?-dnzO~(o*Ba*H|HN~HMZ3U@&{MemG$;cF&S-rWp!?!+^p)+`01J*OHdfXi0>-~jpsj^aah+}f`-yQl!3DL< zA(qT7sjPl*0?V%AFy)$^-s7FA35(c#LHi2oa$zN*set@3d4WY(TCj~9x-G*DJqj-PF&Sw`np z?WwPqO1x5?{(I?Y8f$8byQa>QvY?>~Kio|}?H+dEaQ~I;C)n{@hHQqa#E-Y*x4rS} z(_zU699L2zH-`uVU8Q*II-gX{TYc$(g)b$=mo#{5#)h-Yq(zk&pY${4QX=oBYsMv| z44c$Yt^cz8HNzxyHBD>Dy>mf$>`0lCZnxpts@LF?dvA5yPo6(N8Qy9jc9E$C{30_X zXGsOhtBV$vR3Vd30($I_I8c@z2h=B*SJu^g2PH^L{3(iYKSp7Qbs$8FP%hMYN~(l0 z_1E-@e@$dw&Ky(LFRPkR)21=3uD@#-iWqPhrFE}Hq~YEw>g|F0aYp+p>NT*U26R=0 z&^wq7w`@?C)9-&9k+6wrGhHsJwuXM^yH2pq@GPvPQGAVL#gBD!t0*eMeOZkg>44&- z;D_OB>V!T$Bow*BOpA)DOBQ%gm=owhoe@aulQk|38FVe^cF1ycuWGo}jZBz7UtnD? z=z0_)<|JIle}!EdQ-6IgpgXqf1!+uGnIPSgOrcmxS>=NI%$orY6RHs_@}WSppBq1- zF5p+-5q=?pETaf;O#Ltcg_DIOK^i+2;UT;|Q#7RYR?cmxX{hg)F{W%R2KwA;Pnl2y z?T)42geOj60r3m%fbmoBMP-$4FLXAhUK$IQ64a0JyAbnENokcwxO9yB(y;<_D~Tk2 zGJSC-@IR;qLmO@JR5PEvpaD)AfpB?cqaFlSP z)o_*giF)W0O7PSt4f~Jzhuat7U()EofE!aLC3?!sVT$uH){>H?NGmuhTy8ME45m*h zk-P#mz|%;a!+ZwOQ?HJUN=cTu6p$+-q;93vW#F1j9ZbfG*YL^zhiyo=>+nQQX3I>5 z`j@OU(nic!Rh|+UhNnswppQuqc{t9Y38#0CB3DuEO{}bjt@3}6z=-kHcNByLGK>aJ z*hh}Rt2gi?gXcqQL2Ze*5`%ySC87TvIhSB;R+YFtOWY;yitxX6NH&bvl4{*Hbpgvf zRhaBiU-ei|V`Y7xAzI?~)>Yyk?s%wPgm&m8UE?Y4Yr=g%-trn~w>q&lDRa#2W!K}A zsgDUV!KxYzc^c8+gu6P=Z4H%m9&~7BwHt<0?~#Z$Y&56`+Y9Upr@esy(L+++C#-c# zuz>R<)Gi; z8i;QC%=|p39dwQvInMlihh3Usn>HEbgV5(spDyL)*jz4NJ4>3305<6EHo|l&(1p)5 zd~)%j13&ao>KuH|$Hxp}!&}k@_OWFJ@w~5BemFe|8i(1ir-*5FZ$XnR7&>0@4ir|5chAc z427DJoUN~ zF^ri!Tp~XBSYZQ8h!>YBH|h_%O!9C6|Kvm%iM*x=E21FQGLqL9QxIPq_1-X63bjMf$-`Nh%Jh3gu-y-XpuBi*(flTt09Q5CY|wk*<8vQ{$zmZ1Au%Td@zXicgT) zZ34*Fqf9HxT*}M%qDtan$^tE^JCOHJDpkw`b^vXZ=%faU;~09-m;SV zUX4aOSJ0SP7K?uM4T*X=KHU+Y<%pNF<8v`eq99Tt>@d>TQ-)8DPnQ|RR7O$UH(yu* z_Ez+_<*%WTndakOo%Y%p>f&}b6&i`&II#8{!>7k@h`K4h6^sKbu8al`eJ-8=yx3%E0IAHHZDyhQTci(hLYB<#i)Re&L@+u|kvw#}mkl8y8Mq$gyO1xE4q4f$@CwlvK~=hUj+&z1*@m;ZAd%)I{M&;VvbSCv zwZ~^g9c2R_EiTVTz859Z<>BsBdL_I>B%KTcOZ_vj7ZayBHPRcBjO5$40lG zHIS*m9#A|PG8KFt3gI~z*V%ODcPLLrULo>EMe_0@d1c5e9we^`d1Zs-J%YT70qUbF zPF<%oa7qKGG;m4-r!;U%1E(}_N&}}fa7qKGG;m4-r!?@tp#kXA5D5ZZdMi;g>*j2$R=3)rDXcGw+ZDrB*AkWZuDfM+e;b~ihUF#&4@M6LS&N#ep zjwL#~@ZtlzEL`v~zJe2m6jrG>FHTKy@aP%@pKYTXYsV0rwe(yjHeIt^husR$J1>* zeUGPK@bm~zhg{6-=jkOpy^5!kcsh%xcshxv zvv^w0(}g@;#?yy+`Xo=ENv~j#@BXK#3@2-FifP#Fx{}&jPaQqImWqqH_Ub==q!e=rJW5<6YI6M>occt#HA6B#Bq}Y$hm-hdge7(oienhCqj}7*#HF6N_N0g!05V0LPc{6sLesm)7(#lJd^fODPbbziN z2cVK8NBl5+`j3nVub;yOh5Dbihct{s+Te2F>?d76&<{HT2N`j2Qa=s%Gfrvt)1F@b ze^+(0zG6uF~i$5d#Qk}Mm>-))5`u-Wp z`}w4~JZ|WGR$_&pJ*DrTuf(5G!ZKDT&f@xRJEiZRp}n6^+IfXuW({8#zkN#IKVOMI zqeSW@I*xfz;VPc315J4bByKJJ==$dNLTe`V z-v7URH+ka3Yox@PrFapvLAp9+LQ3k`D;pS;c2`<@O6vF|9m@_OMlnt4TjSvY^ei3J z8{kR?E=su7slRrBjsEL&qHMmw*}n*7M+x^n5T?&6(kbrd8;0FHWrrvG7Y3}?m-s}u ze5+33IB(^Hvx~TO8Um{=}yA!bd#h> zLKe1WYr>k)2_&J(0$r{MxQ#mMsH6Yiugf4Wy64{8U7gT?mLwj$ zcH)h=(d@b2i^gAU2E}tY9pKpnJK3I@7%uTHBDQ@yj=t-0fXov)Bnh64n1#fL$&ow_ z&thIo`sXCfj(HSo6(kS^*}v}_3c z3w}-Ns-^(h7PtoL;?L3uRHeyZo>1Z zq@n6kAQuNdr`E#+*&n!_E%}1n5V(&>fglG1hloTJH$wLYE~VLpqT8XsY><*hs9P}3 zBbH|deA}thVbD6?LfbA){ubXV%}&}+d;KTVt5*`V{0r!mt7y6H?@OOtN5j)F$@56F zDIO&gqv@x@R_4nmutZz>KMwX8R0km+^ZDuy+bLlHe`@j|o9YE(N4~1gQJtyjN9OmMp_?BFI+`UxfD+D-`Yelq%Y z3>b)rdfPFlfnovQgq49PLC>3TWas@H(DQl-;63z~Hz8pV=w3S^geQl=)hST-Z?YPa zGZAujYSC&W6TgM0uc;#`#`1p(>8T<;0V-9Y*+-rfYOxc5SgEXxop{2ScEtuzGFIX# zY1)2KC>>JMu4H^?oJof3gyBRo+#(DdHtkk21os92W@~L0BM8%iCw>cY{!xgH;Htg= zE&me`5?pmL{9s$zO(5-5!EuUsK-k7B;R^`gcxkW$d6+5ujT>9i0x!Zv!3m0sC_C_{ zen2KlODOOINCwFZbOm#8ih2?rDhM3Lv<6Xd7D!3pTk3kO%E5}O4D>{#1gEJfK&k?p zsoZgDI}l`!5RhuM8%RxHBE3CJwZPcLflXb2%vPu^kIXTHx*e~!0%;C>hY$>&pe{w9 z4S_pgRB)d99v|BR9}t-@(sl-_F(ttTsse0#0=e*IaG_ccqOUmqaL z6n%R`;0h}gV;hA)BwaUXPJQO$;2{X7x-4AI`1kPfgR`Jer zfjX@G;0E;`T3!si#6dTz&w#ucI7oHc)Dsx={lJ+Ny;=PM$ftoz$+kuHfy=+n@zaT0 zWlAiSGn#g9Q_sLITjk_a_$lgp$V*c>^EfNpRWf8`s~kQv5AM)&q;k5j->IrA+VWIR z9nJy4)6}uR3RKSZ_z0e^`rxBP<$TK7I78inNvl*jPjV6M6f3J#&OXFKaF>i;t#W=x z1I`qzM&%sj#GNJBVwLkD?fHdZD^<=3$SuLM1#4D0S8*xrmNl|L<+LJ;1>4mOxOpQp>URxli^q8a9W;CK!dwxO8Z>4Gh9~2 z+TWq8&uo!)p0%a_*=%P%PgwcmKT(SR^R1PV2=)AnOuK(8Oott&w>mJT zYtEbcg0uZ2;ryJO)I5N|i-|Kf^67B!+$=QNb`KndGun}pG7Wn!RXA^guupZ^>s_v%}xxGX_g;<`ya`N|))AKf$*cG<{wC#Ivlul0q&^8c4 z3;TM$;~Lxls<6K0uzu^ZihQ0QtmJ=yb8+yRixD8XOV%y7g4fCkB^SSnwu0ADeLELF zg|>p%vs0Si!&QF+;cWj4jO}02ppbtGSME)O`}tqx0zGg&;5`4AOh-3Uw*vn|jL=&M zm-sK|Qn;0HrN4|WKS&j-{2w7{2XChiHC76TxdCKze(>6E__&hx)dIUg{<3;eH;|L?Z|F7f+03m+1$ z^j}I(d_=g)-yQ(`F=M^je<#)cg!CFe|Ij!1DP=D9=QHR2gYZiKnK%YNBi!u&BjtZi zc!U21&iEIExB1snpMNsyclt*Y{*p5H_O8ch1=2VtS%PE-{z>h!6%%48FqbVkiis~T@G+zWa}{%2LEuf0g54C?dP(3u zvh@(rm4VBN^b~1TR@Nxyi(YCS)?$zBU~g4Uz3!$z`WU@dQ}ez?udNUs?5CI&)2ytw zdIRgPI-^gIg5Usk6+2%}mOQlvA6DQXy*E%X!P$W}W{pAWI(80>1#fV$`aO{Bz*Nro z5OoKJ3I(p^vKy-229jrGwb84?)MT({Md+{L>Q{hm`w))ZxWm$$_^ec(Pb2Ej()b!1 zx#N=Xw6bT_Fx?4G{jAl?r=SM zv)-dcBNR`wS&QN7;7Bz`dMscM%Hz%|*28GyqgRR*&mCEFaR`l^XVOfs9eZFb%d{grf|@JySY2&attqtmQmSja74nb&CcUz|3hSnE@AGnE&t{4ahCyuk!NiGJ_^mFQEC2N zIFEE0Bs+14I14Rlmm@)77XdK&T$qrT?&s#w&!17I2d+miD|ba*tCill2Pk%K3*nTp zfYWkO*279qB~P~hujrJX#zY;0xR)R;S*po|oKg=6)+8NnaQNS3x6^6FLS9CnHUbfD1xbVsT~!6`c$Z zQGAvc(~;3dh6;X(gq)F~mVj>ClX1*=8grg0JD`l^n50Y})y?<~oRgU>;E8Z>reDBk zAT=`};MX*_vw$mMOs4oCqc8N$Ocm;AWZ=v+0r|t~%%FguP_BeuMsF4n1D|s3nUyWeUamc}vgO*--!lM*<#0`;S)MDfD6(>o2Y*IBT%XmQ z3o9dx1(6*PFoM9y&Y-}IKJc9Z4}l>^7{n=mLbGg^ey}oLf|#7XkAl)WJ2#~>qEa-; zO^s=io2oV00y}flIH%BL3VP%Q-{e4_!I`<~GJG}GSGS=8-U!XQX9}1NSM(?ra29xb zln6K&3$Le7z*1!KUY!L@LWG5S3789GLm`54qvLewmT6_o=9(%x8A!%d1+DuC?*p*A z&p-iBL4xQzPr#F4T;KTuW?^RfE)Z}T=CS7k9~E#u2JSyvz_YQEx}*^FeBK^ezFSiUU` z`Or6CQ=ynu(1TP@cSL-;oyS~h55rM1Khvd7`j>*!^5kJY(tXDY(@ZikQW{()ISUD= zt%rI5Yd?L0jM@{UrcV@b0H!^?O2GG^clsm&n~}a(WUGL? zN!=!(9}7Klih!-qA+lY-QP3f>LqOKVh@2|muQ02T(*!&jkq|juz(?8b3;`d8n<6^} zycg>@vP;0<(~2_%oQY|VoF(94tc%Dm1pEgKjGQgtDR5L|w}74d06a&)UYOv>9syZC zJkl;;DH3+%TmeU8ZX)Lih<}c=BIgU(fVqfVAmID(dgMX@@22V(3D}oDyI8>6F&B|b z1iTGtAabdI&6InYfY(!p~;QovPMbdj3`+)0K5 z0 z{x`u7SUjXM6_F`yMA(%eJNvY4G))woA!V-?DRf!U`{ONE= z;pkUI6;^2}EdLEwO`892n2W*+ipcigg^4RXhH!WPZG^`W4*7W|E}U>P`sVr1!?AE8 z;R62zC{b9Y)jA1FtZrXJga0=ln|8|Zwm2kEH7g+y=#}cmbzet|r2ru@ZMxJWI zEB)t?XFB0#e+6ARlkf)r*O;%udDLf{pO<8X3kaX;zn}0T^6d0aq01KY1h~h4CtbXR zJQw@z5a6Zc*$+3JjFCMB{Xyzak9nftDI9t)J!3zJqof)-O6McGNzU|wg6AqUdETTu zS5Q<&Z+d6z*V6wjwmpm$S$N2Y0LJ|<9J^3NPXOs;vNG;NBozMkd%zi|U;>LyxgONM zR5fD>CVIrh53#%93|QI55^=GkLIv^7*jz_60yvS#|~mM{!A5w zf$&Jl7>wsYqn^sfY-P!=Dv?XreJ`B8Opi+aGT$X{Aa zbj+DZmtp)h7mi-oq{8=r;0gZ#|9Hc^Ka>!j0f~v>(?CxOPs5__6rPW;u)}+>7=7VE z7$iCTIR^KKOF$2VOF-!yzMFQ3*J0Sy@M0*K7EZu1819VLF5v=33{D5pJ}??_r&&@MYkIkM@8$`Wky@ZNRKWN2rTvCKIVtCjh~9{h81;TB`a0>oydJbRMt-{zT~)L7(U-RRK~L(kabh;3&3T zh-1p7=y4(da|bVLF+wZ4O|1t>_7P_SkSlj{`JR9iMf4PPA=|kIZ)H0(K|IqG#7lIC zx*u)b-f;zR*=R0-&VN&(~mGfD|9v; zc7eK{ogQ;L9l@_zcVk|nd)4#c5`8`;yQCmtXQ-!>27P#GHt;;G6iB&4kEwCbR`4$R#i{fQo?+i=dt?dtbr{nTacsJyH(E9*q0 zgXo>=J8)?;kwUbY!X8syREr&$b*#Y~RP;gB1jPd1)A~mRwgTf6kf+sYXt*JeP2?GM z7LbF1r3{oms&*jv27-)&KPhJRLxFA7`4x3JkS7AXTNeF`x*EuHfg9-Lzp5L7ycpnV zGy1W*708=`G-w_DLfr-A{lNKT`&RLS=hMJ9G|S_82*}qtXK~VP&*MNWJk}$Y;&~bv z?&e-YEYNT++79)tB-bB=arOJ)m+Cr+f|y z-!65O^x!JAZ5m|TmkP{PhP^3{S%2V?9pV{=)iJOWOh}9Fk?;(B$%b~%QdnmN_=G2V zuDH)ePhKI|b{39V=Oe*HcYBy@YW_|Bdl;f= zgBq7v(Q73~WCwm3M#J?UF3?bbkCUPY1j*AX{YJDX2n6WwheW>;or``e0xAQ4x> z6c*Izo{IF~(}a5oj1#hDl}%U%GNu>DZkuq>sfRXLo^%92S=DL6a26S^B?IRpU7~QV zgIe7yRoRUDfaWfN?n^oDg!SIpn?#hm%0|5`j7}X5`ZU04B&t)?> zNZL@Q@3NVsSF53PdD$#>s8K^xsP=5?x>yb4J<+ncY+0!+&&dgBn#VCSCd2<_E4T~H zpNNI3GMKG>UhAC1K6}cAsI5pYxM=k4) zxFM4V7+Duvk*T1oy`FD`%lX20qzDn9J~q_-x5xFJ;PhQ$MfTysLRET4g5>jt94<~r zgTvLyiZmk#)mh&8!nN&QBWolqJ=fvFb?T#?z;&^AyKuevsOBm>m0mmO-9s)AY9c-%V- za{ufdhUqS9veb(r>IZNX?=$R!t1N23_M?-x}d?I0o%U0F3;A2EFYX z^o|%bLcySqjX)aI7p!pxIX`8wihnRnG-#ACi9x}r)|v(_GJSt!P!GdIgSKcUpVt*# znja_Q=mt3)jzMb%!Jx~B!k|xFgFY35Si7>g+z6yWUvxC+4#!Ky#~CIXCUki=bvbG~!Ft^h7>AS5|Rj9N%b%@7@@n zDzl2uisPB*@RXe8xZpa|OiMO7&DoCe#lMZ~eTma-JH1O^kLxW;4t09M;@gbl_CWOK z(ihHu`D{J1a62asn4}6y)C?YKyvhIcT9xI-EtjG5m>Ppn$;#&pavhN3_YHCzkjO9{ zJ#C2}at#Xtr_L*$?Yc#Jn-OO5_?J{K=Cd#DNARXobDJ`PELb0pQfOb8+IfK@TwF#g-fW9O3$?Rj+Wit{Meu(xH;W~+2poC1vc z_0VeURXoX$m(=8ZJN_Zj-^VEv9E8vJ0^iUZ_yp+L|mNvFVCy!D-VIh zx!-u`TyOimTeWw;_0YTC_8YcFYN)~Y9_r|AKd^NIc1k}O?BLeO{Zzz~sq?f~MsDIv zD%o;*+wa^Oxy%pTW16e<+<3ijVQl#O-Vi$D%HX~6_@_d80U6ZDo-=2e}}f>8pCK$!tN2W z7$q%`c9)?GQzta=#go8xmr3Ndxr-*lF!9^-#d!vyslPHk+N;q+_Ig;nu09px4eK+8 z4Wn0VGPJj1wBj=itrz5f6Qh@2V(5cG*V70NHd-}+la#&!+_rQgZP-vMPyZG8WxGQ@djb6-KYD+$J%v9 zeTlVhT;&pF`oOD5L>l8pBJDE(5@{Z5e02yb=NgF<=jW}lNE=@4MA|hHX|LE$q+KgL zG9Xw-+I502CN7b7y+qmp;5v=na73it7>l%KYH*-Kq}_Z(q}?LY-AKDtuHj>mcF;7N zNW0CWBdsz%(k>bbAK=9~P!D)c8HwOvr1_HpFw&X`Fw$NJP)4>6k>>PoW`7-NyumM# z_7nj|S~P&JAICcWK@?R+G*r5(tZV5 zPNcnS=)zQNq%qPybeWEbv~h+>BCW&n0WcMoTugq#%Ph#)O-^i z)%rBGK!V59K2D#d7Mj)?eU@5eVkZ${Wb{$|Xe%v^YbACT$gYC#Uzt^zKoIAS7oa1bO>{?@0zGS4x zoY$C9dEG&A?Z0Mrm0Ua>rdyQ@41>&gyFobPnV7R9X8bABtY`e)R5&43ye5Cz1PfyO zzM0r5>KRX!8Gi!)&@-MUJrp|Ygc%h4SmH9{T|C_PT>@Osc)FQ!v+v808Rxz)F%9Eq znJG5=zHE=)_Z@mP1a;uj$phq5_I)1cbQG8FiFc{f(Hb|~^u!r9M|YHc-wB2+{<_H^ z^w)UJ6dxfEnP%;;-rUt*Ex&Sc>aoaw5I;KrMGeyZ`L7lgVI~MUu>TGelmu-N%1UL9t3}8I`9K3fc`k& zZa(&%F2HYA$YYjg>f)|IZc%Cd@yHR(Q-6~P{8srL*#z)KF6ajQpeW;Q-@kc6qor<> zR&V>B&5<#>hYbBI_MS=nPSst8zG3sRcXKh_B}1#Jzo$#?mZ5W@ z;?!(jFWe)-_5s%cdavS--hAFP6QCnPxSZuHTYp#yxKfN9qL$C9Au`UP&C{+z^c~5F z_K6RG{sCl9yd6MS+$_gkzl;7ej zeZz=6W|AyL*}ct5Z2!gP^Aot+N#^FX^kc(hPJ!G%TE)j3B6Gv@5Pp<}TymP$Q*wjj z7GBBlJ#dGsadCbWda6_VYJel~=`&|Exh6V7KcoKJTq( zJ9e7wFC|yL9^da9gHZ3+9lc9V*Zr)@grOSJn?%Nse!w748QcSQiVRD6q`-wLJ4HOK z;v)twD07V__ti%`_zza~%{JED_q@yH)VX>CvBQ_zd)0W{9gyqc?-}NIVo^+&%i$l$ z{}E6f<(33iDN=@xa;2U_&(nbGDEF8sf38hOxz|MbuiJE#CrPNiLw1#bMvv)#M4ZaFQ-F-GW6se3_S<8&r}x?_8#yEpD?O)GxVEnQ~Sb+Dno=f zLQbSV`)3;Ado~}N3xeuugbze+P+1~8f)p_IEehAaf_)8mG-^`onsw z60DNlDorvvsI%LVjH*&+w;nxBBANsM0)JP;4_toh4{OyJN8mB8K*??nqlp*0LsV_yDE9mXk>Oa*5I+u zRmpEn3-e&3A+pNla%+PYdk!={|`X928rX@RFYv{t1HBxt$>wlL?a%=4R zpBEo|$yEt2G5F%<@qXQdrlSj_!@qF~Rm?NI&$^PWiYKvs=j-hlt@uSl;|+jsWAxGw z4Bfdbk-L=U27pY~v1XdnJn;2swb_}hS$)i8@dqC=SyuxnBbz%}rOl?Z*#>JXp#O%Q zT0LOQ{KfQ6RX{v?pk2UOL5#y9k}I;$Jm&mhFis% z#^;oJQ3t-AF{#6?ib;m=3C-hagKu#9JOrC@`g{Z*H3{we+mgfuPqbN;(@Y1h;)-hw z@;B3aD^8nV#W578<0^qK#%Cs?PL`y-OB1ci9;R8g2BmSlV;tUQhqw56)66q>gEZsf zG5*Xg>y~?)%6h}fvvFl(+!&>2#*HDmyy(ccy77E?w6nIRxwlsZ#hIKZUb*MPyNz<@ zPPAq`g7b`ewgL7juE&|{iQ%+7lfU&<&+?Q!YZo`g&+-Qcv+>CIIsPQXn(<6m5cf#O z6{h1J0i>gR%UpNdBOO~;H(N@+aDM!3{JA{#%@oaBQ4}}l3rO?clPKc2pR8$VwY>KU zJ~!Y;+piMIKlgitllOjstrHS)P3-B#KHlF^g++^6fZq&!mP%yPXHa>bHJ^q(H-Mwz zZ~2h>_`>6H_*!z}VaSv$IlV6c-Zvf750?ki02~G|Vt1bIf|I7aCQU@^320s1AG39y zX_X{-qd}bPz8M`Vzcxcl;>$A}%#Mu^#O%lwu+Gqh={E*pcHHkW$@Isv<5Px7vg6AJ zp{c(&JtQ(99D#vNrd44L)cU{cO0+8a7}~#LwBme2V^$181EXT;7(;imqCC=OR$Oo^ z77^clFIlh@rz^h8S{D{#Ys7c(MYUjHI=<+GcaMq%PhnE{E>2rm?t>@s4*pp!cn8ks zJ2Gv-^Oz962aR2@xf31|;vIQp!7*4Td>5@-(1h9Ld&GnVk70)RjvTjOJIv!doVK8N z0OaG{_tFL3p$y-VHy3F0q1%ZIWJ5s&yPgs?FmuDs{^apVB5byHxgoVLubmRM7wfF>VZmRLFum~;Kg%da5 zW#?J_;oAiK0389Kq`J~7oPIsn`y5#Y>KeG`f4RDneJLCjS7OO$-vfeu00*bU5(!5s zvBa-SAEm?+dpaeS2sjyMEZu=|+!9OAz|td1EYV(5Vu_$rVu_Y_EV1OX z7lY=rDINM5w=Ys7I7oBHbymB(N#5MA&CP2zdz~0*A5-n%ZF|oWCRKFdFO^U9w^c0|UqWZu{v79WSrSI4&JtVPI!jzhZk;8zxOJ9@{ERwFJ#lJs>nv3R zaqBFd0mQAdbRH15&e8_lvvcb#Z2{ueSvnPnTW9GC^f`K+B}#MaER6%=)>%3kh+Aig zzp^-boh6$6GwLiIf;6|z5*^~!S^7I#+&W8N0dearJ&i%#I!iACaqBD%z{S2>XNgmS z8waM&(jPdJZk;7BdZ84Csk4*{45ctkoh6wgx6V>GwB@;VmSzAeaO*4$09Nv|>MU_| zeCgWMSt5onT^|e~4$o4=7Q1zph~Z1urp^+v=Ks3RlF#PaHFcKWga=KDr3cZ~XFq~t z7fRtfod*zaz(-tlC7-UY#6}ih5iZRdS6#`ct1FQy&yHs@oL6FFtg9>e?2gryWYD(1mVZ=}59XvFb`b`=_ca`E1Ja#HuU#>}x2BRH?2+Af1A(xavwi`>?|#)sPo&Nsw;JrY1m`cm3;P8$P>X*T}jy817?Kz z*(bFH;;Jk8Y^kn91z3C}PM!avx)L@)9EzQubUP7(d=lfTEBS1xu0-m$E-SHgo&)|< zbtS6rR99jrr@9j1|5ROxDgFPjx)O)^kJXj1{_z||{4O<>#N`5Uju{DPj;pTZvt!kj zT;BNVO1>kiE5+t8uDX)X?pR$(nB)u`S6#_xcdV`?O#i#pmF~a{I@Og3JJpp4AFa9) z>HpK!l_=Axu0+_Wu0(jhsjfuWsjfu$|E}svw5Z_!RaZg+`k$??6O)ED3B%ED@FxOQ!(%86}oja@Z-c z#4b*WCGO>%5=(@olFqNe5G%2C8`(>G=_EXlb( zR$_^5ro<9!PB|r(DlmyoiKUCUoAB8b@03^~?37qq0oZ3BhvWZUC6+kDQbAxcZp7i@ zzXBgli6xR)VreHnO^GF0_NK&=TumOW#L_(Qn-WVDZ%Ql)s7owSpeeB=pf0f_sxZi3 zMzdT;yCs%>#}2wwl8C9sI_Wl4KwTQ&nNmqY)uoaMN~t8DeFhSeDU~D){N?aaAH&|1N)oCrl_a1pl_a1p zl_a1pl_X#;7DYi1f~Hi`KQV;QX5n+m;7ofudFu_VE|ruEs?QGLC}*!ymr2eh!ilR5 z2+h!?k_7At2c}OHP?t&)P?t&)P?t&)P?t&)P?t&)P?t&)P?t*LaHdp}fVxzYfVxzY zfVxzYfc?s()RsfLcy z`G{`P#ZM>zJWXgarIIMhluCM3`rpd7hp{IqJoJJPSaPWgMf@E=I+>hONlyYcrINVM z>&w%oDV4r$0iV9~Ij;2%+NBG&LlFmgxQz}V>no>#Ufs9f~ zKAY{PRFWh^T`K8TP<{4l9Gy~0?*TTYlC&IMDv9KdrILK{rINl9Id|b`N+od&DV3Cs zhl6JKkYEEgrIJXNjZR5idmh_$sU(t3sU(7?6jC3|M#4^&r9ONYD3iidz|#;u zZXmvN7Je12Yr^G#r-%QI*1GTo;9MSF4*IfiFL161uLgZWcnIj@!l(Af*8{>|fd9Dg z>p0E~FGOo?m_KEo9iEO}$A?>SJT~mZadr519G8ZVMX%=YE^tl>N>pblT;?#M15wKh5 z>2)A(ohOzTiPd>J0Ib_RIHdD{I9P-ar_K}qy`l^0l;$vS6y#q9c4nb?CR68WJleW7x&pZ1HJ3ngpa}hlG0W6> z+Qd%hyPYm2LmZc>^RyofU6w-NogFycI!||Cf8x}6dWyUh@qx=LW4kDPL4s4~={vA> z>mTZ%tf(QqIzm>Ksq-`({N1YEP93?7g{IEad~k_Ao5;?4DcUOsEI_0|pKbVX>pa!t zy2PpTB$rA~ohSX|&D43~PSmaQ#Al>#ohKr0ohKr0ohKsy-8xVF+Z4CX6KCG7^F+kd zd3ps*Zk;F6-8xTPi*B8#Cb&IT=V=}66vvzgAYI%ce&j8A?HH%d(-jWqz2xlZY`4x6 znVdRL&p91Fc00)QyHo1cdHTlT%*Iah2+?L9O`WIjaLMA-c`5*Nx8r`Ihv{bOJk2BP zCYQCN^66B_t@Cs`xWpm%c4UvKF8ap~%sQsd6YC&3b)E>jb)I;8%dPX23B;}Q)D4JR z=czXkx6ac5Aa0$fp+MX^PlZ6-I!~oQ+&WJcK-@Y{(c&t@Ff+V0ct%>O8U7 z5FS;UI!~*Ch1@z%vb>a2=V=GpD&0CyyMR^wM4hK|&=#xnbOEr%Zk?ygfZ>s&sq;k3 zn(@_nA3ph(q?mzl>>pxKGEe$!W<=cQlOvOx&t?w7Emt|RPMgnWBEW4tEG3`KT;{X6 zq5FpM52cPMxQ1qU~QmaqB#7#~IP7^AwkIuJq`XW9mFHM2}kMi3`-N z^F+k0^F%G&I#1N^=yjeXWmry~r_ZUZsq;k55_aR@vv~}%*wg5BQtT;(UpILOMDbyz zm2evwSmLRyry^ZSJeBnlSV}yVPxu?iKKm0KyZwNJ&N|fC@_Y`wd?I(AQthd{DihlQ zpWO||ZXhh- zHlsYQ+Ee)q4w5#6)t<^twWlF;eEBSPs8K^wU~c*BWN z@&*0@kk3Ma zqXop|Y;xrcl97H5IU~-Ys)G{v=NJE^+}Wnv^D((aQs!Y6#CzakrNpiR!^iUiUuy54 zD8828&s7P8hk-fLw62x-6%frL5+kRYZd&D8S~N7r^LDTlx6%*uqyoa5Wm@n9w*FLm z-yY%&Fwe^YIN>9>Jkg^COw(Y>h>Ngq0nX+B0`QV37}av`Lg@J%z|uolocqwHY(qwwYBX5sA)C895n(>lF2fzC%MnW{aSfB`BK zE$cgXfJ>#*z34=Nmx;jbL$$!Y($;Ah62-1Kfzg}e5Rsl1lU{&5$3sTi9U^Th`j5i* ztuR*q7n3qz)-Oiz#9Z=0OS{1*8P9L(YNE6gkiBW zz1En#axLIb(pkR5KKeW}=Ru65A5DCYR!PQ>TVwgvXSLhIG}1gtcea!Dr`r1w{4;_@ zVZqS?;$^T(=YexP%me#JI(cG`w7*&m{A&;R+Xa8Vl!=PJSvc3jgghXMZ;Sr|30qy^ zaf3k;+kVkt`DjEDOG2tEJ=`iBK&C3W6{@MLJojMy0+?5F0D{!j@+0|Kz)P!>z{pTV z8n{Ww8$d3<0EaK8XiN5h1r;p5LVBn;2>M>UE3U=i5xlMPlQCf_Q#I=$vfk<$1A!G6 zg8px^o&{ha-qE`Z*Cst{jI!w-F6tThy%TOnG=lqAhNq_%h#ysWE~2zuXyekT*aF~G zN_)m|?hFeUPO2V;~x6NdH=i|?+-2ot_A zOo;V__BS`fw$e1XZy<76$qn>tj+Zn426)Rk*-h12ST%h-!24Q0)J^76t_1Ks7#0GU z$kwW^AfE!PstlldqgPcSkjFu_CY1v^2yz3-#{l_W$kPCQ3i1w;X95|) z{4!}FpaR>nE+u&>kOm<)0%`{NHj-O_TrK1+fDVHE1j##qyiD@s-GDv?x#}DML7fi9 zj5fA?4_-`IV*88ubu)rW6)?*cSYzfHhHQ9;$#FBMgCRG;GaXX#(o{I{V4`Jx3^u)3 z%e-8yMHpxN1(P)=CPQ9($Rx8;DNNUbY0PzoX|o$Q!8z|i>&n8UjRhkH(W=vEzh;h^_t|nzB9BE zeT-`ARnKBdXTHnFWJ|uA13;(8y9_&?TF6T!dGB`EIjpc>2OxUba2}B!ce;A$^e9R! znyIC^QWLyEq^qe^BH6p$xb2Fm%)M%6w8gL=A==?{GVkc6CfN8(tTDUH3Z4SN@u76o zWHRJft)rVY#8t+fvcI=r!TeC$({oop_&<; zc$>(ys+r1(Zy@6HDtWA8IoWHCYgiq13a23Pb0LNa&;g8ok zzXA4Xa*=Yp$~V!|%$1N?Q6PVU9Zh68 zeigY_5j+e>mhkWB)IRCpsRBNi0yMq|cx4|WB@xJcbV340d7RmpKA+dE+&OxX;o^_? zDyt0A9mr1I#pfO31oxPErg^V!qIs3ghBg3Q?u?_KV(1aj|KNz~<%lZ0DsChw7B>p~ z5gY)YF68u{F6OC1Y~W_n{L{t!nYC+i;iB7St~$v9>D@D=Fn~4dfjsn{D9ARmy!gj% ze5t*GQO27SMch|VafEntqT-Ut$L=J8-?mPa!M_B?JzSLty9klTi!Y7c5A2w!(AMDI z_7B@e$;NS#;(Eu+)>Yl0Nuc|Fm{GD z-jl2(Ys>?N=Ok&z<%;ubLa@;eYr@ILYc;sN*kU^}j>2h`;#WjgFYKyEq1@qe#Y^mk zP_CRd-8EZFVbwErqWZ?M=X#n!dtP?z`5wUN7YsY?vBr!uSPEoQ4Xn5c9~zp$t?_9_ zPwYJ8brKa>iwSHK?+(Q&w=&!COoiBaB1~KrEj4trhxK`@&ofeRQ`dNC4*Hy96j%>g z8;t@>K)=P6h2$)b9sN(!yc*4wFB$oqi8l=Imb4BtQR$tdholbQYaXBXPA61G=bL6p zkF)=n2RDAls}^+umhth7itm^c=Ew^E?f_)fV7dOgfdE#j2C&uym@BJQ!0P4bUjsM& zKxtO>dZ5#Rc9CUoRi92A^Ep`KGmg8ES*}Q9)(NYaia&Q24Q5Xz>K*I20+X%T_8H;L9`K+Yt>XMf#MfFNRIgjooomezSQ1#(KwHjG@p~S7yc;x zvtinR1%gG)q-g#9>rrg@_DwWjJ12&o>8TdTY$s*dBW$ZsK~mNKg;v++2rv(rqU z9gsN|+>jf)J!y?T*EDlYT8nt2PK7ztAwy@8Y`)aqHYlfu3UjC@A$&3a&02+jfJ3vY zF*5u!fc~i_LN)TJ#-+rtd_GYyQ@nWk|I{PmHc zb9z>3;gVctyK3>i6?`>f0Bo0k`%E(fn(ZHa#mnl&{LeC#L8D)QU5h0OGY+6@2^i5|_&kfK z%5tl25F~rbIG~sJR*sh%5M@*uU?Z9+PIoCYUE<0{TGsz&!v~KLo(r{^I(2 zp6gx$4`1W-hjlV;%2cgF!+;YIIGw=Hmt6mw^kGfnNpU|0`D{F-oeYH1%6X1ULnpr8 zM5`LA_JX0}70@SxK9R(#_kiu52jCw7Cd~)%4S>Rh0K7;$PZJ0LSicBB7J!!lH1+|o z8O(f58VIB^fX9~5xoQxQnp$KL2GU^+!&jRAu+*O>j)ex?Dc%Vl*;&0vfICHN3EoXS z;r?!(=HCOM>hLG*G7sdwV5q1D&_G~3fVTilBv;jPU=x=DXaev$feipW%b?l&^R%=_ zF=>xM7M2w@pfY`X)+fP}KJQn~3XHBXMv>WN^lYWWeO|RkOLFw|FVOfQ>iN)%#~3{q zVpP%ddIHolw+`BoOZ2P;CVDyCb zasyTz3jzOd96DjUVbx<`f|)Lzr4vG!slTs<1Ua+ZPJ?)6c^yD>ubGODXO<1;8E;rk z^wIGPHCjzcH$F;9TwR$f0MJKM2+&830HUi5x4sTBH;v_PffSU%4%}$_E)09k7}f@Z z2O9c6K>w$quLu2mLuZkj?D27q=(OHyvAXwFR;AA@KW?P>jDn76LDy@Hd8Ajx=%e2? z^iDX*SB^IdQ=ip_3Hn4CmV9vnv_Pz9$N4WC=i}sym5W#9ohEJYl9@Nm(8cGS<;+oK zSEFJV4Ec&~PG~>K*wEJ>PZ^bX?vACs`qM{QR*O-jfaBm*H39E4(9jNlSDtQKyXE1F zHfTX-?02G<8FIWl$!HKzS!NO-FK8A^AMw|iyA6G}p>36my{$-sSZ~e_pZg8uJ}e$% zaSxaz);bT%X%lY`O*hHg8lA96G`E5i@wtX&diWkRcFnyw6Z1W2#F|*up0yA62NmxM zi_lu0uKz43p*69hKP^t}oz{i><}a_UZ>U=qYHbQFU0YY%S{JHsY^`f-&0AW#bY)%J z%KBw>1KXOKTkD$|Ya2qiHm}_{bZBVxx|Y_^%G&jHAzXmeH=_Tt`qg!fE#PViEm_wZ zTG7-R%CpvXEXuQXdKk$D?+Sbi6bf>0bCD=WsCDqZfiY$38|z!^Ya8k}byR4eR>5lR zyRUZ$iqtnYV>D4`2{}Vejdh`xx)o5S)mj%hQg#^CGj9hav zCL8T5DZ@$mBaRMzwa?mo_BS z2p88jG&C*6q%E`7J&lFTcjY#X%Yc_zw zcb8$-4f= zHDOF9(%s^VL4K}5uU>BH`BHnQ!obA$lgU8v@fDVpgsYmx?6_%uvKmBUSs$vqRSz{-fhu8nws6&h^oY!9#^gUx(RS)EsoFj&QnomY-D@R!GQ zURiA5`?5N(EHUuOffQU2iZb5zS5LNa%R!J+SHV!+b`a$FH9#5+GIbMlu6Ff&!fwW8 zzOdCHr=~uoxT{I+kDQ9E z^wSMbB36azXQd7Wb4CVwS*bfXz$Hy0Uj~?L1mP?Y{H39DfJT?zae&(mk23)JS*g2Y z1H8ilcEEDJXu(acAS-yjp>u#|U3$j>_8OiJ1AG!2V7wOtjK|5+86W@xO{E3L8#)Kz zmq;|d;{Xc{Plo}j!E9V~9S2x)Qfz?Lt{^MeV(1*8-KBRN;AF$oVSw9W19VEj0HKp( z13Y5{L9yU)L+1dWx%7?$lo}pq02-9q1yYR#D>=Z^*jVw!ImidWtcKv*hRy-_S9qG< zae%)Y9%lgTFRav-*Z_w)fOU$2#OraW*A9Z6i~8s?f~sIA14s(w}PYnLa0jdI(P6rZcFA za?S#hC^eD=eT1l#NcP+a)_)W2i=mX~0ieZjTN29<^Nymm23R($4h!aOw~}(C$gqcB zk@mLRNs?u60+#?^IzR_XZcpYU|GVfiY&)(j4l(0)re?E%cDt!F(BQvB;GJ7-Bt2=n z2W_3|YjNjM(A=REBH&6mJa?pvo(r`52`N5Tz$chapK}@V=8`vf;Oxyd}U(`k}R&#@76STzv&E&7 zx%)){xS|XEVRio_jw)Ti7M%4p*s1j-{Ar~ybi$U*;cX;?PByMDM{Zih^3= zh5U%*3J4nk3~7FEK@5LUW-Jk6&QHMkF=4Y6>1PF3uh3(Wb)II$l@4ahWlU-TBmUDn zIN!K`0$5vIR;$-Zadc-SCwnk#b9oNp)~)ZkG={6Ded5v>kD4~rr7;RM zZ6s-)1kVkaK*>fK$xai0uxf=8krGuub?-1-Z2;f11P8nll>StE$N2H(WNgif6x$O@ z@2o2CjSZSXZ#;2(d7mK+sD9F3x)WPtSq7s8@3+QMBLrHHmfmZz<l|kMVMl-5OV>T|M;iS$ILIX!rQb#e z%ie`0pPB&+(kCy~2F=C5TnGCjS(aFSYkX#7!7RtQV7`7=D1{#r29)(h9sk9*^&kR-bqAVWi? zV-{;)Z9NCBaExK)+H7)Jk#Y_F$axx0E$NH2IxVPNe?U{Y z1fix4P*5}toozIY6d7)&j!2n7m~bQG4Z>4Zk7)*>kFCPL4#D3(;9XdVMIrAkgYhTP zco#Xhr#u#+qV=FP;O(&pjd%>5_*PSj`lD+%yge47qK81e2XBx3@Sx~afSUv*H1cX3 z_*Ort!oOp(`7Y{;egTm{sOUVj@50+-5h|>~fp077!3h@CU9iPfU&B8so-f!5H)sy+ z&@O5>o_f4H==ESxgB}7)<*Ay^#m1M~+km6@abY0%$Fm`MisAZPfs$x$@tlb6D)btT zA=8W-Jn~X|U&W{QQVI!CNblXarogHdq!EZ$7H)`a_XF*eZzOnNhb>R;XW8!f054NX8y#pa2J|1={HoBO!FOfy$_G;z1mpS&B^qUF|(p7%2Nv z2upa<2+M~s%g1CTd21LUAJ3R_@(dI5ai`H7g8c4;d<+`NC0z*lxM`F9n@tE;U#fpf zy5F19EhVVBbV|#@8tH0zQ%3ot{z>Ujq@U&Ym!lm@ge_7M4FwzWq!B{tlnNfq9zV$8 zD;9Ecy7H#@{gk9Y^>=o~mVi58(8NjUL_0f!Vg)NTDV;r1dqYk>8q?Czb#x9gPmz{^ zP8B&Mbjef$ax`~(qKukhSfGh!$?Od?yM}SdGGm}OE^S4-L}-`Dv2xP?Wo* z52#%?4nwx?>B`SOVo?tuG*1Nd%yO(1w3iVeW8uM(cvO(yrk9B511;>BkNRd>GzXgj z0yxhksB%A@>I@mqgQV`zXg~sn>)3|BHb5JcALu$>R2dXgg|-YfCP>?mc!zqjEo-Pz zm4-vqFw@td`7!5`xgdt(=Xf$(=o5*_BvauqRgt63ACSW}R+)X_SaD_%80loVg+--C z_^}wfbfTL_4zpr?CZ%(xi+62TXdf;sqvXu(%JO8lkA?hFrvcwk^Q+le>>t&%3&!RQ z2R6C~y3IgSjeo9KT}D29m#)YfE+bcHcJ{0UXmJ?aDc7J8GIkgoPGvbk+Kh{|8OJ2L zvr99^X<_G$&8EEZKkjC%pAcs~%$S&Hl&KP}K$#S$Hs)3BrwhjNaerpCD%L*s98 zoOQ-EGK~VrO1trm-guO{RU6&>F_VmL(~U`(cA7LJo-rf~cP2PYfZ$$Ej2+(+Gj^db z9(ekZ%;v=VKzFZ+8ILcZ#niT}6Wrw@YK8novR3B7A=dn1#vpeIE#N{Kozs{icEjd{ z9SE>`QH+2mbTNbx)rhS2Kbiq`;`2F)@c(5%9o^U`{2ye%BlBD=2RaEa@5nUy-^qX- z6CU%QllYnc%+h2aU-kD z8Ceh~(Fx)tIzf==bkcIpisp1OFD8q#Ha`x9r3;dEikDP(gaG68Mec;eC3YixaZEP~ zuQfh^LPVG7Wr0%{z7Z8SSCa9U=?RiWOPlLD4&x;1<>+gzh$U(>FrJ1EPv=RgKK}f5 zqMmPBF)0 z;~Yac3%(+3WQ%eLCZok{#pyeL910}^v3}ig%eHeIiq@1EN)uD{s!xL$E2!P9_flQr z(oee6fJ+)iM5ak1sagLx6q<#N zRi@_vW0gVkVuCo#z*t}0r-`-E<%2tFtMk&3m~`y6W3q9s`BBj1JGmDdZJIxv97n9P zh<2ShmlVbpGVU3~5NEeDq9e)dwA~HVk5)btpR9b>6l|J0oP4~-1d~D%Vl(MnC?e*L zT>W}dX_6hc`g|Rglik^VAHvR(XO!yOl+LJwp|CD3yj4E_(~j?!4Raz>PRegtyL4DX zedD?fLpPL+Tr{$1Xk=KimA|fKZGL^@(uQ@*>hcX;Gp%S`m%pUGwZ+O`*|fSYe^dR6 z+Wgwa+NDj6%W4br>l)X$(UCk`m(@2#hONY$Sc?{|Ue|yxWbyl7 zi)x$eLG=m;zW)U{0oj#`Yi4+&zDa`SSJ9*VZnrTY`|&XiG!=(vE%W@SV4~*+qnQ@Bx-}eDg_lwX4?FEsrNPZ^U|QTFb?X zb=17Er4>=25!lwUt~s{IbU>|L*QlG9)~&BKiv-KXU5Z>FjNq`ZNB-qsf;f3>CEo4noY-?Q$DGra!o-sGCO-bQ=5ceN)^dHYZHUbxvF{WDpA z*Jh9Y0*{7Ol07Wx;0ecjdoB0Y*`XctZxP~hIeIhp@zc9{vV5f1c$?mlhR>4D^(Nbs z=iYpacTKC^f9_nN+GO*5mSnFmOmBSdJ>a#sd(-xLn_rmxFR!oFyJp|{TTkC{!Xod6 zVP`&b#VJ$0A#KKeZT4qBw0rk?_e|Mkr(bc#gWl)bJcHDfIo=Dsu&|)jyL-zXZ}vOh zl`nXcTfOUBz1j9KZ(~i<1=qn1-V=wNId_LQu>CpvaI|do>}k=t=gy7JO+fiDsQ+x+ z9z4({beq>RKzRqgL(iGU%i7>PH0}4=t-G%9cD8HmWhcyymU;cd>@w=H%6{U{-l2A@ z_ey1V@uod-bW^;4I=RMM-s&B@$KK!_T?4)U?%g-V+gn?Ayr-{%vAgWvSDbOW_py^b zVFgcnFG=9UeQ@$;_Nm&P-WGcg-RYT{ZvXys??VZ8$ea4kPY=4|$Uz@S_^Cm^*klI% zy9|2qPY!xsA_jf*-vq3O)@EMzvidoq~-Nk~x8Cy)R> z_uRNVw}2o#pBR*f3lea7_*}pp6$BSl0t!Cgli&ALop0Scx4T39J%9ht=Rc6HZ=HMU z)TvWdr%u(qeY-t*?uq_~wzf;#pOv=f{GWKcJNzGew(%d_nw*hHy#SMbU~9MktOqN; zd+Q~;cELu!@%A--VFx;E$Gk@EkDp%4&iZgNIX&S&z9qF?Z1t`LZ1uQ5k~}@(zr4j+ zOiTS?OA5vO>l3M`#D4z!6Uns++mTnUxHPp*di2UwscTp4>A>_!zHw@5>J9#+eKIIU z@R6+QuKJK!t-?r30qPquk8WO&A_0AoI-W`Pn9jP<@Blh_t`~4&4cf8ObiYBUNn9l&8 zzL{r!>iL1EZ&{GKGW8#o6Mr&h;jmQs9l!le$um#x@=r*fnedmT_F-u6_Ip$NlCMET zRzQT8^k!p0PyD}%(Iocy|5-5g|0u4EEea)3nT--MvUK!y_vIvPLac`P$wp;ZVYt|U zTpE=jYGYI;0a~S~4A)##23$ER!z!jx8Q3%y!wRbiG`F}3=lG}$l<`p+2w>N!44E)0 zLna=TffS9(kR2D5AwMoELwGn4B1#z2Fkdo4Egaxe_OvpAWGyIIokQo=1Vb2hY%0MZz z*8$bVY8c05Vzu$nnH?yYj$j@h!yvIS9?pA`{63I>AdxsO^`iT0ZO@*V!9cH_3& ztG*%WGx(4tXOF*WfAZ{;7muCe-*s(j#%T=;ESFsFuDWx<>-@h>!~&e!>2I8m#Jc|N zh9`CTKfhPN?|v)H{>J?yclw9#2cA%6KmB62Ke@QepM(X|vV8x}Z%HI}YIpMWQy;#5 z!Mq*+nv%FJS(ThU1*xB!QZ@SKUVqX~WD-@PaI(rT?DnVZdpgGE@7h`vl7GOPFX`Nm z;?u86uAP$HnV3?QTsIZhp57!Y^1mF0MI5;wZR$2_pznY5uuIIM?fXC8l-j4ne>qI* z^XIuMHV6>zqN5#f5<=P~jwWjR3zYdsQn?WPj~rH9m%4uqRCr=jYCAOg#9_q|F;iuuch}MCJpsdJLI!LK9InoiyBXHqlp*|poBzT2 z@MDe-L~4_EODu78`!2bVdW}Eh4!^KHHRd0++rI*_%^h23_`mnMQzx0u+YN0#R+YQ} zF9jXj?a%J^AM_x8zrU>A|LWEIR|O9snqPcZYNKD=?%#NA>Zz$$_@ijJ{|yga;(u;y z>Zzb(onMnG-g^Iz$48dGcb32E`_txbzi8&n)JSTksF18mRad3%+|}WacKeg&U!Izo z`b+9|GbH_i)b7VqJ9Z#u9Bucj=VO<1$%2u2{*j2IMR_|qWs{Pu`iTFv38~%Ce=6eD zf%*Q_OSt{;PulI-cyB|DjCvFdkoH{p8J!`+g)h94yug})l2`q_Q0j4o!_sz{}KOjgy0?i z16#$l6ob@Eb46ldUp^+8ZQSAFW{%oN^B4!QuPSAw+);&yxEj^|U^9k^6!$+lF(oCE zRR^l>ejKfQ?`E;YGZWhl7UiB1(UDeD$T?Ev99bo;J2EF9GmSxy^mmF-A^3k3 zVFvBo3>B4TZ7h7E>duKp{~<3WBj?Df_A37fY_U(?J+E<1FfVn(P<>zOO8xhl|3qT9|HPJdq-I+TckjPvt871i0rNvzqUsLwX~Oz>Fy`k>eLE&5XC6uHgF%3d*zhaT`M@&_7E62p?S_O8Zr<*{9-A%Glq&yE6X5$(kD!+19)FSy z6VqS*UiQ^KOa}k2n^6sT4)wHMq55fiwk?5s*{{KlX6j`SJ1$up)@rDi72JT6&# zA_{)_S~2iflKWD-QX@U7skhpxYR?|_d^vU2y|VBm{GR0O#01~K`>wV9l<@Z_F{PAd z=Yr6-JDTm|g#WE%yNa#)@gwsZAJ{Q{L63iV$_#f5)>NZe1hIIXHKhuB{TXl?{sx>p zq;~jc?Z>`pr~k@??&I>k{u@#=ruvVK9Y|fcEVZ#iP6$q#;lJ?`zdl))NcJ6Ug-eEB)i{NX^{4 z;1TQ+Q&(PiLF1)^@ISE2V_Tj=t3R7K(6I7DZ@=q;xA_82tUgFQ#y=H-b|JLp68s60pU2fef&1vF*RlCOW1`c6WD5g0(+IDCnsMe+s^U_*)zH?sFAa=@o||6 zz|IspPR+iZO!!l)v5-vhrxufQrzR&sI3hGY7VJOgC;Vd%1kZfY|H#%m{1dAEzDxWQ z#4RPN_E&u~NS!usNAj$x4UO~ezcO{&=Xcy+=((k_qie^m=b zKV6!-Fu82nWt}VhH(->t_oNmeQadNnaeqT<=DJiza@Fz4zNu+{bNl7~Xwjdv-5lq8 z{_*ou>z3Y?Tz_1$Z*uCzWX4^S!^X2|YID6~a@1L~Y z@7?!<)N!evy_a5&Q8wGZxZ8hA`^p6kZ-a$#Mg>RH*zh)(BIysMswLef9~a7paX@x!OD$;V@HZFz zqi~p$T%SnZlAM~mb8oOFxngRr@q--=h}rfg*G?^9$Q5=JQm5_3`A=a-L-L%d^X{Ct zBZslIcjSA?H_632@oz!ub*Y&Rcg{=Al*{ulV9q#p$JNPmPQ-}85;emgk>o(?=|XBp z!>$HMT?LbO>}}k;%isKzf7Cvx{jo%1UTUOa$8}HRUuvrTXaqMZ`=&o~iGM|MLt@_3 z^L8}sXy}-ioIW+zk$ls!jWGAq^O9#{X84;^Q&RmoQ5hYVoH@1dLPLP_r1{D9(~|j< zlJ86$d;cTb7F?Hl`ntx(Z9R>vFMM-qYGc0Vx{dQ1WJGpG?BdJp{-wH-bzOiH3oiJ|x@f~{$!F3xO*34{d zT!saC*ZtQ)ShD6=ImAsaoi?02d)k83TT?S<)_i_P{=&uuZ+`QI%VgBO2}hY|pTISK z=D_veO|~5CPfN~EEPHhI@;5Kvmb@u(tX%Y`Vkv%1EOAU~XR>)}atd6?JK`nO5Afkh6}if_a-l0nE9s%%Z-50vIx8 z=PzdoUM6xvx4xt(S&YgKyT;%p(O7B&{OO4k<(Ih@W74XWxR`5chWsBL2cE~~;*w|KBJsSqQ@XCyAL{oQH*!50xc5s8YGTo6` zzG~&tk`C9tLxZ#-nGEi%_vMEMy-wR|+&mk}VZh|G13B%}mN8ViX-zgih})ewvf-U* zSEm_0mL@xA(1L3ZZUUgX+p3z$g4<}4GnEDCr_QQcsAeJivf)H zsA$Cvt5`24f*ESv(xvM$atb(h8tU%O6v3Kl&8|nG&J9cYvprfAB$myr!s}QD&>V9) z)37QdciA(;{UZ>9_knCc5!3TlIq9&LWviN7Wx&f+7&1cu^>*dV`$No#!%gc7W|T<< zof=>lYwp>2;KIwK$1cXXX1#eJAU7~v#1JzU@iKDeY>yH>bv%nO@D|}sl^CHocOA-S zy81Ta#WLEo73R*iXPRBF2+g;-Qq`@nP+Awunpv9uTJRT;% z-FJ{}gBK1BkF|_uhgZtD&}MKDV?O2cnKbU-o4MSFiM2_3Ip5dWD~HV4Jmvu8FO(i` z&S7rIYjc*PN2Tw@(Y4^QhXyEUW~1|VI9jKFls-51!3-0}Wv$a9_p!w;`P@KOs?B>F z`|t#dj^-8Sz}=1s;}T?~TX!JpwB|>07}F#0&TQopSS)QVX4)Asxp^97xPQ!pN`1v* zc0CGZR_F5ThYC6@yYoYv@tPGh+%-B=+&m~#Uq)nG%i`umiFga?Sa z3Z`|f4J#HeTD-iqUaG-pz&t?rN#9_ZHKpY_nw3L?oAL01u`M2Sk+(AE`i-9wt6^?> z*}~XhX9g>Vb>cV($8!=Hd+Vg&_4M_dQ*zfx6vHGO-!yAGb}5_=M%I9-5+iIdn;+{A zhj2aSTrZ^eb+VMWM7yX#YO6e&aj{s@cuLz37n`S0$PYJB7{+9=u1Ur#euoK8O-4!C z2aY`ZexnT-uO5E-Dq}`w9>0PFA5iA`q>Nb@3cP_JH%(-j#OhGQkS~->WE(&9Nay7? zdgy&{nI2lUxIWz?PEl;8Awqg=A?_L2#uY}^o2NCXDaVLA+Q$ogYG!29IdWG%#!fl9Ba@EBc0uruv1;8`3LQEn2M3E_7I)%x!VvFdC+1 zWM)h7WE>!BuwJOeyc((zrZfF%`=f1eJ%)VOa=aI+-}VYxXl^ETccY;wEVBo3pZJ%d z{Nhc z7wz*zwooS~urc|?9j%K!IQ!u->2n#7a)S~7W8N7T3IAdypAosTVWD)7+RpTf^|=*R!)`FQi6&n;R79lFpOY9g*JZn=yDyK%nHYC{pPP)( zrUxq>blfO!_rmI6gS?FKp7J4)Ir?EFH;B5-LIMIEsS95V(j4=FCQ?0?!e|mt02^fb z5SXF7H7>fP-&j`epTYM?#EliD2fcs=(5}sHS=MnF?F<>&8ZPHeoidor{kKdJ56Ve= z3O6`tU7s)_D$PkXnEpjmz}yFhQB~tLcGD9H%MAlcI`PZ0h2EAN-d|V5tbz0)EMk$x zN}~J1uq?>kW+b$%iGDXCUnapq&dowIo{np*2+cIq{8&|Gcvof=I z<)ViAl_p3n6L6b`^}?R8q{tS}p}j8N&Go!!xDq4LYSoE&s4$wvTX@71_GAZat!laS zm~8rNQ_(FN{6$n=h`aMN)Q(MvSx$32(hR{VA?ktbE_8cYf-14Lqxa{Z`(sOA#)2h!py?X6Z%OL z#V^XF($VDl#Iu7>Mgy2T%^HRq{McZ3O3Ptr(~sd&h#n31LgD__Y`|;~P{5qT?#}qj zLs>UvXFfE5d51uxANotwvsvo4nq$e@r0bRVPi%^ui!Y-{S>^P=Kn8~N9fzxp4vlW(guR(~|@J*b`)X3gT1;bE9V9tZeN_FK)3@OP+e{vO7m;CQp{y#Yk<8PPWXE z`Hm$rpWFkvo+4+M%1;$hA7uFKRSfI{G$&&ZQ8NWtQkHgROxLT6uJX0*4j zm@DA9a76M2ObM?j`#g70A^Qf-bo0=rL?*JX8}}h0J49L6X}p?QWpd(z(ZM~sc*yWw zxMcaN`c}{v_v19AfUUoo6;MrnCQ=M0;WR=={4YbPfPJ;m}Ph_BGP0SbWV%f|1l_8e&>+hLgBN zMA>@EHMIUp8KHJA$?(lI6k0~s7iA}x!AKDg)STOhQ%xK?n&A!k*n{AaeLHDv*d=Eq zQu9E52nEE+M2-y*g7uhirF_u*XEs^LT0bMy`& zBJCaOH}iIleV`m0e0WTwFFQjcD$~_gMRN*+o-?OL=n3O;?C@A^P88ZC<}(G%zQ>qP zW@EJXW$kAK-5`evRuCot#NkOKiFrp{?5N7<2obkeqqev+jmxs@`^~`t z4rd^&s~_8!LF<%lJ2Y_D2FOacT8{G69m(909he+Y*ux^}9>ZZ(m73D_OEe3yQIn7u z9D@Rt48_|T*@e0xk}u*o!;QUCJ0SXzbv~#bl3>A#SR~<*|2;c0-YT0GqIW& zpP1>h$AlR3Svj^b3x^$EQp-ZbLl@%*9ptINYAX>;(lBh}wOO|5;kc0*UXarwoK}EY zPFwmqqiYfL5EB=ofrO2hIXq*(k~Z44!*tGuQL|Bo_cvGFmbRfFr#@o9h8{!)xyAOpR(2?6nFu!-ww(p@yr|tk zH>@;!XEw)My0Wc$VQp>o`7`FssHvVlP#hV|EilVNUnjD17P1q zdOcD%Zj=s)r;hVfnj2ReOnO|;cb2%{Ms zxkV34?RuM$UNHWvzPY)+BU7Jlka}f)%hpx)BQkRB-6eZpN!JD)jSYiE`YW+Qt=+!N zvI&lwg*eAXgK#tYT{U-yKwJ*G(001%WyJb8M@CGop%m@2Xvm%!aJPc>$M{C9GZ+=J zkGGx{`f#}y8$@P>^JiQIBvcmfnRcrIf_UEa;+6nDC9<-K|BwaS_)8Oaz*zK$9&YZS zYY+k8*a)}Jdf*yzc^neR*TJGekL-TnSC?ZJz@$BW)+3AE5mA^&eQ`(-36kucA@RiF zX4!CaCC6C6cGl`C$ET)mAi$mUich(C)vD#_4-INF*aq31jU20j-jyB0UAA7c`HlwR zQgrpIhDByWI*Opk?sTMn^c?Ii+bFqfgHdOfM~%SDNY+-`MGNs8PHB3KKSOauxB18! zOVt0$xu9Bb=0L8iFKfT|xv0Lieo=G7>czD)Hw?jvWM|?VU+57W5M~kao9NU;NunC` z!A?aEveh5z4Gj(tIm2>YE4oBYKJHeBn4nF3CR>QOaL_alLpFZeW0o*e!W>#j??WY> zlB+PWAbQk}#j*jo#8o&bMgbjy1G=WcV#Y&??1jtfn?)p|LfjF?K-c&PeV(`XA2B4% z;jV3A`I;~|F%p;5$*$3Si%;w|uzIMkEA39A?b(ObS88r*YgpK5cIA?5u5V~D=b+fo z<`AjbiL8xP|6|KSdo$eZvTfGxaZx^Q1lXG?2y=%qaWcpvFzFu|7{qcJ5w2c>^V5VX zwvmw;6W&F^p1s@d*^A~jRHm;h;?cYef)K1q=5CV4BQnRZ7r~CZZ?ieTGsDLOBT*-g z+VCAeD2V$u_P9-N8V=*jt!e={kM%>@d>5P&Y(0XFoFLq=5 zE^}RuxwKVssM3cYiqJ*XtR97t&Q9GP@}?;oF>6CKp0caba@=Jd8McE0Gn8jB8HDI^ z6IC1vCck|IAkq>Y))G2mO2rqH@R&txh)p9Z9~X#o8-}57zrkyVjQG^?XS-Ol@`lz{gp=W6J$wj8jqy=pOrHH{h=lk6J!wuZFP5V*cq)6!!6ma>qC~85vT9k|)7i>>7@e4A_FzjW5vu80p5)_O z(kO?{vumYHIq0gs8>P&F1UKLqQ}Q(%z3CwhUsk_-2|`pX`PJcoUouj}0TyG;%tm8& zPH5Z0Rh{tZd7QwaW~^j=o#u{BWRANY0u2X#SO`6w2jDyczNRpWgA~2#W6T?drTWMS ze6~D2Va6LK3`}FE1G+>!*_Cx>AYmC#FK=y@uY+~sZj#-DOHWG(X`Zs-H-y|lMr@2o z=`f5Mz|$$PAL_#jVU`hdBQaz^oHaNu&_Y}jDbI`EP`3{e1~|Gb>Prhx2MGyDi-4Q-I3|mk0BSMZ>q?jt?7M~G`U2F!*!1y`bXONA!#KxKWf)0Ic zjIiYGVVev#%5it*NeJyA)yUz435+D7k@!I7R0EFmDO{e&5k?7TMFf6@$DaEM(Vmu3 zG9pcLU?zEUcg`GodtJEeYD7pdE!U#Klf9QGvP28HiqQp^Gl=08*3)L5=V3kKddxpl zkD5cGFa)dihNwoZ>Ol~v%RG?#HFiW4JPG(S*-iG*jD=w{(DC-3sMUXkRe>s*^;~pC zpHi7Nd*f8H@CQq z$%A^G@^Ftl?19byrs%kK;}s`iU_>()nUjAr)&L`84NIl%v2h`_RBsH2TbMMYn}7KA zPdr_xPwCn9PIzc7IVVzVynhBW*w`vuy(p5#Gq)_nwP7gKI$q(DmSXEo*?8Qb|B!AO zXaP)r1kJcfAoJNct^T2&K1exRo&{Mj`y3B4t7mP<C}&==G_GYs=E);Dqv*}z0lba)$+Dht))?x6Dcr#U?&NfP?pPh}+!+bo zCif}xa_qxn2JlQyu$U4f{WI!7fKS7nAGBCRK5ykcYNj5(1mGmqqSGuJg8vvT0o zkEzY1iWV@3)&F13u+xr{r~L6sQB9LxO<*bE=RZxc_*P95UYf8P?_86=gtzu=g}>K# zoP+d&!h5GVE(^lXFKel@Aw2 zX@rkDQb{RZ`j%r!6Q1mpE-A%Jr#hB2;gyb5Qi_)b9ZQ;UpCgr&;-v>1OPcVvj?_YU zz>(VG271R-OL@PhFwjEn%}ZRbrmt4GJtV~ieOL$2ec>P?oaAJ$<{)}YNV-C0ZgjF& zSCaiDl`Vsq_~EcV!mz#;))zLfg@g7uS6|CX3NH*vG{R*esZ`Mx+U2&8>pcoT9g^bB zVIE&dG=B-NjqpoKn&7>VXe*hhdeOs+B+XyKYb9)Oq&C8+tk#OM@Zw38)kH`E5KkDD z)l^ZI`%wa~neYZDzJ(CKgkb)pjm-gbUkDiqqm?$7RC@O$TjJj|g%t@m9vRk9YDE^@ z=L)tGo}=c3EW)U)QY%u{SDmaTLJEKt38S(~tw>p4ce0uZKk4*ZOZX`zSt~Nf8sMDO zo-R3G3yf-tHv=!zE(u;IuX)Aa8cn0udX#pR8Z~;&mI^;vN0-yLB`z9Hni6-bWRSkB z@P{F3mDbme!?S-mbg1mNdpvSt*-gsiEo?Emg(YfM{dclI5yn^ctmqvCK8L z?f5bVtX3MeipqR7-C&ToO8-dX#!caSyVB0Kw1FYU8*hkg^{>?;hRDoIuAAEkX&u8C zH;W<0?eQdD#Up>e(iD@H@MAiK%DYh({EaKvO4xyyP01f*5k_T|IsnRg(aCBeqyRVo z!ljtyMRr*I7H(s5Ygz*Xu zG2VDXWUEio78)WmZ*YCuMo8-zzPMQoF>a4)yq!+|&eRl>mhb}_P?dK8EO@Rf*h=`% zj?_jNl~w8hC~K}O*hEMHZ~%l+S)~quvTB{IX2S0|QVZeFm1G?NbMMxc%Sccdt+X_b zV!#bpf4NAMA%r(mreo{`A@+PrfgjX#h59>m=UQ|+vA;wESriAZmioPb? zm9Vt(ea(VDbOl=pPtt;rMHrP;+SioztdrG5NCD_;!leQq~uptY*UBIK9>qKBpvWMdnyTnI-c`Es)hzHlD-}IQi=+X`)-W*AdQC zl9VL=kC4rpF8pML&5;U;yIh57!kZkag*yBxBvFSWBRUWt=X6LDhHABh3Ry(8W<=B? z-tTC00xH3plir+wm^F(fh*_oZ#p^sx<guFK?udReRN+?swHr*86;C6m(gc&{UDet9N@Mc%A zm2k~bwkhyZgi%?gUW&5b=VUbzQUJUZVN_PBm!hl>Ia$qwzjAu5BmAwBte0YrH6*-s z@lKhf1@6;SW(ted$#0h7PiVwHQW_i{;XTKMY{cJD8rTRMCx>jr^5ehq2R6blD@iIP z#_K=Hc7ejMNhOCB2W!#^ktPwp#Wg8S$fjZ15uOlfCGnZAm1)8yjZatPl@9=n0?q{ zDl|Jz!p~_38Q#B-3`k<%jlSW@#3#Bz+Dz4s)=?-;Emfm*y*QW}Pe;@sKEhR*CN#g) zjD{h?{lqop^q@|hEQXQFEpRg^9U8}?gX#8UL^opAo%Cv5SFvuf2Ql>!KB~c4E3~*C zuQ^CpeLtcH@sE^-)sm1smh?tF;xp8gFvp&2n7Pn(^Qg!6Q156PPILV{>Xo)(G~NX7 zU~RZ7(gxzshixFFItedB_^i{mCLsP@I6??%aFmSe3Q}CxgQ@z_h^oX?Cfe;(ModOb zEe^HY-+Yj6zb>K&G4)7zmlLw;a%zCB)EJk?jd3s~Mk7iPzb))HLMnrpgs_~e{#w{? zgtQn+ma1A3MIB7l{)no?R7S!+;&O^n3*n2aQ7&Ypx2IP}GjDZAp9$j9y!#tzqk; z*Vra5?l9-v+6ZYq%NVztrIzVqs1^3E zMchYpgfn{6jrLMSUPk5EGrgKPsIvc>5PRhwtXagj2F#u4^k^Yu&6bh5%R@%OEfIxV z1LkfG83|t+$H?4qu4OHRY^<%?S}6QvdmnaOmtUye7g&ogS+RS+j_5a(cYz zZqcqLJXvpELLA{$&KRo+M;)n&@H$6YO?ZPNH4(D8SV#%kR;wCwVarRRn+w(Q?XCyX zgbz5#(1Gh!Rz4kaDt z4SQ81-l+izp9pzBt?;vPUfTF2H%;4SsBJ^Dw_T;A8=cup9#o{{C$$FOoAgS+%rj1t z_TMW!+I4;V8kM@yNpF93lhinqux88(4l&k?XoisWVyF?m;IwKbfQ!UE5B$mxlGO@>11jrm9~*>&!xR z)#;_}XlFalb9J{9veS=_Nr~&miTE#0L=z`nPUnyeTjFhQiZu~lucRX}Y-gRh(L25B zDvVuH-yUAw8eS#!6QMXl&TEK^6eykIB~~tNNxHNpHDy|oj<>}3h7$o*{HRkgU8+P$ zm{O`lJ55J{=>1(PlWhg5RCdYSnWi(XR#S<>u2g(=ghu=qrC~l1PH;|nEn$}| zNh)}^;`c{r#2<3B+90hrhS$0q#A!m_JBA*=)T_^FDm4%{iA=<`Q2cgBM zUb(*exc1pMT;HW90`##g*zO@7-F4 z^)vr(N-uM8!=fAw2UF=(9Sd}W2QvwIIi#pk+;nligVp~|t^a+RqW-2zCL!-_%FFth zf3MQZn1n?+^be*Ip9*-&Jr&SKc(bHnq~_d;*e%;j1uF0SZah;8iDP_ z1S4dl?M9ioQD=fyLbe)Z3E67&8X?~dRmrC$HEa~sWs&xd!I2suYYR6VeXPq?YHnvnf!8JVNm zaI8eg$zvIrqg}yR$*7z^*!d80{@}chFly9xXVh)ZsO^L_qv+A@<_+BsDj}p9EhBTM zIa{?8a%MmqLR!|AWsdF#jD$3pWn_*v10x{~W*M0a%~x_BxmU;HacTmIw~1fx%#tRg z9dXQ0Q+^~;lK7smBq0rol7vgtNKz*;O@lp28GEqhCph=8mhgNfi7;Y12^ffQqRx3? zBj$t#8zJo}Y=>aQ=hTXSb~8Io=)0NSMEKe;MtJhtvo+K1_|k0B8LpLSLbePW&iF)V zhzUx1!zUdAPj#GUzB&#La4|ugXM*yOx5_)PX=j`XdGc{Z!41McUK=pS(&A}je|KmF z4{3^^@9JADPl+<0Q=YGCdKBX8f_KVDg%JBuW#d$Z=xvom$9c{`TB zPhJ~qptpLerDQb4fG1@CW2&&f+uTrCN4Un3+Tyaes?4NgDcO{kvXr<6vV`g!0!v9n zH>xi03)@K;P)+Yz^#ox<$a6D}*_+Z4Q_eg6 zz2g7s#IHLP@eeEe*EOwzw;9F(HS=sq8E3(wWGSaTd7e=57z*0f(D?7rIQk=+!iTlr zYA*-x(iQI0nnthvRcWVcKpVXlO1Fr0Fi^L?n>p5uMzTVN)s~StssTpU&am1tGDn@t zn2pa^z8p69ZZ6Ux_y2rb`3&JbN|JEB^eGY+oTFUgEr>60!qS8VM`{^Sct=Q5#uJr`V zh210JT8$jU@T~!J?^G%n3I8^Zk-4=RS%Z=AQ6*W4%$@2Up;=AHR#---zkdv(u_D}hQbyQjyO`&KPlX+sWgiC zex>=|Ra*mQ{yi*9_!UQLj*ssrRPoDx@ zDNX1*6)qtD%*!#gUo`qm`BpuF_^Tm1|a$tMp=5Lla@NhSIoW zu@>Na^AMExRoYXLca@fbXYh#Go-%6!U*q<;dt{)E@ChZ^zG3e9kdg3}ag59z>ss4J zc(Rhp&+<}LTFa@@8c>zrJ5^caYgZH@Xc# z8{s{Uw3hG}O0t82x#M4B=Qi$35Kh;18C#%j8=3jwsg}={W^S)i4J$MIPK&ZK_cf&& zR%VVnJ<4iv_{sKt(lw@)@ZXdKX)lkkEf3A zr>pgh8qBmileQ3^5!H+tnik1riu)~Xa`kI1x6p`dqGN5gW)^6QJ|bjW%s~5!UjO{-tUyC#<{>4bFya0v z-|scEqukJCBO zRd}S2sXbUDNbw3mir06l3ZU*~W(0fb4cc>)yv%K>rrRdm_EeIj=k?MNa%tC%8_mUy zs?fEXR>{TNo=q}O6Gh~)Pmk8_87Gg#IePJMFRY9Gd1bY2;JAV${+C{i7JQQqfVrBE zE4WOF93|r}($R6WjJp^&^IOzFYcz$K2^Xk&;49dC4porYuotz_C`d@{z(v@q3WJM~ zx`T@ve=p=BWFx@!cc^}^YfTH`Uma;R;VVk=y^V>0nMtnIt%S!Z$?(w>b6vYv6TVeR zn4ULoWacI%%H0&~af5@{?SGue7DA2;9LNzS)kvlmX3kV1Yhj+F2;~VmpkPr#_7b=V z*-McQyV~_m3t`xYs}ItL{}UD@oT%#yY9u^JPjjNPBgFR*pXmB4O?Z+ct-f9%r?ix0 z57H$IjlIGK35_kHVnU0L!MEuUden`swS;s%N$)&YB?pB(xS^drq;nr6_V83^ueIwG zZg!+J;TA_)_cMjm6rC1Vv{n`5q{iskvyqnIv^K=JDGV`g3qy?8!ERz*c*?meVCFtI z*-$YvKXqnH6F%oi>j-H{RJwCxyhlwXY{;{&=hK9of2J_oM&p}iFf-ltM4IpcC7Hs^ ze8ustBYcB)fZ=0?6ApKr2|3|Ry_{<7NJFHdy)XdVvoStJ3}@W>SE=DPIs>K&f8a># zXsP)+{N#4H!sf7DaoH5`YZWA8`qg-Fqv}MxyEKk9E(wO(oCgVxh|lEgcJltZrn&wbmE4{ z44nWydb|?aXsT<)YQlw%bU7jGg(~q{L1LXpIG1uc;W2JxZ%G8qP&a&aCaxBZyG(VE zuf)Z9Ue@kC+3A@k`;{T-I)$O?=~C58{M=d9wl;=Vbf`cxbeOqH ziP$z1e$tUHr&st!NFw~6lLm8ou|HOMuPe8jqNs71$-)+g(|0vtsFU+lk2%@v2!HBG zX+l=#d)Jn#887UI|DrI|v!p+Nq%71Z>0MM}w>S?x&&+b3hkR@$ zmPta2FWX=)A$PpaGk-Xu+jN_SCA!?MJRjFoj_8Q@Dy_^BUAIm%F78l`kg*m#wXHM! zWT}v4ZFte5QkIrD%g40tPr2EjCcM{?N=iM}Z{-{{VAG^G;t{XfLfRIO_-ETl+vpK5 z-A>vTkNEXFNZaTUU*Tw5JmP;>nlCqUY?1Jj${p^MD|*Dyf+f377JNxtm&BKXiK7LZ z!n!PC!DC&Gd5`#GS6iMWRu$q-H4&abx$o*K=01F#<&;J>GegT*R_3DGl-!M+p>pPE zS}t^bMVm-1%LI*amZRmE^Y#g~r`UUdC4>LPR~PvQh&qZ>KQ4{9GlAH&c5 zMx`gbX2RcIYeis1c<^dw*n<i>H^|7ttKGuDujH*1CBO3ApU%jd^G{_ zgN{}k5I^i_bpi3W9Bo!W{9{L}35b8`Xte?HH}v3I#sTqn9Bo!W{6j~p35b90Xte?H zvyN655WnPTvjXC)boM8`*#Yry9SyHtQv8ym)ds}JsD}`piKjT)?11>Ijy5MC{)408 zK^etU)U(LACqCcN>H^|+N1GK8?{c)+0r4J3n-dVzTgg6)_?R*F`K*9=zN6Iy#G{T@ z8xX(U(dq)?e{r-~0r9^&T1`OwO-HK@h*RoC#iqnZJKC&(_$)`89S{#XT1`NFm80Rc zGK#7AvgiRqnX1`yLbC%w>9(-TWsVtQgJOH8jTe8lv+!beOGD}2PSb>6ln zAg0%qvc&Ydq91X+^T4$M@eW6;3yA5Rr7Urs^U8ROxng>3(T|uOTSh1GdtKSOfSBG{ z>L8|f7J0<sLafSBG?$`aF?O1;GN zrt&In;zylltqq9()6wbz;_o`ztbq7Ojy5|WrgxS)i2v;PW(CCmBVHNdBc>M@KH`TWe8lwF!bkjb$A^_z zF}=9(5g+ZmZQVgV`r-&5F}=9RBfiz~)ds}$?7~M(&u+&3-|6KqQ*SR(8u3RQ4KE{8 ze4nG$1;mdz+N^;1caAnYAfBMhhl~T_H#ypQ0dc*fogW-*$vI=IU2;Tc;+>9G8xVim z(dq)?Z#&wofcUqLHaj4mrfZq#Oni=`ofiL->dpafm!(MjYj% zkGe9^M@^aNqqa=+QTM-yK4z7PK5AX`!Dyg-WB?DYt3RUwX&W)40Vzw&m0!F%F;{+B zDY-&&A0uUnxe`hp#9Rr5kC-c=@DXz#BYebM35AcCE1~cab0ri$;)Ab%v)rDBtDcl4 z=Bg)k5OdWNKH`I~cC+06g{zs2HR6AAt5a?8h+?i>I751NGc!Nbizy-E94g`VJ6v!? zEb*9IKx%_$6mucLD&p17i(Z;yz#+EMbk?}s&QpZ?Z#xd5Fag4JL(J?lAu|GOS z9;?bAB<7K-Oip4Rr^@6c=25DI48)wc(a@HgZYfKAfjja(FCgYosqD6i+0{vJZpr>w z__$v3P*eC~O?BVp+p6E`lRFW_U)K|9oSpC4sxQe3%kx`+L#|}eBW6iCif;{zSj2)q zcLlHVh@%CYDis`X1@j&;3qsdtw^Y=X)g--5B~_y2&dRL}OOBN*`DAIG_o|jGImyd! z)uVrssPxe|;NEhOTLWK?iNgGhi-5!0uZJ+9*?ytr{-7zE_(Z9q`951>RMRRLNABM< zf4$N%KWA>RtP3>-iSTDif+FuxuNNvZr`zJ&G!=h&W3Qz>psCOvRovmKTuZq4O3M!U zFK%I`!4a=X1f!aXX!7C~X5Oh|8HYrlQ}`uKrQ%Rki#BU(SUlzZX=?3Db6@gGIBC7m zO6{dqYLCkw=&_V9Xev5z)4T*Q$WCouG8wWK;Cg zG{=&Sv+UAw{xMA>Q;wuh=_MG5@Y70yZMO!@aDGKZkVJ#L+N>?t$L#NiwBL{NHk&z# zS3ll9;q-vyNA)t+zi%mgynG3ge;k&0QsJ-5mmvA~VTl6@|60BT$({~SXf;9MQB23T zgXAe;iD?R7SH1+vZwN~aE8H9@k%wnG{ftw@RwQo?OK=F^qeI#HyF(es^pu;OQCbM| zZoFkm$99IpoK<675z^1#5h+5d0xm+1rZV5eoUGpREax?cJ>n}>tE9KlqgoUtwD{lk zotE-;P2tc^I`>c0&V5Zl%I$h3;pwytMVz6R8XSw-JY@OFdW(;7ES-cbb)<*-UOacg zmbPDAcuGJj6~k;={^tqdb<(NU_|jp;7ELYDuq)28M_HcIR3aVXBb*B;^E^ap!D88h zC4tWKs>1TTHit9O9xzBZ_x(MH@B>Wb4D}_>$(Bh(=)26aKY9rhnifS(- zs${|1ma)Qrs;O9y_`kIc=4ldUrs#kWM&etPR`#)z(i)0oYiMN+DOW?$BaYTkVvbtX zbAhG_Z-DSF*MwcyU-JXaJRS1IEqj7qdW~br5>g}rl^YL5)Sy~)MAM1hpN>8EbwRww zKb~nRC!S?fu`BT zk$7>Gk(m_{M&gxG#<-PAbuO0GxwTa1qE}Am)>55|UOAmxOLZ=K<#aBw^EPe7J2VwL z6Yq4if=9edX=Qu+3caZ77=9>t#3L?7gc0lW0{9KnGM zNnvD`r^AS&1xvPl7pdAyHI*h2vlhH9!F>y=P5H_;iNY>%!ipa8N=F;w(x(xhq?@diz$NyM`qt>6*QQCit1RqMrbUBRM9Jj2mOJmRRZk|t3YnWag@ z(SjvS(u%ypBQ+9FcC>;=e5}&SHtAX4UOi%hO@}qzsi~HANzY*xeq7UZ%I0@!x<}Iw zO>#xX0hS-s@}rvG{(mX&d~NxtrqxH->iDvPVa@aX2IG`LhUG8R@=G=C`oEOdca*Jx zKO*_$(KcUwj7_i5^gd0W&@`>RwM~1Azfx*Tg`e@?Km13_f34~}L(}+A7KWNwbn6tG zzE{%^YC37ErPpZseobFE-qP=#X4Cj*px1w+>F+gVxj6sQgr)HJfu5~W4joT5nqH*o z6`Jw~lOEGNf2?X$$1Q&$Dqi2Tq@~1vvE#m1sW+|lOwjqkU$Wtkl$@aT@&`Ds()=z> z?jjHa*A^mI+@G`&F61)B25 zVVX71Uwnw`_g=leM^pZ~0)L?)&i@`QuhHu}H04hh@aGTW z{PA%V?^phW8~ZU{zWPW@xmeRBn(_@OasC_BztrgbZSwcK|#rJOiFUpf42%&%Rc{*CWs^X2J$9devM&UcfRdzYqs zlkvTp=i7?;X5+ZL52~JbYkIQ!t(nJJ&WkkVyBy=?`PN>(N%n7Ao^Oq<#DA0W?b7sJ zn#TDn@$qH8AJK9juTZ`cAK&WxSuOYZ3gs*DeMR{n*7WNY%2(p!dwzeQ<$hA3d?mhT zl>fJy{-HwoN_>afU+>WM`6HV0?W=sxYWyXuzOHk8y=uHapHMkJ(exLZ{#w)LHRZck za4r#hk(=D3vovU$q59;-YG<{UlZ)^GkO*xO_O7%enx*S_ zod161dq~qCYx+w~`IgdOYyNkd((lFP%{jqRey3^rB%9|u@#6edCtAuCnm(whr{&`O z59#$sH2tQg-_!Kxn*K`Dc=>xywiLdl2j^z`u5yrB^D? zFFW(g=6tC{T;BfISvh~y^ibueUu(wWYIMod0`z{YRSqRMSfQY{wq8%ZD}Pm+a#5`9-!N)qA6+{O((v|Dueg z@Xj&Mka%7c=Ra89#~TFqCrUYd93A|zEZz4?-(~3Ggkpl%t8qwFn(s~Y{#)n&czTt$ zzl^{2y~Dg-ov-8hQ{Lp)L@5^Imr1?-)#d1idry|7AK^WCW;y;NJ=!sP)%T9_TCa~% zEcU&ly(`MnkMUl7S2_O4-rRSWqaW*~@2oDl;CoZN_OkS;-W6r($9d_y%E^Hfy}KM8 zi_AUc=qGr49x6vK@0|P!PS5w$NqOV{e6LI=g2ngLX^m$GU7bYVvN36>BBMl3Z?hQ=xs{x zjnUUCy(>nS`v;=u6*2mipqFd!2A3yA%Bg-mK;-WizzMH({3jr=oc_n-+)4DTQ8`~yS+4`X=(AtPt!|gSmDqO+KM(K7 zU#t8N#`qV4F8q(Z!3rKjUjEvY-%~qN{x+4LKF9K1qx@GH{%KzIebwdyXHK>9KX$I= zpN+iy?NI*7%CFl9?=#9jdAjBMn)1I(`7ese|E2OjIdeS!$Cdwz82_I@m-gOWXZe3N z(WakQ{wF_IZ3r?ST1j|C!SFsDC5=Pk?L7($=**uFW25mej^D7n84BS z?{=$zNK#JDWbk`l={k6Zk#l^>dUS!IPxEeks5&Cq^A;;TzD_k8dg(fN8R)w887Jqb zuwWr3UoQ&F311T3iN#y=-&XnoSxe-|E6F*@?PsP+P^=kAm`-@bRTl6 z%d~fjkzab@>|5&sL!Sxq|$OIM*xJ?^7$#XH=le`C&OZ4WL&?8>Bc_fnT13EGOsY z3iMAK`jOs)OReA!slWYl1^(|Uf3@y2Sl6#A@c$R+<@7%d>wt`BKBBQf>pc$si87hnP(2fsDnU>Q?~+w{{Fz%B2#_%(Le^q*o`g5*3lPdVJlPl1tSD@E}Ue3;~h9Bc(ht-2T z8!E_oPX+oXE687~@^`Mc@+YYMyUcuq-%vZuvzNWkD*uzJ5B>$tU!MdboPZ-X_Eh_!2gR1 z^gmaiAA#{+&JOa1@N)DwfnH9Z`U?DOD$x5Y&~L6l|0w9j&R?sJn!)>I1^&+(I`%bl zEqe>rOZoec3j9B-K>sh5!x!QF0WL)TCc?j!v%^X1_q^D8S7Ye5Z1i(E=#vqb#>Y=f z1^MeM$Qi0Y-)`uo=Oi~);J?k#u}{odEw6?@lE3$pu6~LB%InDu=JMHMF4H+QIGF1! z_6-dd>SlCbP~+vZgIz-dnV#&xKsGaDIGZodsPVEhX3eRa>kT8BHOuRi+y&l9w`Avd z8}pK$=XFU^W0E-+7^XRMz0Q2GP#o#*p3&)bSUU%K_SjK9e8ORlT zhq?-umLD1^_6@=wx$XkA?aUZXp=QdqJ~LG8&E*e`wRdQwkn0`l?>Z#GI-~}N@&noa zzognxH|+P}Mm^yaJEN1Kt`G>nluSnNYvx`*;bkkR{?Ll`)PTrty?EoQw!cB8R8 zIWqa&aNcFQa@qd=A&j%$T<3-iCLG4{rc77gKyI)gGtn`2g3o0jJXz8L%Em_jC zxHZ#SfAR9g!q~cEp)1r?7|QfQ{r;SnS=6zzenrDVFFm$+Q)i}OU>K%ZY{q_itfiRm z8|=x{7cB6Gheo?Ex@hV0hKmSolPk8d3fpmh<_j5%Ima2TbEQASYHW0?j_<`rhD z&d!`STPs|h?H|e6`Y`ov$xLVOhEV3L8N$0FH;`G9&q3i#&%n@Nrci{3g9XxK3s;$5 z#MPmJ)kA$sMyZ_lye?1@igK${~mz5-+5NE*&nJ z_%$=W>!fn6yP?pUAIW)`CKw)}>UE~!V%OHOVYYG*^WIG9`sU{Pj!b>J!RsIDiJA)v zI>{?@qtZs>S2r0kQfo_wRz^;HN!o8={a{y1-{u@HWtQ{}cC9ep1vzQ%5}C|JFVoPv z0v?8r9fL#MdTSQ05f&~K^PSmZxo!{#L8;3{M>P7PLH8BnP+MxX-hZy zaJ+$*)S-RZ;R0NEW}!?h57qbQ2FG(5h0Cg&9;?syWYXDuA=f&B)d#I_%#Do><>64F zLH`h1BIb{%%RzxPY;gR@fRl*zM5c37R>n)VA0{2IlWU)qZj1n}raRl$9}W(2SXwa}6R8;U8V+YL2R%kirx_2F8-CIZvqf1!dn9;v)<{YX(e+j)zHLxY&& zm@A8N{a`h+0A=8)tZ)DnibKQV6R;MHGeCOGE*H9L*b(U=sOyr-p~%AiTsH4X5Lb3> zQ!bc@Rb~p6Rsa{dWMCpVuKaitl^bKVvh2}MmO8Uc&aC+1EHQ)Ct~4v`I*;Y8b9gMQ zsdlz!LQ)z`9GZ>3GG27nMwcRQMqzBAm|c&wn73&!(cn-qH=}28WX5_K;nVxNJd^9q z7J9uIU1NhNX;X|$b7dof02w9GOa|BTx&EvaP&C|M^k#_Jy%~TrdWHasxlQ<#*c(hk zrs*?sy&7%zc0mYpwx~K9Q-T;ZWC!{>QTmtxZ zHtXke08DSk`C#-c{rS8A({Csz^=Cay+i?AFk@9@5faxABFZzoN>aS$4 z3qY<>9-lv8`U%ya^|L(NKd9v|R)Ks@f$4p^QI_&@h7~Wr2{f?@&9uK78T>Il!Dij3 z=o{=5LCX z=kpd!b@-9p^p(p1+@|H(e^*>@uQI(oc1RPi|NXJ@e13!JA~9~{FK+)&;+n9qem>{H z^gca5p$OWI`FjE65m|Nu=rcP*njn@Vb-vcmJkxC*cHrH86e%{1mzR8O zh#JP49!tGK57AEpRa_y*E!Pk6ZybB3pnmKev*quV#X0g9uV36-x$+C&YRg{|D-g?i F{|CErc#Qx6 literal 0 HcmV?d00001 diff --git a/full_network/korona_full_network.py b/full_network/korona_full_network.py new file mode 100644 index 0000000..1827e40 --- /dev/null +++ b/full_network/korona_full_network.py @@ -0,0 +1,702 @@ +#!/usr/bin/env python2 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 14 01:34:41 2020 + +@author: ziga +""" + +import matplotlib.pyplot as plt +# import networkx as nx +import numpy as np +from scipy import stats +import datetime +import os +import sys + + +ens = True +run = int(sys.argv[1]) + +### QUESTIONS ### + +# How long is it going to last without overshooting the capacity? +# What measures are needed to contain the virus? +# Does eliminating the superspreaders slow the virus spread +# How long does it take for superspreaders to ge + + +### PARAMETERS ### + +isolate_from_family = False +day_isolate = 20 # 20 = 1.april + + +# healthcare capacity +res_num = 75 +beds = 590 + +# number of nodes +N= int(2.045795*10**6) + +# currently infected (12.marec) +Ni = 200#127*2 + + +if not ens: + doubling_time = 3.5 #days + + hr = 0.0637 # hospitalization ratio + # infection fatality ratio mean + dr = 0.0116 # death ratio, need of intensive care dr/(dr+cr)=0.53 + cr = 0.88679*dr # critical ratio, need of intensive care, + sr = hr-dr-cr # severe ratio + + # dr/cr = death rate of critically ill on intensive care + dr_cr = dr/(dr+cr) + dr_cr_wc = 0.9 # death rate of critically ill without intensive care + dr_sr_wc = 0.1 # death rate of severly ill without oxygen + + # incubation period distribution + ipp1,ipp2,ipp3 = 0.54351624 ,-0.09672665 ,4.396788397 + # inc = stats.lognorm.rvs(ipp1,ipp2,ipp3) + + # home to recovery + ihh1, ihh2, ihh3 = 0.60720317, 0., 6. + + # illness onset to hospitalisation distribution + ioh1,ioh2,ioh3 = 1.3991675382957414, 0.05548227771906829, 1.44454597312056 + + # hospital admission to death distribution + ihd1,ihd2,ihd3 = 0.71678881, 0.21484379, 6.485309 + + # hospital admission to leave (severe) distribution + ihls1, ihls2, ihls3 = 0.60720317, 0., 12.99999999 + + # hospital admission to leave (non-severe) distribution + ihln1, ihln2, ihln3 = 0.60720317, 0., 11.00000008 + + # icu_length_of_stay = 8. + # hospital_length_of_stay = 11. + + + start_of_inf = 3. # after infection + end_of_inf = 2. # after illness onset + mean_incubation_period = stats.lognorm.mean(ipp1,ipp2,ipp3) + infectious_period = mean_incubation_period - start_of_inf + end_of_inf # 10 days + + R0 = 2.68 + sar_family = 0.35 + sar_family_daily = 1.-np.exp(np.log(1-sar_family)/infectious_period) + sar_other = 0.1596 + sar_other_daily = 1.-np.exp(np.log(1-sar_other)/infectious_period) + + asymptomatic_ratio = 0.4 + factor = 2.**(stats.lognorm.mean(ipp1,ipp2,ipp3)/doubling_time)/(1.-asymptomatic_ratio) + Ni = int(Ni * factor) + +else: + doubling_time =np.random.normal(3.5,0.5) + + hr = stats.lognorm.rvs(0.57615709, 2.17310146, 4.19689622)/100. + dr = stats.lognorm.rvs(0.44779353, 0.15813357)/100. + cr = 0.88679*dr + sr = hr-dr-cr # severe ratio + + # dr/cr = death rate of critically ill on intensive care + dr_cr = dr/(dr+cr) + dr_cr_wc = np.random.normal(0.9,0.05) # death rate of critically ill without intensive care + dr_sr_wc = np.random.normal(0.1,0.05) # death rate of severly ill without oxygen + + # incubation period distribution + ipp1,ipp2,ipp3 = 0.54351624 ,-0.09672665 ,4.396788397 + # inc = stats.lognorm.rvs(ipp1,ipp2,ipp3) + + # illness onset to hospitalisation distribution + ioh1,ioh2,ioh3 = 1.3991675382957414, 0.05548227771906829, 1.44454597312056 + + # home to recovery + ihh1, ihh2, ihh3 = 0.60720317, 0., 6. + + # hospital admission to death distribution + ihd1,ihd2,ihd3 = 0.71678881, 0.21484379, 6.485309 + + # hospital admission to leave (critical) distribution + ihls1, ihls2, ihls3 = 0.60720317, 0., stats.lognorm.rvs(0.42140269,8.99999444, 4.00000571) + + # hospital admission to leave (severe) distribution + ihln1, ihln2, ihln3 = 0.60720317, 0., stats.lognorm.rvs(0.42142559, 9.00007222, 1.99992887) + + # icu_length_of_stay = np.random.normal(8.,2.) + # hospital_length_of_stay = np.random.normal(11,1.75) + + start_of_inf = np.random.normal(3.,0.5) # after infection + end_of_inf = np.random.normal(2.5,0.5) # after illness onset + mean_incubation_period = stats.lognorm.mean(ipp1,ipp2,ipp3) + infectious_period = mean_incubation_period - start_of_inf + end_of_inf + + R0 = stats.lognorm.rvs(0.35535982, 1.14371067, 1.53628849) + sar_family = np.random.normal(0.35,0.0425) + sar_family_daily = 1.-np.exp(np.log(1-sar_family)/infectious_period) + sar_other = stats.lognorm.rvs(0.3460311057053344, 0.04198595886371728, 0.11765118249339841) + sar_other_daily = 1.-np.exp(np.log(1-sar_other)/infectious_period) + + asymptomatic_ratio = np.random.normal(0.4,0.10) + factor = 2.**(stats.lognorm.mean(ipp1,ipp2,ipp3)/doubling_time)/(1.-asymptomatic_ratio) + Ni = int(Ni * factor) + +print "INITIAL CONDITIONS (model data)" +print "Number of positively tested as of March 12: ",Ni +print "Percentage of asymptomatic: ",asymptomatic_ratio*100. +print "" + +print "TRANSMISSION DYNAMICS PARAMETERS" + +print "Basic reproduction number R0: ", R0 +print "Doubling time in the initial uncontrolled phase: ",doubling_time +print "Multiplication factor for initial number of infected: ",factor +print "Secondary attack rate - household contacts: ", sar_family +print "Secondary attack rate - household contacts,daily: ",sar_family_daily +print "Secondary attack rate - outer contacts: ", sar_other +print "Secondary attack rate - outer contacts,daily: ",sar_other_daily +print "Mean/median/std incubation period: ", stats.lognorm.mean(ipp1,ipp2,ipp3),\ + stats.lognorm.median(ipp1,ipp2,ipp3), stats.lognorm.std(ipp1,ipp2,ipp3) + +print "Infectious period: ",infectious_period + + +print "" +print "CLINICAL PARAMETERS" +print "Hospitalization ratio: ", hr +print "Includes" +print "Fatality ratio (if all are hospitalised): ", dr +print "Critical ratio without fatal: ", cr, " and with fatal: ", cr+dr +print "Severe ratio: ", sr + +print "Fatality ratio of critically ill: ", dr_cr +print "Fatality ratio of critically ill without intensive care: ", dr_cr_wc +print "Fatality ratio of severely ill without hospitalisation: ", dr_sr_wc + +print "Mean/median/std illness onset to hospitalisation: ", stats.lognorm.mean(ioh1,ioh2,ioh3),\ + stats.lognorm.median(ioh1,ioh2,ioh3), stats.lognorm.std(ioh1,ioh2,ioh3) + +print "Mean/median/std hospitalisation admission to death: ", stats.lognorm.mean(ihd1,ihd2,ihd3),\ + stats.lognorm.median(ihd1,ihd2,ihd3), stats.lognorm.std(ihd1,ihd2,ihd3) +print "Mean/median/std hospitalisation admission to leave (severe): ", stats.lognorm.mean(ihls1, ihls2, ihls3),\ + stats.lognorm.median(ihls1, ihls2, ihls3), stats.lognorm.std(ihls1, ihls2, ihls3) +print "Mean/median/std hospitalisation admission to leave (non-severe): ", stats.lognorm.mean(ihln1, ihln2, ihln3),\ + stats.lognorm.median(ihln1, ihln2, ihln3), stats.lognorm.std(ihln1, ihln2, ihln3) + + +status_active = np.zeros(N,dtype=np.bool) +status_immune = np.zeros(N,dtype=np.bool) #immune and recovered +status_susceptible = np.zeros(N,dtype=np.bool); status_susceptible[:] = 1 +status_infectious = np.zeros(N,dtype=np.bool) +status_till_eof = np.zeros(N,dtype=np.bool) +status_incubation = np.zeros(N,dtype=np.bool) +status_onset = np.zeros(N,dtype=np.bool) +status_home = np.zeros(N,dtype=np.bool) +status_hospitalized_severe = np.zeros(N,dtype=np.bool) +status_hospitalized_critical = np.zeros(N,dtype=np.bool) +status_hospitalized_dead = np.zeros(N,dtype=np.bool) +status_dead = np.zeros(N,dtype=np.bool) + +incubation_period = np.zeros(N,dtype=np.float16) + 20000. +infectious_period_start = np.zeros(N,dtype=np.float16) + 20000. +infectious_period_end = np.zeros(N,dtype=np.float16) + 20000. +onset_period = np.zeros(N,dtype=np.float16) + 20000. +home_recovery_period = np.zeros(N,dtype=np.float16) + 20000. +hospitalization_period_severe = np.zeros(N,dtype=np.float16) + 20000. +hospitalization_period_critical = np.zeros(N,dtype=np.float16) + 20000. +hospitalization_period_dead = np.zeros(N,dtype=np.float16) + 20000. + +#### DEFINE INITIAL CONDITION #### +# randomly generate initial condition +rand_p = np.random.randint(0,N,Ni) +days_infected = np.random.exponential(4.9,Ni) +status_active[rand_p] = 1 +status_susceptible[rand_p] = 0 + +# assign incubation +inc = stats.lognorm.rvs(ipp1,ipp2,ipp3,size=Ni) +Linc = (days_infected - inc < 0.5) +incubation_period[rand_p] = Linc *days_infected +status_incubation[rand_p[Linc]] = 1 + +# assign infectious +Linfs = (days_infected - start_of_inf < 0.5) +infectious_period_start[rand_p] = Linfs *days_infected + +Linfe = (days_infected - (inc + end_of_inf) < 0.5) +infectious_period_end[rand_p] = Linfe *days_infected + +status_infectious[rand_p[Linfe & ~Linfs]] = 1 + +status_till_eof[rand_p[Linfe]] = 1 +status_immune[rand_p[~Linfe]] = 1 + +# assign illness onset to hospitalization +ons = stats.lognorm.rvs(ioh1,ioh2,ioh3,size=Ni) +Lons = (days_infected - (inc + ons) < 0.5) +onset_period[rand_p] = (Lons & ~Linc)*days_infected +status_onset[rand_p[Lons & ~Linc]] = 1 + +inc_ons = rand_p[Lons & ~Linc] +Nincons = inc_ons.shape[0] + +# choose which go to hospital, which stay home +boze = np.random.choice(2,Nincons,p=[(1-hr),hr],replace=True) + +#---> stay at home or asymptomatic +inc_home = inc_ons[boze==0] +home_recovery_period[inc_home] = stats.lognorm.rvs(ihh1,ihh2,ihh3,size=len(inc_home)) +status_home[inc_home] = 1 + +#---> go to hospital +inc_hosp = inc_ons[boze==1] +n_hosp = inc_hosp.shape[0] + + +# seperate hospitalized into +pd = dr/hr +pc = cr/hr +ps = sr/hr +bog = np.random.choice(3,n_hosp,p=[pd,pc,ps],replace=True) + +hosp_p_dead = inc_hosp[bog==0] +hospitalization_period_dead[hosp_p_dead] = stats.lognorm.rvs(ihd1,ihd2,ihd3,size=len(hosp_p_dead)) # starting the death count +status_hospitalized_dead[hosp_p_dead] = 1 + +hosp_p_critical = inc_hosp[bog==1] +hospitalization_period_critical[hosp_p_critical] = stats.lognorm.rvs(ihls1, ihls2, ihls3,size=len(hosp_p_critical)) # starting critical hospitalization count +status_hospitalized_critical[hosp_p_critical] = 1 + +hosp_p_severe = inc_hosp[bog==2] +hospitalization_period_severe[hosp_p_severe] = stats.lognorm.rvs(ihln1, ihln2, ihln3,size=len(hosp_p_severe)) # starting severe hospitalization count +status_hospitalized_severe[hosp_p_severe] = 1 + + +# now remove people from home recovery +inc_home = np.where(home_recovery_period < 0.5)[0] +status_home[inc_home] = 0 +status_active[inc_home] = 0 + +# now remove living/deceased from hospitals +inc_hosp_dead = np.where(hospitalization_period_dead < 0.5)[0] +status_hospitalized_dead[inc_hosp_dead] = 0 +status_dead[inc_hosp_dead] = 1 +status_active[inc_hosp_dead] = 0 + +inc_hosp_critical = np.where(hospitalization_period_critical < 0.5)[0] +status_hospitalized_critical[inc_hosp_critical] = 0 +status_immune[inc_hosp_critical] = 1 +status_active[inc_hosp_critical] = 0 + +inc_hosp_severe = np.where(hospitalization_period_severe < 0.5)[0] +status_hospitalized_severe[inc_hosp_severe] = 0 +status_immune[inc_hosp_severe] = 1 +status_active[inc_hosp_severe] = 0 + + +Nactive = np.sum(status_active) +Nincubation = np.sum(status_incubation) +Ninfectious = np.sum(status_infectious) +Ntill_eof = np.sum(status_till_eof) +Nhospitalized = np.sum(status_hospitalized_dead+status_hospitalized_critical+status_hospitalized_severe) +Nicu = np.sum(status_hospitalized_dead+status_hospitalized_critical) +Nimmune = np.sum(status_immune) +Nsusceptible = np.sum(status_susceptible) +Ndead = np.sum(status_dead) + +Nsymptoms = int((Nactive-Nincubation)*(1.-asymptomatic_ratio)) + +print "" +print "INITIAL STATE" +date0 = datetime.datetime(2020,3,12) +print "Day: {0}".format(date0.strftime("%B %d")) +print "Number of active (from infection to recovery): ",Nactive +print "Number of infectious people: ", Ninfectious +print "Number of people in incubation phase: ", Nincubation +print "Number of people with symptoms to date (proxy for tested cases): ", Nsymptoms +print "Number of hospitalized: ", Nhospitalized +print "Number of people in intensive care: ", Nicu +print "Number of dead: ",Ndead + +print "Number of immune: ", Nimmune +print "Number of susceptible: ", Nsusceptible +print "Nall", Ntill_eof + Nimmune + Nsusceptible + + +# print "Number" + +# households in Slovenia: https://www.stat.si/StatWeb/News/Index/7725 +# these numbers are hardcoded in fortran +h1 = 269898 # 1 person +h2 = 209573 # 2 person +h3 = 152959 # 3 person +h4 = 122195 # 4 person +h5 = 43327 # 5 person +h6 = 17398 # 6 person +h7 = 6073 # 7 person +h8 = 3195 # 8 person + +# elderly care +ec_size = 20000 # persons in elderly care centers +ec_centers = 100 # elderly care centers (100 for simplicity, actually 102) +pp_center = int(ec_size/ec_centers) # people per center +group_size=25 # number of poeple in one group +gp_center = int(pp_center/group_size) + + +#%% + +# GENERATE NETWORK + +# first add households as clusters where disease can spread infinitely +maxc_family = 25 +maxc_other = 450 +connections_family = np.zeros((N,maxc_family),dtype=np.int32,order='F') +connections_other = np.zeros((N,maxc_other),dtype=np.int32,order='F') +connection_family_max = np.zeros(N,dtype=np.int32) + +#%% + +print "Generating social network..." +print "Family/care clusters..." + +import generate_connections2 + +# call fotran function +connections_family,connection_family_max = generate_connections2.household(connections_family) + +print "Outer contacts..." + +#%% +# we assume Gamma probability distribution + +k = 0.3 +theta=22.5 ## 22.5=normal +rands = np.random.gamma(k,theta,N) +random_sample = np.random.random(N) +rands_input = (rands - rands.astype(np.int32)>random_sample) + rands.astype(np.int32) +rands_input_sorted = np.argsort(rands_input) + +connections_other,connection_other_max = generate_connections2.others(connections_other,\ + rands_input,rands_input_sorted,k+1.,theta) + +#%% +# tp = np.zeros(450,dtype=np.int32) +# for i in range(450): +# tp[i] = (connection_other_max==i).sum() + +# plt.figure(1) +# plt.loglog(np.arange(0,450),tp,'ko-') +# plt.xscale('log') +# plt.xlim([0.9,450]) +# plt.ylim([1,2*10**6]) +# plt.xlabel(r"Number of contacts $x$") +# plt.ylabel(r"Number of people with $x$ contacts") +# ticks_num= [1,2,3,4,5,6,8,11,16,21,31,51,71,101] +# str_ticks_num= [str(tck-1) for tck in ticks_num] +# plt.xticks(ticks_num,str_ticks_num) +# plt.grid() +# plt.title("Contacts outside quarantined household/care clusters") +# props = dict(boxstyle='round', facecolor='red', alpha=0.4) +# plt.text(1.1,2000,"Average number of outer contacts\n per person per day: {:5.2f}".format(\ +# np.sum(connection_other_max)/N),bbox=props) +# plt.text(1.1,100,"Average number of family contacts\n per person per day: {:5.2f}".format(\ +# np.sum(1.*connection_family_max)/N-1),bbox=props) +# # plt.text(0.34,2,"{0} persons have 0 contact per day\n{1} persons have 1 contact per day\n{2} persons have 5-6 contacts per day".format(\ +# # int(tp[0]),int(tp[1]),int(tp[5])),bbox=props) +# plt.savefig("kontakti_new{:02d}.png".format(run),dpi=300) + + +#%% SIMULATE VIRUS SPREAD + +Nt = 120 +tab_days = np.zeros(Nt+1) +tab_active = np.zeros(Nt+1) +tab_infectious = np.zeros(Nt+1) +tab_incubation = np.zeros(Nt+1) +tab_symptoms = np.zeros(Nt+1) +tab_hospitalized = np.zeros(Nt+1) +tab_icu = np.zeros(Nt+1) +tab_dead = np.zeros(Nt+1) +tab_immune = np.zeros(Nt+1) +tab_susceptible = np.zeros(Nt+1) + +r0_save = np.zeros(Nt+1) + +# intersting statistics +day_infected = np.zeros(N,dtype=np.int16) +r0_table = np.zeros(N,dtype=np.float16) + +# start simulation +print "Simulate virus spread over network" + +day = 0 + +tab_active[day] = Nactive +tab_infectious[day] = Ninfectious +tab_incubation[day] = Nincubation +tab_symptoms[day] = int((Nactive-Nincubation)*(1.-asymptomatic_ratio)) +tab_hospitalized[day] = Nhospitalized +tab_icu[day] = Nicu +tab_dead[day] = Ndead +tab_immune[day] = Nimmune +tab_susceptible[day] = Nsusceptible + +all_isolated = np.array([],dtype=np.int) +while day < Nt: + + status_susceptible_old = np.copy(status_susceptible) + status_infectious_old = np.copy(status_infectious) + + # remove 1 day from incubation period, infectious_period start/end + incubation_period[incubation_period>-0.5] -= 1 + infectious_period_start[infectious_period_start>-0.5] -= 1 + infectious_period_end[infectious_period_end>-0.5] -= 1 + onset_period[onset_period>-0.5] -= 1 + home_recovery_period[home_recovery_period>-0.5] -= 1 + hospitalization_period_dead[hospitalization_period_dead>-0.5] -= 1 + hospitalization_period_critical[hospitalization_period_critical>-0.5] -= 1 + hospitalization_period_severe[hospitalization_period_severe>-0.5] -= 1 + + + # spread the virus + print "spread the virus" + # go over all infectious nodes indices + infectious_ind = (np.where(status_infectious_old == 1))[0] + + # CRITICALLY SLOW PART WITH PYTHON FOR LOOPS !!!! IMPROVE + for i in infectious_ind: + # go through all his susceptible connections + con_other = connections_other[i,:connection_other_max[i]] + con_family = connections_family[i,:connection_family_max[i]] + + for j in con_other: + # infect susceptible connection with probability sar_other + if status_susceptible_old[j] == 1 and np.random.random() < sar_other_daily: + incubation_period[j] = stats.lognorm.rvs(ipp1,ipp2,ipp3) + infectious_period_start[j] = start_of_inf + infectious_period_end[j] = incubation_period[j] + end_of_inf + onset_period[j] = incubation_period[j] + stats.lognorm.rvs(ioh1,ioh2,ioh3) + + # update status + status_active[j] = 1 + status_till_eof[j] = 1 + status_susceptible[j] = 0 + status_incubation[j] = 1 + + # compute other statistics + day_infected[j] = day + r0_table[i] += 1 # compute r0 + + for j in con_family: + # infect susceptible connection with probability sar_family + if status_susceptible_old[j] == 1 and np.random.random() < sar_family_daily: + incubation_period[j] = stats.lognorm.rvs(ipp1,ipp2,ipp3) + infectious_period_start[j] = start_of_inf + infectious_period_end[j] = incubation_period[j] + end_of_inf + onset_period[j] = incubation_period[j] + stats.lognorm.rvs(ioh1,ioh2,ioh3) + + # update status + status_active[j] = 1 + status_till_eof[j] = 1 + status_susceptible[j] = 0 + status_incubation[j] = 1 + + # compute other statistics + day_infected[j] = day + r0_table[i] += 1 # compute r0 + + print "check illness development" + # check the illness development + # where incubation period < 0.5 --> illness onset + inc_ind = np.where((-0.5 < incubation_period) & (incubation_period < 0.5))[0] + status_incubation[inc_ind] = 0 + status_onset[inc_ind] = 1 + + if isolate_from_family: + if day >= day_isolate: + inc_choice = np.random.choice(inc_ind,int((1-asymptomatic_ratio)*len(inc_ind)),replace=False) + #print inc_choice + all_isolated = np.concatenate((all_isolated,inc_choice)) + connection_family_max[all_isolated] = 0 + connection_other_max[all_isolated] = 0 + + # where infectiousnees period start < 0.5 --> status = infectious + inc_inf_start = np.where((-0.5 < infectious_period_start) & (infectious_period_start < 0.5))[0] + status_infectious[inc_inf_start] = 1 + + # where infectiousness period end < 0.5 --> status=immune/recovered, no longer infectious + inc_inf_end = np.where((-0.5 < infectious_period_end) & (infectious_period_end < 0.5))[0] + status_infectious[inc_inf_end] = 0 + status_till_eof[inc_inf_end] = 0 + status_immune[inc_inf_end] = 1 + + # when onset < 0.5 --> status = hospitalised, no longer onset + inc_ons = np.where((-0.5 < onset_period) & (onset_period < 0.5))[0] + status_onset[inc_ons] = 0 + Nincons = inc_ons.shape[0] + + # choose which go to hospital, which stay home + boze = np.random.choice(2,Nincons,p=[(1-hr),hr],replace=True) + + #---> stay at home or asymptomatic + inc_home = inc_ons[boze==0] + home_recovery_period[inc_home] = stats.lognorm.rvs(ihh1, ihh2, ihh3,len(inc_home)) + status_home[inc_home] = 1 + + #---> go to hospital + inc_hosp = inc_ons[boze==1] + n_hosp = inc_hosp.shape[0] + + # # assign hospitalisation from illness onset + # n_hosp = int(np.round((hr*Nincons),0)) + + # # indices of all hospitalised nodes + # inc_hosp = np.random.choice(inc_ons, size = n_hosp, replace=False ) + + # seperate hospitalized into + pd = dr/hr #death + pc = cr/hr #critical + ps = sr/hr #severe + bog = np.random.choice(3,n_hosp,p=[pd,pc,ps],replace=True) + + inc_hosp_dead = inc_hosp[bog==0] + hospitalization_period_dead[inc_hosp_dead] = stats.lognorm.rvs(ihd1,ihd2,ihd3,size=len(inc_hosp_dead)) # starting the death count + status_hospitalized_dead[inc_hosp_dead] = 1 + + inc_hosp_critical = inc_hosp[bog==1] + hospitalization_period_critical[inc_hosp_critical] = stats.lognorm.rvs(ihls1, ihls2, ihls3,size=len(inc_hosp_critical)) # starting critical hospitalization count + status_hospitalized_critical[inc_hosp_critical] = 1 + + inc_hosp_severe = inc_hosp[bog==2] + hospitalization_period_severe[inc_hosp_severe] = stats.lognorm.rvs(ihln1, ihln2, ihln3,size=len(inc_hosp_severe)) # starting severe hospitalization count + status_hospitalized_severe[inc_hosp_severe] = 1 + + # now remove people from home recovery + inc_home = np.where((-0.5 < home_recovery_period) & (home_recovery_period < 0.5))[0] + status_home[inc_home] = 0 + status_active[inc_home] = 0 + + # now remove living/deceased from hospitals + inc_hosp_dead = np.where((-0.5 < hospitalization_period_dead) & (hospitalization_period_dead < 0.5))[0] + status_hospitalized_dead[inc_hosp_dead] = 0 + status_dead[inc_hosp_dead] = 1 + status_active[inc_hosp_dead] = 0 + + inc_hosp_critical = np.where((-0.5 < hospitalization_period_critical) & (hospitalization_period_critical < 0.5))[0] + status_hospitalized_critical[inc_hosp_critical] = 0 + status_immune[inc_hosp_critical] = 1 + status_active[inc_hosp_critical] = 0 + + inc_hosp_severe = np.where((-0.5 < hospitalization_period_severe) & (hospitalization_period_severe < 0.5))[0] + status_hospitalized_severe[inc_hosp_severe] = 0 + status_immune[inc_hosp_severe] = 1 + status_active[inc_hosp_severe] = 0 + + + ############################################### + # here, one can impose additional measures, disconnect some modes + # + # to remove all contacts of node[i], type rands[i] = 0 + # + # + # + # + # + ############################################### + + day += 1 + date = date0 + datetime.timedelta(day) + + # compute statistics + Nactive = np.sum(status_active) + Nincubation = np.sum(status_incubation) + Ninfectious = np.sum(status_infectious) + Ntill_eof = np.sum(status_till_eof) + Nhospitalized = np.sum(status_hospitalized_dead+status_hospitalized_critical+status_hospitalized_severe) + Nicu = np.sum(status_hospitalized_dead+status_hospitalized_critical) + Nimmune = np.sum(status_immune) + Nsusceptible = np.sum(status_susceptible) + Ndead = np.sum(status_dead) + + Nsymptoms = Nsymptoms + int(inc_ind.shape[0]*(1.-asymptomatic_ratio)) + + print "\nDay: {0} (+{1})".format(date.strftime("%B %d"),day) + print "Number of active (from infection to recovery): ",Nactive + print "Number of infectious people: ", Ninfectious + print "Number of people in incubation phase: ", Nincubation + print "Number of people with symptoms to date (should be proxy for tested cases): ",Nsymptoms + print "Number of hospitalized: ", Nhospitalized + print "Number of people in intensive care: ", Nicu + print "Number of dead: ",Ndead + + print "Number of immune/recovered: ", Nimmune + print "Number of susceptible: ", Nsusceptible + print "Nall", Ntill_eof + Nimmune + Nsusceptible + + + # add statistics to table + tab_days[day] = day + tab_active[day] = Nactive + tab_infectious[day] = Ninfectious + tab_incubation[day] = Nincubation + tab_symptoms[day] = Nsymptoms + tab_hospitalized[day] = Nhospitalized + tab_icu[day] = Nicu + tab_dead[day] = Ndead + tab_immune[day] = Nimmune + tab_susceptible[day] = Nsusceptible + + # compute other statistics + r0_save[day] = (r0_table[r0_table>0]).mean() + # print "r0: ",r0_save[day] + + # rewire outer contacts + print "Rewiring outer contacts" + + # nasvet 14. marca + if day == 2: + rands = rands/np.random.uniform(1.5,2.5) # /2. + + # ukrepi 16. marca + if day == 4: + rands = rands/np.random.uniform(3.5,4.5) # /8. + + # ukrepi 20.marca + if day == 8: + rands = rands/np.random.uniform(1.5,3.) # /2. + + random_sample = np.random.random(N) + rands_input = (rands - rands.astype(np.int32)>random_sample) + rands.astype(np.int32) + rands_input_sorted = np.argsort(rands_input) + + connections_other,connection_other_max = generate_connections2.others(connections_other,\ + rands_input,rands_input_sorted,k+1.,theta) + if isolate_from_family: + if day >= day_isolate+1: + connection_other_max[all_isolated] = 0 + + print "" + + +print "Simulation finished" +print "" +print "Saving fields" +# save fields +np.savetxt("./2020_04_03/tab_days_{:03d}.txt".format(run),tab_days,fmt='%8d') +np.savetxt("./2020_04_03/tab_active_{:03d}.txt".format(run),tab_active,fmt='%8d') +np.savetxt("./2020_04_03/tab_infectious_{:03d}.txt".format(run),tab_infectious,fmt='%8d') +np.savetxt("./2020_04_03/tab_incubation_{:03d}.txt".format(run),tab_incubation,fmt='%8d') +np.savetxt("./2020_04_03/tab_symptoms_{:03d}.txt".format(run),tab_symptoms,fmt='%8d') +np.savetxt("./2020_04_03/tab_hospitalized_{:03d}.txt".format(run),tab_hospitalized,fmt='%8d') +np.savetxt("./2020_04_03/tab_icu_{:03d}.txt".format(run), tab_icu,fmt='%8d') +np.savetxt("./2020_04_03/tab_dead_{:03d}.txt".format(run), tab_dead,fmt='%8d') +np.savetxt("./2020_04_03/tab_immune_{:03d}.txt".format(run), tab_immune,fmt='%8d') +np.savetxt("./2020_04_03/tab_susceptible_{:03d}.txt".format(run), tab_susceptible,fmt='%8d') + +# np.savetxt("./save/day_infected_{:03d}.txt".format(run),day_infected) +# np.savetxt("./save/rands_input_{:03d}.txt".format(run),rands)