From 821e18754b229b6a30d1cb66b58ab50b9d75c762 Mon Sep 17 00:00:00 2001 From: Nissim Lebovits <111617674+nlebovits@users.noreply.github.com> Date: Sun, 3 Dec 2023 19:56:27 -0500 Subject: [PATCH] Built site for gh-pages --- .nojekyll | 1 + assets/permits_animation.gif | Bin 0 -> 59933 bytes final.html | 9670 ++++++++++++++++++++++++++++++++++ index.html | 9670 ++++++++++++++++++++++++++++++++++ 4 files changed, 19341 insertions(+) create mode 100644 .nojekyll create mode 100644 assets/permits_animation.gif create mode 100644 final.html create mode 100644 index.html diff --git a/.nojekyll b/.nojekyll new file mode 100644 index 0000000..8477f44 --- /dev/null +++ b/.nojekyll @@ -0,0 +1 @@ +f50c4fc0 \ No newline at end of file diff --git a/assets/permits_animation.gif b/assets/permits_animation.gif new file mode 100644 index 0000000000000000000000000000000000000000..911fe1f4dafbb762e2e3115fdaf31331712bbb7f GIT binary patch literal 59933 zcmdSBWl&sg+i2OrEhG@!0>L3@aQ9##1b252?rx1cjk~+M2X~jo-6g@Tr}I4TcVy0) zsrfZiv#WN~`|f-1<{#8rtFEt75?tK+Q!phkzJUK6WDrz2aV14jVR>;@MrPPI&}Fh$ zL{oDUMh8oC0NDc^94mC`^WWWn?+X%pqXdCKcoo8DXJ?^n5O7Qk7@|F7NEQ?EazNP# zftWx~2=oN=0C;!+CSd5h7@im&0Ho2R(FA&j34kXCVA9C*&%3`1&}}i$Ezl#z1Yi=w zV*-Zm1VdsPV*orR5O5>(bCa0Hzt4hhXaqwyHJU&-o0veqz@#xI#v}#|Hfe0+f!<)k z13hD6c%Zw$U`Qhv{P&^IAB`~(F!=9R0RRY)ZeKge&jis!`S!zBJXur4j?xtxRZHd0 zqi612o7JwejA@gH_3q;~Z5x-(o^n;-i~Dx($%`k5ub+QFU{G*KC@?HMA~GsECN?fU zAu%Z#l#-g3o{^cAos*lFUr<<7TvA$AUQt<9T~k|E-_Y39+|t_C-qG3B-P7CmtAAi{ zXn16FY_y0*TtxwXBsySIOEcyxSndUk$sdG-6x_08?w{X-rl z8y*}Sm7EY4CQN^qFKoOlY4xTgEY2%Z{ZdXjQlF;kLneBOpL zPm#uT2gTtUj$Av&b)yoAh}^dcQS2ZrGySkrH`+pg z6;UA|uC4YTBztQiwmOpNyO;?lDWN^uh%u9(Ct7&=s5}HFgGa*zku^z5WQrX&63Z@G zZvYhEe2~&oInt)WpU26Rb(EER8UDc_CGiTu3V zjMZ=kwr92RBI06@oiZ{`o9e-!Ki?UeJX{|G-&_F3Vt2U)%dbcfl2I2MOFp;V=Z zxe*-4hk1Jdz)}7nK##TngkpG9nBrVIT9_^EZt@0QPuV}nHvQO~E+@E@s$`(!%}Hr- z{V`o?U2iOXMgEe7rG0)cvw{eJrgAqd2goS^mR#(Y@Q=5#E2n~@WJlXS$Y_T*1%3BW z%IYx~3FY&lJiELr7mMKUxY zS!gDp;V&`Vc!N}PArEJAcE>Hp4Mh%57qW+)9v5fNHtm?%$FQ6Pk`{p+S5g8W_Zk!Hg^*^2K*0(&5Pw94HfCn7} zUq`000Q1z*kT=x|u(EoPb^UHwOi{o0>0eQ#zO$?mQ(2z%aeUg7Rrn|_H!b(7R1=iFL{eKlsv6HeN0IOmXYrxhfCDBVfs zIEe?rD)Jkvm5tOI50Tc(K_GmLCP%18gBv%yLkf)ujQYSdLG;ex5q>(8r9R#WYfl`U zJ)LB4kp=m6X~qkT6FYqVK`EwG=+43HpX!eBnNz!jKudbA=0UL##`inE`tN9XeY>Q_ zxxDX+8tjn?3dPzqr7CSe`Fw;DAuqIJ>MOjjG47bKdgICDHPXQqtvA_yLKX^4N(=qn z{f3cfrRvHGi)1lz6$Y+ldQ7xvEmT`oUxlowx$#G`O4;IktvB@G!9S|(jiY{!?J?e= z%zxCMN{X5L%tY=HG>p!HA1p3&EdMBLgcBRe_9AV`vZOI%-&3DZQ7&ukfxNnZi-&|% zX&1t)wZ1aZWcv-_XEMFkjvP!&7dwJO%E_vv5NB&Y=-j^SDfGl38E=Ks#!6h}M_A@u z4GJxU$!=<)XWAY&W8JCDz7MPgi%*{)O(`N(@F>##%wC$(a_q75w5C7R)akuX^p9ITFZ;bYY@*}#!79pW(odMzJ+I6Q1ynxQs*On8z#Md0*pynjP|p^BLVj? zu;V0B7`PkDN(%pd!ak@-*4W4l@68UpwTWuoa+Y&r5~zPfH^Xz;?{KCgmeL``Ptz!# z;%DQ;^H)EiIN7DrMmC;;iflhLHU6xmNoE}uQlb|#L!!+Y(1c?lyZlUMT5~79GXzbI zoi_!-s~Z%suKB9_mV+q48)Qbm#KSrcYrIe`jG1_Qh0?2!6Oi>|X{o@UKMSF1Rc9ku zdVp6+f}6U=Niq}ef5EfnDT>O0oM-+k~K3Y}lvQU2`GhiMnVIv+PnvzZ&-NY>9y#(*cJr;BG zDe3`KBRLVF!4z#WW=w*Xx~RVGaQ?IFVoV()8O~UvuVpk3J%i)U{tzjvrT3jRBwY-V zpg)_ozM7i`^8+IR$MoalpQ+Fr4Ij^wS`%}y_O}SQrI&NpL66lnW|!7?{ zci$1oFa16~%cI)5>XDLMhJAom8 zi5cQYx#*2_sEtzWw;D#1-xPuqXR@E*oaj6vKQYMn3en${%pZQp?dzbU zjf%Gtr~N08zceCd^)CUIL%-H)-vt$CBSdGkIAz@kR#zy`|9 z`%o_!RSnj#tWavAe#&A{m`ks6(NfrNY#L@`irT|4M+Hf7`|vs>8ptBWr}m|AZx}uf z%!pp3uvv1lfo46J4M1K z5ba7T?H3EV0DR_Rq8PEG*p(2iuA?rO6i&=U?W@BQU8C>S!Yr^5t3klm;vWF1Z-JaK z8EP?b^VBexGRSU>aKkbC6f|3b(&0tnLqf5aSTr|*F?d0|Z+^uFIMN`)$N4$&4G6{2 zaMMUC$I+%ly!#brX%xq_9H+?~*(Vgwd=zJY6T=0Be}YBhQ;!e)9lrq=+r*yGA)mmH zOe2|+fX6Mkwv-@;AM+JIaaJsGyBSp#SzNs&@jyuf`I=VuK@6hr7XP3G|7eD43QQ89 zRmEqAx06D(<&OSIt9X15r$&qFvYdodlH`HUCxMLWqn>2`@IBv5UJ>{`#4QRqESEfF z5Lto}4Fvg>B&Q95G>7Su)e~cBL7Qg!#O5ej$f<7li9#)rm&&Of_Nf<&G^Mn1(HzOB zuxT|ZsdX)>{&*NH88(l_&g`iALShBRrp>Fu<@w9WJhSDI|r zw5CwD1*wQ3d^$no$ob`rRM?D)-?0GZ^x@&~9jVOUY%~YCnYbyLe`v$cjxsT&Xr#y! zmrFpWDctudS^k&;uMb%Og6zVf$ZhhJ?Gj`pjp&DBe#psnB>H3ad1%B@K<21g4w7{E zbzaWEP|iJV_wg$)cb*j z6`&%NVT%ABYVWZ89~zknMme}B)IrC^`gCz&7CBxTCBetQP`V-q(h@kf64SWvNh9%R zG^H@DRItif>h2jirKNtURJF&fS&!iu>SgfCrI1#4#wL^!`^BU1XTN;%!bYaI9 zvvRHbQ@X|y$LhBfG*HcI2?K;Atz81Ys~gk_thH6qF6TAeh0Kg#(@-=qi3 zaI$Q=d(60BZgl4@a6QTNerl4GX=Qdv18Uv(0@I#--aF7UJ zm{@xPi|emnDZg~(17%m!d%I`0O0P%^hYnjFyveWWez`P|&l5I2uX6BY)8?gnA-GbDL3a3zMU zO^0gelCJuPSPq7qvoqEyhY`OHyO5?&1r37>hh0fCn)in55k^X+vvQwR@~tcitq>1r zGh9l4l#Pv4$Rc*ajC!R1fYiyFH&~6L#g8yHj2bPBc7Enal2z*EGydf{)J2|blQ<^Q zIA&X&4d1Ff^_)9{HV8*IB1NHqj&7cYr2fH8ycZ6o28F ze8q>LEIWBa6l~_1^gvW4W7YSP9xsnJ^+x*_>`Nki1-!8A6!JHC6SWDCr76$lK2_T3 zw zeD6qaCwf>kU0%GRSPCFss_+B_@+TEOgTnb|qWIe2mzUzsq7z} zI#$@dz@}+GTEETtSmnHEj*oP##u0Q6^DhF4*E-)Ov^+=8SCrBdt+2JNtyY9sua%_` zSJ!4V5TY)85BGf`sz@Q85*T$=`c_ZOyK%ktla#+Mi)ijqd*s2ocCc_3AU6uLzR|Zj zn+~jt3f^R`sE4=@Y`*PW2V1QqB5yH!Y+xS$z_!`KW!(BelIbzJg?zp_7m+~5xXmup zaBjE_0&YwBZXLpI7awhNoo=z5k15FP@LBDAtW2qTF6Z%{&TXy#95sX`&?(HgORl+F z?72&kwkyrJdhsEXz-C<@{nsGF1cA;Tj>jH^Zcn2!UHg1bcO85*HvC;~{~jeLide*e zvEE`Lk(sBja)#UHrI3{5z&B#Immtx7V$3eHE6GURr}H4_WLEYiAEI<<)H>f~^gUvt zK*zf1Blb}MGXRu%w8tO>nkaEFKJwS2Ppv%4Jnq|dXFKCkssU;FK1q%^oje~I zqCE45KW)=FE0Wdgsyu^*KaUwdTd=WMs{A$~ckc1tp+!gPsPcSd-DLn>cm?(1B+B>&3HxW$^u3r}XiTz~vLo+)QP`K-nc+bP;;i<$J;LD0BcqbTW3= zcj0nXCztsG+hIiT+Pxeze8Ss|8CElu>2Do^(Uj` zkCSXS*2N#@ch`iB*QP{&_-(K8rGBSYTysXxJDgmLLgs{-Zf?=KxovM`n69L%ZqCte zzFpivUOl_PjJHeDH`=zhrjSzOs$12`Tl0%sTfv?i#ydryo;jO4CC_5Fs=M!2bJ$0B znirwHjQ4_FGeOq(w1Tt1s{57asi2GdfavL1rU#b|mkjAY`nC^-(e5d_C+U!z2+p_P-D?%iG0NPhGxe`SJ6o z_2d->s{*Ew?TqzM><@?dy4gP;*p#rQYUi%!YB1-KzBpWCbEr5cH99+Eb8@OWHQY?% zaH+Y1nqAd!xz*!TmnZJ=c{JQpSIQm;c%!&^D_5Qf`LwhW+S)zAtJa7{{fzs;R4$|1>{ixbc0(>l6iD*hZ4OD^Lp+7U2_lVQ z{=GmrOQpe`U0SK+YmM{Gu@e1wiBeY*`x9%qTE6LW;~GfiOry!+2s4n^;q!8<$3taX z8f*4yX8^L|6Dlk7VsAM92i|9vi}8VY+Cg_KmfH1^G;YgQEtWqs6M0gEFXhbjhcjjB zd0RM}*J*QqG!`3)sR&P(*k{*oUwBbD*4lg~k`Rg9>}h@-;0sm~@!0SGPA>}gUc0+? zoFCV}Ti0mgKHuw1{z&|~>iTi%=id-wK09U%n4K7=;C8XqVF#VFvRPn$1Wk&SlXY> z#c@PEV+Vi+VkROema8negg=K_I7x!|h6a_{R_%r;9%wO1NrMzIr4vAGdy+l5?e^)k z-%DnwKX+h#9c9R{+8<>tLslHaAJ?kGV=LLMrRLb<8yHvkji;ExeFIl0C$kyhFezky z^d(Nv924t@I+E8f`9w@F#{_Zw{xA5F3@>6c{MXq^$uj(F6pG(yxvypwFp!lM<#^c> z=Y{3JDUYidJDn->5_~gX(3&j1iP99D!B%5+ZBSO$^5?$|*9TH!FPqfJJ1?)NXWXKg z9nm;rTQn5(u5efL8-wkVPn;>6fAJ?H+c$`+YB?<$yJ)(k1aRuO0rQd}8XgDzE*f5j zJDj>cg0RiHeoW-e>H!}`m(+tajqcS#J_f6)0u2j5%Hd1{ONx=y`%6mEoCwPbvE!87 zhVg`A+(wBNMh}X~y1@?$DQ)>JN@;KUTZ}SX_8)#^`5^H8$RVa|HOa#j;V~|7GU8Dy ziuL2sD@3AB2LS)rnTpp_mY#)gu zbO|sruIJh|(-WkWeE3Otv=ZTNSDIo1V#=&H*6f|^eFE<`J-QA8?y7-F)h z>m%g(W+Myo`U3gr`vz+hfUG72OY+_iS!ojne@O^Y;JrUa!X_N8nlP%)`vCm8O$6>G z;dkEe1Ie(qkfhW^K6Fwb{q<__|9CY)03IL~#sww_s-~b)4QkP#_6+LSphoSlt@`V` zc%T9eD%+r%4eHFGY7J`C{^Qc1A`Pm*{_cVg^JUO^xEzp6E;a{Je={WrDt_t*O$4tNMDDLTaiUPy7J z)k&w_=SfHF4cE#G&kIO*_J4evg6bzqyu5<{(YN{TXWO%6BW(VcZ);zIBONon_WX}; z)1?G_$0KwL_>XThfBOLe%LS1l?qA;)`Q8)>kdpJSZxf=%K!*iY{^Q%C5-?FI3W{3) z@ofru1azR(!N0z3+NEYRZ04_T%lzI!KQ_Ml*SA@_xE`Av{Pk@zE`Sr)E2wXK{NnU! z>j?mh!bUx?(cp*lfsTU8B#D+Em0UUk2-5mKr( zTw5@j%xpN4E>%}Jp2p>HygFQ0G?^t7fI=)?Up$>B8BeD@QeQG#B%ddpA>B|qU#3=P zu{P3BwpgXt7e*}8SiW3mHd_kyZ57}qyPc5?nWoCMHrGGLYokq78&L6wN+R1_z10`` zfnH~3A`}`TZvz!np2j{sP!v>~(vM#o6?`%2CyWSRio>>G+@$ z9Bu~w5x@T!glAy76@>5HOcp|HaGw)M7RHrpMj7n-$&5CsdE1o1{ocToC6$XjlAXsj z&y<^^nLL^w;hsD~7=?=>M!e94B2-GG*~sjh`yEBRupCe$PK9;JFiB%!XwQPh8`&hx zps~;(#bhOPAB0tj3=Xu}-ZxCQhmkz6eO!~$@^?QkG0N&2Mx@C+eJg$F9sH=GpBu(H zSezv%uu`0#pki@U05T~(Dol4iJ}Sx%raLaqPqH{JDK0EME-i07J}#>spgSqAU$8i- zXx=Y9scgSJKB?+Ppg*ntg=2YIGelK(T06#Za#}YfL4Q_1r($^qX;?BTJ8N8ZJ~?aJ zL?;9o!z~eXx9o^WbhqwJgq^oU*U;OxRd!l3wFAd_Z9DF6)6VrC5!y?;UNLd1yGt-- ztMw2+MwIrVz%Eqx6^XQ8>V8nvD(%O(OJwQIw@16u`4~V{I|LuE^}97E2aGbpXxPr$ zk=4=G+m)F;b~%!JgX%lZeL{pc$?vj?JSOTROE4&m%JAoxjOJ?nEJ<>N`q*b?FSao* zRsJudac?U87t-w_t8{fwP*WxaS| z4Xgf&&iB-xs2RA%_1j+VhIorA)O*y*;_UUmcsN%v-hPs}f^2>gvQgh}$)4ceD42Z< zI8t@Va9<7a>wuZNi7Wfrte?`JcMM<3Z*=(ayz|pl|I_Q%0rj@Pr{ClRU3#rk)bBh3 zajRVUZ__NhI((M9Xr8aP!Y6qi*3vF`9<4U+z1$?qy1f0aJZ#&7j@oQnz4lJJ1R-Nj z(@%DkydC4ws_<-^FV2ZV@UN3fZ`;-37Dd|dBW(4s$5!BG?VAYFqm_vhHW8*3I&>Po z2RI^US3M^8H10B~CVUV^T~-j3jJV{Kh0R2@JVYY?glJA1u?6-F?d$Ilx933HvLjH% zrLqwT65QRwh`bN|fLY{YN4xq###NkVpx-Kj+yKLrOB_OK6B!}XPKNlhz_y;)XJdgS zMzqF-#x)TeC0w&Ze0BJN|2#ID_iE?k8{h}w+}Ic|wp~&T+7IF+ak1R7sW@a0AEXpw zEKCcBTa`DjFIDjZ-N~-Y zq*VI?Mq^}QTfi*n*s$CqT(g(BX_LA2qPk!+h^X(DW;U){sR3P=*erJRlx za-nX--B?p)W#zW9dA!2LWfZ*rc-z$d7V6t% zR<{Z6nuoqs+9!{$?(^KWOhr~Y7Rjs~YuvR$mc}ccn?~2p-IJqGUd|VOc|xcvak{VeBGDHW_V$eM2BZri=!^Adw{yc@R)uF; zig9(i40yrQ|AUU#Yn^x*h+MAnCT+(7Kvfo;@rDiE!|B)4+ZvT0=>}YxOT7zoMWm`| zI|6gA$Vf?a;ceB1#9Xz7@NccyY^?%6epqfJ5%7=oh|pT5bVHJfI&uC%ZH)EEHFElv zKHgN`mr?a{xP0tQ(uSs9BT3FEyu5bYPaYFiPP{Bri7yME(iS*I7pq}>>oBHHja`r2 ztSB1lMEi~@$YA{@?Y&h~mruV7?=)2@{D#2fK6x@qLt+!6yxDZf*Mg)xxTl2>>KVT9 zn^k;5LCutKH2C`M_)+AsN=1@At7o+lw|-^5N(T$uXwA4gST(ST!$A`j+fp6&X|boX z&QM?bRBAD?ic7MwHo4+N+oRPNS*xUeEyC19b!&y{?OkJ?tc>W+=69v{fM)uzxu5ZE ztICCb*_nof@LW4ddb1o(r&3fx zggX&7qbw`aMf;8t!9R~}SraFWS?S@i?(R|MWhMX@{-RrpwvL>p>kd|~CMTDxHU3hw zjt@j;V`MEW)Fogh)YMDjv9>L2jtiV>jOHjN4>UcKd-07?$;kz_pAF~_k4G5?(uiJ% zMdMoJJNzX@F3Y0Y~3HfG=rkDeBtlmC3ZQYm70WYT?l#L-2V zAHKiHc-0EcSh-lho%pi|n**NmzPS~DIoMD%i!B&j=Ddh*;W%SKeLm|2GS+(Has|FN!IUg~>cI^XpYm674FF*S zkXW^)U2QaqVU(AGbU@Ve*ufk^R#Zbl?+1M$2*GL~2ifMJ_8))*Y+oi*Ez3c8t3&(W ze!l#8T4}T1(w8A?erm5Gf&D?wR=%OvxDE-b{(XC)mdwGcxS`7-fpmusBG8f##IPwW5s!z#1-=j|OcXJE#5h#O7x*X~rirWBDjEJW8CZ1~ zo@pH6%b}AK8j)lS9LJ5|qY5qGizpP2JiIaNng=2W=?4~jjyDBVg}_x(Ymqg(U)BrH zA^Ku*=tg{ntsC-QE!O!G>>|bNav9{I?i+a<6rDpIWv52iYwG)^1>W5(*eld1n=NGD zHF(M|9DT^*h|9A{EeicUVyZdN%{2!5$c5{c@nf?uN=q#9QS3dL3odQk4amL!Fp>tD zi0ErPK&D(5AKkA}7{ke4+>CU$zCT+MbDo`aQ z*&R6zavw~w9MCA0Ty5?w)RI;{9OHHmDhp0+-!t(OPoFIY+8m|YrGR>Ab@v-W^Wih5 z&EsdKGP+zr=gq0J#c}v887{Qmel6+Gl8U=2rX9;Mm_@ex!x>H3(E>tgqg*0F1D=^J znIfja^=`3{h>|dji121-2>&0+^+Q9q=jYearE&GASK zeNP7u>CK88RwmToC!;76u#l;})ZB#po~C6(kzS)@M}fgVd#-cPNGPMOC*W9QPJ}H!GDePe~O`6-yS?V>?xFmT^-f z)v;mK567mDE16rR+1#!*OZXuvBQ?aD;nyg&c4pyW8nG1NKJVyCY0GM!mTF;)nWBQZFP3^<#vR1)*4Ns zyvDACP42u+)}xKK;Z3G(O;S%)ful|4c*+c#rL#%3pEc8zMoUWYTIzXPNS|tyU0UM9 zB2vL|9K$tfPc68Ft>-AML2m9@Pa{ zl|J!y7abWO@RcH^ccaMWqVd6_KX<=BEx;t|M$+g(RO`fx@H=Yj892@%NpHU)>}|5B zrx?q3_UP>^>O}}@{_L4oiw*O?=Oa)o^4}22zkCFWLjDqu{{WD`lmv=Kptu7HP5v#5 zq0j`%M4$^$Vgf}Ze*{%vJ&XUzYyf#z8H!>py1#Y) zgxK?McWfSnAoUf{H2qI^+y=t#pMD7k3VDEX5Tbj_ zk6kDTuNYcG)j6V4EnLd3k~O*g;WnWF`G#W{09q^aKO6*0R%Y|tsNP=;g1#4YNNNP@ zAGo6j#UKvaL69_i0L`fsl!KTMNdRsDpAM$~U*jOH4JV7$hNGF^+8WQm%?>B)<84hB z8{Gk@ALZJcuXcvw>2)XCTmBqO<;i5pb+q1`EY(?VOmwulT)_0A-pdIj-T&DdEn^CL zea?e9+4(B9{gU8(HMQ+xK)zaHQ3C;yBSKFDNPKOzJfU(yzbc;6x%;ssDynC7$)(*@@OS;LBg z{Y2=B^W(!piwgn_!uEsH>{m*PlnZGjTyvwOE%Jnd$#hCRD^ri~4tu1F%7Yn^B|I8y zfn`p&H$3}Q^9(drNxu||wrajp@n*yovsxC`O$ifL)Xzzqo;3`b5Uw>IIGvn1Y(g1I z^KKF$gt7dPqlhsmxAEk>?d~TXlfyNNnNj{7RZ3N*sLqP5J%sf>uoM4W2qA^$Ahjyfy7BpFG|`6Ep@ePG zb8DPyHccX)uP*h^q&9+e+KBaM3pSez#txEZIS8Wl+$*h=)BG15FTG3+!wU07ohFpq z1%v2sjjkqXW_Bx9O=q`YyFrG#RmVl^yEWH?io12syR*9uUqr_H%|Kk6`>jyw%KPmI z&hz`7n6Hcvy9uf`4|^cf%7^`Qm-B~%>=4Gs!~A5M$0J{O-$%qy$4>z+i5Kpqsb>ng` zd;Ei}%1^f_t zOl0yyHyVJ5+>5l1#ROqH>fc^&3!5_An(L0rTuLvYor_1hi_{bIBrLJl!v}}1`EFz2 z4+x^Ui-wJ$Mq<|uiBIU)gs19%!s90&U_rO_H{r{BGuzybo*NsPU`|4mZ2F5^8X^(p z+)hsN2MPb^swUdw+0gKrW=I&JHYU7;{1djrkT~mauc+$@9CB_+X*rX)B&<#B1y=$& zh4XNm2x1ztJBe$zI5y*SBE&qqk!hjd@nKd4=PI-dYlATpvnA4tG9(=SP5IzL1wK1~2w?xo0;2Q8?p=RjJ|w=~oNJAEk=p zup|}k;^ihc=rlK_kIj$sf^hbQimsWT12t4iU*{rE^W`}6J|B<4TqZoY$0Dh=tREJ;j#PD^^-H6#8lkl%@ABk!(mToaUUZ+6xUF;N#$OD; zpYf=P%dJT*KXSGm)ZH`>1pVJ8A5ssKkyD@<`Fw{{a zA8YUkM64@F%zz z-BJx#9kO=RXJ|HZzG!MW8r^~SO#$w4qd3bPoE@<<=k#{5=n#!y%7tT9mq!c!j*#T! zbpL*vT0n27)p;fHNjD{|=JL7KKCAYkh#>ZGLz-%G!7)?(L`Zw)>YgMarU*kqDpXL1XI2fy#|T%y}i}y=Z%u z7=c??cF(&ed&MB*(S9GZO&%Qy;i4w*DQwR#2kz$Y^&pYBzRT1)ql50KYR-NKzNgAw zQX#%pz8X_5dZbHEm|Xsv93CGAcmZ4hW;5@WWKEY~S85zD03ZN18-bz7FZsabOQF;E z_W^vx8a&MoZxikm6^I1D5M87em$frEp-4=HNW z1kp582-4W81S=aa!CMcxWGM%W(dgrFxk(HJ|2TxVFAkt<3`v0x;Z87iZ4NO}gAe2a zpcFe87n{8e4JdF9HQ<5|0fho+;J&J9rs9THBLKB%;NwAVVKhO;imp3@779baphIAw zBruW7X(R-GWzf7B4=6eWEawX7op)M5bn0iVe^7#SiYM7j!l8>z zbg<8Ht}{MTIJbF_0&RebgTG3Qvr1tcL~~Hjh}-A1$j~e$$8C|A(ri4+cM>+ zBOc2c*Om^MxZNAw>t% zAVoeKc<368rRQb2`yTI{gj{vtx3RC>pQzDY~?hH4U=SzO;k9A{{1GQRxHTtF#T_&CsF+QiGM zI=ze5;*2U|M?IrZ-2!V7aTFMOVxDOsnPw>;TuUBuHEeMShje>nH!0z)e-t~WZk?m9 zQ*@p9$dyg6llOhnVP#2no$&&C3>8L z)9R>Kl>c=l?+z4XEAA-wm@ncUqF|x)nJ0&eGmmveeZft}{HQ<(m{{@1*?Lm~g(60aLf_^*-_oL2^Z7n~8Te%pWj@N@y zD+&b!Hn(KOu#FZehRw8JDowK002v>ie;6Tjz6=;v&dCu@QiA9kR5@jl|Nb^}1f>uL zukxt0a$CBp%U!KAzc?JF8fBpB#-dWKL1HPjf?T3HwZw&DF#s90>g_=Fi!?&ADjVcJ z59LXp7ZkH&fq+}acC!>motm;|2BbcTSNK$GC{auKWI)Fq`MSv45>mrp8Glt;MW|Bq zjkK;o8c5Gexq%oK+EVjDrIP=ou3XBewIG1QBXQ-S)+wcy6kN|dR9pTifoWgjI#4Mf zB`!x_FQ+bz+S>rJ<-FENJ%57>mA19jP)K&a`@@ zCF7(eX#@MfcGm&x z6RVj8=_(*75lx7v-N-~ezfX{{JoRxt_k4(8gQAd+yeK4SU8F?4xxU?qD!r>Ty)?3V zbXL7NF+H@;y|~ZFP=L}b-p3x%!(?#s*ZkaX#MdFN+01L&Kb$yVK~z3H zl4$d+X-70DGa?-7t|j8x@6z67&nIvgSOkMZBOQq5>4l~piZi=#w3@(Hm>FR2z;M>>3-3I2@ZC>+b(N{z5-C zFWb396eFcPUiDdX&1$$DZ6bhqVs~tGp9twNePRP`V!?C#>~x$BesZdP?2qguqw(Y& z(fGsa5J~@JXT&6gD4mLOil=#mPHhT-K>_J)fWdPE8iPN<%@j(<)Q6XN%%CaK(kA?h zX~MUYSlSIA*QV(xr>HU-DKlo!$0teGBI#esS^B4OJK$OP$M;TWJ~JdmdCsYB0W-mw_1rQ)42yq3<|6 z>!N=Sr>5`PPMciOSfArgxsnAJ@JQY(BefvGSYwJRxBU%p{a;`N(g4ve^WXARVWzYTa((FmBcEZ~rOl ztfShA5ZYmv+or4Rs4v{{@7}3V*-D6MXF1>DtZbt|*_O=QWzuPpt=yG+-Tfw4rBK-g zr@5!{nyJQ!tkk)up1G|ffW$$v{{?;Dh;jekYhSQ(-J){~c5R<`ebSEP09oL`41LSl zd*?k4-2X}*gQ5KHU-Afzte}0?zvK~$9idqkw4VAq-T?fcp$C5}t-n0d1d17PZp^>4cMzeSMH_zOxc|F2rDe~IM(Uta%`Q$6OGgh1r!K$`!1a*F-V%z(x2e~ zKunG(LWMW*@5BcE^!L=`i~me)XiWzokkM@aGqJ%V`#tvgKNA~%WY^t4zq~@Gg@7=U z)!97(sWCw!P;%NA;{I7E9ZF8cKNA96E&NY&p#>Ro?=%R>H!y}qvD9jD6BY|;O66}r z27H@PJoprmyhqHJYBDL zXYZyNcb{EuF0ayUKMHid@Sd-2Y)n!?JYNBsTe@CvaJUA1=kbe4@sOyRHv_cU?nwf% z45syh@GMyLf{8S?Ai5!>(pUPS6dX1BK#I&zgHXE4A^mVh10=%;>?ldYkcTzUCl5rC zB61hL;H%Z>wVc81STSkN99&^4TElobkz(Bh1v#k#e@sR6o*mOWQ(_XkR6{ZIb~ zb8i{ccH?jR0>!Nq_u>RC?$SbWDFuqVI}{J@5ZocS6Wra3JB3obSaGMtH8(u}z4zI3 z_dfFCo|*f5lbOuqT_%&YzH694}n`!MU0F8i?qhaA5QT}Z*2v`tx4 zb|kaI*J^YJNvz1P4Y6DU6mL1gH7GgGlRwdD?#?317u+= zr!dph_Oxxkqyl={e*9vZo)58{bzZ`{1zoqh6=&U#*GFePa44+j-w=uH&U;ZAE6@8d z`H#>0ab#FIIY*TB(rO8}=q-obwMYNxwfS|XR?&2X8jUg}Mc3;v*JNpovuDy9Pq5}H za_exfvvN-eeWLxw%XecJ11GMFS3f;~qt7$lLvzybV~yvekyBoh!fjCPNVREBlPtrH zQ`cp6vfj{^f^Wt=$<2MvI=qQ*#oqgh@7Grs-j+o-Q#bxKpk@>Qy05|&|7swHn7~FT zx0}EwNVrMhcMR>7z;*%=ui#D!hMV9n7}_Yfm2LYhxR?K(M`*wJ!Byy>th7<+s4DnU z=&-JW2P%BhROl*v+UDITeAcCT`7qymBqnk`DC{coXVkD!TE8OBwSL6QT%XeE1&>v3?>rV(H zaqzg){%68e2n(g(kaB7QPfe*1kp#NYq|cq9Kg6M99Bt@RK$49k#YXrw#(Q?1v}18d zkOoQ@_U3OAMOU2!)(7i2l>t-2B4k->R9rYjSvvTVZxs@xT(Qb3MRWw#G-QH_HvYDK zbo_DwS^Btk@%RXg^m4g*zFfp_$?55{zcr|3*gF=Q718T6VQM+oqf-`|v+WTEXnF0& zwIdF#W;MuR*t4cXM0DM7N0E$ZT!@Tkr4zpcM2{d8CX(M zQ|Hc_CGTe&*N(ioWlWGm3}#+BP(tkDL@(cF5S8Q1qf(~MUH-KP_+=`kVE~_pFd6GU ztu+>0&7F^~Z^eTZ%KK3azhGD7)7>KA<1bVN4j% z3%fD}-8)#d)3qcfq)f??^oJd?Yw_TA8MDv5T6l>Xq|8%Utf^MEXe}t6)6Iq#PDi~A zH@I8|mLMxDpDUAYs+hfH)+*$b7maNqPGUVE%J7{3P()yA5@D&k$}HC!w3*}L)~`D$ zB=xvMfa6xRPY`>%@SZO5{oDkCB*qa+#_!tLX#F_T-1yNE1e1EZ%1U#x%yDi?=LWBI zdV+gYoheEe^dhz@>l67U$wZ|lTLF7~y>ZRO-=!x0x`{gWnPJO*E{-}TY>uXxi!;nF z5+UT(bPV^|SrQsjJ3bPaBsTkH<#<%8bn8NHHB{TAhR z?n_;%J+7qHgMmF3%1LgFwnbq5+QKoQbbhnYb)jSPi{6oH)~ceNCjsNUO|<>Y2BBd= zA7Asw7#9%`&PD4m>oce89gCU7MC-`=RF*`}{B3&SyHSlWJP8P?&u!WSyUo=AYEWIc}CdPAO(cXmctFa_W-Wt}!C8@tktO zPaNMPibAS57ddr6VOw!bNoS4lTPYx6v$VTwWsP76DDthaQ1-wWUH}9@pi(Da%Kth~ zX0tB@mU13`)}{S%LE+N*o9s}7C=cs{d17e`)l^1NhE#GnDEX3lFIZD2p zwg=>#i+9~yI<&(V^Y~l1; zmI$^B5apD%*!w+1_jinXJOB|__T8(~XD&~=0BW*Up_!AKU}*64B))xfTY!By&BW$& zya2DT@KT-=Nm>c!B9~^=!u@p`Woqe`KFuCD$?0=+f zsqK-M3@8(TY@o2k9E(J{@0munpD8=`-a{RGXE`qW$42*wq4fC4r%;=c6|XY(cfjSt zDwI0hr*y08*nH8c(PE6=c@lVv*SG0YN7LJy zo9~P|x8rdPzcYHSR%so7l4bjymcIqs3sdCXfdp@~PkrY((?CUUOYa)G7W@$WpD2C1 zN!JZ$<14yPZ8`iiFIA_8fp;Fh+;2wMJCKLen;-wcT~Tn6#H=NY|9 zE?qoAeX~%+?FmI~Z%}>=BFhK+KQdV06N}MmQu3ENQVIAsT&F)98?&c1lESV+2qjJi#-B*3H>RHrNUh5ik-QlMZ4G0omVv zqT&zoHUnu5h=gfH7J5d0ON;btjm#&GkO_CNMWw7viGseD3e_8k46FrJHnPHwdd76o ztv^G(t|K4k9R&nr##*DZteo{r92n-5W@Tg2uMnLzf&fi1lc+3KaIsYMvA^dz$lLo!+cKZdP~#W@aTO&o+LkDXvLSa?N>{Ra6#H?1WE372k>{FX(=?P(fOHm# z@wceE>7m$t>DtrfnWpV0~M> z4k`lYGNrs0$N&naWLUwzv&qVGqVglApLi1U4l>-(bj^lSYDYwBhjfL* zGTTP^x=K@&>FtJP9QuAp$RnxLGpRuvKv~ivXGOB(OR7*cAMJwF&BZLo-Du{kRJ=Qp z@18dEu;wh2=|)@X%90hp5=`tS+P({}AIL#-6RMI*sIW#ZSjxFv%FsB_-bBf)Xvo+n zwsKzlQbU_ZdKUo#VWK)+b zl@b+w;Bw2ww}p)6$=$lYKa4-O&gUG?AV(JiHF|dO=3{~bM7&a5w|s4QHlO?i_|=pcUghygm$vh)xpp^t9u zNZ&u$(ELbfSIXKP>amXMckmh|ytPp?Yx7MTd4(J0Fd7w)8kmI}6jvIlGaKc}1Yt3< zzN+hgHa-Cwt-z?(F?rum+g$m>RuD#8f^BVzTw7U&P#Q}cU3uG! zZK_Z)g53VJ2rP_h1KYN*dbX6=*1Z{SuUg@+8E;Nv>5!QfXtM3djp|sG@@VpI>dx$_ z?&#?BZtZ6&GqvojCF>ZiNc!5}S>Dn9UAI`2wkyZGW6ryoF05-fs&l0yUt+6kT{mIV zwu>sPOWC{ow{Z85@w#JSFs!`7N7i$4)P2d)d2QRmbJ1PV)bm%k=b0rQPO+0bwFg!E z8=G9mL?#OQaWix`^xGtJYY_i8la8u_m2dcX?I~MbKRM!cD!NTy^`fx$nlbgwOLa#d z^fFcUF&_7Y%XLpGWOIr1-)r~pjH~gorgMns2|Vfv*tPN#49t4>Gx4=IADCuW=(x2E zjAr;fxPi zvkqBeY8Yn}*E3hy?__?)H};G!{W_65W0T^mIP4~pDIq+9tTTe3lf+31iPQ^jAFEX5 z&Grl^0A-aH-{0l~CE)!Rm&*ttz&H~w8z={o;?d4_Zv0!f9!GAnO}_7aj0B+{QDlKt2z4{$>bdj8qt^aK*%*=CIW z^F7D&X*`lsy%3d_YMe=eDMFO2*zj$Sy`nv-q>-8_{p2r$NB6v^6tZV<&$ zj&W%(NZ=3*Ng1bB@HKrxx2jqvMdqIp)*oCzGV{yX_ixS&qtc?%!BO&=2CS%s6J|ED_ z4G_iz%JpM%GSVRBi8j_1`O2ual|)nW6$Q=J_tw?br$a18w9lkkIeq^p{lNch@BB~Q z9Y=lm$4qrNcnm%KUL-OR{5~|6N4$PS;YW@EEQLzE!8f`-oI?b*ig?4Q-gbB+)KMaM zqx6~Ncw;0LkGSJ(9TT_{oa0%zlT5m4h6w_DtQ2{Y*T*idz35D?largvTQjn-f6Q@; z{3oteLegwkKjl}fuBPQ4+PCKQ@UVDdjl$Wk7kN{CufALLL#|=rxR~E77Duce5&wka z$ai7kIN_}m{$Jsk909sMY61eAVZ?s7x&s!io`6_!2X5jdf7Q0lH2(x3eU{@C5W1)5 zqYT8$OpK5PlcZznshJN}pw4UVgclh0R z*gTN#aMucTtSBxD@F`+i@07vge@*{kR$3CI%|WE{iV z0;<0z5p0VoZ;Xv4h(+u{)mAibQc$HdWH~|Vl_vO`yS!|aGaujnHQ{LoNJFjhtD-G< z2uNokL%fF;?*lJBV6A$hl>dP~a^iOf|GcR`DeO{6$$-@Tq7ro+&L+dwNgHK`?DB*t1qzISikP$|_ zj^P?UshQ!J8GFP)(@iqz6v_Q&TBMY|yi8ecqCPtrzdJ2W#{*BbUSafn1^$dj;hhjSNjy zx}LC2%*ffc>S~C){*m%U`9kR5xg38iJMlY4EpNW2+D>6sld5m?gQObOu~HI}w7Q|~ zs_r@7U3TVZ71}=cjSYENRaSpYre-t(GfT|-c-1-(m+J`88d*+kLvdD?(XPv+A@*hs z4CACf=bJ;3u;D2o{0p!yA)@T56VomKPs#;>F{Ns(bM&coQ6eGeObR=-pf|!%Idg;Z zlD2#JN_`;Zc&93fxhJTC*pgndInp9tFZ;G$hk;gHrtBEp`iBK5q4h}^B+(I z0bIP8UO%OMO>zCJ2$WPDg_yQBa=KfjEpzk9Z!1gm6#QD{m6U97*U;+MjTS?=nl9`GjA-+pGav;B!{I{NUjfMZz5{fuzuTo>|f5V0Ua!Z4MRbjj#?e8bl; z25Q-awJMgGp9==2w0B+X0APUF4YsiYrNZr4cB)LsTf)>oA=HFl z37c&uf1_k!$*t2brLh<+r+;Y2;=Oe2cj+_G+3@+OqC(@nNBdbza9zRyhSYqoWe^lNy^cay^h`~=Zfz8+9BAyaesOH@k16pvCS-=6n9C!s*9%arknU&GSu15wj*{NH~e zSG%-N{-)WPf^5+Er%= zv|DQ?SoB;KE5DaN?!KIB@ifN6N5wGr42Rb5nnO<%+X-VK|G;46i>7V)lP_?kNC;dX z_sYrU`^ZYaF@Mjy!i#mdCw*H^rv~?F$sL`ieGeSXPVu086^6W)cb|{3_OR}9`^_(M z#UG_Lp$ut;MnW08x7vAA_Z_d@IQ!OyHHgR`EPbF?Eu!9xh){m7O-}zyP1kKnFOLBq z^dKPSh`(OAXqcuGnw2-0D*$oPeBw$oEOB%6Wy_37aS2FYZ^VQ71lA#*wN}>apU)n-y;-4^P@F-LMukvHHJSe zLd}YK77!w}8zY8Kn%V4JLmWtGDZ4S^NXZqc!3eV3iB_hE)7%w8nv+yl3#f#|JuJF^ zz>&LvIaECHk6MVxqn?33BCLwyhvpDs2}6(jft`e;s;wV!t@R3$j0r*FX6he-RcQUWE_3O-mV~Z5%+u>VUQD|*(9CE^3BmwsLT&zopK%Rs(H*-<6K%iCl z4spVFWIX^_pTs&zZq)qOzI;8wD`*o^B7d8QT3ND!U@|z>x%yR#7?XLxUPusraB8}o zrkrs_u>;`BUtJqyJ1T1rO2H_LVxmp^T$1W`kmw=!YAx(l7h*#E&ooyibYPC{6AOXgJUV}^Fgv9<6 zpfJC(%nN}GAiU7>m2k6KkMVO6v!!zggf^+wolCgDC`AG6x@3Qc{@`*F@ zN!#TtbCBI)!{*g)cLl?S&%C zg<{MMhX;8lqxm05bFV=8zvS{|GV-4e@)aVBu9oDbY;t-b@8#OFn}|L~T@+2Z8f`Qb zhY2ddCvz;nDxoe_zBetgi45*(EOGD>?Y}H>zUM3DDRoCL1-q7dd6k4Umik4;t6Y`_ z>EsFXK*GF=DO@38%O#Qbr5!&ZL0$^+hmbg-G7r+yRH0J+j55devaIqFllwA@NTmYw zlHcm(C0k_>`BH___VUW*@@jOx6Z{HTEwFKn=iapf4yfGTp66CxVRls(gj5#m21(l~ z?I|lOGcUgzEm0OSdMgi+AT1%1ubPX5s4NzJ6e`-m{V)A6oPYLn{-YlTi)j2yzrKVq z{)uS3L@{7!3Ro-yHm3vAudvVsZ2jh?C-zT_12&`si&((Y886I<7up0&*TT{hFm3yS zny`35Ou#fQEVcnFhr!BVFu24&>5Bhe>Hc@>_kaH$j7t&3F(U{EmmV#lCGfxGQZz`Q z(CTQTuK#B)MF<)$GAfx2ku;%F=>OnSc;TU-;wU4L!sTZEJC`CO`!z*aVb)> zn$?8Mr5L0Jkx9fL6OR9jOJPfkf{Y9Ydb|7&mtq)=W^H|q5`F*A3zx#r8JU1# znj9Jb{2BqBT(*`f8M*(MqnD^RuRk1zmRZ=03IR_RwN&n9Cr3sR|AU@L(E0i>Y$vDO zOmc<>FYy_!`ei0Zi3!zo75>|K$xF&X^moPk0y&W1K+wxhjzOvWu?_FbP7W{4Dkfcm z_5Zz{od1KA!@u{#TF8ysuX5q*4cC>qJKRq<#J%1cigwl=um|84KVv=}c9<5CYZXC{tyB+hrkb|T9RC@n3TOR0LYHRucSJfTD=!9aqtqGlOtJZQJmr4z?PO{9!ruvj*(I~@p@e9z z)4bS9JfYIi1(P3w>l2MXl*N75zN?0N6i0n=^e~^_{XE$;C*?27`@+ULJ#ZDuJDpSnvN)Z;V`)P^7e&3;zL!5_ zxE?QU*%U{c7SlqVs8c)~HU5Nzuy<>1Dj+=GnSJgQy*mzfd!joU7Vq8;@;TESj&Om> zJo(RPZSC@--kuElM)kj?f@>8=WQg?-lr@B#v=qm_i+hJznTJFj@d4c?o*Vm-%8DJ( zH-%FhGQ35FPL@-x_e z4nDNxqk!Q7S7=?V=R_gp&B(CW{hzo1%pxijU^go+*E`5l_lqi5bhRA@tEt?$&hp<>S57+i zm9!E4g1C@wL=2ALgA`-jS>|Yw-FKE~#RkmW^g}+|-1&kN8p;=RGs)Y$Zxbiw9U5q- zl9Qjinci5~$fazv?h#Idm2=H%vriitMEq$cfkWKjYeaOhr@V2Wzcg8193|qK_1^>M zKV~9Ne#Drw{1(E@6L}^Hkz}kJl;?3RNDZ+1xC5Gr)=2gvKiMZk-Uww0s>!2q>xa{j zSET_r`apLTt;w|H)!12o7Ca?elUj|zd1B@UNu(=ehHE3E}5=aFxm<+%1g z41OXKCpYJ6wewRytxvNQ=!}m9Vx|fm?Z$dpFaY;L(xa?1JI*Zx@ZMGG?uUwHt?ta zs%3>^jBP%4CY(J-v*GN`)G2zCqMEOJV{we8%{#Ki$=xzoMq!7hnzVW?70{4CXJ_(O zrGKW~oTxZzbzlU@KWAHWs`BiyiXtO?^l^vd%98$=5G99NfmM?h0B zLTJqi&|$TPmLllPi0v3?LcPai*w*o2f0+&ZvXAQTN|=#@RHU1|iDb#Gmrn_GE6~;X zq1oELNW`XHox6O@D-*m}HKm_nn7{UA14uF-dueGwileG2Ic!Yn(QYc5u71;R{XO2f z+a>U1O#SjLvLiJBpK{;gw&UF-IqH>srfZx`vyOc?$~b?smpA_-w6S_P6{D|9pZD*Z zFsyY{WXd>{g;$5X<@c5?extc-rpF{1W`L}R-9S2^VcOJ#^MW^| zh~6_P#`oMD;`7}y1P6sZzFt|K4-5#q-hxQJrye7>`)~;hlN&CQqHmuDH9JLqo?UY? zuA)P5zrbzBNgU*EZ67frFLAmrh0`n#`?nzlTrmc`+V|*3rbChSsVWYDM#~`e>^WoY z2jH~&6AtU+YpH2Z8*k?a#^V~jMfD&A1kuWRyw@;f)cQhxsd}*K!G>!5L}$PDA&BE| zz;Th6gkC_oJtfb2Zv0Shj&E(FpH0uxNX#P1>Uc)E3@smu-gsY8IV zjPAMX(b@5evM9HflwYhJQ9k^(g>w z*Dpp6?^E``3d}@`F1K#&iwxHq$W;7D2=&q5ZWTvRLwAG2M@J%5onrElL+d_9Rufyy zUq*IW8G-ksw$ekh%wm6uPcKoSKbk{|<9$9EtLr=;}i#-wQL z6jh^#ToU#j3uge7NUiW*7*#4S?oO96tqpd-kN7_7@uoWDpA_@T@O zp{a1=OCcYr_ZXkqpR23^w4e)swVaMK{7})>=fL0i;UHJj2dQHCh~j6pAap|gskVU4 zG{rgil3pZ5lS2w1N{;bB{zqo*9|LMiHpTrO;FSXKSLRY1A#YQg(m#Uc6CW+Oyo{8+ zd=#Z&|IfA87V{_|5N07QScJl&q%<^ACy5kt>nhr^y0dT zb>)0(UI}T1#3EI?l7fRJ=&eHO^JR5MB~*03^EnctVAGt?S*RqvbE4tSP7o7a+EXY&Da zuViZ<%0rG`)wzQ!@H6ZByzupu?IaU%; zqw2YY8$yR_1)?m4$B{Re8a}nx2|P3iSJa9gmA&t%Wz%ikJK*~OgG&B)i?{z(vBQKq zOqRbKnqP!EtN;PyMZ*0{rN2n{mjM^pM$60VKe`;YjoS

