From 9203003f1a9d0237259da83ae81569e395521c50 Mon Sep 17 00:00:00 2001 From: Francesco Minnocci Date: Fri, 4 Aug 2023 00:00:17 +0200 Subject: [PATCH] First working prototype --- README.md | 25 ++++++++++++ hc.jl | 103 ++++++++++++++++++++++++++++++++++++++++++++++++++ solutions.png | Bin 0 -> 23589 bytes 3 files changed, 128 insertions(+) create mode 100644 README.md create mode 100644 hc.jl create mode 100644 solutions.png diff --git a/README.md b/README.md new file mode 100644 index 0000000..62f5345 --- /dev/null +++ b/README.md @@ -0,0 +1,25 @@ +# Homotopy Continuation in Julia + +## Implemented + +- Total-degree Homotopy +- Roots of unity start system + +## TODO + +- Projective coordinates +- Parallelization +- Endgames(?) + +## Example system + +$$ +\begin{{align*}} +x^2 + y^2 - 4 &= 0 \\ +xy - 1 &= 0 \\ +\end{{align*}} +$$ + +Plot of our approximate solutions: + +![](solution.png) diff --git a/hc.jl b/hc.jl new file mode 100644 index 0000000..e569d73 --- /dev/null +++ b/hc.jl @@ -0,0 +1,103 @@ +using LinearAlgebra +using TypedPolynomials +using Plots + +# Define start system based on total degree +function start_system(F) + degrees = [maxdegree(p) for p in F] + G = [x_i^d - 1 for (d, x_i) in zip(degrees, variables(F))] + r = [[exp(2im*pi/d)^k for k=0:d-1] for d in degrees] + roots = collect(Iterators.product(r...)) + return (G, roots) +end + +# Define homotopy function +function homotopy(F, G) + γ = cis(2π * rand()) + function H(t) + return [(1 - t) * f + γ * t * g for (f, g) in zip(F, G)] + end + return H +end + +# Euler-Newton predictor-corrector +function en_step(H, x, t, step_size) + + # Predictor step + vars = variables(H(t)) + # Jacobian of H evaluated at (x,t) + JH = [jh(vars=>x) for jh in differentiate(H(t), vars)] + # ∂H/∂t is the same as γG-F=H(1)-H(0) for our choice of homotopy + Δx = JH \ -[gg(vars=>x) for gg in H(1)-H(0)] + xp = x .+ Δx * step_size + + # Corrector step + for _ in 1:5 + JH = [jh(vars=>xp) for jh in differentiate(H(t+step_size), vars)] + Δx = JH \ -[h(vars=>xp) for h in H(t+step_size)] + xp = xp .+ Δx + if LinearAlgebra.norm(Δx) < 1e-6 + break + end + end + + return xp +end + +# Adaptive step size +function adapt_step(H, x, t, step, m) + Δ = LinearAlgebra.norm([h(variables(H(t))=>x) for h in H(t)]) + if Δ > 0.1 + step = 0.5 * step + elseif Δ < 0.001 + m+=1 + if (m == 5) + step = 2 * step + m = 0 + end + end + + return (m, step) +end + +# Main homotopy continuation loop +function solve(F, maxsteps=10000) + (G, roots) = start_system(F) + H=homotopy(F,G) + solutions = [] + + for r in roots + t = 1.0 + step_size = 0.1 + x0 = r + m = 0 + + while t > 0 && maxsteps > 0 + x0 = en_step(H, x0, t, step_size) + (m, step_size) = adapt_step(H, x0, t, step_size, m) + t -= step_size + maxsteps -= 1 + end + push!(solutions, x0) + end + + return solutions +end + +function plot_real(solutions, F) + p=plot(xlim = (-3, 3), ylim = (-3, 3), aspect_ratio = :equal) + contour!(-3:0.1:3, -3:0.1:3, (x,y)->F[1](variables(F)=>[x,y]), levels=[0], color=:cyan) + contour!(-3:0.1:3, -3:0.1:3, (x,y)->F[2](variables(F)=>[x,y]), levels=[0], color=:green) + scatter!([real(sol[1]) for sol in solutions], [real(sol[2]) for sol in solutions], color = "red", label = "Solutions") + + png("solutions") +end + +# Input polynomial system +@polyvar x y +F = [x*y - 1, x^2 + y^2 - 4] + +sF = solve(F) + +# Plotting the system and the real solutions +plot_real(sF, F) diff --git a/solutions.png b/solutions.png new file mode 100644 index 0000000000000000000000000000000000000000..19d771313eaab4755aec63dd7911ff041bd1099a GIT binary patch literal 23589 zcmXtg1zc3`^YyZ{bazW5AR*n2q;!MQC5<#FN=busE8X3!g0yru(k0#XK74=wclF~6 z?7jD?nK^Uj3}I?2av0A@pFtoH3_oGzWJxCla7uU&p zUwl(jQ+IdwT6d(Yi_6=LP;ALi{t?nuR%d5tOH0f9n{y3K&9jRO6?u8n|NYQuJNSn0 ze0wy{YB+P*?~-DV8azxoY6^u+QRV67I9>j~?{x)T;5>97q7eoM2SZR$QFmsl)>;Tp zPfs=8ym_Ooja2EhD(UL#T19_nJ~=rlOd0JS-W}U87-m>*q#SwGVdix!P(zE2!o$NO z{qp213Ac8sE=_ps2C@U`J(g$|6&013d2bB3{Z~qno3eqK+1W@cv8A!GxSlmIIw>tJ ztw9Rb(d>m(py z7yOu$oxQTTNk~lm;cr)%CqYhXYHC(iR!B%lat$wH(!hqx4}0#$RtrZ|;Z+x-V6ttx z-udAa0u~kp28#}*-1wv(aO)&&XA27l#1#h@lf(yT6|ttV1^`&SghMaExn|M zgGRkK9FR>31U(vxT4y)iMVx}cBH5>9Y}{`7C)Lr|_{G`T zAO^o_?v;r8FFwyy4XViF09Bx+)8?CQ>btK(ED&2+^fGuuK!(w5ieV3LHygTMdA#zq#W zslWfjM?|zQ7-WR___gMJap~#lo-b13;>McY_b$;BGcx@4rplT_@?;|gomOCs(zNe$ zX`lbar=(me(kN&!c6c|GF2KUV;xJb$;IbZi^(l-bLyjrTxZQVWxvj;@`{v88?b06O zJ;FyD0)mIpQH_tE8bV(E^R4Oc?;pw#YH*l;#n0c^+IqJf@Yojc*d7}@+G5GU#Png< zIu;Cmb=AG#1`z?Fx6>HmxDt*MtiyK|dwaH!&LW?}@^VKjtIMCBVAE_5=IWBzjl55P zG`Vh-n{*f}D<8mu;OXR|zJ`U_jORyk;=LwjkyOAZB)p)dd*|0iSo#BPzpSuQ$jS6y zv=}fB-@8CmR8$)qo3OAjLv|wA%Fj=G#@=^V$3H(Iv$C*wU+nL!u8yum0#l;owM{54 zW(h1TWCFJJ?S+AHqZ7@G7Y?0YF#e?SqM@U&AAfMqVxkiBC#R$|*3;AD(>FF=nys-C z^4M<#6SlC}I9}`7+uMU@Vq;rq^SO*?RQ(zmiHuHcYGyWil@*U(4W)$l6;EXkmD5e* zvQWGa@%HvEP|5IqG+;MuMRf>1ZfIyIq#iS7Q zbvimZaiNU9{Q3JkX&H`U?qE-}-?E?DS9dQjSy|aUBpyD#8RWwU&?3Ua z^=mAZ!S?P>miYPkB_<|fQ3`+9=udRto2q)ymW~Kl!!0T*N^HBGd?WLF)7CbU!WEn_ z;7?wwsEGRAoB>w>ZUQEhl#qZ=Osv}?9)d{$ZZVSJY9AjLCmT&^bnD&|MX@K#{qp5Y z9v+13!a{3v^Ms@%D+dP=_gzg_*HcOn@2lI}Xi5?NO)(dLVWAcFk$87>TKIp9Eth*0 z9ng{S@!Ruh6F*;HVCvP?)pK=rz}-H1+?M2)mX`8z_UN9HMbuv#F7KDy+g!JXii?YZ z8w8$gX>Klcqk+a-GdQmBk5!gYtA~TP#YP$jYG8!$5-wkZD)G0@Rer&z*SEyP41VXg zORZkuYWJS(DMyHnnmlm$?jhPqiERZLa4SZ6lVS{!ovRa>a3|9y0r!Uik3!CWyTHdwZuCu47j*d=;^IA7pJF@eG-d8^Q zbkC(*3sLAN(hocMIJs_4S`ORG3!Svo3K9MmmnbVKMZI4H`##+rwYj~WoRN_c7Z;b6 z_1@If)WSl`XDYV7{^I&{OG86rpueA-oZPm-d~RmuIX(Sk`%6yFk)9q#5s}5o$vR*? z7yD=+(p_C$jg5^ZahO17YOJEe!=FJYC^FE9*m;Qm{P|PnsDlv5s#9KARK(1}Vrg!^ z(BfH9T-^Wjr-YPL-j3MI6cC@`GvQN95HyL%ww zIQGcG!XmPb7{j2BXy=<1ySVVvGf`P2W9_C&XVdsv32NfvaEn^JUn?k}?xCqyzyity-6?&CO?&z@d&64y1nigD;I1xSaMMqZ+ytR~|&+k5#PHsxX#K@=(>|}5-f|b?Qu-IdI&sy8{QjOIxaBG@! za#p|$z^Qwy^aoOg(;gQSbF>yMhJ=9Nd;FUSc*N2*?2~Y`Bm7mE`V*Kz z#9pek(YW7YR7$S2n`Y5^7%6ba7w>| zQGv?_P)S^0U;h&yyHT^-Y_$c5mIiI!b#--5q4jriD%>2;O+NsI{NySJ(;y<4nVA{2 zcvz%{#^?yb5d}$rEtti|!&}?lwzjrL5m!|mo}CpC6r5XFXmt!aIXNNYb96ggYP&E4 z{xB~$H(xoeo1t9yDUWi0fkMsAJw7v&Y(^}(Tk=9uGMH9KXdVX3{+E-JW9H%E@d*(! z@#1t7S{(v%C=mQj3=J22ULl%WS*fOQ$w)~_DJbkTEc@3~R$^PVs%n5~bnN{2_}#)^ z|E4y(HM_EvvqQK8{O4ltBxwVZsU?bsf-uQDX_085C`_*e<*S+TAJM z-F^CZo3xN*yJ(lBgpAb*U(|n>ke(hjiGxE$!MU>TV5%%FH<$BqR&j3l`knl{P2fH- zFfed1mh$AHZ?FJJ2iCMyR|z2LKoZ9urM8X^D>E~%E<}8h3+xcUk-xsOdKIEPuI6!LSMBkjf0mYPH{} zd2;xbO3HdWA2?;Y3lewkk-26QJhm&;xRKqIiUK6tXdL$ws;I*A71gF{KV6M&T zB>@2e2L}fh7FHD^9j%tOkq6yg_d)vnvA}Ud){@H5>sN31Twb}=kJY)1cy2o^w$zu{ z3LV}nDk_R@Y6K8;u-xl;8x2sldAc3{{#_(fG}SycBs;_|%U&*9Zk%u2T->~~vqOW& zm^m;!UF^42V}+F|LN?p*UP)2$lj^Z_ru*U2{hCfJfEWN@I4-sBY;JywiHXU~{I|D< zwgRv(AOVbxjfE*I0K@=iU+n%&1E3}lABkS*KkvWV2dO769v-*-Oy$VPh)RZF)sUn# z`>o8y#@B(QF%F}l#`NU0?YHZHE%&FUz8_q@TFZz`=S`5VLoN9h_c5-nysoa?{hpEa zIjkjp{>IsP)`D-d-0mkS`7t(D{XmL}Dl3wVuiae%cp&ApR|?alt8(Q< z^7uJc>}BjCm)i?{3GuO5bP@d#b7^&lrc*gm1DFoVb5s8&?F4BoH_@)ZRF3Hp?_EQAW*enkw+D#uk8eJ5Wo~5D z4^*ByuKw8=30hvy07>c6dFSbY4pftCA1naI zNLaPqG6lE`8yi7VM8I3e-P}C?t5c+{^_$u3r(ei9s5xuVUn@3skol3nx=?OpAHoE= zXvRV=a9)uf(BVhRT9o5VNQ-_o^nADY`z!)}9|?SF@=vZo-q0l8Go1T$z&b@!34QZ= zd?59_?0k4IW?{h~Y`i|({CBnJgLj!~m!^7=EWF^{82+}?Vo7E*fbTg%ykGD--6$zd ztfRmWfq6%VdF&PmIT@?Cp6^DVJd0QrW=SlmzE@O;c`Wclu3}!uB-eT09$qGzv0o%9YAFnjf6>_SQ3z3VaM+N^u|u zuE27-TdZ&qSorsC@e*U+({IcjxV(-G|8^nYXfb{j13wfmF|eE#2v7-W#o%F&EE+SZ zR)Y^>H=Co~{$%s1wH6t#$(tdbfp=IhEOQ*HW;~te_-ybN86Vsp7`Hky8kXkte^=nj zNlSjmCZLHW^@S7Tc{$;pvDaF2%Fw7Pf6cP7=V4;UZ0C2CQ7wsj1V=0CVutu5R)y}p zQuDFSVaEmiL5-nJ6Ro9;jeC(b9WL$XLX+J9Catt4 zt_mr&MMUo?XD6xn(GqwbB7vhblj2X<`UkQUUb==bY`OCQTF8~UHUaq5wCCPhhy_NTr=3F#!6^R^zV z{l28#*5bG)czd=_q{R^&AB}AzB_sd5H?ov0lS~-vITgRyzdH^i|JV{bNwG}P6}WWV zZZ%={_bzFDs2X@+j1*J+17GS}pnabERoNZ-0M*O+whWW)VTJwVF_M~wg z5Gr*(Vih(e0eVKjDs<>qlczmawO7;1?t`5Xm2O`ICrK6G)VX3?s;XC0iF z)Dk|c4EdWgI7;ZHG&g_r0NM_zi|JJmNr9!3=&Ka=LD?s)i_>sdeV5td{qO5C2{(;p zT+egxq(oY5$UgQ4epZT~$Gfq}liibgvN?tg`Y;z3NnOC3y>;K)c`+QBnVFwS)NNEb zQg2Ej*Q^Xz52bChFZ@7Ke(v$!>~*loO#fW1Qaqz&=bQxyN!kQCeV&LAq>Dr_G zBcj_A%M2grMs3XIx>K@~^-OTXetW@DqXiuQE-przzaaII-ji7U1g1E*J)2J9I#roM zv&Ua83wK3+WnehY>_U*_VIsI!qi59ME^@nF+L%u-T?2_4L|#YbztI;bzYEct|Pq2pD+m)z*U;95VC*eRh00pJR`$`aSy*ksAAw z0>v3cZns=nbG_}ql%~HKDcCf5w9kJvu=-CYywgyL6y2~}9T1E#tVC{=_V`_Q$vMpx zO-kjyyIQovm_LnT+a)S?9SBMQPP^)F-mq}Izw+vPLga7(0RapGDd1L;7h8#9TUluV zj{5?9@A#(l)_ciK1z<$Ct+Ic<4n>Ou?Cq2dNBa*Eg$DSahyvIXzfuQy!Q0#0lH%g_ zrkOu!e0E6!9{b5xQi@P8aU}p-Onm$fey>D$czcxA^d=1l62d46HS9{P<)WenS^7v_ ztC>m>I!G4B^kWZLha5eenz!#3x=+`J(ri+zhG?Vxp@dKlJ*Bd4 z4I!=*c?rFrjh$=OWJL&q{&$Vy)Y|~^7?>O^G}SaWHzy@cC$j2vkt<08D_)&(aUCeO zD;bIsuq!qL2QYmXNf$y!@~1qv(8YQ^XK*XUMAo;D>EoT?sYz8;6$ar`8ljw^?~(fCIlgTsFDnWMu3a9Apy}1vLJ{y;Mqavi)-V;_$F4z5qQ*#+p*< z?n0hxe6(^D%Gd1skAqDRk)RFO@K>lclZ`ZjtlR9A%xgWnO)GAaA2C9JkS536GOBQ# zuP^hvy`ZxH3CJzLVFI*mRuB)>4F0ZA4*U)_rIYCu+~wxK(i&{Y6`~U&1GMyFe79r= z?AU(okNEoVS4OWTaP0xc1Uyj1#Y5 z*FEb;n{%+|t2k{vcIZa)T-UFq+Kv~2-w%|BZO{=C8g1ICmR}zpN0U>`b~i>k$$g8- za?J~MqoM%?h~a$o6A8BorfYl-JL9MsvA5QHIkdJCWBB^8Vc7G)121%tIoRElliR`*V!Il&S!PGZgufWP#V8O&Uz>FfebrC@vZl+(byK<^iRv?-b|5ahQ$!|e4nnO@g zB6ZNYE62l9*(a-44a$QFqIn6iR}v;8uUD&7SI$d15FCJ8Rbfi@3aR77bvFU9n^EjDj zeD(~gJ%5&*omIMXjh$=SX9fFcOz3$;JZ{ZBr{N7 z4D%Y5w*^)WIm=(nZoiBMaZ5wNNTTc5yh)-5ZmR~OOFzZP9wFqz<>@Zu6a_ta+{+0@ zru*7+fs!T9kMIdIMs8*8pGDEA2t^k>?8v}J?~w1+ z>8P}itKaVDRI&P$*%wX5xf70A8=5j!B8Hl0Nb5;>4abPos0=GIE16|{j(G$#6ty2p z5db+HU50Emrle>rO(hyJvU!!t*j)DcB}E67=Lhn|UrF+qERrX1R&Rzu?03>S)B*!$ z*n#@OAt=Cxe!Xw)P}kNp_V*G1gJ_kMkzd-<+o)Ng`w0%@xYnda#%G>MVR4eHFaNEv zJ8X*JCF@_;zhHR~0EKkOcY9_ZU=Uz}wv5H0_?(hx2F(zGutsIBh`U@(z*?@rV9O>0WW-54e47 zW3pB6K1I)$K*^!?eZkMq59%b;)Pe5U3dcQ(V*0jqb9q;@8822%fj@z?Eh7KzKiI()Re4A~-Z0trD3SsmoRG0VjFfS;gLCAEa~W+}OFj zf`Mm!{r#z+kaO@SU#C>p%);Wg;AlwOhvQuJfdI1OgRWP2X{@{)%*;6I-?41%EfuJZ zQs$up_yi$Eh2$?@)RmZd;4RY|eey?6p8<3#0=7(*g)n?nzT^tUclw4>wbuKDWGGIj|Mxjg`Hjwz=!1 zF-qxyP7)B@r}7FReBkcqTyVe2zfv~C`juasadLbZSDwWk>-6-5($alUG+bT%DHvKf z2%v>T1LB}TOlVK&?yfcq3w+QEaVl97|Gi^_3Zpit5tKKfa0vyzHah06=NP%CMzOlP zIm!a8l&a?n)rxE2WQMCzq@f?1??dSIr_npS*}Vi=?_DT#2pQ%>?5R?YMJ4uPETiN|k zOEvwAJJjim`FY+1o}1Z=cND>#rdWnhZ@6n%1w^j3reVyHch3DRw@{|!c7*vtdf};J zP8nUG5#8J@xF`KryRQPh3HzjExuWzFwxvEq9Y$cvsm69OyXBrr^bP-_a(gv|h+eu7 zI0ROM_l^KqIxgpm0NbE9T0tMdJ%TufdfNK&-WO$RBpK_oEZLv}eijm##LJN+!;!Y^ zB?>+<^f63Ef7!dxp}*f6lkGOy2Cubn2B)y!AXdJO%=}U_5;PW3j&lvBb46w2K5{ zh#Up^6iDl=T%4E7z2tD0rSpJ44df1PmFPuvcPA$@q7k|U)Nw!Po#2u&Y7y(98Na!9 zz-=2J4g@<+eVgZe#tDdsDoi@s|Cu#7F1EI|lClAozXtFfuanqQQd0cRCv`xr@$Vld z9P7xK6}R0T7Jn^{5@WCgVR;C5sHIS7P(={e=&*IM+3x$j4v^uDs?y(r;+K@9wB;Pb6y&%_wuwRg*+-5E1B_49u_v1FMK1^4o(b5l@EMh z#Pd9T)c*@$Z$b;={6In3b#u^bx&$1zfcx|5zIevA)>hdlavKW^0uqwhlB5MchJ0R0 zugl79GkfiQlpdU?5$40&BR@d%A#mdLOj|BoZq?&8LcW)kOpLKafGMP*p@|^kXa@{g zS_LTTH8eNRjgR9&k_HaOe!O;btDWJlmXl#RoLbfuTXJX(&waK>E&c2a2^WqPMxgsF zA;Lbl&w;Tb+URk0qt^!@6%0yYT_YnAF|p<0;rENpU&zvHYHH5U&UWjONQjBets>5( zqJndk^IXcpP3qPzq9T40=^@j+ks{DVP(tvB*auMs4#Odu4RiJSkMJ>eG?_RpYY%JB zyc@^{T=zsk;g28gD=;|73;o)eCReM2cnTqpbbjZ*3Cxd8zeB>q6#+B4*xBE&2&Oi= zmB_Cd>~bJ}0!KtgL}o-bB@vqN64?VX9jFIWT+xA<6t1x{Gez2QJ4)?;0SPmjyXh>C z&&B@D?Ug%VXAHX+7uzZ;j{p~R&XHzeYYRjuCV-)v4uf8RPB=vd_h zfB6C`NE7ud0O2A+gFxX66rCa?|2*8^_IxEd1uh>H8iBu;E*+(ZC&;;LP7An{`1uLM zNU&sM4xA%cGwA4u>48m>f@>p^fUD`EyHoj)?V8Qm9WPg+T?T4-A#H7bfGu@&{0BDI zF3Y6F<3-l+aEpDeL~lr-QDsCf?$$D`KR(E$I#-b_Asde)!w4yL1LNm z+_WkX|D~s#f^uN)rWj6`)`&{Tabh-1yqT#TWseWv?|&)F%fo|nz`{f0mT<~6#T?4E z_iP^-7vcxQI)cN5Xho7ghcQ$0ywt=$MkxXt8Vz_2a0e*2SQ;`^qV8|te1BR#jO!!0FU88}FEz(b=;dLN zY{HAh9qlX*fy!-9D=@5&GebmSq${YEWFt^^_XBUW*CjvS7Pz~&3A&la#q@H53tYK2 z%vDP=frXrSki+1w;N<^JvM%=9w#B_UpcprSb`%$y^9vaaY-r_#$*be^Cy`40Ntxgl zDJwz=`@ZJV^muoBxQ2=+{54XDgf!KTB|U|wSx-8cendg}z+iV4TWCR9a&_2jFH_NsqQqCx-`nEX-STn`Q!K)IOH0~4$abenYV>Hi?^!}(-iSbl-7 z&^V#VLp|%kr>q2bwiDM2U53PJ0()BD98aTUpnTfYu!Xdi ziuQVa-8>s;7F^b(Pd>3c(m%*$PzKD18iNzj-Z;tnfQS&Hs!J1Uo(+Z%IQy&?zac?H znylW*8D3FFW4ATC4)w@M7~RUxy+_(Wd}=&vOkL==z~6h14v55Uc80Xaz`*2vRazrr z(f06zeQCLsYd{!wv?*uN{LM{=Rf$NtuUUaK79=Om67!54`X-2H1##q=HDWUSX5j4J zBX16`d>9C29O10tpqRi``$89)W_?_MM38=+;lqy4%Rv}?DWrrilLL_wK01hp8Y2{W zEC_vtIU8$m^SGVRNAn~Wj6epApeB}^f?4W&YDDa7UlrkbTHhqgK)RKy?o;0tX0@9W z`cF#hYa?MiC$um5M*Bv>+@DFNcFpz*yyS}o>PmHpQqU$8ZjX9BKzPJRaNEWdb+wA$enIy zt+JerT|GUC0^}9MWtfOu23#U1BlVsx{W=i|NK(Cm7)+oGWK_)rdW`F>MznMT8&O=WMx7DtRh&n*k^Ybka(15#RIZ*Kh;*TZY zL!j$%aBvVH9|GYT9ty?9#VsveYKTr4y-@d4MbEJGjz^ql6O+EYaj=AiD^%Tn&zZ{4{*w{36bciIxfu53ysK0n1^*FUkMIsct@OXm7 zSi2Hz5U9Xqh2b%N2@dW8)w$olo%z#qq(Y#&&2Baz-iH&gW(s>%15Fwz6*@O90E__y zF_+t;pWJPc3QU67UoLI||@^8`q91taBB} z0t)_n!LPV9)yY^pI)V|R0_x<2m{`E-ug~O|J(H8?fJOmz@vl-f;xbo?^76evX#~VC z4~HF(H&<5x4S+&e6(&88)*Xw5j0Ny}I8Y134Ml(YJ^Zns%TKSzKU-kEZ97IrNA&;= z^Fd01R?rU$VPR#RE;DGboh$}q-s4dx7Qh=oVfT(gF&V(|UUYiIAG@ixAg(Lr&>{qa z(jO2kER2uWT$nw5>fO6pS=@3({~y&*U?6(g-?PuEH-7~FbA25g@+uN44(3u{|0?Nr z_gzqCrm?T@T??WA$+`kKNvWx%(LDwYjsTM3FC5hGT4=djoCC0|i)<4Nggu@VD2!E1 zH6NH<@X@k}6sN6P{#T=)V(y(+$AJbmH+usq(qiY6noOkcy@j+8fMByjS|$_|Uo=PiVr#~(oLl9H10E&7b6;xolWj-4wgE!q`S6tZJ0LAz>AMFSxw zCce740pd?eac}P%kV49SCH^Uxz&VPuC2ctpYa+c3GQ30rR_(LB>BFPdPEdj5wj7M9 zIz`zDb*1HIG0_ILwwg=p`x}TY=;=3rO@Na$Jsw|I5_{Yz3^Y~MK~pGUBkDL0Vn94J z=i}K&NgzYe%t8>0;l_@QBRD6=L{ytl7vBGrgh``38$hemss$F$^<>K1#k-qH8^E=N zVy`D{>=kcU126ZCAjfb-YxIXfBhGA$$+wCE7CWC~R09xlBviGc!NcsCm{w&q@2Y)`wB~pegTtO~T(rB9xbdIC;3p6PO%SaixIX8zWyDu2lK-rCXL!uO zy4Dn6QG9?n8m*f#4zPIy2dEtuk|AS&l7+C9Q8Yo7BJO|lmuec6IJHmA*3vlUS>HCvHA)9 zPgk52ZmY8%XoHdH0mR9nPgbka|N1S!5DW@~%5GczHI-)i`3rz_YhfEM`wyu zi=YpE0pS+^^&Wlg?K*B2-tg0nz&q6nISfEX-X}8RWaVc}tFXn3GNl!sBjR3!YKUJp z!N@aCwEk zJ3%{1)bDsgf4X^JhGdWJDFQtUftn&E!%;lm`trtLeuowT4v{DZ`W@07h+&F01HA|o z8BsiYw3=|F_(A*p4T~sad!2+GVifcV)vMuz6Jk|DOCN5(LwXS?0itMu>8()4VCX{2 zf~ErLMo261E-bYGeo9k)|IA%pvWroaF%yKxTydFR$kzF8o@yH-Yv3)6f*N%O`8k|T zAmR!_Ba}T5ZKT@fwBBy?-yvB5p98vNa#&yJ_254D6!Z5%VB(<7Ar`2CcXw|N8kT{4 zp-Sa>7|`H-?%%5~(Lh{4#0~sl${PbofMX9NTfzJQxe(8Lh1*2hlrQ$0gN@^R%M7*+ zeg=+6Nhce)p^~^bBoGykq12?~0mPQu%R^8}!&S#oWdAWe39K{4cjXl}vNR+j&>lu| z1ZfTQ2s{Wu(p+b$=e6;nKUJhk6U~iO39Z~n*u+mD;NKRi#H3Yz%I7ZjXVcTkF)$*5 zpo*A`tRO#MUw`a}7iliu?vF6h6|_O9WY8@P&lJ@Z?g_bii&H%z*AEv_ID*OUp+ljB zcGi-?G=d_{ZF)ZJArI!fzd!&4l?Tx60kO8T6A}hRBSXHt;13x1=7w+(gky5~f0 z1d#R{93SL$_W0i8gbBxZjztJP(GCwAtzRy8u28bOX@+FKu-}cQ&gJ9fr|rExpyc5m zq^l}Ku?WE#q1wiR6e7J1dJS`2K@)G$<2Hd2RZhFL)?F$d{q+CBkbnR2qoCXVRgeT~ zEP4_*JI8V2@vsxw+N<)8st3HKCxP~34$!!pn?HvDnF6RaG;aYV`nG9GLoRKCvM`M7 z2+0x(UD5Ztp0xZkW!&((b(6sI(LJ1P;Vv&y#<(e0IZv zq1M;)H!Xag@Aq4T3qF(^G(-Unu1H?_gDsg_1mO$d?sqc*4C3SFw)0k*1txpTFTR&b z$`NF`b+r0&cGi8%^ojugtDwhz#LO~ZwgCF9v^NGLGJ^|je!FAh$`SNCpY((V8ZiU zzy+ymYMOU#@qHLW;p`h|f$M_y1onyJa1s+k-R##D!h=IY#i{wvk0x4&PIZxE74aC9 zw&xoT!L5M8_3PKKpxmhj602WhV_XtBW_625_UY-3fs8AtY>KAtwqks$1>> zQgG4L!j)yKrM93pm2pC=G7z4o@jWKIC_nl1(A9?q2sH<*Ew!ablXf1eD|M5AlXQV| zq8j8)IA@qeHo+0EDdrtkBjdl;;8xsY>f8q9_fLuXQ>v4OMG?L|+Bpm}ZO`u>u}sf= zp%k#x4umrgWoQOa*FXc9!V3Ke#ufSuBq3CQC)NjuAdu& z$iI}7jkL`})g`KR!n4@JH;@~sUr+(eBb(0@;){rYa%-#}JkcOzcIni!>V_o=u9Iad zK$^q4yFBgySTWKs?DgA;k03fD@9IpxJ^W96b9GTV*_1J(Pa&w^=YOp^%gu2K4)wzQ)Fzk z1o$Wz549B2Dhk}i^LtrZJ@m?nMK!w*aPoRlS(h(R0T>rpFVM2or;Gz_6+1zQ^M* zH!wSDcn6U#jeq4$_H3q$a z3lbL+zp{Fdz2kQG(-|UC<_#)V>u3ZtUK_V6%N1?y^^Li@_S`Ok6`%^pw&Ig5YmbK) z_8c3<|KW1$I%8}6w@g?vD7@*Ex|Jf$VCj0>6eVw+BRN7|!kDrZ2VLQs3$k8_D zo*ImH0hE8}LG#WN%gLa01^kPim6@f)G-BG{S3ksF6P>nt!|@n3^nJK}2o<%? z-8kOVt1idNiK>0a*$i;Nqd{gu0#uusB0e{lhr@y>?_7n3BIU!d7e8ArNx@x_5hK+i zzv&VaG1~%L{3KS|Z%GE%6@Y?604Prn9`JiUG+mcM2jSY`n<2-BHya5=9YzmjEf$_n z#b>Z-KS1f}3dsHPf3$7XZAk6Au^7U?L9VyZg%qL-o^>jqXLd-rinoAgu%dUYK&YuP=<1o|QfF zMy@e;I(&A^YB|*lQpAA~`U0YXeZ7JU_x5c-m@>nV`gTh?a%WRY4<;wTBqgwbaU zlweDDvHvh>PVl#p>cA-d)LWFBn>#qD1k~FiB5iyoe%B`!_mgL!{SJ`&VJ|bTN$3>nJn+gk&D{U}P2 zxv?=9(0jzoi&0^NgNOGHv{RhVfig(SIzBwhFCCDa{$6$Jj-v3nzj4P5kqU{P+n!T8 zPrum+1nqMXkEl96!>z$qa17IInd1qnYO z?1zTRkLo}Y94_|#1E`vg)mW)Yh;ZG_Zms><<2f|`*(T+QRTX(GWQzF}VqE=4B6t@_ zbwEj0`&{q0VCWqwr>Xcsz6!`0g@jrdRhQbY896vAawWUPLf_!TdA@n0rhnM+Fo0)5 z_>G9`>rc2>&<)r$x)gR!G0ub%zu;;QyO_SULFJHSOmgxEkbGu}ECWps=vwmh^aQ;h zKrg$#-n&mxx3O(Xa)K&Yzt6wmeAN}e@C6IYidYI21-UiI0mlDDAT)-p15|d@xc}IZ z?oo5S1>)KB3eEB#KPL0!upn}vjmfQWm5>4;zmA8Sl7D6(-3L5A5xb%FyLaQJ3PeOn z02hKTyshCZWCR4IEKxsj2ijYPC$$w-<@a0Mib!J+LXZWKxxVJWkwXQcj)AP5PeNpO~c&C`_7fhpG`l$9um)xUf~cU%a0>AZOjj(Ced9)6CZFj z@4S!?1^IfMEU~qi)7QP^pJty*uC(mqmEd^-#b&=yQ#k;cC7=_~nqPqE8Pfn#?Z}$! zJ?K04HIf42!&h67y}%!CXeKt5-2PidAq5whi?6p{WA{DDQ9?Rb?I{}1f32Ybv{0@5 z7IhGDJHZ~6D2alP+hy}74)Rm-4-_K@g(wPnHcU3;~{2xW;-{KJ(o^@lKNCm)+`D-?G+xWlMHhccdH>xL6_miJUd+k)E=sR zf!GiAAKp_<8>l-jLA1>(tWaKl-_~)^(sEeSZ-xG@7P%R2fP{hiy#`$&iL3S3hL15` z5D3|?r{4g%3SD+I?l*SQVX?9@)eK4Y zhJaOPSassxyDgi3GfGMz8qQ(r@$vRCkWL)8WJf`SwW9oPYe@`80*6%bjmI$&Z39>) z4Q3~N53D64eJ#QJE;ntG+8_a#q4s$Ebp7&DJFVsKFc_Y9u{V%aoq!Z844a9W>42(T z>=mmA&^VTIDZ$4DF=>zEf%=vkSZJ%Bbcg?>F_NS~;V~#7U*cq?)w*^j#$AU7&2Yt0 zuYxY{tMO8aZx)LVRgf}jr|K-ha(08|l*}tA_Ou*K0V1Jx(QDrLmQ)0|9T@$wD%)lK z!U7k!i{<*ZJfjrjlpZK(JgM}+pdO=c`Iz~bx7F0#mFtqNbPz1zI>NoQn;LfKia)EP z4M&TpzL5V8cu^cKclz_Spq_z}0pt)2_5eNh&DkC}XYt%*tHgI5vX9LdH*h<3B zPxAz+c-Rp$@GU^}WTM=0S-u2Rni1}{E3Zvc_zMzrpg}A18q5Dw=csr^`<2GOx|2(a zQhCyqpOyPb%)U2Lersu_@S)z964v4 z{!gYV{ct>hfAJI;Yuv*OXhp7@6Gu)A*W*K zsAMNT!^;7ux5`f%H5?jR2^o1}>vszb9&bvUzBP^O7XS%cs5PNR_6sLz>AZpFiQEat z`DUX4m;3JTWT2CvGD<2G1&B9aY2s^r^Dt?CHxk9%lMjT)O{L9B_%&bQ>m%M8rZ8Zm z05ECYP_)R|V}=ui!$8VViibn1?{Qoq;_#PuBv8L^*R9_q6Fz4pw;t4)K3t<;SpA*{Ni#O>g^rq?hsTC9iGp31RY5}_g zatHgC9ET)rBX&fST*3q|a1&)!)yogc5;?1)D*!nHDyP_3zN+earyvV*^r010VEiu^J>OsaK+ootq-A!~a&JJNsRpYVgGK|yCQ8hT zIQ^AIL2N5IYTa=%p9wz4Kcw(byTL&7L?(snH+{8o!$v1F(9Q=GEw0+yQ$RBU8c;yX zPo0V9_Q)zwq|DCF{>5MdNtAT8+si{2*ONa$k9v1;^M{`|5A?@{!Mg-Kz$9{l<6k7+ z^#kN3+~8GwJn>JdgS*=uJ*3l&A{<3SU0uDO?5VZ!wz#I|6sR=?-FFlI*F=Bmn*pE= zDBH|cXlZJKP9BRo6q|9TSpIOs-@w;X*z;Z;&IinwRi~|pn)z>C?<($Q39n$UAe(-} zrc)FM5_eh|3Tm?2)Zf*mr0|}Y|CQ!WJNA5PYU-Q{bdnjh`;x#w>Y-hxKWSy zJH{d0q3i4HyvJL%<(WyfPgvS}pRhdR2x^mwhm78A@BWB(}Y>8n0brUwYo3DW83? zEG1<&W130r*4Y{oG?3gSQ8QIX_+r_spt5FXD?B&up9>j;D7PR1S>dOGB7)3*A_Yke zD@GIT)nvA-sb0>nPs$R$e7I^b60m~3q5M8Iud-oc?2|07$Ds5+VFlqy%~qOFdBaQX zyLx9`wx{fg1?2uIxP@W}zHRA)% z*%f}?`FVD4hJ~gVu?MmDU3005xRZs%Zdh;CZx0Oj?j%Rr-$We{o!C#IuZ}uy7!9|p zpF}8LQn*{4DExa#Y%%Zq{#mo}R;lQD(5JDod+m)jc=G#8M?7>icT;!!8u{jDQhWAx zHYXeBN$T3_=L?}Dgrt3{yr0g&1RO{Z+u)LfDSny9{zi##YO>BAFP8R|^uqDaxD{1C zc|Pfb>k1?3e+~zD#00O*(mxz!Qwy)1lU)eZ3-h`oq|4)bmB&ujwb?#t8uv@yNRF6^fY+Z53pj&mnW<2e`}(^;Juyx4#b?ooWFY3x z8mI}5n|J|8&Y;IiXHdoiQNpD1B!cymCF`$}cH9Wp! zyAdPNd6a#8e9v1oG`H)aNZvO%rPi&{ZV6G^_^{-#Z?ne}CBM*hx&3!O-?UqNAJaFuHhB zH`LK#7Q6I@B#a|7r;=MP9+NBi+iq@dcKF3)t^K?Dj2|sup8brDS=m{!axhQy#CqAYi>4j+@A?(QRt<*|D6VMPIf1Od$w<8CG( zo}}u1%`vUh@vd95GMFC!G*deh_9cwzL1dYhR@Phfw=X-4pVj4zYa0n%MJJLdS#(>7 zMS*8s*;2W1@nKYToq#hTSG1-A*DtEjz|T#eZax*XeedvcwBv`xI6O7*;t!z1-;>xu zNi!d0+{;f+P9`#J_iSiknqPY!z6Q(<2ylblgsbz3)kk^_?TrLmMHyfto`=u)Jh~Iwi;`5L3`IYLpvaD>w4}E)9L}!ue;pf5NE~%Pu zyCfzy%IoAP>-{{O5|JU38<`FlY5?8CIa1PC1um*kBOoCZP&etH;i$c%(J4-ZFp#h< zGc9W=Dv~adr1M}hWirV#YYO5JAI^_-oNMZo>f%e}h(3o}z{lJ7L%a-APzX6e80Fsc zE9%@aiH_umFafXZR2C<$Ib_A&oe-H~r{P>1oW=dtGC`1AJ{W@3#0U%lf{{hzS!a-- zphYB-bGIGh@?v0L{d4f{0`~@I?Pq9-a4?48I8o+hi#Yy&f5Se=v~Jc%6BMov4yT`8 zRHRba;QZ)JjKbjD#TU!V%iG)Awd$vrRI{Mb!z}>zkENbNS5Uq#Rl!U5hWb}WZ}0HA zgoK3PpdhyEHtCz>uK z8F-)J;PiSPw%d!}Zbz{&CQtfYhU=erLtO;!+R@AS?Sv&C3k!?)a!}5!10I+lHtF@z z+|Lp~yriddKrUMeuP80;Kc}g!Esw)VZgyeu8ONg2;X4MQUsJ`S&l;u>r1v_lPK>*T ztF?YV-Rgqr?O5p9*&9KjqdS?Im)F+Rl(C9F%M91BCqcm%L!b*F5OBb0SKDM*g=rUp z5LQK5ne~R(A;@69nN{TG{)R3=`d~f_eqD(W&x6=fePfri=IhI#glXG=%{2`n?MCo{<@btO?Kh{Y@O&!