ac1BOufS-{fmr;aUx+G zFR%s$>=@l_@gme;pp^e6=3f;4|3?q3mk{6ts5C+!_+OADf1yYo2mQP&6)v#^jl|aBi(9&L^t<)SZHMy*J`RcYfwsm?#DZBYv?S7BUkPc^ewk_^dW-GkLSiQPD(yP!No?0!x zKHaFWS*PswnfbFlSy2k|>+;S$Ss(wT?c42pKYn&J{LZHbnmF^sxDkLzXtfcD$^h92 z!sOfE2*#22l!L?7*vAZoL(Iz$BXl2uj)ze>?{7v_g%Hz#03{{`Ag1Q@QP}w!KugaN zB(N3ZV>I$HhEFJ;HuimN-6pO`4&XNq2XYx55hLzrvlt2dw(Vq%%82bGEs3RpICXY3 z^Hd`Sx$Ojf4FOVo;}WzIT%}Yd;{e@qroAL(uD#vh;S+(q91#m*I7g4)U`t#d!ZA~) zKt5DJei$FRi64kL0+JM~abNm1fg$1`C6y0dzXVlcxvcOikU7sbH&%zKsN~**!L?kc z-MZqTtKGV?NLjwZZ(Ryxn_qd#yzUwOEM!^RS1n^#H&k+GU4P0pZr7llie3f$Q7>oT zG;cs=P~+vGYu3DGBEz2XTg(gU(AHf{1fww6+_HD{=P$hNJXltE=X|}aP2Y5Xs-V{c zv7fQ)L0C|5^3o($aOyh|!<4SXAj14(ixHl0Ge~Hdl`&v4GWG{d=}*o<(iW@7HFSY~ za4PoX_8;cdvaero?-^FOh#M`v%Ux zi@l4G5Pxns%Hf~lzuk;lzdTuD5fcq|5MFyepXOfc`g3I4!+m`=o-6eHXQc;Pb~l?l z|6C~2afcc3r1o}JFGUQVFb258z5cVAsOnzMDJX*F6C!OmESa$9lZ!F^YVSgHg(5y= zVf@!5xw@c~nTiKXNG{47&Q|Pf6|XJsk7&jaF7ngb*I3EMuTjiCh~dXaaMBs!TsVnK zt>OexAr_o@3cN@DZI9F2fHQ&PDaC~QJxu$f3HFVO)D-VbBrn|-j$B9!5p!LrA}ej= z5UMzZGW3r$w|X+gtMtBythzWl5gMvq7a6s-1P#MqUnt*b%S!frH8C?Wq#?K)dDjT? z3p&Q7=%8Yt{yU?z@{r#v1RABkpTx`+A>ZXI3$0J7kDAgpBf^D*^h$kA_SKnw%upL; zdrUxWHZ$)p+8&^JW{5FJH^UQhANPOb!XBoE%vwyV$jZZo){DP`ug8@3jaM=g{nZ{p zby=SsjZ3;qRMcA5%y=NL8b-`z^h^n;Mo+>edp+4)pu&2xNdF}Js?&^rwk=vx+Ag^J z7i=xRO)0Sa0<2p|FT&>DeVdh(lu#5UT9G;yg5Z*KZ?k81KKdENfTv2J79$zgJkIH= zQUpyU!CmIk4o&O3L6>Kw7^R?J438stW>X}wkH{gMMx_e6;se>FqLjaH$Jy3`_ z{8>6h7p;(eNTO&xmsF&Vsk2%x*RlLY5Z5ip>Xu%6ynU#Foui6JpXEa$aadVzQ?ax8 z5uv`#Ldg(s4Fb?ae_j5IDFPo2?Y6DSp#1k6xVh{ow+bVY@o6n*g8GD%V*@75xdYc^ zOynmFRk@^kx8nk`w^rgp_ z_K($d1;5Zsbhv0wQiioUrcbJDZ7VUdGB^IxX4krFBWoP`&_29_BUd%P>X?$&EYk$7 z^4JN;fT*HNH?c9g^LiTOH1l_@o>Wl$dN9BvX=%9_B=kd3s2}Sy?>c@u#IyOSWrsWR zO}vUaWGQIljG0)U`tM0_=*lPmudM?2N*o$QtTjE$jNNGVwO(k#D63+|?O-8g-z{7d zq1M|$j_xEUi{n(%h0USyZAY*Z-)3{|Ep|11O?;aA%xksVkq4!EXBov}fY9A&CT~Hq zxYs5%5%8V%HBEe=2tDP+V5Nz)i-UdW0+&C@gu>9PWYJ%GBHh42hU?EfRUCWib0jhr zVJgwL1*DoehvhOAd8hChd7G^Wa4{dg1u=*lc+28J(ieNpC^27YS8!mY0YiCyfJ@Ot>&X={?5CjagC3G`M{n&e3N?CitE$ky%yDbKJw{d zAkE+B_wYRCbe!LM&9m2T!_&PfhtEt?Ze%e=;B{V{?;p|E9!p7Y2^3`fPX@c6nIqPp z&iDs=cl=!o7Yxx)Zv5|B(4L0;pN4Ix7=d~yey<=HNVmXS{?C|J_Ru$r0f@H_NsIy8 zsR69~euSRRn}h-65vJSMzLd9tB~)A5 zwkses!?W_GZ4GY83mV5>Xoh0^df=4ui7O zgEa^vukswXB;7B=Bgx)4P=zqN61!F}z7Z2}chHJ5Pp~~oih{N;l3WA__bd|ZNrd&I zMw`NJ)>ee1&C$)DL#GyrU>HfulIWkZrt_XLj?|%$`0(X)^k3;gYqv4#0E_t%C zNk*|!o`wgI*pY}>SSjJVv7V4@oDsF@O+=g_yutBE+?a736g7SV5Pu;Xk4K|48xi03 zO}Aqs_SI;7X{~A{TmlMHLQ<`=cM1AZL;?>k1ErvfHEyC)9V6WmVW@MW?F@h!L>L<+ zwHe_+v6Sfdo5_YQdSW+6s|0WH?0hAW7Uh`2#4K3oUu-npeIILm3h> zY2_VWkaRp?lBwVeEF&ET=_Y-n})_CvQ%)ngfTURGlzfsy;fptRtGap;eF|| zb$)L~ltPg=zy_&J2-N*IJKrmZJ~B;1EP$=enAyCqa!8e@&=Bm_Ktx zN?~;}v%GFdIsA5Q7K1&aEoiqaLp``g;h-A5*{o5fgn|qdf17=Jkf(G}?lPi($H_w6 z4s1e6IK8f0VroFr{~F#A>QwHHCBb9;C3TXFJ5kP*hgOE+PEZ> zah>7qq*+7yPK{7(xc+*P5(k)Eb>38IJ8r9tzrT^2*Z=#~(Au8@1BlRbl$JCMZ zs31~fc-wtPcew}IT1yg?W0Zeo9F(MxaLo^gjv0XY>l?O+4~`;g26ykTaknI$Btq84 z-*kP)E8X>*z1#a8Y8@2J6HOKi#hHgrPb+;xoFeVrK14Cr0_ zvfT$V{gS^RjmV5=ER{xHb?~g;UUl}fJPz=Sl-`W@nqqu?OIUFz)1x!dW7p8FP+4{1 z*1aAzU|alk&-ok2a2G+RvP`somY0Ep-Z#x(F266u7%R=iF$b~==vm1P_X>xQABOzx zoU?g*BGe#MA_D9pY+gwtULsK!Xfz9OT! zJR;8Uwz+?p#!I|WX)ufBNRhY7UCbW+Tr2BinkPXx6nv~Gn~MB>dyS&Ml4n86J6Tg2 z?UfDDzK_EY9!$?-=3x3ttdC)Wq!Mr!mEYCKsi08>Fu)M8>h~^=>G?Wpk*H~-;~ss+ zzC&%`U1j`@$P7|;_XhxJv1I^VlJCXM_m9NW0cby-q_gO{#u0@34(`X9(fV1)e$d!g zaI(!RI*r_$Ao<$`Qr5mx-B%X0jUXt~np_k;19X1TM6WowK$?a6#mTk-_o z5uB<@Bl9JqlNBY!xM|%A53HXKyeLGgvo))W=H<=Aljw84q@*; z?I_^VDAq_U_0!Y*5@QBaLe_2{SF+ox}I5Ag(c-hwxX2@Z8Bt$LMqKH{R3 zRUbm9DhDjj57;H!lKm{P?2j^^cefRS? z;y63yy`tPC%-JPw6Mhss0|uAa*=(D!pUe6EeEM}hN4;*Rnequ6Qe$`;S-y_6U% zC$167xMS}13pVSscsZ`nsp<0xvVQi{1)M&iZbhZDU8>htdLsjNQ*%;|m#~_IY2y`C z*>C$3^oo+=+K)X9^W@Up@w$;9Tm0j8<>XN3>9q#t%>j}(7QeJc&0lxNuMyLY)tJ7$ zCd-yyNB+4Jf&Nu?jaNjwpfQF=kd-~ z@yCuy{~}4cCO7`X=ivTQSzE8bnegZ9hUx|rSM?Q!TGR2*8=`+ml8aigX=%-0_T@O` zuM_6ZW0stCsO~&7J62#MNs(68CJ1@v(gA_c{lT%$Bq^nlrtaelNiuicaqWaWAmDO) zjgX?X`u*q z=If(2+odIyJd`H=ib9n!>cujuRc7gEp<1i}R;`3`S<;Z23so%Q)+3Hfh6rph!uub z(v)jPt1gKd1eb!a8X^ors1?N-zLVv`R8tFX{bn?4H$quxSuIXtwq`NRLX8?NLR$Qe zGBMflRyoGBN$AhZg=QYrHQMKiDl=$FJ1@z8l%8K;RCL2U#s85}%Eo;jL4m=2tBE)! zfUd<#@_nAWUY1wxhoXciUFSW%*kK*&pm219zPzO0K*EB2_pIYQDlnE}RL)^wZb8)u z9a?F3%{O-Om)0iJoSN51V}iPs(aF^H66LoQL(pV#<>O{&mP^P^GilBI(sW7fgD}#1 znXFKHe$93e3Pyv5La|ZyS2L8$hVy1*-bU9<<|AT)yB7v<+j&26Sx9xgCVA5ns>yk~ zG)iQtNbo0a%B~JPd&<1;`6XaIm@evJS8Vx`y=6&$|86g4cS3bv<9n9N433Dln_OH- z<!!P$hZ?=+wQoV*?#Bcnw18|-OnpU z*OARh>EJnk&bR+~Sg|tkAYQZC-V|8dJi>Er^D^NDKJOXfWV!AqVxYQ>C;Mzw9c6}O zz19{8G{0Vqw}DEMzuUT;uLpm@z5kvR_~JP;tpn`5uQeg?_S(qi|MME2+HLdp@EHks z^90Ps{J0|kgMAMFgO21k+@k}I<=TDAhUDL1uLDWt8up@u6p*i~3o1#5B03`ldZy~a zNVJ5&Cm{tbbL+vWaprw)LJG#y&_ggwNyS=13c+pD`{Lj#i2p~{R{%2aYmn-XPhV;P zc)Y|Y@qk{kvg%M=Phzxu@m?zC+At(d{cjpABJ@c&p*p0y=!2jgQde~Z3@z5}oavxq z?5LB~!XD0iJ`WdpOceLeO+0u|hIBPH8aTcW^ykMA|Ndo;#E!M#50b>l7~Pz=x1%%> z5mvV-H;R2?U^eM|jhJ-WSAAgK0mU$eB$jzp!XC>Zl|{iYbzz;yJBrcl-Vvp$xlQ7M z{Wok>iDC`qUrB8ThuD@cQYTOjNvanG43Vw8hH*4s{mqMRg;@~HWLR9ya6Yrtw(?k6 z*Z#a6Ad5adKRNGz-e7W4b}DFnRqrqj3C;V7 zyVd))S{&+d<2&-2WDq-3#@C>dA7tEeuu2VK7=xywyY;iQ&cped!ZP9_W3qAIkh3W+ za|M;!=A@by0CoYTB6l3~q<^N02nuy%++aC>W(nOFN5K70_sC43o@X{V-haPG{_A}E zXBPVxWCB4%@-M;!0ZS0Z{P#lnU+Vc6ZvK4=Vw#{dMa(}q6J&P-5l_$wGU##BU#JkgW47#@cDoG<1gKmMKi}lz?%Fq0*OY6uV=sdKksr~f6}KkmY=;q&Eg%G+ZE(z2xQTU%wuu5n#|L$@NL3Gm&85%yR3IY1d zHz6EMXgKK_6wlJmuI_)VYf;|KRY}7wUq?h5mw%I?+tg0Sz)cxm0hl-otUJ*K&F{gjYJCC;3T)MeeWbGc|! zrsNoH^H!3DXpTx)JFqP6s7ZsQ&9_gAy7g)aJ~&HjB;B~;@o{}WRTLkts2j3-y#V9p zljcpIY0l~(;4=#?6V?|8IlVemlXsEYZ#olILxcvszrVqO=5a;AL;7(=KReQS&75RXXVtj*@-_}>tdeCDdd-VtOXOUdW!vj5l6lAf zv$X534`@G~_fDMmycZzJ2eTO?=FNtepfY9B7-(DSmWhcG_i*eey1}(xP;K^lP*R*# z*jGVt*_KjcfDe2EnQXD&GFzglYH;hG{8{uLZlA8n+aK-TADvs#!1CAY=(4vz z9?I|FD_t|lH|o&5JMT$Ff(EG-ew+G)kX}x$u-|6p2Gw?P8^+-}8`*DThe0sw3p!@UeXg2H3z)BmUka z?^B9V&H ziLAQ#OsiSBuzt}Qn$eWXpE<;e^6YXCLa4C5OIBFjCZPX3<+TTJVLM>~Tz+{3d=7vIsJpq}BVw9yyElkopJm=;^?m7D=E|tnarbL;t;$IklqJ zz7l}B{%-0vLJ`~IY8CrlSf*CGaa#*kr3Eu;#-=<8@ebBB6Tl(+o>zwEu5{W3I>8W> zR_BAC{keyr7!Kt{frZw}wbap&CWPd`bR^e-VLWC^F@Zec_PZE&_|MxP}G z4i*%?vz49eQY!&N%gV_eaTu+rRel^UR$Eh>s6PyTH&IYROF}PS9#PVe)>~LANUX4F zM3V@romt#aQm|p8(W(5bVicvsW>t8IdRK-BCRtaLpHddQ0Xx(~-Gm+3NlS&lhJ!ox zscyZNmSLZp&!%op%=GZWkcy}E7U0}aKxV{Fw41w7rp$rk1wq+*r#`8#tmja1nDvU0 z2LS)tw025o>T^!EzjQAoNMda_s>XX@EY>Q|r!WdhNMhhDq0zjF=pWzV-lA z)tXL}TJV3WZlPnqHO5)nI{V>*lz!g5Fi`$UM!d`VdER=`V%y8s#R)e40;`Ozv=iK< zXZR8)uvlToYvDNMr~h>){)#T}&vLUB*kNFw%gJi|n zJ3ecNkuIZS?SeiA0|#QS^p%AKPi530#xFQ)%x2j)rVsbeVPH{HF$Os07F+hi62TDh zZ`1@ZB5p%@gQDT_s|Iv7Zk^Xymw-Y`HKr%bCL?y1B(Jq>jjohQt0qyPXuC9hcW#rn z=roJbPyHcdlEJEilF^x;A^mXflP|x|xa!_VZ4q56r;YCOuW?WNrPceBh8k>v5YhzV zF6xE2PWK)nUZ5X-+I;a;WxCd-VQ?8YpUzR`@y(|*b)WF8rbR~&(EGyJqaTaS#sPhv zoo2mfX^m=?=Bjst=7mZVi~Ye)wE+>@W)wK1j$0UcOa+$|^G_MNg{<}7(hHnofSH}| z6iYo|*UJ4E3p?M6+nvnX*?=i4GXu3PwRnw9d+FPDzN~;ltDjEMlh>$ zYYk4%{Z!-^0EKr{c$Eqhh^pWI@a_BMM+jrnqBj8x%bY#$jqPyIst&vZx0>R0 zg44#)Qt;eBxxYhYX7esDgUgt@WK}Lj;faz$)g16=f#y0+qKp9CLM}Pb$ydQkbpJy* z9ySUs<+$TqAnHfFf5xNEEX_F&a?4spdKx)R6k!lDdGy@ghr z;oIBzWtTNV@AYUdtOqaEy8_HU<*-BTF_K|`fvU__aZ$T*Hwy*`5h?U{n$_!-&^msE{x4X3wloz#I0xEE3aFh z>(?BJ&C3xCkL-E@ufxKaH~Xdcll-rLuHu_t_uV@$Kis^slX4 zY-xwT)IOF`fo6o>N#LqbzO&k$dt02U-nzsP=D%^v?dQVe=d&UN&F!5L!A z9|6!x-~MoL@Raz*>{!~`TVFInV**e_VVFfQ{ky@;-2G*UR(ug#P^4QpaOG7{G1<2; z7~J-xg(%omTY}P$xzYU@+f<;eD53WoDGfPgWs-y5o}lEQ1CJB@^>Re*pY^|r*pj%0 z(Ck~fQ%dWI+w`JS-a^T*Bn1R=LR%FF24Ol|?8qp!7)Jr@V=zOHP0aI5pg&Z|5?ajV zTl9XNgeJr2k~D{XHFYa|_A@p0DlQ2lQVt~wv3obsYB4bfvb*|}nud!)tGIj%BU%EG z_F1KH0Qvyo$%Seq!CrF4L?j26D<&b6C4g>Hd!a$c4X7~fC1ZzXfZZ`Fh(dn zXs&P=w?H_vn3oWPr4}9J&@T;7s^F#i`_8clsj;UyinW-&C-ol0Tye~ffH$YO)TFqx z9gtVf>UqwV5i6c170{?cQ!eJ(?HXTV%Y)w%e zCew;a#ubx-Vn+rZtJEs>V>G5BFL7cek-9VzQSYgV zu9hO1x;I?LDA^!1*|JpDX)syyG?`29m;Xz=x@x$LnpFrF#vM$8CoE<(Y`i}>y(3d9 z78T>`UW!*}%*MVd=Spb!DSr-DnmKIhs6raYZ~)ayDj)>ig(~^v&h_?L3wkImzm&Q` zEj>lerSmjZr!M`Wn6|??EpRk`MmS^hTXNuwZEdO1SZhYuY6dG+20CX-;&8@-Sy*PO zsq#W*9A`X#b7p)%It-LZH*APRs?Who=+FxzQeI}2YvRjX*88z7EoSygD0G{+)g!Ed z7QnNfGTVGdlMRppERA9hR%nsP5j?VIOUVJ=wdNc(sAx83=K;J?)@;rWB??0oG1Z}G zTXKsbnDDT*LhMNi!*b9+=aROC!+y)rBc{XQ%v?OoW03U1iK7kQ%TQ5;BT36=R3{8^ zQ>Ar}z@{#^dtjkVO61#vA1Wz89V|FEN|iG}tfVY7K`Kxf&D==Hl7dTAAIp~Gru~g# zK0%f?3`mA9lf5|35rGQVISsRH%lfPyu8Zyc%2^CGRvc?u>@sFYbQa$!9xHC1lm$~_ z@>;~sZLUY0WeHZ=M_Cl;Ucfgf$yQNU&d=MqTD3SqA@D3eg%;rk)4j zF1lC^cs40JhK;LI_rp#tT~jR|H!TICjK;Nbs|LLChP{@%V4MB`l>MMCH%Rm945`pt z%N|g#SRP;qUe@k)FPXuv^fgoHW42o=D;ixZT^mc-fUDYu%Th=y*r%>M9Lt#xE5!;f z=1?y>AIsE#uA+vs#G|Yd`Bs&nUL6S@5u92r8VXQOu0|9uY=NsWF;g5T4x$}!Q)0JM zTysw!s{zMJK@E2glC0vDw^$fq&izaP;ai*4md6h;W(>+e#|nA{HFH20%H_rOw`H2p zGV*4O-=uPBQR*2jpg`xz!u_%9uMA5`wZ<(ZOzST8*wr+9_5ACQ0yGVqj@H(04T$Lt zEGYFmOL5L(VV}Sn^~W0(Eoy-AkiUvU8r90HC3zb83?$=O!3sV%>G3ofXf$m}h@vP4 zv*(#D3X{r(LqXJtI3LtwqPo87HT#7(!@Wv+taBd0CWF_s$m7({um}b7ID5bmhT$|s zNM+=pNqO@4`${z>TD0PTH^MIIkoEECzqL*r7jcM(M{+0rvS^qzvTp-y*5=if2iioHPg~EDQ)vf^!1f`_3!YT-b)-ThlGdC`wzL#@j!oC2b0e z$4Zq)9m`L&byiZhUM#ehy-A(83r)La<=JzMBo|t=|Ga|l8bR+SGW54J?Sc;KA~56q zcy7C}P?NnQuU(JHW9wJJG6B} z^h4oXv}(rNbWnG*z2p#k!Z^eCMh})y@QlU(_n9l`GOWHT&jy*!;q0c}o931fzP@&g%O6J9@qy0XBfh$%_ zD^dE>T;wFy%Ay&_mIuhrXjP&Wtq0N94ZGF0N}G!D5tQiOR&H=Y-38BNIo2jMRv`p(5x_(Smw(cg9qFNB>%T;+}9p zeS`#lhezal$pXSOv;ELRh9RUN=(1r(?m{&7eCng46JcX2GQCBIcZLey4$5i#qwIRHNG<+hzt(}dwmU)-~D|2q^!UK=avthiIm~McSujz?}f+};iBy@Bku*Tkf zlBAOW+)G2)ejWo9T45XP7iV7rsau6!82W-sBe5J`sSEqktmb(2#9~y^5?@A=x(KnbRX762Qa){F=X;x0 zhHBBmz`kdW$K>(|XPfhcJ8T71+|Ly`qZN^MD4og`5@!P*4^x>v#aKL3XX(zwOeiQV z0xM~NLSTJ*rB{pQ?DJEW%Vd=!G+}Y4k29;_B<}Lhl{J;I5t_?^TCZw(1;XUa$Uln{ z+kxZe;SGq-GvD%MzTPfE>TV2+B}8^wkMl)kURa1Q%xZ7oo%*!wOH--ny7&*1QM3*F4bIIWZ5fZZc5^@-O^lx&M< zP6DY@dcY!jaW8C`B_k(2cVMc*D#9}nzU*kd?(kD*^Kz80zef@x5@3GGfXDodHx&U%93rSVgX*24wbj6T%2koPLx774;8LVqBN&m@rv`vZ4I-uRTAeo z%nn4`jwaTt=imx$5cdJ=K|*o6ViF2?qE)?gOZ-iSr9*SZfV#+?U`)56zqHRznYxGZ3*&_4hmr)Gvp@h<@5yuVbr%QAr@`=bR z(&;-B;uoD`HRnB8vfmAsQ_ZRQ=wm6o zi$x1Y*EQ{znZs&U=Do>IUeG$Xr_EcT36+H_-pE(C=h9BdmL|xEdaT#G?9(IU+p>YE za(vr!WZOIZVN=#Z_ed&t0^3_oFb*9MRwCI4?ArpkLjiuMA=1}(zpnfCL@$B84YR!4 z56aiKW$)Wr+P`huuJwQqJ6AQ<5;sraPw8Y&QKGlb*}n26PaZ%spO2@QPs^+H&!MU7 zQ`XNo$P#Un?vB;Zx^C{T_%GgHUuxHndTk$*2%gAA&I8+CfEf211;3l1U&q>)dWB!D zwO?GrA8zp7bA}&du3s&xUy8A>b!;jT+t2H=YxgW($3B~!9XvGay!8;gd$K=nTfYrg zPmT(_>~DC#uf20!ixa#(3&_YTufMSq9FTc-rB*ewxNf?R*KdFQWBT$){p(|vfT_8Q zdVhFyrQ^@V^jj`Vjmc3-Mj}fVL8Y0|0ElkFeG=KPN`JYQ|r*JJUBvOmcRywX?sGZ)>7YLQ!<*ZYf#d z80IKz%;y|Hbn_mKP01-;$S1@)*ovg6J?dqPxV^YUN1`HbwK@XK+@=9dGP+9M+d zN3%oA+AB7X22SJ3b4rrijH$~z7B_}gjAwf1<_8{sSI7FVF!!zJKXv>}?r;qXgl-V9 z-&kFdkog>>W4Ov8tqHmV_bx6aD3J6-BXPLw10O1P`DKVu)S7UcqV#H2 zU(r$_1Lm@&GI`t{D_~i0rFoFBR(V6?!+wG2COaZ<{+MVcS5%5!As;e9p;Dujnk>{t zkESdfkQTd!Zl&IM%RhAgI7VI2#4=t5qMOx>sdAHDR;(Cc$ze`+-`G7lQ92`hn_-}- zbm8Mh08}bXykft5TO;w>(fe9ub*GXglG&xntr;|aN%#lX89F8n6`L(LMGi2~6g1;Y zWqznjSkA`VU~p)2xs)z~Ia44~C$Jq0_C>2Scd>Jo9!k@PCAu1q-)}FECcC70EjVuo zE&P&>4{vjOxl_ygmHnV9)3QC76aI^!2MoTBf#+?=v>Nwgg9Rd!2!M*z2ebK|q0tR# zmyR52kb0l^gcL+KwVxrZ2jG0>Y-?D9XN-UqNLrq?X+{8;a9SzGK7DxXQ*eYKe6RbU z;k;WWJ5C5O`H>NVL`{3*uNDaa##Qheso`TZeUx$9(9qjz>6VKKhU=JxxC!QT+*xFY z6F0_*-XT`hDb7uC=E>+DUSk0%|AhEy{?*&dDf`>NISB#8$+#(MJ#=3iF1DzRH)En8 zRBvmQJ=H~6EiTv^T8cQ+jlZpdkAwG-Cd_CO?X={@%AK;Pv`M~hl`|>LBFnSsUXS7r zPUULVq;Z*M4W0E~_jHTw>qU!_JF0TY^(MPz>x1#<(aL%=k7%7!Lf7lQ=2IL~SHYxO zOWvf9>6;O?Qw<9C&&ZxJ!$ip$ck(CN)X_0mO-VOX(#@~umBKDOVtA4Kc#i`=WyXsV zrA$2TR-AsloZZRSe>{-*Z)Y|>r+pHAdrG&NUc1`MOwM(Wt!&5&SH3X%FwH3_+mU+fnjU35qobimN65Pw4jmJhWVOG;ykdNTeW5s@*^Q z$wCmXq|G(@EIsQ9| z`yXbYznBfg*g%c3aBWbw+uG9 zwXp{>23c?Iv>u&bTwaOxo?ibo=Q90XI(q_jkv1B8K5#%?q&i>Nyc*`eT_k1cG`_#) z+`;gF2;E8_<89@U?kTPiUFQVJ{8%Fy7)4x38IcYjdXk)MBunQn}J}E!S#yZT!*hTC&q?Ho-oXfntBa;>wYIb@JkXHCJi@ zG}z>hxLEHWKmEc7?fUyhv7h{GA> zcl9G!$4>OE*^W#HB1cb8_P_FKKkdhGLQ);jmLO9l#U2u!rp9S8t>(x75-LTHdi`A@ z6`-umm6N1syhNs~VC=S+q#32Q7)~;QX&huuU#bz`UVvqyX0!x*oN0JiTAbn;b$FbP z3r&3-Si53&64!O#e3A#nKXz;%CnITA_+8Iz%`N58@qn}AVX3ru(+~TsgkPjYx+E_z zZ9a9SX05E8eUme?qJB;iH6^+R7P5-utc?b&=j5vL{S~vbOcX95zb*kwQ$FaY_V$0NBgWvw0gtvk6hMID_Nn)9m ziISli%Wi>lS`%SNm8SB6fxrdvW}&Nq_OVmbrt^N)w7AlF*0SEPoY{7`Qoe6U(dtQJ zi@*NvZ&sdThYkvNayu}<@QM{cu{FCF#qwl*5H}5=JFK3RwwzD#DQZoPbg@narptuT zyBikJxVj%%S6^wIS1m((TvVh>dmK7#^Ji|g(TsXrH`#sZxJU*WaOb*lrd0#`Xlzn; zrlPH_?h~nwS zYT%#ur)$u)0zwYhXYn3zY+2kZo*lKBV!?)rNK{rG0+|v<5nP`8h;tc3xK@b}p#}tJ zdXWH3$i(m#ra#E7Yq;ngV^D&g{i)(^0{f)!&~uuSp7}8x?E7t$3?SbihX?F!A>%qne{S*%=hijAdEpcmkdcdG zFEBuXQTqy8r>99gpO1!xB}y~2gp2MEL$i$Rl*>T@5tSSlWQDZ>(-aGH}O(UM^@_VfIp2C_t)#F~nzAkBm) z*f>x)QR>&Zwnlf^@z5^{@zBWkl9cV{A`B!+5wv14T%koNr<@9#o6ocEpm}Alol3FA z3l_55#|;(0v_O`PMpikNog04r>a#Yl$vszkC@jQA5F%EseOspLMHF?>I=)bFTdo^) zL^q1vU8e5D*wgSu9JWk)v0xEO*NIH=_%+rqi7`x1{HJ!jq};Nqe1r_5vEpf-YQ$tk zP9qBqZY%Y&DIBp1T`6rvQP}i^c~W%--hrWptHN~ILCHL?>5cxX|43O<%@XSw%Nz09 zMekzPKt7$Jl}6wq_iT+4(wT(&y6gh~VguEtrLGOUx&^UIIMBcg9l;NNv%8P5W*X7o z!p3#2m%3TA=l#rv^>w}FY`mREa6?D9kTSOr!)5k^-WBRR=zR*6rytNfj_&##;@Xq+aXQeU-ksna?j2tyYx zQ*@9XnaLXntBUQt&_M+a4x#aLZmJ_uja2)Z4f$tZ3btOOvlfM~16dkDv(<=DFw>Wa zihX>K_y~f%4Hg~TgZA%SW32e}iGBn4`)Me#nt$vi<2-q(KrNXcFjH=*;YZ~4YJJvO zOn_Nf6Jo;SvB8TQz+8q&>t*X?P!*k>)`aP8m1RdII*@fP#lw{JwRKz{{|+V&(v&0^ zD<^J?4*A;hR0@={Ek$&4j@B?s{7B|^XL2Y-KRV!O=a3KiXb=Sb_xi=zeaR}N#v3C1 zT$AI=fcB!juMgq{TcADkn-XPpai?@YIa2fFm!=v6q)Vmbm}PCp+s4#PEA{b}CEW*S zPm9Y_gMQT|c5^2PtO#rUc%0-zy1ciY%0Pu_oG$K_YnbcIJqzN8&6OqaB088GUwF-h zbAGjxmC2ew1HV14dk)>q^DSBn$6bUvPJh8ETP+8i4j%VMyAZ*fFmmodatix0X|K!_ zJC-9GS2y56W?(8w0$;ZIUbriL4O*M*F^?{pF))05K+d4lO`O4{;|R-%0CcLIGNB^3 zNbB^jh`29CS&2luWm}Bed9K2l9NbT@fjTJKicaq_KARONxYZ12uwa*J;x6AB-hm4d z4cPe~pb|-VdbSX|*9Se#eg&h#%G=y01G`{3_S?wPRZ7!z*PatRserrT@#+12T%mc=xemIT#AtU`>gXNhKU7*R)sQay9o0i zUxft!sU5jD% zXP52eFG|Z~-&v!jdG^sc$JDp@aAFe16(-+qsvEP_!pvie8UA$ypYVSe<@fVq+k>V#+H{b z=)t{#n;K(c@Ue3as_dpngboh5^f4&F~kyaNfMRz5-Y)N%Tf|shu~Fj69ZK#^aZ3fFs8XC<=;KE@iQ=N}=CRKWQt(PrjNHr?XXCS|1bX8r($%5?{;9vzQfpu0 z$odlvLX_WT)1Dul(qNOT%+du%wCVDr^`SIGO4Cm$9S2quCt(PCU0Ej0w1>EuyVPvE z>N4_2BjJl;@&NIJT$#hI$;h8MRz@?{QuU5s7+jv9b|ez(xnz4@vPxglWw>0YTcKKD zIF>ql-p(HulpL5|k zawA|rW42`;Lx=xW=)&boxw43+Z!3>b5Q5Q zy!y+U0NesDde>cIJaM zQaBoedyElz%jKrOR@i69dR^#kR#tTLoPDN}Yx7FSB3bA_o$lmL=K@!xp`Ps_nL8_y z%>$74OXJ`=^EJiJ4S{3#qAvDTpOHFTUJ)z7PR>273@^LH=4Tot;l(`e7V)2 zB`q)|Odu;QOv+5Sb2wz0x|~6w6q3~qf-p^Wz<%@$Haf8c>YkqEJ_XmzUxYgL`g7H1 zLsFiC%0<&;A?zptv}!>|VoAjq*EN!(*OC+NYP&SM116DYNe90A>NEAc(z}AmU|evG z*jNdQuVeX?{WJvX)~$LXq@!ioD7EinECDTLjxJ%B*kD4eg<~lUkgwIK>BWxff`~jI zJg1C<1GY|Hdz9)WJ0>xb8q4*gzVb&kg}aAxS-p;VCLd2|-)go-TD_-Y12!iUe|y7m zTJ->T15lbK7_wjPe7?e>p+w(ZKSUzbvzTmYtq7D#YHdzxZlT-f41GSWYIkdrFn2dU zufQ6Q26>8BAASkFHM%Xw_F*JWz%;0uFq#+nXudZ4@_6~LTN9q!L74>u!khS?%>k`V zvb3Q|0M;~nl$!9>lt&XApd$}`|3S|ZK_*ew8?ezoZJJ?;rCt*p5c=Z-wd%`_Ih zEp?&IDZEB;mc#44`Aq`U2(XMei4S>|R2ioA?fSGEeQs!3uTnmuC4sU_ZA^-_wC+A| zyeeVEnL*0^xPkq*)JZ>&Kz7F-xdPGlwgIG!eqqZZ7zyLAj>#1~WKSObc36rURZkV2wjhp=}M?0_VWfCaU>ujy3BE9{@f0o`Ir47&H`JmyH)6>!J;*Q^#H}@ zGoV>f#YTV*!_0YkHdN+I$s-FUKgM(QV`v4}%Ja=cBob+dE!;+v%}mldI(EaxfLR^v z^@<(Cupy4L!*mp#O&;wd=;N(t89F}#_Q|qRU>P87{n@!HTQf2%5KNJdT}jy~A!jCf zJu}e1)el-iLV(d_CQbCh5nm=#K0BqP3Ma7Uq>So>w%#WLG14}< zr$lN()65@#me+!)buZ^x$cj zhP17-v#dpwo&wjhf>#TVbeg&3ZPCEZcW4~_^3kS)>wvMYin}~XL?^Z{Hm9f7q*^oG z)Dl~;*v>^bdw4%r0aa*_IIJW%i|0jZwhqBCr~|W!kEB?&NV=dXIE!vo3RLKmAoiGl zd(x!8oJ@GO`QV(F?Sz!)TU5)8^&ApY+4NOIlqVTjWObO`K`pRu(O}JlGLT+Ux^GrB zT1*Bk_8aSKL@vibE-%L`OkytS^DR9D&%N}C(_AjQA&O3Bxb0Q&hJUV$;)QgfTL~NT zn+aYryIf*EUR_Eg#q!GaLl|&7g>ffZi|w4Pwi1cPTbX)aEznwn9v3(*T)T#c4a=RU zL)=J(7R~Kku%O#ebZT#t){vx)iOwOmZ1ifX{3*}3St7nFXEBHHS|7bWIyGX~t<*a~ z5!f!gu*A1z2VXg{#v%}0u=SDIiM&ZupvPmh)xW51{t+r{mO68*47~c-5dLedbI)W?u(}g`c|zM;#}3_ zT+bt0sdNaE_ngnUC}BjiRpjoN3Pd~~IL8;k+qW&cs#&5xfu!1cVpb51mut}x5#*68 z8&CU!S@{%f2TU>hVR(YY@H=9dpLrtz-@n+(PJtis9V$j0BnFSaeL!a=NV{JiATA$r zcC8Na367;?mUYW97i{}}$=&`lAL$h`PTSw`5`LXz2&QPOC^y|p}bHHPo`pb!q1$2IorKWe)?5BayJ z{I_a5e@yS5)L7wl(mQ>zX!Bp_rnA19^xA^%zE0QLa}&pRk#>OuM60U|=2by1zPAi# z-i+)Xd3|ue0dGu;X)!46sPQ{fxT~5fq6aRo;dOTK32>vn{I=9cWZl@{%i3vn{4J$@ zPAjk!*6`bl{;m^QI%Q~Gw?^m-#JZ2w86U&qh;ic4waqr3 z-FCi_Aq`c1eJ3|IV*As>0;sA$N=W=DLH&Dq+%<9215o`C)*b8xq$m8kpdi53ecljS z+8-GEY_Wa16U5++|B`LPp8IR1>=>&^fOrvq!_$g>A}g^2B5iMcxBchQrmIw{+?%FP z5B=?l=c!gxMWTazIh3a%?@z3OZUtnkYn2FLg8KK53omzvyA5oQsN74~KXSdZ!A$rA zWpTGF@)py;55M%6=d7nF>5I$j#nZe$A9O|pGFm`LxbX&^OF~8fiCS8KZ5;S|6BzC| zi`F_v);Z7in*bvvE=u?i33I9Dt>%;pEuoiuq~bS+)i~&k4I41oHe}N>zvl{irTQxl zUi$3$^G6zgZD(;Qy8|3dr6s+kKU(#7uUC}dW5>3wkETPKYh&NA3JE!?xF-e;)c+Kk6_`vo! zn>J0!nR!3iltpm!WcJLYkN-lAdkQHTXL3W>7YNsR z38DK7En6SqVOGJu1f_=&(&$09>4sz^hS&xe)CNRAoW^Az+zz!Ke7hx;Zd?x%EBT#oIOLh0|sus0B9u_T1zETh$j?}4jr z`NC}M5*f`)&!fsrN|4ABPj8i`Alc8(uJVy#fEd$WY%*VE710=zl|tN|^;5X{a;Gok zv&2VwGuZ)$?i;bxJJns;{%bg`2V;7bsp4WSgUSMT0?$N`c zDRC*AfD+%I_lWhk%k!HLfrX7?0zW7q@xYez9vOR%A0l^43_mkv%CtHHzbjF2ny*T~ z8y{q1LJ(FMCo!JGCsAa|q1^?dP`b91cpuctg}^WlKN6)#hO-s}JC35~T|S(s1CnSF z?A1|ORyZn@STXJwQe>v7n1e5ZYbhiNYUXYzKCH5#22#q_T;!>qpqLH=3FO-W4#QIOxX&%iB%h(sbM0G6rK9Iw&&_n5 zO69A5eKe<)vcJI5)~CRh>(}N1t6sP2;ILGz>Q{*$X)WkAUkMsJu!Rd+*3!Z;xPq)> zOMYdd;LubYkDb$&m%4qQDyMA;M{Bx!wMeCAxbU;#KaYErldSvGZma*vr$30uRg}hx zU2GJ9(vOv#G~1(fLSrWY%Oqymft_BltbpbxSCgWy#5*&}ljZ{guK=+=NhlG!6^$r7@ z?pr~DbV%?0!ow#2g?{%QI9WV;kTee9Iwfp0In=fC$G0A(pw{HRNb8(}_pN8K^0;Wd zW`4FxqTp{jd6_EV`8{9ZHN4XeAy5&5{0;biNp;#0TI?T>ii2RHyZ+OHHbT>%{EdE2EH zxaSJ)xr35$EXW}~X>@%{K-aY_0AnvpX*{c|-tjHK6t{+>3t5VQyF`dwGa8>PLE>P6 z$$t{Y%#fyo$dj@O7CFim%4anO(T-AhSANF;^2wG{)Y|fu$@ACb-oF#;Kbdo2|s0yJvt9E z=);Ex4Z;E_lUNIJ(Ty_%vp^M|8|%=Zc$#Y7!MYb6zT%`@E-E&Aze{Twe3Xw6T6NW-kt(i!}=bcOQQ~ zyfr}vUKy2`PPJ~UXebu6mxSOBb5l@T zUNip9CK5|`dBVTmM8@XWxCUoT@I5_ouIuhcR}+~_3N8rzy(IrhZhmdSKKYI>9P9Xx{#j56gp3F5q${8bn-K` zYWgrweHg~jEO2q^S?-)Yg#1ilTtD|vmN5qDjiCvbM)kjTw+vBh-$dkU>|<=Fba9bi zM7h216AcrM(33OAZR#GNQM8UqsNRNV(i)KnVht#C)ghH37|=SjeiprAPN~N|V2oe= ztV7R|=FkJc#O0c_uJR)pK|1DenC^3waR|BgG{QPh9rwJt%K{>pl5#lAxY9S|jz)9| zhw!etMLgE$6Kon;OJa|&mp2qr08U9Q+~+gCH54;!oXURX>Mi)=nDwi2ORn5~D98Gt zOawIiFf}lL}$QD{}&(*dJJOW|Dy1Jt-Jps?|+u4|Ey5|x8NIO;r(yy z8-&sS;_vmpW7YqkcB!E6_8;-CI4ySc$p3?Q_rIp8GS%zX|8puQC^!TV8W#R{n)(m; zPE1NpNli=7$ov<47yO;3mX?+On~bfiZ)p6FRL*}+Q~yonj89BX{RQ9i|1nKn|C@~6 z+5HQ?|DC3S!1wjdKdGF@r{|YRLOd8KC=6EF8L^t2Mo1tEIHhl{e{Tp~p_rOmO z1n*1$P%5Ve779}n_cXW1|1-WqSx4cHpdT~|Lz8D~U6HgZlh~A{Er@#q6)=J#GF)## z)6_gs{|QE<()r>Pdg21u#xf;wXcXE6TTm)T80M4kvt-i+0Tk{Z2}a$f%5`oQ3`>~T zrYbEqJA#Og@n)NpR3y=Xd#M)2Z-Gd-U)U#F?BbEWp-}&9nXNq>2P2pK(|n_2G?DWo z9ACPv;dH)4pV>wU zLz#)LaWO2}s;W=D_OfAfO2RlwEH5o8uT`_ zC-U>DXV1xaU6-k;YpvzTTTIfmhkQKa+lLym>3(109A~kvb*!JaU$xi&zPag!uyME^nZK}~ zM&p>apR5d^?Yz_VvHp0F-7e)D+i;m?ndJ}JPD$UkvFXw?lE>$IKCh2{JTem_@X^DJ z&*blW8V9mB>Gx;>Ve&g*l&)PNE~Gvshx9<+7P|JvCqtk=cPEos`a9QV;?RRqIdlJ4 z0M8mQ>&P92SN#okT{|ng78i_oW}vcge(k>-&Pny;Ot32Fsxy3j>Mp?>2QZo zde;wkFt~oMITqNX9an z@r-Cp;|T+&M%|37jk}>C9BI?QIllkKes*k)0Q0yS`0)`nv-=|;11ZSTxUP_+DIX#e zslG*u#*dB+O(P>INjy>#G?qjI8!s6V=4n!hMIa$3-LXdof-+a5yoVwW=tFFYGL`2z z-TvMsxzd*Q_#h>k`JPo)laSVQrszCZ4=q*#oAH?>69{?D zU1oEhS4`*NED=xltnipBYyv#}S<6yUF&+>k=rY?$%NO1Yq3_(D1UK1GhqCS*OO)tD z6L(7JoD!c!0N*zc=cYuV5oM2>I5McRBz?{=)_hLHB8{FhXPRNO{F^3p=OkvIy$O3_eIlB zFp3jDtKbn(x*f7A)2J-6ryVL9iY!I3Kd+o?)Bbu_o0iqC`t-tD+4@UPgeszefNE5! z7E)%?l{RdlAVHxEn82*hBBa%#*D*0@>+jnFyt7}EEau@7d@S=`# z4;EZ~pD<&yuH$)9MrVWzg`+6PH*`6(@{O^(G?W@vC1Ta_rI=jg4&zjFV-oW}ap+sT zv2(&Hc`=LeJgO#{6>x3U+lHgOw(bU+Qb0X(S2o)z4gUGdlN)r9=e*oDQwqdxCNx>_ znkp)C_?kIpFpU3sENDY(`d2gknO$bwXqQ^Lxo=eS;y&6UMT0BVM#HpGY<(*PD>T-! zp0#xcgkd~Wx3pi?^-rJ+tKt57UOpZ*J>{w5tCp9%p{DGyo$~2c;}~w#+A^58$!1nx zSg(6FEw-WZ>Q2iwYr$p`F_}>AWv=(zu#W4tpHgRczgr4#)bJ$kfW(G7ybKVJc8AX?WBoqv&?Dz8 zOk1w;G^AXn#ud4=Uk>w_3%sY`R`M4HFaRW1J7tsRTG9{b^D*ekIuu_cN zdY+}uo4)_^r$YhKS(y6N3tj1}e=y;~lrz`iOsb4K*-lXqJJrWrcC*tD?OovY+SSS` znrGMwaxXgF>ke{s;=Rmj-+QTN74(2lA@G8ybjXT(`1d8g(0eD9vMrp+zrJhTlh10{{=c->lcoRPN>8e<_4{kj)L&?ge8(sK< zY`ZYC&9rUz{@NVxIlX=D_jXr3>xlGJ@r#D;>35l?!Y^_CjBoUqYu^l0oA#wcMCQ4% zT-fMsy!5kb{va_Nb3Vg&*URhsmcoABZfAYxGID-ZVbzmqCKP{45P7%uPtdm{wL*OX zNPho{*MAryeip<$-G>BY*H9K%Dc|C3b|!!o_<^b-g3wb(07wLK_fF&EQR-KMY```O zcrpW)a03E)MD;&|!z;)qgD6LXD;Pn9hJBR=c?|enJjggeMuM?`eKYfIwD(D7WrTDl z06NHd=hjk8crDDv3KhpOHP{84M^vX4Y2(s@TDXOX(}T~VcdD>~Ix~bWrhirCgx*JZ z&$lYtb~!*(DsYGjc8D}nIBsK@Q9l=Zr*wt)Mu*dtbMwblC}=T1Scj8XSRb|_MOJq) z$aj#qb*(UnPg91GCV|ydBC+yWjL3#R=7neYgonsNO+tz?w<1ABc%qnP%;1P3)HVN= z=yWxrZ7$b|a1sqr_(PjQi}<%KgvW-*7!JRfK*U%F%=lg;5{>Yu422j!J=h0q*oeJ1j2K9U21pCB7(NkLj&iV$>`0I9D2Ln=f?-$(muQLO_zLfcJTOQ? z+emT&nTh|F3;XChd(%z$s0Ig#kPn58r4Ww8^NtZ2gas*%3$&1iK#@0gLJiqPxwvp1 z>5mh+kpt)k@;F2QQ#Sv2aTYm&$U^A8yg<9ozm^QauKqzl^ z34^V$buKiTni7exStAwnlUFHywV9cX_nM|cX^mMzpg3eYhjpk)W!)y6#dw){@QUf9 zP%l=T$R?f68Jnf4b;pT#FeiX&nQpiV2;<3+huEE_6lBYpZ{L}8)+wG87zpv%JS#Y# z>DgGKQ=WU)W2MHO_UQ+BC^(gab^&Uf|COKr>1JmonEr{L7g@|q3Cy^D)WWx88{9Cp5dvWi5Q~p*^iY&E(v;~Y_udQ zx`r$oelMDBsJNc1mzLA$JNo%aPZ@|G+JNiFn;hz(QR$>Q(<3wWF*$ZMgcW%Htq- z`kB32tCU1ftE!oFNuzKTR7#pk*0rq6Dlu01Tv69_>J+SS7o=}Bn>56&sn<*5T1#!p zc&kuHo@%AGl%c44eyH{d3>B#sw>h?_r+p`&#=5SGm#O0Fruk~FvYMFuT7hkNd5B7- zoq2(inyT69m?h(1#ut--i8b?Ps1GX%4C`M_TC5d|HCZ{Z8Jh`Kd95GIDi}MmuvMh( z=%t!RY$B_&moQu{tFRV(jWH{;GmEJv8>K~rvkj}WnxK9)>yFl#K(*>xa0rwW#}3fWclxI6!uqGB8RnPU5@Hs-BDxGZEw zt&dQyYkRUDOPZ7esc9=;bX&J#*|q&-h;0Z);nt;%z-RF)vqQ6>Y$uxAYPgJmxc!&3 zWUHcg_d$IdxoCB{Z0bn55{Ql}w})GJkF%=^x0gE#e7dKne2{F9Wx51-WPymVyoI!_ zJ6IZIy9Qf9lFMQp>mtJlyPlzq|= z#eXYWT1)`$%BcFge}<^TP#0HVSH`S!#);>Lt9Wie0Rh{!Z4anAyX*3*o5jFyoT+_;7?98}^mqzo6{X)&v>^`bIiP_7|aquH!Y@^4V zQ@2N?x`w{!YzHkm#>~f8t>=Shw~s|T&x0U2H>_5A?6m;w&wGH&6%)`z3TA>^&<5SY zDf%!9t<44YjE(%zj9{wARK3)Od_NeL7@g5cc1(AypDf$Z^BdBUHPUZ{p%v^(xw^tC zeY4MVoJkvx7N}kiw!q>%SsM9IU(-yieJ*q`OZC<{5on!M#o#}Z-O^rn; zT}%CS$a>GzGJ6)A%~MSWpZ3%90)_ug4WFdvi$(p_&wPd~z12Z2)JxhR1VF$ov(sw5 z!}D3pHA&V&JvJAr#WqXPQ+(G?jGH<-*1!3|P+C#o2iS`T&p1uk&xW3dEkwmsHYfaw zXAD)z+{}_qr{l(XS}NEq49>o##iZqsp6!x@?Fg}0*>XdYU9v|Hv&31v*@cL~vN+0% z@Y>2b&allt4#UI%aK;2sj=C+9$=wJFSl87l+aA)x2xHt!9FWW1i#bfPrMlhFtty!` z+|+&C{LI@Juv}D}*CS(wv*7ZB# z|Lk#o9NVpz;RCSY9|uYvE-O{N%^lRXcR0Exesrg&;?@^`MP}TJsLdfx<3sI!1b*XJ zx~dd2bvur3Z3WmLY=N-Ulch!EB~HtxJ%35A;oI_hTnRVC?Mw)D-A5_qutw#FP`KcD z<)E15E57A{N1QV?*XKmTu^X9~=m^|GS!xcH4x{DMl7$_OY0PDft9?dxZd|gW2uO{) zC5@d$gM@G{%&i1gg^oste&;W|){73k5-kH$m%h~wFFR-3qWx7u|H@(Dcc!sR%BNdMO_~&Zq?akO@);6-CiH_xd z?e7jiM2Vc=?&?7d?}IMy$zIs?{(#qqoC+rG`~I-;u0^D->sn0%Ah2{EG3`ePSM2KSy`unCsOhs zFE@BtB=HR3F+cOp{fqW=gB`!}GV4P-kEe!;S7=>%GsI6W9j)1k^!TaIzMZB_pOj92 z%JC-kQ*TxH94-G6xb!E5zc4x8w(RlZ>BbW)_J8X0WiQyLvXFaIb9a?)Dvgx(xNg4~ z_dxH<6@5jH=k6dZv^?3ec&BTAFUO|r(T6hlymj7UYwCz^(Uxbo%X#()ZCd5Z(QE7W z^a%RMT#L>)nwsC}oNrsQ?)6yneE zZ`n@ulc$gR;22mx|M_YE(v3Us??q$SV!5${1BPqw~y%RSd^3eUAY2>6l`|Y=lLagERui7qRjR0PyGm-t&xBJKmZU9iACcPnN%*B zP3IFDl}`VuRqGX-)o!_6?=uhJfBiX>vjOJ}fA6PbT z^>qLG6lzo?Ka;izDssujZCYXails~fP;8JR=8ZL+)mgD+GrN5pd2(CiQ19N1iU`QXP2F)m zGEsSU?c0mvjN4mi63YXvUth}|eR}1!qJwuGwY+(}>gm@vt62V=@9w4djBnn)00oqa zkpT=maDY?z!;dVg27D00vkFAe!0U!u5GhiqTZFI(J^YX-R4A-)I^*_p&=l_iK@tDN z7G1O`m=O~!ai1j8!bL@&;ChkA9s$Z`#C{-R&@JFhV^7B&eO!{sY2XP`Es}_YCdY2H zY7)yVy?T;DWAan2t7oe0LKr2rJQGb)fP^w6WcvGMOfLRQlg>IhOwEl0o48ROB)7S{KsQnXW}GE#_EVx8r*Swxjx z2v|ZdCC1igwcYmEX{ikr7-Uz>mfLc7mG0U-y$v_rNvky%UQ_#P7t~n6y>tK7c=g>^ z7LVFBwv_q)AeY~Q&t2DEdKD2ETZ9dMIKoP;i1tK>J(1W-iy_{)zBt>wS1^SU!MNdf zIX;;_j~g;MWP(R#_)nKno;kbOiqyDUikB6c=AO61&Qh6!7Izn+KlT~vh+nB#Rh|dV zH)P3@#X0G!MRoI0@ffWbYF6^m4CtAtwwi3(Ry&nsr!xk5vapL9yX?3j6Bdd??L7oz zqnnnS@9l8anqlCsVOt7J)b1N`AO8k?Ya3I6dkKgVpIpE#$plfEF%S#89aL`@1?R3upEXOaM(WWH|?+$;R_r^NR;y2kzPo4k9$ffQbdHSvd z9yYycG_mc~g(Nxo>O)Udqvht9Re0KQpWgUStsg(V))^N-Rr16u~UaJ-?7E-OT8AP8Bwlg33*HEz=&3AaGB7Tz!_;~`-206{}1=*e-; z!{HDqlEGN`<52ONp%Rn8rk0tHh*7Kvb}+{%_w{gvN^s%{yTU{%esNJW^y16Jr$q!w z%P3_8;~M+q!xxIegiPqt&A{P9Hs0|eW6Xsfg)m12Wyp?r6lDKo0#-ZKJ;9Hs;^80} zNsTL#Q8*YpUzCmjL~&Fqj*Zl0FG|P35-yMuU96%2y=Svco)S(O6yo|oX~;`jrg_S$ z5h`iPi%+f6fxU7i_@Z~g@Tswu!TdxZak)NE38y<_X^kg^*~~q zn<;wZRzML-;EeJ-yQ3otvsuny$tq0Xw9GYE#!Pd*vpNfWr8#s%PJBcYp81@MJbBqm zO}I0K5Y*>DBlDbj3P46PVjVNhAx?rmbQl}L!xnwT1kAA$q8Swvd-8O~xSa4~Bm}5N zNxF@TnX(9$WXrEGh!=;J^r0hV-uKQ|z!6sJr7+DWK6C%cKLFZuoHWd-PR(h%KcKXf z=#%LpLq=4m=8UIxqiIc->Q5sWsH#uJW<_yx!0D`#gj>|APM6mSV4@L!DNW>9VG0Yb z9>$ZfwA)(QN>X&%@F;N05nbyF(F^HW(+$F>G^&6ZTwyNKgu)p{rG(eEU!fQn!y?{f zRK3s}^#&9?Djsow7dBwl#TO|lmT`uaN@4#l6viap@mM2?AsmY>u|u{jDBW{pB|lkp zvD7hplMH1D-{rtTWU-a!C*57#rMFJzG9#DENwHMn$73GieZPF>scl)!z2ojP-5lrF z>c^{vth1djJ4-6vh0n_MU^KuSXzEV2pnEnnYV+J?MK2n+eoAzsAx+CUrxMYX?lW{T zO=(S^lGCZIrk6jB$}$t$c%{xJH*%|LSEK(ieW!l)t3AW!S*yC%p}uvcx7+Gn^V&J( z{WYWwBJ5!=x+CsAHkXUN&1C}{9=d)uvn%84X{))~m%g@@v#srID;ZBn{dSytJ?L>` zyWHYFH@9=0>2*Ij-R6EbxZ^GFdGEO0_P)1@Ijin{54pX5{QR^a)Tv%|t6v@K zS=aj3x!!fJe;w>$7yH=BUUsve9qm(TSNq!8-gdXY9qw_L``qbXce~#m?|Ikz-ud2l zzyBTZffxMX314`_A0F|ESN!4`-+0GA9`ccw{NyQLdCOlO^O@KD<~iSa&wn2Dp%;C{ Hg8%?KWrL-> literal 0 HcmV?d00001 diff --git a/final.html b/final.html new file mode 100644 index 0000000..352077f --- /dev/null +++ b/final.html @@ -0,0 +1,9670 @@ + + + + + + + + + + +Predicting New Construction in Philadelphia + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

+ +
+ +
+
+

Predicting New Construction in Philadelphia

+

MUSA 508 Final Project

+
+ + + +
+ +
+
Author
+
+

Laura Frances and Nissim Lebovits

+
+
+ +
+
Published
+
+

December 3, 2023

+
+
+ + +
+ + +
+ +
+

1 Summary

+
+
+

2 Introduction

+
+
+Show the code +
required_packages <- c("tidyverse", "sf", "acs", "tidycensus", "sfdep", "kableExtra", "conflicted",
+                       "gganimate", "tmap", "gifski", "transformr", "ggpubr", "randomForest", "janitor",
+                       'igraph')
+install_and_load_packages(required_packages)
+
+source("utils/viz_utils.R")
+
+select <- dplyr::select
+filter <- dplyr::filter
+lag <- dplyr::lag
+
+options(scipen = 999, tigris_use_cache = TRUE, tigris_class = 'sf')
+
+crs <- 'epsg:2272'
+
+
+
    +

  1. Predict development pressure: how do we define “a lot of development”?

    2. +

  2. Define affordability burden: how do we define “affordability burden”? – % change year over year in population that is experience rate burden (will probably see extreme tipping points), growing population, % change in area incomes

    3. +

  3. Identify problem zoning

    4. +

  4. Calculate number of connected parcels

    5. +

  5. Predict development pressure at the block level

    6. +

  6. Identify not burdened areas

    7. +

  7. Identify problem zoning

    8. +

  8. Calcualte number of connected parcels

    9. +

  9. Advocate for upzoning in parcels where there is local development pressure, no affordability burden, problem zoning, and high number of connected parcels

    10. +
+
    +
  • transit
    • +
  • zoning (OSM)
    • +
  • land use (OSM)
    • +
+
+
+Show the code +
urls <- c(
+  phl = "https://opendata.arcgis.com/datasets/405ec3da942d4e20869d4e1449a2be48_0.geojson",
+  phl_bgs = "https://opendata.arcgis.com/datasets/2f982bada233478ea0100528227febce_0.geojson",
+  nbhoods = "https://raw.githubusercontent.com/opendataphilly/open-geo-data/master/Neighborhoods_Philadelphia/Neighborhoods_Philadelphia.geojson",
+  roads = 'https://opendata.arcgis.com/datasets/261eeb49dfd44ccb8a4b6a0af830fdc8_0.geojson',
+  historic_districts = "https://phl.carto.com/api/v2/sql?q=SELECT+*+FROM+historicdistricts_local&filename=historicdistricts_local&format=geojson&skipfields=cartodb_id",
+  zoning = "https://opendata.arcgis.com/datasets/0bdb0b5f13774c03abf8dc2f1aa01693_0.geojson",
+  overlays = "https://opendata.arcgis.com/datasets/04fd29a8c022471994900cb0fd791bfc_0.geojson",
+  opa = "data/opa_properties.geojson",
+  building_permits = building_permits_path,
+  acs14 = acs_vars14_path,
+  acs19 = acs_vars19_path,
+  trolley_stops = "data/Trolley_Stations.geojson",
+  subway_stops = "data/Highspeed_Stations.geojson"
+)
+
+suppressMessages({
+  invisible(
+    imap(urls, ~ assign(.y, phl_spat_read(.x), envir = .GlobalEnv))
+  )
+})
+
+phl_bgs <- phl_bgs %>%
+              select(GEOID10)
+
+broad_and_market <- roads %>% filter(ST_NAME %in% c('BROAD',"MARKET") | SEG_ID %in% c(440370, 421347,421338,421337,422413,423051,440403,440402,440391,440380))
+
+subway_stops <- subway_stops[phl, ]
+transit_stops <- st_union(trolley_stops, subway_stops)
+
+historic_districts <- historic_districts %>%
+                        mutate(hist_dist = name) %>%
+                        select(hist_dist)
+
+nbhoods <- nbhoods %>%
+            select(mapname)
+
+overlays <- overlays %>% clean_names()
+overlays$overlay_symbol <- gsub("/", "", overlays$overlay_symbol) 
+overlays <- overlays %>%
+              mutate(overlay_symbol = ifelse(overlay_symbol == "[NA]", "Other", overlay_symbol))
+                        
+building_permits <- building_permits %>%
+                      filter(permittype %in% c("RESIDENTIAL BUILDING", "BP_ADDITION", "BP_NEWCNST"))
+
+acs_reg_vars <- c(
+  "GEOID",
+  "Total_Pop", 
+  "Med_Inc",
+  "Percent_Nonwhite",
+  "Percent_Renters",
+  "Rent_Burden",
+  "Ext_Rent_Burden"
+  )
+
+acs14_reg_vars <- phl_spat_read(acs_vars14_path) %>%
+                    st_drop_geometry() %>%
+                    select(acs_reg_vars) %>%
+                    clean_names() %>%
+                    rename(GEOID10 = geoid)
+
+acs19_reg_vars <- phl_spat_read(acs_vars19_path) %>%
+                    st_drop_geometry() %>%
+                    select(acs_reg_vars) %>%
+                    clean_names() %>%
+                    rename(GEOID10 = geoid)
+
+acs_tot <- acs14_reg_vars %>%
+             left_join(acs19_reg_vars, by = "GEOID10", suffix = c("_14", "_19")) %>%
+             mutate(pct_med_inc_change = med_inc_19 - med_inc_14 / med_inc_14,
+                    pct_tot_pop_change = total_pop_19 - total_pop_14 / total_pop_14,
+                    pct_renters_change = percent_renters_19 - percent_renters_14,
+                    pct_rent_burdenchange = rent_burden_19 - rent_burden_14)
+
+
+
+
+Show the code +
ggplot(phl_bgs) +
+  geom_sf() +
+  theme_void()
+
+
+

+
+
+
+
+Show the code +
phl_bgs <- phl_spat_read("https://opendata.arcgis.com/datasets/2f982bada233478ea0100528227febce_0.geojson") %>%
+              select(GEOID10)
+
+# Create a complete grid of GEOID10 and year
+geoid_years <- expand.grid(GEOID10 = unique(phl_bgs$GEOID10),
+                           year = unique(building_permits$year))
+
+
+
+# Joining your existing data with the complete grid
+permits_bg <- st_join(phl_bgs, building_permits) %>%
+              group_by(GEOID10, year) %>%
+              summarize(permits_count = sum(permits_count, na.rm = TRUE)) %>%
+              st_drop_geometry() %>%
+              right_join(geoid_years, by = c("GEOID10", "year")) %>%
+              replace_na(list(permits_count = 0)) %>%
+              left_join(phl_bgs, by = "GEOID10") %>%
+              st_as_sf() 
+
+
+### spat + temp lags----------------------------------
+suppressMessages(
+permits_bg <- permits_bg %>%
+              group_by(year) %>%
+              mutate(nb = st_knn(geometry, 5),
+                         wt = st_weights(nb),
+                      lag_spatial = st_lag(permits_count, nb, wt)) %>% # calculate spat lag
+              ungroup()%>%
+              arrange(GEOID10, year) %>% # calculate time lag
+               mutate(
+                 lag_1_year = lag(permits_count, 1),
+                 lag_2_years = lag(permits_count, 2),
+                 lag_3_years = lag(permits_count, 3),
+                 lag_4_years = lag(permits_count, 4),
+                 lag_5_years = lag(permits_count, 5),
+                 lag_6_years = lag(permits_count, 6),
+                 lag_7_years = lag(permits_count, 7),
+                 lag_8_years = lag(permits_count, 8),
+                 lag_9_years = lag(permits_count, 9),
+                 dist_to_2022 = 2022 - year
+               )
+)
+
+### distance to transit---------------------------------
+phl_bgs <- phl_bgs %>%
+  rowwise() %>%
+  mutate(
+    dist_to_transit = as.numeric(min(st_distance(st_point_on_surface(geometry), transit_stops$geometry)))
+  ) %>%
+  ungroup()
+
+### historic dists---------------------------------
+hist_dists_x_bg <- st_join(phl_bgs, historic_districts) %>%
+              mutate(hist_dist = ifelse(is.na(hist_dist), 0, 1))
+
+hist_dists_x_bg <- st_join(phl_bgs, historic_districts) %>%
+                      st_drop_geometry() %>%
+                      mutate(hist_dist_present = 1) %>%
+                      group_by(GEOID10, hist_dist) %>%
+                      summarize(hist_dist_present = max(hist_dist_present, na.rm = TRUE), .groups = 'drop') %>%
+                      pivot_wider(names_from = "hist_dist", values_from = "hist_dist_present", 
+                                  names_prefix = "hist_dist_", values_fill = list("hist_dist_present" = 0))
+
+phl_bgs <- left_join(phl_bgs, hist_dists_x_bg, by = "GEOID10")
+
+### hoods---------------------------------              
+intersection <- st_intersection(phl_bgs, nbhoods)
+intersection$intersection_area <- st_area(intersection)
+max_intersection <- intersection %>%
+  group_by(GEOID10) %>%
+  filter(intersection_area == max(intersection_area)) %>%
+  ungroup() %>%
+  select(GEOID10, mapname) %>%
+  st_drop_geometry()
+
+bgs_w_hood <- left_join(phl_bgs, max_intersection, by = c("GEOID10" = "GEOID10")) %>%
+                st_drop_geometry()
+
+phl_bgs <- left_join(phl_bgs, bgs_w_hood, by = "GEOID10")
+
+### overlays---------------------------------
+overlays_x_bg <- st_join(phl_bgs, overlays %>% select(overlay_symbol)) %>%
+                      st_drop_geometry() %>%
+                      mutate(overlay_present = 1) %>%
+                      group_by(GEOID10, overlay_symbol) %>%
+                      summarize(overlay_present = max(overlay_present, na.rm = TRUE), .groups = 'drop') %>%
+                      pivot_wider(names_from = "overlay_symbol", values_from = "overlay_present", 
+                                  names_prefix = "overlay_", values_fill = list("overlay_present" = 0))
+
+phl_bgs <- left_join(phl_bgs, overlays_x_bg, by = "GEOID10")
+
+
+### join back together----------------------
+permits_bg <- left_join(permits_bg,
+                     st_drop_geometry(phl_bgs),
+                     by = "GEOID10")
+
+### acs vars--------------------------------------------
+
+permits_pre_2019 <- filter(permits_bg, year < 2019)
+permits_post_2018 <- filter(permits_bg, year >= 2019)
+
+
+permits_joined_pre_2019 <- left_join(permits_pre_2019, acs14_reg_vars, by = "GEOID10")
+permits_joined_post_2018 <- left_join(permits_post_2018, acs19_reg_vars, by = "GEOID10")
+
+
+permits_bg <- bind_rows(permits_joined_pre_2019, permits_joined_post_2018)
+
+### clean-----------------------------------------------------
+permits_bg <- permits_bg %>%
+                select(-c(nb, wt, GEOID10)) %>%
+                clean_names()
+
+
+
+
+Show the code +
tm <- tm_shape(permits_bg %>% filter(year != 2012)) +
+        tm_polygons(col = "permits_count", border.alpha = 0, palette = 'viridis', style = "fisher", colorNA = "lightgrey") +
+  tm_facets(along = "year") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE) 
+
+suppressMessages(
+tmap_animation(tm, "assets/permits_animation.gif", delay = 50)
+)
+
+
+
+
+

+
Philadelphia Building Permits, 2013 - 2022
+
+
+
+
+

3 Methods

+
+

3.1 Data

+
+
+Show the code +
ggplot(building_permits, aes(x = as.factor(year))) +
+  geom_bar() +
+  labs(title = "Permits per Year")
+
+
+

+
+
+Show the code +
ggplot(permits_bg %>% st_drop_geometry, aes(x = permits_count)) +
+  geom_histogram() +
+  labs(title = "Permits per Block Group per Year",
+       subtitle = "Log-Transformed") +
+  scale_x_log10() +
+  facet_wrap(~year)
+
+
+

+
+
+Show the code +
# 
+# tm_shape(permits_bg) +
+#   tm_polygons("spat_lag_permits", border.alpha = 0, palette = 'viridis', style = "fisher", colorNA = "lightgrey") +
+#   tm_facets(along = "year") +
+# tm_shape(broad_and_market) +
+#   tm_lines(col = "lightgrey") +
+#   tm_layout(frame = FALSE)
+
+
+
+

3.1.1 Feature Engineering

+
+
+Show the code +
permits_bg_long <- permits_bg %>%
+                    st_drop_geometry() %>%
+                    pivot_longer(
+                      cols = c(starts_with("lag"), dist_to_2022),
+                      names_to = "Variable",
+                      values_to = "Value"
+                    )
+
+
+ggscatter(permits_bg_long, x = "permits_count", y = "Value", facet.by = "Variable",
+   add = "reg.line",
+   add.params = list(color = "blue", fill = "lightgray"),
+   conf.int = TRUE
+   ) + stat_cor(method = "pearson", p.accuracy = 0.001, r.accuracy = 0.01)
+
+
+