r;l2|TMN2Ube%CLq}`usjUUlb0Ld1dG_;zLWri9~XCJcGd2 zyotu^W5_fwz#_;cURebGQw%ZN##|?lav_%^@ez98MP}3bFCu8#o3Ve4Y{xa=O0^1n3mp9|Yacn<+WO2*ve(-;iY5Wj=(q z$ejFp1X;|ZOw9kiWtFtK$cYHS5yaC;PsUPgz;n&1(RSYdLap1?xYAoKJRtcwS!K z&W_X2MK%?1qseC(q|8ivY_#(foOg&%OJk>}uL1=eaGsn}b#-)(4wi~xWxsL%SCaF} zkA-8|NcnQZ)+ktG6w)f=x^(HTo+~jrtmKx3#UB7Q!iu^CtUCG=Q|t)m-pdDOQ~%=< zRH^sL@yyI|!89X?@Gc1n4R0oO95~nY17`&$G~}90A3~y{L~4cZqh{DXo1p4nMT3R&@3tv}+e!%1sGzp;85{04p)YKFJ zhCJoi_+^Js6!Au1}$gE2QZuX*AH4e61~)8Hc`BLj2H zTUOxZEuJrNrjc_|1vs6l=>llt^VN8vjHT6zA_}68N9TtsY4M=u8^N@>(n=p7n^9bB zO(9fTRW%%`$Q3Ji7g`S%e0WT~Mg1?&sH1i3oC-@ymTJ}?K70tQn}KueH@rNCi*)w) z_d)De7989sMRG~ydn=%Uwh>a|2`z4}QbE@bzWGpkr*3 z=bGSEa-TLWMrChT2m6f2hBO$;Wy2L{?g zs|#jP;;Oa|4wBx-l4$TsSU|Z%=T}u#)oQ)adopwLrC+~{$DwMF1@x!{#_Lh$rlzGO zC4D{m=H}l|CjA5qwXr)uyE!e`)G0YQI)0;_22dqv7PR#A*g+Y|$hZ&1!AWN^nD;AY z$A$#RP6M5e=-Q?Xq?pTCtRX~nN!O3z4>NlbOTm|oMX*@rH@^h=`7NG5pOjUc-DJBf!4tbjV4FJ${D!&3k!_OgTPgQFM7mWb$Hh@p_ zL2>!$$OwQ2=V$2nq^9Zvf#6k;TbrkUo(@YdFS!CC>{)q08;vQ0V|p$o(9!4U!=G zlozxx7C=KnZX?l zVjP(MR@}eJ6e*w8&sskG{3=LHoo#GvfD=n^|1q6p)zs4RLB2+9c?N>`P(b2CKAKC4 zYZ>dqH;1P3(vZe0UbFg9T~Z=ks@Wy3Av9>V4SpmKC_RFm;#58w6vWhWXEQrH8|p{A zLPBX+lPD!_L)PEN*P|qWG755%>6MiRyE(<>?KYyK>gv5ufn?yW0pttZAl{|xA7-zw z02)aGx@Akt!p25#P0d%+di8@76RFYBHgg@c@rrY+tFVaB%X}Q>7aXB4_52msMF3?3 zD%)pQw3L*(fEuO6v#0|eD;F0RNR1FeB33{4zTOooRoRk(dBf?*5zg@!Vy~e5un-K> z*54~F!3GCvPw$;k6Rp+-)>d06NcH#k10J(>VD9V_UEPUOY3}0U6CItbtSoR5o9|Uk z6ag&6V-NU8kxJa*k&%#G1($F(>?7C25v$1tfw*>16oFAt%Vgos-MitM{6o53FO!mh zQK^}IPt2%|a;U~FqyM>*l9HsPK>?JbvY}ke7hl$c=j|)- z2?byv!MJwvMs;z$uiqgIXh4xD;f(CCaL-`VP6U3Xw&>Fk06CHJ<-5GcBFD<^&yA3Sdkd{W5&;_f>I)9ev_fW(3gWN$|a zdS0&ZkN|eARE6O5j!2s=ZaMftl^l{*bK@p|t1UztY2Xa#+H0BRWlji*h!S7FP7$*Y z3IXA}7Z~(TONx$pw;}f@L8^?6i@QK|W{IBymwj&`C>;73$oe1E4v81>o|Y4H+6zx!*j!kjEMKGPWffGU2n8GL##O16U34+uI!A$4h5hDB-Mr?@B@PBy? l#(~JxH5lvt|Bvr@?=0Ump1*oI6}aEX9XS