+
+
+
+
+
+

3.2 Models

+
+
+
+

4 Results

+
+

4.1 OLS Regression

+

To begin, we run a simple regression incorporating three engineered groups of features: space lag, time lag, and distance to 2022. We include this last variable because of a Philadelphia tax abatement policy that led to a significant increase in residential development in the years immediately before 2022. We will use this as a baseline model to compare to our more complex models.

+
+
+Show the code +
permits_train <- filter(permits_bg %>% select(-mapname), year < 2022)
+permits_test <- filter(permits_bg %>% select(-mapname), year == 2022)
+
+reg <- lm(permits_count ~ ., data = st_drop_geometry(permits_train))
+
+predictions <- predict(reg, permits_test)
+predictions <- cbind(permits_test, predictions)
+
+predictions <- predictions %>%
+                  mutate(abs_error = abs(permits_count - predictions),
+                         pct_error = abs_error / permits_count)
+
+ggplot(predictions, aes(x = permits_count, y = predictions)) +
+  geom_point() +
+  labs(title = "Predicted vs. Actual Permits",
+       subtitle = "2022") +
+  geom_smooth(method = "lm", se = FALSE)
+
+
+

+
+
+Show the code +
mae <- paste0("MAE: ", round(mean(predictions$abs_error, na.rm = TRUE), 2))
+
+tm_shape(predictions) +
+        tm_polygons(col = "abs_error", border.alpha = 0, palette = 'viridis', style = "fisher", colorNA = "lightgrey") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+

+
+
+
+
+

4.2 Poisson Model

+

Given that we are dealing with a skewed distribution of count data

+
+
+Show the code +
poisson_reg <- glm(permits_count ~ ., 
+                   family = poisson(link = "log"),
+                   data = st_drop_geometry(permits_train))
+
+poisson_predictions <- predict(poisson_reg, permits_test, type = "response")
+poisson_predictions <- cbind(permits_test, poisson_predictions)
+
+poisson_predictions <- poisson_predictions %>%
+                           mutate(abs_error = abs(permits_count - poisson_predictions),
+                                  pct_error = abs_error / permits_count)
+
+ggplot(poisson_predictions, aes(x = permits_count, y = poisson_predictions)) +
+  geom_point() +
+  labs(title = "Predicted vs. Actual Permits",
+       subtitle = "2022") +
+  geom_smooth(method = "lm", se = FALSE)
+
+
+

+
+
+Show the code +
poisson_mae <- paste0("MAE: ", round(mean(poisson_predictions$abs_error, na.rm = TRUE), 2))
+
+tm_shape(poisson_predictions) +
+        tm_polygons(col = "abs_error", border.alpha = 0, palette = 'viridis', style = "fisher", colorNA = "lightgrey") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+

+
+
+

We find that our OLS model has an MAE of only MAE: 2.44–not bad for such a simple model! Still, it struggles most in the areas where we most need it to succeed, so we will try to introduce better variables and apply a more complex model to improve our predictions.

+
+
+

4.3 Random Forest Regression

+
+
+Show the code +
rf <- randomForest(permits_count ~ ., 
+                   data = st_drop_geometry(permits_train),
+                   importance = TRUE, 
+                   na.action = na.omit)
+
+rf_predictions <- predict(rf, permits_test)
+rf_predictions <- cbind(permits_test, rf_predictions)
+rf_predictions <- rf_predictions %>%
+                  mutate(abs_error = abs(permits_count - rf_predictions),
+                         pct_error = abs_error / (permits_count + 0.0001))
+
+ggplot(rf_predictions, aes(x = permits_count, y = rf_predictions)) +
+  geom_point() +
+  labs(title = "Predicted vs. Actual Permits",
+       subtitle = "2022") +
+  geom_smooth(method = "lm", se = FALSE)
+
+
+

+
+
+Show the code +
rf_mae <- paste0("MAE: ", round(mean(rf_predictions$abs_error, na.rm = TRUE), 2))
+
+tm_shape(rf_predictions) +
+        tm_polygons(col = "abs_error", border.alpha = 0, palette = 'viridis', style = "fisher", colorNA = "lightgrey") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+

+
+
+
+
+
+

5 Discussion

+
+

5.1 Accuracy

+

Predominately, our model overpredicts, which is better than underpredicting, as it facilitates new development.

+
+
+Show the code +
nbins <- as.integer(sqrt(nrow(rf_predictions)))
+vline <- mean(rf_predictions$abs_error, na.rm = TRUE)
+
+ggplot(rf_predictions, aes(x = abs_error)) +
+  geom_histogram(bins = nbins) +
+  geom_vline(aes(xintercept = vline))
+
+
+

+
+
+
+
+

5.2 Generalizabiltiy

+
+
+Show the code +
tm_shape(acs19) +
+        tm_polygons(col = "Percent_Nonwhite", border.alpha = 0, palette = 'viridis', style = "fisher", colorNA = "lightgrey") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+

+
+
+Show the code +
rf_predictions <- rf_predictions %>%
+                      mutate(race_comp = case_when(
+                        percent_nonwhite >= .50 ~ "Majority Non-White",
+                        TRUE ~ "Majority White"
+                      ))
+
+ggplot(rf_predictions, aes(y = abs_error, color = race_comp)) +
+  geom_boxplot(fill = NA)
+
+
+

+
+
+

How does this generalize across neighborhoods?

+

How does this generalize across council districts?

+
+
+

5.3 Assessing Upzoning Needs

+

We can identify conflict between projected development and current zoning.

+

Look at zoning that is industrial or residential single family in areas that our model suggests are high development risk for 2023:

+
+
+Show the code +
filtered_zoning <- zoning %>%
+                     filter(str_detect(CODE, "RS") | str_detect(CODE, "I"))
+
+
+tm_shape(filtered_zoning) +
+        tm_polygons(col = "CODE", border.alpha = 0, colorNA = "lightgrey") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+

+
+
+

We can extract development predictions at the block level to these parcels and then visualize them by highest need.

+
+
+Show the code +
filtered_zoning <- st_join(
+  filtered_zoning,
+  rf_predictions %>% select(rf_predictions)
+)
+
+tm_shape(filtered_zoning) +
+        tm_polygons(col = "rf_predictions", border.alpha = 0, colorNA = "lightgrey", palette = "viridis", style = "fisher") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+

+
+
+
+
+Show the code +
tmap_mode('view')
+
+filtered_zoning %>%
+  filter(rf_predictions > 10) %>%
+tm_shape() +
+        tm_polygons(col = "CODE", border.alpha = 0, colorNA = "lightgrey",
+                    popup.vars = c('rf_predictions', 'CODE')) +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+ +
+ +
+
+

Furthermore, we can identify properties with high potential for assemblage, which suggests the ability to accomodate high-density, multi-unit housing.

+
+
+Show the code +
nbs <- filtered_zoning %>% 
+  mutate(nb = st_contiguity(geometry))
+
+# Create edge list while handling cases with no neighbors
+edge_list <- tibble::tibble(id = 1:length(nbs$nb), nbs = nbs$nb) %>% 
+  tidyr::unnest(nbs) %>% 
+  filter(nbs != 0)
+
+# Create a graph with a node for each row in filtered_zoning
+g <- make_empty_graph(n = nrow(filtered_zoning))
+V(g)$name <- as.character(1:nrow(filtered_zoning))
+
+# Add edges if they exist
+if (nrow(edge_list) > 0) {
+  edges <- as.matrix(edge_list)
+  g <- add_edges(g, c(t(edges)))
+}
+
+# Calculate the number of contiguous neighbors, handling nodes without neighbors
+n_contiguous <- sapply(V(g)$name, function(node) {
+  if (node %in% edges) {
+    length(neighborhood(g, order = 1, nodes = as.numeric(node))[[1]])
+  } else {
+    1  # Nodes without neighbors count as 1 (themselves)
+  }
+})
+
+filtered_zoning <- filtered_zoning %>%
+                    mutate(n_contig = n_contiguous)
+
+filtered_zoning %>%
+  st_drop_geometry() %>%
+  select(rf_predictions,
+         n_contig,
+         OBJECTID,
+         CODE) %>%
+  filter(rf_predictions > 10,
+         n_contig > 2) %>%
+  kablerize(caption = "Poorly-Zoned Properties with High Development Risk")
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Poorly-Zoned Properties with High Development Risk
rf_predictionsn_contigOBJECTIDCODE
27512.167103604I2
85821.5859761524ICMX
90312.9468361595RSA5
91621.5859781615ICMX
152211.2443742559RSA5
163121.5859732736IRMX
164321.5859762756ICMX
167021.5859772803ICMX
167121.5859732804IRMX
186534.3893033128IRMX
209912.9468363492RSA5
225111.4208043744IRMX
268512.9468364533RSA5
316810.5881735506RSA5
343921.5859766067ICMX
355412.9468366289RSA5
361621.5859736405RSA5
383612.9468376869ICMX
384634.3893036901ICMX
386512.9468366943RSA5
415529.6927037646ICMX
417811.2138337704IRMX
459312.9468368805RSA5
469723.0356339094RSA5
475216.0827359244IRMX
480010.5881739371ICMX
491021.5859739662RSA5
521135.68783310411ICMX
521235.68783310412RSA5
521335.68783310413ICMX
534518.21887310760IRMX
550139.00537311150ICMX
550635.68783311161RSA5
561611.24437311449RSA3
599611.00043412339I2
620611.00043312808RSA5
625811.00043312932ICMX
631216.08273613058ICMX
641223.99983313314IRMX
670322.58767313979I3
670511.00043313981RSA3
672013.26530314019RSA5
6793.211.00043814223I2
679422.58767314224ICMX
680822.58767314257I2
685511.00043314373RSA5
686411.00043314401RSA5
686511.00043314402RSA5
696513.26530314649ICMX
701412.94683814748ICMX
719211.24437315167RSA2
744639.00537315720I2
748711.00043615820RSA5
763713.26530316181RSA5
788539.00537316719RSA5
810515.75340317170ICMX
813213.26530317217ICMX
822015.48503317410RSA5
854313.56143318033RSD3
907513.56143319078RSA3
936211.24437319593RSA3
960721.58597420075ICMX
980212.16710320444RSA5
985413.56143420536RSA2
1037713.26530321529ICMX
1065113.56143322004RSD1
1287729.69270425778RSA5
1292911.24437425868RSA1
1316413.56143326249RSD1
1349410.85073326759RSA5
1359911.24437326935RSA3
1381711.24437327284RSA2
1389211.24437327394RSD3
1411013.26530427776I2
1415816.06680327871IRMX
14719.221.585972828949I2
14719.312.946832828949I2
1472012.94683628950RSA5
+ + +
+
+Show the code +
tmap_mode('view')
+
+filtered_zoning %>%
+  select(rf_predictions,
+         n_contig,
+         OBJECTID,
+         CODE) %>%
+  filter(rf_predictions > 10,
+         n_contig > 2) %>%
+tm_shape() +
+        tm_polygons(col = "rf_predictions", border.alpha = 0, colorNA = "lightgrey", palette = "viridis", style = "fisher",
+                    popup.vars = c('rf_predictions', 'n_contig', 'CODE'), alpha = 0.5) +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+ +
+ +
+
+
+
+
+

6 Conclusion

+
+
+

7 Appendices

+
+ +
+ + +
+ + + + \ No newline at end of file diff --git a/index.html b/index.html new file mode 100644 index 0000000..352077f --- /dev/null +++ b/index.html @@ -0,0 +1,9670 @@ + + + + + + + + + + +Predicting New Construction in Philadelphia + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ +
+ +
+
+

Predicting New Construction in Philadelphia

+

MUSA 508 Final Project

+
+ + + +
+ +
+
Author
+
+

Laura Frances and Nissim Lebovits

+
+
+ +
+
Published
+
+

December 3, 2023

+
+
+ + +
+ + +
+ +
+

1 Summary

+
+
+

2 Introduction

+
+
+Show the code +
required_packages <- c("tidyverse", "sf", "acs", "tidycensus", "sfdep", "kableExtra", "conflicted",
+                       "gganimate", "tmap", "gifski", "transformr", "ggpubr", "randomForest", "janitor",
+                       'igraph')
+install_and_load_packages(required_packages)
+
+source("utils/viz_utils.R")
+
+select <- dplyr::select
+filter <- dplyr::filter
+lag <- dplyr::lag
+
+options(scipen = 999, tigris_use_cache = TRUE, tigris_class = 'sf')
+
+crs <- 'epsg:2272'
+
+
+
    +

  1. Predict development pressure: how do we define “a lot of development”?

    2. +

  2. Define affordability burden: how do we define “affordability burden”? – % change year over year in population that is experience rate burden (will probably see extreme tipping points), growing population, % change in area incomes

    3. +

  3. Identify problem zoning

    4. +

  4. Calculate number of connected parcels

    5. +

  5. Predict development pressure at the block level

    6. +

  6. Identify not burdened areas

    7. +

  7. Identify problem zoning

    8. +

  8. Calcualte number of connected parcels

    9. +

  9. Advocate for upzoning in parcels where there is local development pressure, no affordability burden, problem zoning, and high number of connected parcels

    10. +
+
    +
  • transit
    • +
  • zoning (OSM)
    • +
  • land use (OSM)
    • +
+
+
+Show the code +
urls <- c(
+  phl = "https://opendata.arcgis.com/datasets/405ec3da942d4e20869d4e1449a2be48_0.geojson",
+  phl_bgs = "https://opendata.arcgis.com/datasets/2f982bada233478ea0100528227febce_0.geojson",
+  nbhoods = "https://raw.githubusercontent.com/opendataphilly/open-geo-data/master/Neighborhoods_Philadelphia/Neighborhoods_Philadelphia.geojson",
+  roads = 'https://opendata.arcgis.com/datasets/261eeb49dfd44ccb8a4b6a0af830fdc8_0.geojson',
+  historic_districts = "https://phl.carto.com/api/v2/sql?q=SELECT+*+FROM+historicdistricts_local&filename=historicdistricts_local&format=geojson&skipfields=cartodb_id",
+  zoning = "https://opendata.arcgis.com/datasets/0bdb0b5f13774c03abf8dc2f1aa01693_0.geojson",
+  overlays = "https://opendata.arcgis.com/datasets/04fd29a8c022471994900cb0fd791bfc_0.geojson",
+  opa = "data/opa_properties.geojson",
+  building_permits = building_permits_path,
+  acs14 = acs_vars14_path,
+  acs19 = acs_vars19_path,
+  trolley_stops = "data/Trolley_Stations.geojson",
+  subway_stops = "data/Highspeed_Stations.geojson"
+)
+
+suppressMessages({
+  invisible(
+    imap(urls, ~ assign(.y, phl_spat_read(.x), envir = .GlobalEnv))
+  )
+})
+
+phl_bgs <- phl_bgs %>%
+              select(GEOID10)
+
+broad_and_market <- roads %>% filter(ST_NAME %in% c('BROAD',"MARKET") | SEG_ID %in% c(440370, 421347,421338,421337,422413,423051,440403,440402,440391,440380))
+
+subway_stops <- subway_stops[phl, ]
+transit_stops <- st_union(trolley_stops, subway_stops)
+
+historic_districts <- historic_districts %>%
+                        mutate(hist_dist = name) %>%
+                        select(hist_dist)
+
+nbhoods <- nbhoods %>%
+            select(mapname)
+
+overlays <- overlays %>% clean_names()
+overlays$overlay_symbol <- gsub("/", "", overlays$overlay_symbol) 
+overlays <- overlays %>%
+              mutate(overlay_symbol = ifelse(overlay_symbol == "[NA]", "Other", overlay_symbol))
+                        
+building_permits <- building_permits %>%
+                      filter(permittype %in% c("RESIDENTIAL BUILDING", "BP_ADDITION", "BP_NEWCNST"))
+
+acs_reg_vars <- c(
+  "GEOID",
+  "Total_Pop", 
+  "Med_Inc",
+  "Percent_Nonwhite",
+  "Percent_Renters",
+  "Rent_Burden",
+  "Ext_Rent_Burden"
+  )
+
+acs14_reg_vars <- phl_spat_read(acs_vars14_path) %>%
+                    st_drop_geometry() %>%
+                    select(acs_reg_vars) %>%
+                    clean_names() %>%
+                    rename(GEOID10 = geoid)
+
+acs19_reg_vars <- phl_spat_read(acs_vars19_path) %>%
+                    st_drop_geometry() %>%
+                    select(acs_reg_vars) %>%
+                    clean_names() %>%
+                    rename(GEOID10 = geoid)
+
+acs_tot <- acs14_reg_vars %>%
+             left_join(acs19_reg_vars, by = "GEOID10", suffix = c("_14", "_19")) %>%
+             mutate(pct_med_inc_change = med_inc_19 - med_inc_14 / med_inc_14,
+                    pct_tot_pop_change = total_pop_19 - total_pop_14 / total_pop_14,
+                    pct_renters_change = percent_renters_19 - percent_renters_14,
+                    pct_rent_burdenchange = rent_burden_19 - rent_burden_14)
+
+
+
+
+Show the code +
ggplot(phl_bgs) +
+  geom_sf() +
+  theme_void()
+
+
+

+
+
+
+
+Show the code +
phl_bgs <- phl_spat_read("https://opendata.arcgis.com/datasets/2f982bada233478ea0100528227febce_0.geojson") %>%
+              select(GEOID10)
+
+# Create a complete grid of GEOID10 and year
+geoid_years <- expand.grid(GEOID10 = unique(phl_bgs$GEOID10),
+                           year = unique(building_permits$year))
+
+
+
+# Joining your existing data with the complete grid
+permits_bg <- st_join(phl_bgs, building_permits) %>%
+              group_by(GEOID10, year) %>%
+              summarize(permits_count = sum(permits_count, na.rm = TRUE)) %>%
+              st_drop_geometry() %>%
+              right_join(geoid_years, by = c("GEOID10", "year")) %>%
+              replace_na(list(permits_count = 0)) %>%
+              left_join(phl_bgs, by = "GEOID10") %>%
+              st_as_sf() 
+
+
+### spat + temp lags----------------------------------
+suppressMessages(
+permits_bg <- permits_bg %>%
+              group_by(year) %>%
+              mutate(nb = st_knn(geometry, 5),
+                         wt = st_weights(nb),
+                      lag_spatial = st_lag(permits_count, nb, wt)) %>% # calculate spat lag
+              ungroup()%>%
+              arrange(GEOID10, year) %>% # calculate time lag
+               mutate(
+                 lag_1_year = lag(permits_count, 1),
+                 lag_2_years = lag(permits_count, 2),
+                 lag_3_years = lag(permits_count, 3),
+                 lag_4_years = lag(permits_count, 4),
+                 lag_5_years = lag(permits_count, 5),
+                 lag_6_years = lag(permits_count, 6),
+                 lag_7_years = lag(permits_count, 7),
+                 lag_8_years = lag(permits_count, 8),
+                 lag_9_years = lag(permits_count, 9),
+                 dist_to_2022 = 2022 - year
+               )
+)
+
+### distance to transit---------------------------------
+phl_bgs <- phl_bgs %>%
+  rowwise() %>%
+  mutate(
+    dist_to_transit = as.numeric(min(st_distance(st_point_on_surface(geometry), transit_stops$geometry)))
+  ) %>%
+  ungroup()
+
+### historic dists---------------------------------
+hist_dists_x_bg <- st_join(phl_bgs, historic_districts) %>%
+              mutate(hist_dist = ifelse(is.na(hist_dist), 0, 1))
+
+hist_dists_x_bg <- st_join(phl_bgs, historic_districts) %>%
+                      st_drop_geometry() %>%
+                      mutate(hist_dist_present = 1) %>%
+                      group_by(GEOID10, hist_dist) %>%
+                      summarize(hist_dist_present = max(hist_dist_present, na.rm = TRUE), .groups = 'drop') %>%
+                      pivot_wider(names_from = "hist_dist", values_from = "hist_dist_present", 
+                                  names_prefix = "hist_dist_", values_fill = list("hist_dist_present" = 0))
+
+phl_bgs <- left_join(phl_bgs, hist_dists_x_bg, by = "GEOID10")
+
+### hoods---------------------------------              
+intersection <- st_intersection(phl_bgs, nbhoods)
+intersection$intersection_area <- st_area(intersection)
+max_intersection <- intersection %>%
+  group_by(GEOID10) %>%
+  filter(intersection_area == max(intersection_area)) %>%
+  ungroup() %>%
+  select(GEOID10, mapname) %>%
+  st_drop_geometry()
+
+bgs_w_hood <- left_join(phl_bgs, max_intersection, by = c("GEOID10" = "GEOID10")) %>%
+                st_drop_geometry()
+
+phl_bgs <- left_join(phl_bgs, bgs_w_hood, by = "GEOID10")
+
+### overlays---------------------------------
+overlays_x_bg <- st_join(phl_bgs, overlays %>% select(overlay_symbol)) %>%
+                      st_drop_geometry() %>%
+                      mutate(overlay_present = 1) %>%
+                      group_by(GEOID10, overlay_symbol) %>%
+                      summarize(overlay_present = max(overlay_present, na.rm = TRUE), .groups = 'drop') %>%
+                      pivot_wider(names_from = "overlay_symbol", values_from = "overlay_present", 
+                                  names_prefix = "overlay_", values_fill = list("overlay_present" = 0))
+
+phl_bgs <- left_join(phl_bgs, overlays_x_bg, by = "GEOID10")
+
+
+### join back together----------------------
+permits_bg <- left_join(permits_bg,
+                     st_drop_geometry(phl_bgs),
+                     by = "GEOID10")
+
+### acs vars--------------------------------------------
+
+permits_pre_2019 <- filter(permits_bg, year < 2019)
+permits_post_2018 <- filter(permits_bg, year >= 2019)
+
+
+permits_joined_pre_2019 <- left_join(permits_pre_2019, acs14_reg_vars, by = "GEOID10")
+permits_joined_post_2018 <- left_join(permits_post_2018, acs19_reg_vars, by = "GEOID10")
+
+
+permits_bg <- bind_rows(permits_joined_pre_2019, permits_joined_post_2018)
+
+### clean-----------------------------------------------------
+permits_bg <- permits_bg %>%
+                select(-c(nb, wt, GEOID10)) %>%
+                clean_names()
+
+
+
+
+Show the code +
tm <- tm_shape(permits_bg %>% filter(year != 2012)) +
+        tm_polygons(col = "permits_count", border.alpha = 0, palette = 'viridis', style = "fisher", colorNA = "lightgrey") +
+  tm_facets(along = "year") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE) 
+
+suppressMessages(
+tmap_animation(tm, "assets/permits_animation.gif", delay = 50)
+)
+
+
+
+
+

+
Philadelphia Building Permits, 2013 - 2022
+
+
+
+
+

3 Methods

+
+

3.1 Data

+
+
+Show the code +
ggplot(building_permits, aes(x = as.factor(year))) +
+  geom_bar() +
+  labs(title = "Permits per Year")
+
+
+

+
+
+Show the code +
ggplot(permits_bg %>% st_drop_geometry, aes(x = permits_count)) +
+  geom_histogram() +
+  labs(title = "Permits per Block Group per Year",
+       subtitle = "Log-Transformed") +
+  scale_x_log10() +
+  facet_wrap(~year)
+
+
+

+
+
+Show the code +
# 
+# tm_shape(permits_bg) +
+#   tm_polygons("spat_lag_permits", border.alpha = 0, palette = 'viridis', style = "fisher", colorNA = "lightgrey") +
+#   tm_facets(along = "year") +
+# tm_shape(broad_and_market) +
+#   tm_lines(col = "lightgrey") +
+#   tm_layout(frame = FALSE)
+
+
+
+

3.1.1 Feature Engineering

+
+
+Show the code +
permits_bg_long <- permits_bg %>%
+                    st_drop_geometry() %>%
+                    pivot_longer(
+                      cols = c(starts_with("lag"), dist_to_2022),
+                      names_to = "Variable",
+                      values_to = "Value"
+                    )
+
+
+ggscatter(permits_bg_long, x = "permits_count", y = "Value", facet.by = "Variable",
+   add = "reg.line",
+   add.params = list(color = "blue", fill = "lightgray"),
+   conf.int = TRUE
+   ) + stat_cor(method = "pearson", p.accuracy = 0.001, r.accuracy = 0.01)
+
+
+

+
+
+
+
+
+

3.2 Models

+
+
+
+

4 Results

+
+

4.1 OLS Regression

+

To begin, we run a simple regression incorporating three engineered groups of features: space lag, time lag, and distance to 2022. We include this last variable because of a Philadelphia tax abatement policy that led to a significant increase in residential development in the years immediately before 2022. We will use this as a baseline model to compare to our more complex models.

+
+
+Show the code +
permits_train <- filter(permits_bg %>% select(-mapname), year < 2022)
+permits_test <- filter(permits_bg %>% select(-mapname), year == 2022)
+
+reg <- lm(permits_count ~ ., data = st_drop_geometry(permits_train))
+
+predictions <- predict(reg, permits_test)
+predictions <- cbind(permits_test, predictions)
+
+predictions <- predictions %>%
+                  mutate(abs_error = abs(permits_count - predictions),
+                         pct_error = abs_error / permits_count)
+
+ggplot(predictions, aes(x = permits_count, y = predictions)) +
+  geom_point() +
+  labs(title = "Predicted vs. Actual Permits",
+       subtitle = "2022") +
+  geom_smooth(method = "lm", se = FALSE)
+
+
+

+
+
+Show the code +
mae <- paste0("MAE: ", round(mean(predictions$abs_error, na.rm = TRUE), 2))
+
+tm_shape(predictions) +
+        tm_polygons(col = "abs_error", border.alpha = 0, palette = 'viridis', style = "fisher", colorNA = "lightgrey") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+

+
+
+
+
+

4.2 Poisson Model

+

Given that we are dealing with a skewed distribution of count data

+
+
+Show the code +
poisson_reg <- glm(permits_count ~ ., 
+                   family = poisson(link = "log"),
+                   data = st_drop_geometry(permits_train))
+
+poisson_predictions <- predict(poisson_reg, permits_test, type = "response")
+poisson_predictions <- cbind(permits_test, poisson_predictions)
+
+poisson_predictions <- poisson_predictions %>%
+                           mutate(abs_error = abs(permits_count - poisson_predictions),
+                                  pct_error = abs_error / permits_count)
+
+ggplot(poisson_predictions, aes(x = permits_count, y = poisson_predictions)) +
+  geom_point() +
+  labs(title = "Predicted vs. Actual Permits",
+       subtitle = "2022") +
+  geom_smooth(method = "lm", se = FALSE)
+
+
+

+
+
+Show the code +
poisson_mae <- paste0("MAE: ", round(mean(poisson_predictions$abs_error, na.rm = TRUE), 2))
+
+tm_shape(poisson_predictions) +
+        tm_polygons(col = "abs_error", border.alpha = 0, palette = 'viridis', style = "fisher", colorNA = "lightgrey") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+

+
+
+

We find that our OLS model has an MAE of only MAE: 2.44–not bad for such a simple model! Still, it struggles most in the areas where we most need it to succeed, so we will try to introduce better variables and apply a more complex model to improve our predictions.

+
+
+

4.3 Random Forest Regression

+
+
+Show the code +
rf <- randomForest(permits_count ~ ., 
+                   data = st_drop_geometry(permits_train),
+                   importance = TRUE, 
+                   na.action = na.omit)
+
+rf_predictions <- predict(rf, permits_test)
+rf_predictions <- cbind(permits_test, rf_predictions)
+rf_predictions <- rf_predictions %>%
+                  mutate(abs_error = abs(permits_count - rf_predictions),
+                         pct_error = abs_error / (permits_count + 0.0001))
+
+ggplot(rf_predictions, aes(x = permits_count, y = rf_predictions)) +
+  geom_point() +
+  labs(title = "Predicted vs. Actual Permits",
+       subtitle = "2022") +
+  geom_smooth(method = "lm", se = FALSE)
+
+
+

+
+
+Show the code +
rf_mae <- paste0("MAE: ", round(mean(rf_predictions$abs_error, na.rm = TRUE), 2))
+
+tm_shape(rf_predictions) +
+        tm_polygons(col = "abs_error", border.alpha = 0, palette = 'viridis', style = "fisher", colorNA = "lightgrey") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+

+
+
+
+
+
+

5 Discussion

+
+

5.1 Accuracy

+

Predominately, our model overpredicts, which is better than underpredicting, as it facilitates new development.

+
+
+Show the code +
nbins <- as.integer(sqrt(nrow(rf_predictions)))
+vline <- mean(rf_predictions$abs_error, na.rm = TRUE)
+
+ggplot(rf_predictions, aes(x = abs_error)) +
+  geom_histogram(bins = nbins) +
+  geom_vline(aes(xintercept = vline))
+
+
+

+
+
+
+
+

5.2 Generalizabiltiy

+
+
+Show the code +
tm_shape(acs19) +
+        tm_polygons(col = "Percent_Nonwhite", border.alpha = 0, palette = 'viridis', style = "fisher", colorNA = "lightgrey") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+

+
+
+Show the code +
rf_predictions <- rf_predictions %>%
+                      mutate(race_comp = case_when(
+                        percent_nonwhite >= .50 ~ "Majority Non-White",
+                        TRUE ~ "Majority White"
+                      ))
+
+ggplot(rf_predictions, aes(y = abs_error, color = race_comp)) +
+  geom_boxplot(fill = NA)
+
+
+

+
+
+

How does this generalize across neighborhoods?

+

How does this generalize across council districts?

+
+
+

5.3 Assessing Upzoning Needs

+

We can identify conflict between projected development and current zoning.

+

Look at zoning that is industrial or residential single family in areas that our model suggests are high development risk for 2023:

+
+
+Show the code +
filtered_zoning <- zoning %>%
+                     filter(str_detect(CODE, "RS") | str_detect(CODE, "I"))
+
+
+tm_shape(filtered_zoning) +
+        tm_polygons(col = "CODE", border.alpha = 0, colorNA = "lightgrey") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+

+
+
+

We can extract development predictions at the block level to these parcels and then visualize them by highest need.

+
+
+Show the code +
filtered_zoning <- st_join(
+  filtered_zoning,
+  rf_predictions %>% select(rf_predictions)
+)
+
+tm_shape(filtered_zoning) +
+        tm_polygons(col = "rf_predictions", border.alpha = 0, colorNA = "lightgrey", palette = "viridis", style = "fisher") +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+

+
+
+
+
+Show the code +
tmap_mode('view')
+
+filtered_zoning %>%
+  filter(rf_predictions > 10) %>%
+tm_shape() +
+        tm_polygons(col = "CODE", border.alpha = 0, colorNA = "lightgrey",
+                    popup.vars = c('rf_predictions', 'CODE')) +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+ +
+ +
+
+

Furthermore, we can identify properties with high potential for assemblage, which suggests the ability to accomodate high-density, multi-unit housing.

+
+
+Show the code +
nbs <- filtered_zoning %>% 
+  mutate(nb = st_contiguity(geometry))
+
+# Create edge list while handling cases with no neighbors
+edge_list <- tibble::tibble(id = 1:length(nbs$nb), nbs = nbs$nb) %>% 
+  tidyr::unnest(nbs) %>% 
+  filter(nbs != 0)
+
+# Create a graph with a node for each row in filtered_zoning
+g <- make_empty_graph(n = nrow(filtered_zoning))
+V(g)$name <- as.character(1:nrow(filtered_zoning))
+
+# Add edges if they exist
+if (nrow(edge_list) > 0) {
+  edges <- as.matrix(edge_list)
+  g <- add_edges(g, c(t(edges)))
+}
+
+# Calculate the number of contiguous neighbors, handling nodes without neighbors
+n_contiguous <- sapply(V(g)$name, function(node) {
+  if (node %in% edges) {
+    length(neighborhood(g, order = 1, nodes = as.numeric(node))[[1]])
+  } else {
+    1  # Nodes without neighbors count as 1 (themselves)
+  }
+})
+
+filtered_zoning <- filtered_zoning %>%
+                    mutate(n_contig = n_contiguous)
+
+filtered_zoning %>%
+  st_drop_geometry() %>%
+  select(rf_predictions,
+         n_contig,
+         OBJECTID,
+         CODE) %>%
+  filter(rf_predictions > 10,
+         n_contig > 2) %>%
+  kablerize(caption = "Poorly-Zoned Properties with High Development Risk")
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Poorly-Zoned Properties with High Development Risk
rf_predictionsn_contigOBJECTIDCODE
27512.167103604I2
85821.5859761524ICMX
90312.9468361595RSA5
91621.5859781615ICMX
152211.2443742559RSA5
163121.5859732736IRMX
164321.5859762756ICMX
167021.5859772803ICMX
167121.5859732804IRMX
186534.3893033128IRMX
209912.9468363492RSA5
225111.4208043744IRMX
268512.9468364533RSA5
316810.5881735506RSA5
343921.5859766067ICMX
355412.9468366289RSA5
361621.5859736405RSA5
383612.9468376869ICMX
384634.3893036901ICMX
386512.9468366943RSA5
415529.6927037646ICMX
417811.2138337704IRMX
459312.9468368805RSA5
469723.0356339094RSA5
475216.0827359244IRMX
480010.5881739371ICMX
491021.5859739662RSA5
521135.68783310411ICMX
521235.68783310412RSA5
521335.68783310413ICMX
534518.21887310760IRMX
550139.00537311150ICMX
550635.68783311161RSA5
561611.24437311449RSA3
599611.00043412339I2
620611.00043312808RSA5
625811.00043312932ICMX
631216.08273613058ICMX
641223.99983313314IRMX
670322.58767313979I3
670511.00043313981RSA3
672013.26530314019RSA5
6793.211.00043814223I2
679422.58767314224ICMX
680822.58767314257I2
685511.00043314373RSA5
686411.00043314401RSA5
686511.00043314402RSA5
696513.26530314649ICMX
701412.94683814748ICMX
719211.24437315167RSA2
744639.00537315720I2
748711.00043615820RSA5
763713.26530316181RSA5
788539.00537316719RSA5
810515.75340317170ICMX
813213.26530317217ICMX
822015.48503317410RSA5
854313.56143318033RSD3
907513.56143319078RSA3
936211.24437319593RSA3
960721.58597420075ICMX
980212.16710320444RSA5
985413.56143420536RSA2
1037713.26530321529ICMX
1065113.56143322004RSD1
1287729.69270425778RSA5
1292911.24437425868RSA1
1316413.56143326249RSD1
1349410.85073326759RSA5
1359911.24437326935RSA3
1381711.24437327284RSA2
1389211.24437327394RSD3
1411013.26530427776I2
1415816.06680327871IRMX
14719.221.585972828949I2
14719.312.946832828949I2
1472012.94683628950RSA5
+ + +
+
+Show the code +
tmap_mode('view')
+
+filtered_zoning %>%
+  select(rf_predictions,
+         n_contig,
+         OBJECTID,
+         CODE) %>%
+  filter(rf_predictions > 10,
+         n_contig > 2) %>%
+tm_shape() +
+        tm_polygons(col = "rf_predictions", border.alpha = 0, colorNA = "lightgrey", palette = "viridis", style = "fisher",
+                    popup.vars = c('rf_predictions', 'n_contig', 'CODE'), alpha = 0.5) +
+  tm_shape(broad_and_market) +
+  tm_lines(col = "lightgrey") +
+  tm_layout(frame = FALSE)
+
+
+ +
+ +
+
+
+
+
+

6 Conclusion

+
+
+

7 Appendices

+
+ +
+ + +
+ + + + \ No newline at end of file