From a548cb834d98d8b68cafc3ec571c0324767f5b55 Mon Sep 17 00:00:00 2001 From: Sylvain CALINON <sylvain.calinon@idiap.ch> Date: Tue, 12 Sep 2023 15:39:01 +0200 Subject: [PATCH] eikonal example added --- README.md | 1 + data/sdf_circle.npy | Bin 0 -> 64480 bytes data/sdf_circle01.mat | Bin 0 -> 27305 bytes doc/rcfs.tex | 2 +- matlab/spline2D_eikonal.m | 378 +++++++++++++++++++++++++++++++++++++ python/spline2D_eikonal.py | 320 +++++++++++++++++++++++++++++++ 6 files changed, 700 insertions(+), 1 deletion(-) create mode 100644 data/sdf_circle.npy create mode 100644 data/sdf_circle01.mat create mode 100644 matlab/spline2D_eikonal.m create mode 100644 python/spline2D_eikonal.py diff --git a/README.md b/README.md index 6d73bfd..7af2b9e 100644 --- a/README.md +++ b/README.md @@ -13,6 +13,7 @@ The RCFS website also includes interactive exercises: [https://rcfs.ch](https:// |----------|-------------|----|-----|------|-----| | MP | Movement primitives with various basis functions | ✅ | ✅ | | | | spline2D | Concatenated Bernstein basis functions with constraints to encode a signed distance function (2D inputs, 1D output) | ✅ | ✅ | | | +| spline2D_eikonal | Gauss-Newton optimization of an SDF encoded with concatenated cubic polysplines, by considering unit norm derivatives in the cost function | ✅ | ✅ | | | | IK | Inverse kinematics for a planar manipulator | ✅ | ✅ | ✅ | ✅ | | IK_bimanual | Inverse kinematics with a planar bimanual robot | | ✅ | | | | IK_nullspace | Inverse kinematics with nullspace projection (position and orientation tracking as primary or secondary tasks) | ✅ | ✅ | | | diff --git a/data/sdf_circle.npy b/data/sdf_circle.npy new file mode 100644 index 0000000000000000000000000000000000000000..cbea57cdf66dd557dd2ed8dc7b7acf190158f67f GIT binary patch literal 64480 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1J<FBrxpq`drR8o|f7oT60k_r-bOUx-w)lpC{ z&PXgsRku>m(9}_=RiF%TH83aVmF5;y>LuqFrRwFD=9FY678NB{a>W;=Cg<lBmlTyI zmv9v_KvZ$%r9fm08NC@=3z?i5lM0#J3t55+Sv9;FHKB@9N-7IdxeD1Lm>3us{4EOE zJ3HD7If4o~HN2T47#SECY6`jh{QSKB|Ns9VOnBRu6mlnZ#uoC_7V>I<l~fe+WiVub zZQ;*gZZ8xFDiqZ4X7pzAbzo~N6v}`Yo@P-f3^H6Ks8AGUxJ{uLs^JVJh2lvaP7Emw zIDym>N2w)_v+Q9CApBKuzKHckVe5^X;e3evE;xSzTjm6|%mZ*fMEzkne|rAUY56~o z!TAvVr{Mf~H#+Cu=sXMOL(IDX=U-Cb1tB;eV*X_~KNJkE!ub&UuEF^?LB`L$0p~;P zzXj*V?0Xu$@9AwgAL5?7aQ^*Ix*!DSL)`ZO&Ig5#u=PVYALLjC-bZl$XRvkzAL9PU zaCuXYW)qL*CvZO4JE2eE{IVqVvLy9qa6ZI;2>y<m)7xuKKZnah{QCmV2ZjE${FiV( z#Q(40d~oy%zlQT6;ep`Sf{fq(1}+Z?pSN&6IQXL9!TF%T`lS0F&WEI{5B3my(bL&1 zxHwE6La&00L&9q_oDY%T1?NNJ;{cowQGXcDhlJNLI3J?_6r2z7&sjJhV%`Nf9}=IJ z;CzVrm*IR!d|id}A@*H^^C9tZ1I~xoe+$lsg!gSYAL5?7a6TkF@5A{J_dS5~Vd)QI z?kszVdmq8&Vd)Ph4{`rvxI8TV!Q>(Sc?y?@r9YTF#DCA=^04#=lZW{C1zaAM{$TPD z`W0Lrmi{2(v+N<^@ft1<OMftVNcg;k%fr$iOdb+m@8R;O=?_t^LDVDC87dzpzRDgF zUa<5B6Njb0RrZkhfTcf}dRY2fWe*82So(wMho!$&_7MNT(jUw`So&K9jZawmgP9LY zf2-^v@dZnNF#BNXZ<ReHK49q&W<M<bt%8O(Ed9aU151Ca>>=R^OMfu;!P4I<d(`v? zb1y9Yt+Gc=e=zsM(%&k3)bt1Q4=nwyvPVsSF#o~Q-zt05^at}VEd8yrM@@e)|HIPX zDtpxQ2MZ5a`dej>n*Lzn151Ca>`~JnEWBXp50>6Oz|#>VJu)z$rmI<Sd00Ax$)nQ{ z`CV}Ju=Izj?l4>)mi}P!5dEj%^04#=lZTjh0WJ?IPcFgv5OtT~d|3K}h;N3b_iJ!@ zSo(v>L+rl=mxtu{+i*U_J$K=JSo(wMhq&(nTppJGVDb?6K7z}`(jQD7LO+Jf!_pr_ zd^0pZJ%!7|(jQD75}wcC^04#=lZW{C1zaAM{$TPD|G$FE!_ps29ugj};qtKb2a|_{ z&s(@WEd9acA>s8NE)Pj>u=*F#ZiBc3R*u2yVHh8lPY~%0RXt4IF0}FwrXE)Q?SiHU zSosIj2P^+}p_P9y^I+xQF0}FwW<IR^+Xam;SosID4_5x|LM#7Z_QT4*U1;SW%ssI3 zZx>qm2Xh~+{M&_A{=wV}EB|((m47hz!^*#1XyqTwKd|y|7h3rT^B=7I+l5yC!Tbv= z|8}93e=z^U%D-J`<sU3OVCCN~wDJ!YKCtp{7h3rT3olst2TN}s;OPXC9vK)Q`50Dj z&VtKB(iyCtIRH<`sOn+yZ~(3kCJ)I6hvDLo{CW(|hnRl~&WGfuvv5AdoC|P1q&&F< z=R?fD4ClkjKZw2q(DZ%{E)OXWZov5v`)|Sdko<ld&WE_?E}Rd^&-dYch&>PBd|3Gh z5kCOUFOT5zu<{Qk4{`rvxIC=<gULhu^As)*EB|2fknns4mxq;qFnNf7U%=&I<sVEQ z;{R80d06=elZS-IYq&hD{Da9u!sjhq9#;Or<RRhp9xe|_Z?O6o(JzIB1FW2a)xR)4 ztet|WM^V+o@-ZTv!Q|0tNO&ECrem1-u=@8fy!?adgZKkh{~m_MC#?R3xd&GN9)_2H zF#QnwVD;}|c=-pDhlD?@{yhvY|6uZv@P*aChvDTPOdjGsSp9n#UjD)4A^8PX{~m^y ze=vCn4Xb|-!^=O2ILtq=`u8xr{Da9u{0FOl55vnpm^>`~!T7NB2jfHh537F<!^=OI zJS_de_^|W`<3qv+R{tJ`mwzyMSo(wUVd)RzZp1hUj1O_oENHm_t0!UYTv)n+wX0zA zuzCPduEETM*$d;N(-8Yk!R>?i=PaBLiQfxwJ|sRb!TAvLFT?qe^5QC-53%nWoDYeQ z8*o0v{#$T9B)o6K`4IQqh4Ufdc^}S)xbFd+4=ev5_8x=gmq&1USosH&hq(VSTpm{b z!Q>(Sc?y?@m47gK2>lE$4=ev5;>V!*?FC#Omi}P!5dXh|%fre)m^>srUc=>K=?^9k z37@xcd06=elZS-Yd$>F-{Xy(SjJv`35E|D0gSE?G?LSyK1?&I9<YDayL_G>I2c{lT z&cWOTlSii^{(<#>VdlZwf2W}F3G4sD%!jrAPQlARn7I)9VC}zC(D;D$e_`&2wf|1R z%RiWYh<jk|zf<t?4<-*we=t5Q{lWN<`~qwLor0HtFnL(|gYjYM55|Y3KNugD{$PBF z|6uLEQ}FT+CJ)JPu=d|6wDup&|FHJoDR}t@(+^93Fg`5(!T6Byfwlim!OK6GJS_de z_^|W`(T5o4h4CTkX2I)02p`tZJ`1fkH>1hJ+M$Sc8_YaN_#*06NO&DX(|-!ihv|dS zkodg-mxsjXB{&~q{$)5H5?@#0e29J5;Cx7Y+<@~T{<{U|L&EzuoDXrwT{s^Sp7-H= zi2EME`LOf{(RUV_Umn5bVdWo89^(GTaCunzgULhu^As)*OMfu=>EQAFXK;B~`3IAS z`1b`|9#;Or<RSil1(%1Fe=vDSc)W(o!_ps29zwr`%fre)i1=A(eti#@howJAxFO~h zV0?&sVEqqRxdZEez}jJ(q3ubSJgonL=qJF;gM=@vU5KbhVe%0B5ak@q-7tA{8kYWG z;t=y;{f`Uq@((5ti62=1;{v?=gULhehxI=$z{@|FJj6Y){>KGqc*6Q0F#p2(9~a>5 zKbU@qdtv>L3-Iy}CJ%8xtp9NV-u{Eh!_ps&4@-YAKE!{p{>KG)`3IASr9T)Smi}OT zi2q^zj|=ef4<-*we=t5Q{lWN<@PYL|F2Kt_m^>uE!ulVu^!5RoFE7F8MIhk>tJh%T z_Aowd92M5?-wf9eaX)OF0nsmonFk49M85$NUa)o%%zQ+>3h@uD+<=*P0d616To?^8 z|1w-25?@#0e29J5;Cx7Y+<@~T?!N`+L&EzuoDXr&T{s^Sp7-H=NH{-$^I_#5%zY60 z5nLWt{z1fH;rSRY4=ev*@(}+#h0DXzA50$Nzh`iHSosH&hxqpeTppJGVDb?Ezk<ud z%0HMqBs^Zj<zeX$CJzapw{Uq_`3IASgx7nxJgoeK#T$eV8-IqiqhS4C7$4UEh4nKx zL+fFfJgol<D~Dj?G%)puekvrqVEq`FendMC;vZN&4>J!|?jq6|%>C#zB)+b~-3PG` z*8jZ>FaKcrAmI<||6Ycde=vE7dtm+F%kc6KCJ#%0Fg`5(!T6APfc1Yb!^=OIJcNez ze=oz!KZy8cXn79n|6WGx|HArRu>S95c=-p@4@-YAJ}mvg_z?fY`oEXq<sVEQmi}OT zSo(wUA>jk-|6Ycde=vDi`h)Qy>Fop5Ux;~nSiC{pGYc*c;lt*EVDSnYmq(2A!t_Jp z;{aSg#9r7q0-|3EQx7p8(QknG=PX=5#5_d33W?84aCwONuz0!(7l)Y-OMfsvBtCAy z)j{mP1?NM;`!<{panD^i9}=GT;e3eu9>Dpq^aryK;@(Gad06=elZUwfF<c&2{=wuS z^i#Mzto(zB!@~0!Tpm{b!Q>(SeF2w;r9YTF#Q(40^04v`CJzaZ*Km1Q`h&?s!sjhq z9#;Or<RRhp9xe|{f3SFg@EI6je29Bs<1es&BW(NyR*t~NUtsdE@fTP-3o-5s3E#ue z{x@PA1l2y+_zR++0F#HrFRa~&s7GP)5c6T>z%{r(Ve;rSBtBr{FEIDR#$T?%%RiVt zNch6WU#`K+KbSnkeX#MDYw-3TOdgj0V0>8mgYhBmhmF5ngO`6Wc}P6M#$T?%%RiVr zgocg3T!WW?5b<l!@*Xz+at+@9g~>zW8#exO4PO4i<YDO##)qXp7#|Wou<@5`@bV8P z4@-YAJ|w+;fVvB@4ha@75ckZ2%R~6EdF&g|dTleBJZzo^F|PnK4-&qJabJi#VB=DV zaS)h#h=0z)&4<J<qTc|C&r5K5i1~<m6%t=p;qnmsVDWGRE)KH~;=fyPaY%UIhVvor zxeMn*!t*|y4{_fEI3HI2!OVxa_Yqtkmi}P!5cfZZ%fre)m^{QkPvP>g@((5t@!vDJ zJgoeK$wTNDaCuny2N8#b|0}pWto(z?L&D=VTppJGVDgagc?*|^m47gKNO--6%fr$i zEZiY{1_l@(Lc_-YVB=!2@jqBQ2{!%*lZTD}!TOP~c|4eU*!Ul;oPmv7!}P<(p%CLB z5O=`LgY~l!{REgiB>rISNJKpflZV6)qMU<;3rrrJhJ-h4{R_-Ju<^fJ@bV9)50?I5 zd|3K}@nPu?#)qXp7$4$(*!=G;c=-pDhowIlAC~@Le2D*G<A1l{<sVEQ;$PVK-z|9g z2a|`;u<^fJ@bV8LehXS3z~+B%q0Rro!Us0~cMD$r!Sut@AB+!4Zy%uMBGy@=hC6KC z7c3rN>u?b3LQvJi=BW|$rZ9O(_#)=DAmMckZa>6bh;d&?emV=6hr};pTm%xIm*Dac z^AY_9NPJy|%R}r#)T@yAxB-`k*bfWm+i-E1`yuYR3m1ok=Y2RI;=TuPJ}mvg^h4bH z2rdsxe=vE7`ya#QVdWo89^#*;aCunzgULhu_Y5u%EB|2f5dXe_%fre)m^{S)ui)~q z@((5tp<l!0VdWo0{5Eu4?=4&&R{p`{A>s8NE)Pq8uyBR&85m%Eh<jk;&#-(C8-IrN z17YLOFnQScGi;pf0MuPD^|0}0SUU<f{tVL(8-GTO^TON%iC@_GGh!SBCJ%`}*!VM| zp8%7e4j!L{wJQ<zC`=yWUsygyq%)X2It>X=*!VNdzp(M=yYTiOOdl-$!T7NB2jfHB z4;z2J3orj*^04#=<HOP)j1Tc2Z2b8yy!?a7!_ps&4@-YAKE(g9@#nkn@((5tNk_2p z=ezLo4<-+xVdKwt;pHDh{4TV<fQ>)H){%dJir<Hi|G~l);+|P>c?ci2j_^Lb{fDX^ zw$27NE_(p39}>QZd3#8B9fQk5%ty>?LEH@+H;0*r825$5=Ows)i1~<d5lDPph08<i zL-ZRU@o@t#53wInuR_B6He4R!9#}Zuhl|7f56K@7;Nr0I4<-+B?<2T8Ed9acA?|++ zmxq;qFnNf7p2Fo}=?^9k@!vDJJgoeK$wU160xl0re=vE7|6jr7VdWo89ugj};qtKZ z4<-)@pSN&%SosH&htTig@{sZlw*CsSZvhsL5ck05KVa+HVDlfabsn(!512e`{sT4- z09%IzQxBW}fb|n$>l9%6Ve=obas)Py2Qv>b&j5)}*tj&zeAqZ4VjKkKA4vSb`sIjz z0!$v_UsyX6QIEppA>oTC=OEz;%NH>BA!u0ogXxF37dHR#06zW)lZT~07$27YV0>8m zgYjYM55|Z14>tet0ABvV<YDO##)qXp7$4$)*!;%>c=-pDhowIlAC~@Ld`S4f=06_5 z%RiVrB)nkrAF%ZH0bWl)(jx-{EZiXOnFW`J@Dc0EA>p+dO&+oC5@PQGxI9EXVqFp> zypF-;A?BZg^CA8@3+F@3L(FSI;`0(*9%BAwI3E&USK)k!eTZ=pNPOIY%R}tH1?NM; z`!<{paSx(ig@osQxIDyt58!-QIKpU%dmq8&VdWo89^(GTaCunzgULhu^As)*OMftV zi2t6!<zeL?OdjIj7jSu4`h&?s{Qn9r4=ev*@{sU&4VQ<dKbSlueBQ$4VdWo89ui*f z;qs952Ah9F?8AbE6U05R`8U|Q4cPn}Y#bjp{|1wX&A-9&Ibz)<Bz$4>Z-{kBsP@6; z-w^ZkFnLJ)!sg!)^H?x>i21O25yUtzOdb+HuyI4gI0#H0;$K+58qrUH$wR^y*6u{q zqcC}h`(X2Lh;#;%N2g)w4<-(AKWzT(5xo5elZT~07$27YV0?)GVDoQ};N>4o9+v)K zd|3K}@ge?)&A&Z@mwzyMSo(wUVd)RXhlCGo{_PRG{DaBE(jSZuOMkF%fbba@V0?&s zX2I(T2p_R;1`=MI(c}^9$|3P_04@(vk63pJ39n;td5HOlbx9EaoQ2Cn%tOrEL;P_G zE)OvuF|P%Qud8r*h<(@Kd`NuUfb${t--7cY;e8v<hq&i1oDT`l`*1$QeTaG$mhWNV z332ZuxO$j4jE1=XF<c&2{=wuS{&@<QhowK5Jj8#`;PSBY4<--s?+ds*Ed9acA^v{_ zmxq;qFnLILyoSre(jQD75<YL?^04v`CJza(_i%YgdV{ThM(jI<g#*Mru=!uuxD;&u z7dB4`oBxH$!{&cs>o{QZzcBT%`Cr&LKWzRNrXM!{3+o5M=6_-4!RCKq<pg5h9^!7; zJTqb*3nC7)4>oUu80Uq_L;MRHr$mf{z~mv}3+u-t`Ux<3So(wUVd)RXhowIl9~Lh# zJ~|Cce=u=a`h)Rd=?}(-r9T)Smi}OTSo(wUVd)RXhxi{h|N9tT{=wv7=?}(-r9T)S z5<al`-^cLs4<-*we=t5Q{lWYR;WIG6_z?HZg4Y8OK4M=W#GK7&@`!yikoY(NmxriF ztSg6v*D<&}#QalmKEyv~;e3dBh;>Pj_`C#{hnSC;w}-^nRk%FFKE%8hBtCAy<stSX z#(g2-eH$(haSvi#1QMS2;qnmoJ%ICJ=?@m35ceYLRap9i$wS=#7_J{C4x=Idc?y?@ zr9YTF#DCA=^04#=lZW{C1zaAM{$TPD|G$FE!^%IHJS03`!{uS=4<-)@pSN&%SosH& zhlJOAxI84i!PXxj_SM7u32_f>{Q+#>6>R+hY#jz{{Q*oKu}=#UAF%ZYF!iwY2e5g} zW6*h8n10y$1H?K@n0p}R!PXxj)*-><A@K)We}I^$hsi_Y2R1K^n8$+2L;MSyH$jZ^ z!sH?03maENjDx`BA?}0qyAk~am^>`~!T7NB2jfHB4=cA3=?o^1PQ%h4OdR4r*!qK~ z@b(`}9+v)Kd|3K}@ge?)tv`4QFaKcju=EGx!_ps&4+$UG`h%zN@((5tOMfsvEd9a! z0pT++!1xgN%!1EzL->e&&yeuij3$rR7YK=u18{kWdc?jNNO&EC%R|gR1?NNja~95r zn0Eorhs5V4I3Hp@VqFr%pI71X5c?4G_K^6v0hfo^e+$lsg!gSYAL5?7a6TkF@5A{J z_dS5~VdWnzJR$CV1eb@UKbSnk{g2`DuyO+?5An}axI9c8MnnAf3@#5V|6uYE|Gt3B z!_ps29^(I3aCuny2a|_{$7{GeEd9acA>s2DE)OgJVDgagdJmU}q_+>y^(csQLty@Z zxCgfW1vak;TmJ%EX9ipU0+WZWe}Rod!PdXP)Wg=lz}9ITgWC_W54QdVHjWS5e*`lR zw*Cdy&x7qhf|(Cn{{kxq5bKg4{(`N$LCn)b#9{7-%}XQZv0(C$@P*CWAjWxN^04#= z<HOP)j1NnHFg`5(!T7NB2jj!iAB+!6e=t5QUSWK68kYWG;;{4w<HOP)j1LKa*!q`e z@bV8P4@-YAJ}mvg_>k~{t$%q2FaKcju=EGx!_ptjoe(|)1B?%G&n$R34&futb%2D| zW;A)kzGp~$9DvJ1)Fbu<Lfml-E)Ovuv2O<ApR;gzh<S*0<&gNi1eb@Hk63pJiLa}0 zd5C?8bxDx;xB-`k*pHaEhlKZSxIDx?h<PnYc;1K0L)?cL_l2cDSa?F*ix?Myr9YTF z#Ql%q=E2e*OdjGNM7;`2e=vE7|DM71!^B}U#J?}#^04#=lZW{K6<i)x{=wuS;qe+S z4@-YAc}V!Yh0DXrKbSluyxzm*A?Xda{t9s}63m?t_rTU)!PZT{)?XpcaX?iMTYrVv z=M0mFgfDFU6>J?EZ2c9?eAvD;#6B6Aen|Yn)?Xpkk;CL6=EK%sA=X*K<RS3`TYrUE zhXj*{_!qYR3NcR)lZS*aY~C0#j|G#5xDPh3f*9w8$-~kgj1NnHFh0cnuzoqBp8%7G zr9T)Smi}OTi2q>aIwGCH<k4wZ`h$r>{101y^#Wf0!Q^4-55|Y3KNuepKCtyyFW}`L zOdgj0V0>8mgSi93XJCNwA?}%lwoV6eZVe>7HlxWS&UJvq#{sxJMEzkn9}-^2;CzUE zh<$+&|D1)(L(IDX=R@N25}Xe)|1z8piLa}0KEyu6x=Tnn+<?nN>_@Cif`s>NxIDx? zcj0_Uc;1KeA?|ws=flcBSa?F*`v@)%OMftVi2EPI<zeL?OdjH&r*L^#`h&?s{Pzql z4=Z<I@(};NfXl<gVKl`5ui)~q@((5t36Ixsd06^`$wR{DEnFT}{=wuS;q@La4@qyZ z^Is6>zQNo9p<(O)VC%wQ>;GWmM6mULFnQSeKiE7XZ2cchJ#76StepW{{|D0#TmJ`J zhXLFF3o{S4{tq@z58Jl{Gat7857v)^?F)d}2V4IK%jbx7mk|HL)*&L+Awk4p?t!h- zLCn*`<YDO##)qXp7$27YV0>8mgYjYM55|Y3KNugD{$PAq`h)Rd=?}(-r9T)Smi}OT zSh|4m(P>!vgNeh^AB+!6e=t5Q{lWOK^ata^(jSZuOMftXAbbV}7$2f;7Q7sW@Db;7 zLBeY@nmpp%8c2K`fXhSFBhGbzgx4{+Jj8s&zGsMk&cfv(<{|b4LgMohTpnURV&4oT zzOKUMA@(8Gl|$m=23#IuKVsb_B)o6K<st4ttV@D~=Y6<5#C?c)dszB|g(t+lh<Pno z`h&?s+>aRdg{423Jj6eUaS>ShgULhuhv+xJ(jQD7;$K9)3QK=5d5HgC!TkdhhtZJm zcnz0_r9YTFBz)e&<zeX$CJza(_i%YgdV{S$MVxC1vj^fH*!ol0x>wlxQ^Yx3sOn+s zPZ8(Tz~mv}3tNARIL84d53vun{uH*(47UChW*%()DPkWWOh3eY*!ok%J{g!iBz|D) zPZ8_LVe%0BVe3y3>nvgNknn}A3q-6#g2}_uAB+!6e=t5Q{lWOK^atZZ+z*>~L5%al z<YDO##)qXp7$4$4Sic?7Pk_n8(jSZuOMfsv#Q(5zACb;r^5`@y{lUZ`;R9QL`WoK; zgUQ3vAB+!6e=u_(d<F&>AL5=_@O6d|KH}U#NO)~VlSiD(1&NOXaCwOO!*D(%ypF;7 z5c5yL`4In{h4UfiU4Zi;@p%c(hnRmE&WFU;RX87F-!(WN5+66<e2D#sb>)z7x(%0y zxCgQB5)z*G;qnmoJ%ICJ=?@m35cfWU%fr$iOdjI?$8dRA`3IAS_~$8H9+v)K@(}+$ zgUiFpKbSnkzc1kOu=EF$hxq>$Tpm_#!Q>&~@ft1<6Nk}|@OcZDhn0UYc}RG@hs#6K z+Xv|SWW>3>FmoX8fvtart?P%ae}>KT!PY;+<YDWdVe7zP>z`rjVe6k^<3zCa&oKS4 z_0O<%v}fV&gTyav{WENx1?>C{nE9~v&#-ZL*!deU`(W#zVf{4N`5Q3%Ve6k^`S>>6 zen_~$)@38sSwh5N{)MeuM65%C$-~kgj1NnHFg`5(!T7NB2jj!iAB+!6e=z=Z@cM4p zxG7>B1SStle=t5Q{lWN<_=dG}5%nlc9+v)Kd|0}I@zH5m`h$tX(jSZuOMftR5IzF~ zj1O_oEO<Ez;Umsfg@o5;G<n3igOK!c04@(vk2seL5?;sP@(}Y8=hi^nc@{1YF%NOB z10+5#!Q~<5BlbN*;_E709%3J2UmzqtZouUs_9OPqK*IYrTpr>c#JX}wc;1K0L)?d0 zcL^*1VBraIFJfI1Ed9acA?`=a+r!cyOdjGN#Jm<P{lVlR{zHuW!pc9GJjB0<aS>Sh zgULhu{|fFOSosH&hlB^BUWKJUm^>tW-oo|6#9=ffyxzm*A?Xda{vL5IIZPeIJ+Srn zuyZS6>+ccgsG_Qet-pt@bA_$HhpC6Hzek+Ig{mL6{vNhZ0Ji=fW*%()J>nb(n0`q7 z!Peg+_Bq4kA@KuSe~;J)2$P5S7q<Q$u}=mj4+&q``g_DWa+o|U{lWOK^ata^(jSZu zOMfsv#Qm^!HHdk7m^>`~!T7NB2jfHh2b-5cjPt_eVd)RXhowIlAL4&lzaG&~fXTzs zAB+!6e=t5Id|>qgBAvnH(P>!vL!}`+1_s!;4D4Kd7$4WU(J*;j=c>ZidBM)Pgw6ZG z&K-p5hn)ilYZt)I<$~=Kz;$j7%zW7W3$S$+uyY+?_QCF7fQ`Gu_C3Swhuyyb>zBdy z1;X5eYu^maeYn<@qp$P*Kwur*2LkIvJ`k8k{Xk%z;RAtjxDN#SsUHZm<3125=RQF9 z8-1`Zlqe~bOafmpQduYkyI@2bbis&BP@yd7f)PgW1tW49qMc3*MhpxL_fm`(FTZRr zk&*VLrs<Nsjqdi+iK!Ru?>$U!?+CnLFY)AfkWSw@`+(?bxv#xv>>WH@S28a-WxwIR z_e{1+C+s~Yy;yNT^0<AB#cgHQ+N1U>?#vh5{P(bZM)K7+>*gM^pY!g~&0BE??L!Xi z&E$1HU|(_hwgKOaefAa>Jp0dn+hgx>_0nFOrMvAvq{p_OS-Q(U;Y?3V`mde#5mwIo z;+7z2m^e%wOdrgg2B^6(dtmm$+yQeJ%$+cI!~6mB7tEh9f5XB77A~-GLJzkoAw8E{ zLoeDJ1Z`EGZg9bVhtBhRdo<43Ux-}TsK<Q9{z~DFD`KBc+CSL9|6ld{ar+#%iHFoh zj@g&gXf!f<AAyHg%$k|XwtNTe3--w~b=dB=&)B(<b@!q@_5!mbD$N3R*{h^BaDDOI zVQ;bYM#YNB+w4zB>-R4+*lNFE>1@+O%3JJ3wC~R|&EIU#QN+HvO9Vl~#9``S`fQ-) zz|4i&1G5+A4w$=O;aLE6H_RU}f5H3-^EWIU^q}Da3ny4O!on36ZXwgspI=}*XP>ZP z*0%h;r|nZ#cGUE@owR?GASHh-`MCXxoKqg$;YaPCY(HTrk#yL8OPzr|Tl+!#9VZN3 ze)R9Rmxye(DDvB5zhRe?mW<y{dzk;0K*MK?__@FjOq=XCJfAbS#B_tbfzCCNfV#Ez zJ06@^yvA*{eaBJ7B=eIi?KxgbMD%kaXqY%m9ZcUAs5v{J=EA}kW-rVgFn7V+33E5h zA25Hx{0Z}S3N##G;Q|XMSh&H`CoG&{;nvW+O?_I@NqGM3$@bg0^vMzX4Bl{0<)eq} zTU-Q>+V4GJ&+x=$a?rfJ@O%ezH!QvII0qi~pR~b#22=dg8kRNoGalYp5!12UKA`*} zhw!e&_AlaopPja7f&GEZKeN3i&#^!7_3Hi0Su^bm;!iUrUzrZ4Vd5}#Fnus{VCGJM z+5@u}<_?&<VCBG%#}IeJ`~mY9%%8CQ4GRZYxWK{*7H+V74ok1FaEIw%Q^1-Lz3+&9 zLHK@im9j(dats!quzcg8uE(*oc8mQCR-+Y<-fe)VCy74aKL;I`+83PT_W!YUzWs*I z>^a3Y)9eqlAO6xA-D7{?)r0zz)phm_cXx{w2&LLHYzXug(R1I=kovadsDCnohKj?~ z!Sp?VnFBR90csD-URZjCxyuLYPMEu4{(zO6Fn_}Q?E@{}6QKDJ7EZ8mgO!)Cd=E>% zuy}x(cOy1rg@V)pM0wvI8?n)OyFCxIoO%!w`Xbb7o&AnuM=qT^whSI#0_Pvfad1qA zhZ8(pjEmCjKeX;Qs5O|k{?o+$0=uT|KjA!mf5JQcP3v#Y+Hde@P49!{bKo>g9HtJY z?|~7-9GJPVaD~|mi#M3NUO?RmD_=iA{Q;|wVE%;pdjd4yy?~ZqSD^I(EZkuA4y^oy zl>@N!42u_-`7G1iE}nk9%|68b{j!X>P4=+-0?S9Re8&*<WTSOgll_72piGOxg8dBX zt<1uP6ZSVi!$-l+HMHjU!u=N*F5cd}b?N>EKep)w^sU%`<5O5?bJeQ-D?CmY%`I9D zr(xnSbufJ}a|)p5!t8<D3v~y~T`+fcK=V1wA25Hx{0Z|nEF6|V%Y9fl!NLtzf5Yk_ zSa}L77hvfd7EdtyHbBdJ3vot?E&7Y>Vet*iUkTG6S1(TQgojtgHTyjV*B0;J&=B_d zK=aD|0iUn786>aUZ{l_Cce>%m{Q^d7&Uo(Jv|qxKDPZ!U&HF91t`}8$Y=P4-ahN)o zJ`bolFmnr__Q2u^?ha`A0CT4S)ZH+D!2AXCC(Pd+&~hFYF0gQdg&VA#gSA6o^%AUJ zfRz)ld<si%uyDJe#p%_YJrQ1R!SdA)>ASX1_!q*{mByp3A08&G-_HXLrx2kU^SfbN z_j{-<SzTYZegBK>jMsr`JNNgz(K`0Wb=UqD!@9>g*AX;K9HtJY4`vPr)La9oJ+O2N za|g^_u=4c))ZH+D!2AXCC(PflaDasiESzBB25V2k+C{MX9ac}l>IGQ&22202c!as< zL3aJ-ply@)?}#yBw12T^{{m?K^YO`6u;09Ke*`qVT>8}mxo7Oy&+&7a@b}E!`(Lcv zxN_D0z5BoL`5WEs*}wmnYV*{p<OBPa>K`xW*^HoJ;xKhEeJoIO9H8dH?19+}a|g^_ z7ohHB*aWGkVE%yl3+7LlzhU990op!;g%d2?;O%l~yBC&jVf8(%-h!3CuyO>JPhs&2 zbDsrcON2<y8hAR}<Kh2P|MNC@dis$2P+KlyA0m8u{#PCSad7{$GuNBzz8>1I<>S9R zcG2Pe{ioi@nfyAue*%<-iNn+>LG{7H`3=-uSh|JT3rp`XcWr>iFU;LAf57|&^C!&T z4$$^FEL>pia9Fs(!V}g`hqX6h?Iu|L539#u^#ZJ1f#nZaJj2`z%U>sMYTsP=clZ7e zDWA{I5jn751?q2YXgHm-FbUntesq8JyvTLm<c{s1@v-E))cs@o`@7qVWE_vfX_z=n z9ZVn099Ve5?19+}a|g^_Fn2P@K*|@GKVbfX`4i@E8EAiK4YWT13ny5(!TX8Oeikg< z!P=j&b`&g~!|F9y`3x&(VEF|W?=bhn@|{^l=eM*4hxS+BPBXl|=E(l<_m}Wx)E(P@ zuJzp$XU-G*lj9SnoOU?5Uvlrmi+3eX?Ps|7y?9Ogsr`&l8YT`?2h#^LrwVGW1=Jpx zy)bvc+~ol+PgX$ddze39{(|`v?r&)OAJ*T5^$TI)2CJuG{YF?n3)WwN#Rsfi#Q<%G z!0J<2y#OnBVEG4@4q)MS7Mjjr=?WGOk^DCgl>Rui-=+4N-p%kc`=#=}>V*28-JcP< zb<*y$2pT31QwP(R1P%9IXnKXEOPIYdcfi~Qt8ZZS49p)ef5H3-^Y;R1I|tS;goP6< z++gD&u=)npKZW&&VEqMH`xF*mu=WBhU%=`GSosFar?B(}3%AAIPenegKDnQvIBHAw zgVXSE@_zI3nn?Wl{gP{9XD_w8uz#Lg`k%&51Pv31se^?t%$$W#b0<K<5oRyU9WZyn z>KRx&2<8u%zhM4^`FjDhzqtY0zlMbqEZkt@Kd|%(>(|5TBUnF+0oq@HwO3*7H&}ea z>Qh*G4l9?wL-Q#ty}|rj5jX9MR=_!Uc%9Moi218_asP_XpI<(>aB=@F)>CGU+Lz!o zOdO^T7M|r$bN&}W;v)r`9%1&v+yQeJtiFzcwg+JTfcXpNPnf?Wp#AhU&~XY_IKjdV z*4~DV2f@ZQVEuDgdWO|kuy#MZ{R(TpLE9m)avYYAVC57nzroTC%>T=FD4#N_y0~A0 zNowWlu1otbE#+MCH1hI(0jEEj0cQ|2OdO^Tt`C}C)1di13R)h(?1i}l<}O&f#|P?e zm_K0tg838X?+?)N4Grix5iFcw;RYLzgynBoy#pKPfYtA?ekZKHfTe#}dluF%gS9VU z^&~9a!OAUI{)43>Sa`t0>tUnt*$Y?pPus99t~={0oECtJ!_>j_orb1cFKGUEgO;N( zdtvT?x$6V8o%sgZUWfSu?k{LR9Omx}&~YJHKLr*}uyBL54`J;U*f<+(ya?7lhV|cJ z{ZI*LKMU3`fYo=fb{ed`0joD*<v%Rl!SW+4UBSYq`t+5~T)Aue=WpN68N+xTPQ%1u z>R|d{=14%(p$xSCf!Pal2h3eRpzZT1(0&fgA25Hx{0Z~-1?YHM19V&x7EZ8mgN>8J z#v@?kmaz6DtUQH{hrs&vu>L8mUW1h{uy!x3zJRquVEGnSFTlz-SbBn`H(0nm4O#oL ztKi0d7!4DLse|c*nKK7k-(*7bG0a|=J7DgDxw8Y>KZN-M<}a8(Vg7D_jzcoQ#&e<L zO|Wo-(bNh*SbqWLZdm(`-0*t_4M$SL4`wgS9WZxMJN#hj5SA`S^Dkulk=XnTD+gia z;%NN=8OJ5I{(#lPuzGp4{Q?;eC${|pYX`#my|DQ#So;|^zW^I2hqW_d<Cd`YGpt<- zD_>ykR#-b0*3W{qpJDA@SiK2rhrsGnSiJzFVd5}#Fnwjv_VYhz`x(}5f!Pal2h3eC zcf$IcFn_@O1@kA&->~*GtX}~OCs??_=AB{v64*Q-tp9cZIxhpOk6`2EuznA0+!EF= zf{i1=+N-dBJ*<BU>u16G3$S)Cto;TH7g&7?qhaDObufLfa9aTFzrpGYn7uG}z}#g6 z9Vdjf`(XZn`3vSxn7?8DH(0+G7EZ8mgY|=9{bE=-2pf-u&1b>VGXr$|05(nz>xaYQ z6V{%A_1j_NA+UY`tbYpYXTkaluy!x3{RV4?z-X8_OdU)gEWN?P?E-Y10%kAF9WZyn z@{0v@oEqj2n7?5Dg!vmbjxq;2e+>&KSh&H)VPNAju=F`Peg&C#6(BHv1sey1jSG&B zpF-xTVeJq?<EOB3XxO+kY@8aF4q@ZiuyJizzXVpl!}{AJb^IPS4+2Yv=<`Rg`FhyA zGpt^Qjnl&0nb3I|*!%))oE$dZ3oCbE^$e^YgpG&5`t`8>DXhMP(J*nCI+#9~Ik0*G zW)I9>m^)zZg3lj8*E7KU0rMBkpD=&J`WLWpfrS$+++g!?u=Pl=c{<qo5LkT!>)*lV zoniB*uz5h(d=@M{!sZuX<K(b<23GIE#@}G$NU-q`Sic_DKZVgSahN)oK3Kg4GZ$7b z!0d&&1MV(pxNU&WU&8zW^B2sYFn^=ZU&6u(7H+V4VA#AcY@T>@{v5K7jL`f!Y#jn@ zUBc-44aj;<So?&~`VH7R7}&a)(e*=+byCExAA+sJg00IMUB3odr$_AiHP|{p*t$U2 zdP~^4LD>2?*g7s)dWMw?u<{jFPr>HvVeKYZ`x8D72%XP@^($fX3$Sr=_;@dL+!8i^ z2@{8@gXx2r1M6qO`U^07VeWvr3+7IkyJ7x-`3vSxn7><~>qv2}pM<SLg{|*}ty6`m zhpiigt$%~9<AT*Qu=*OdJ_Oc&gw5B(=Br@+LfAYYY(5J%F9VxjfQ^&G#(QBjOdO^T zt`9md2J4@~`dKi0VeWvr3)T*Sxf|vWn7?5Dg!vm54zO^6g%d2?VC$e^>!M-ndtvDj zmak#uBdnf*)z`4~Nbq(MbR7+>KLhI*!upM{b^*+PuyHe3_`=3<VeJSQ4HJi{gXx2{ z^I_(~?19+}a|f({0COj-9Rl+Q%wI5n!u$;j2Uxhl!U+~`=<D!d;R?&wu<{XB&%o+y zSi1+-4uZ8OVf`7HJ7N86Sbc^*?+hD%gpK3E@)vAg0Y<~bVd`M|(8rNr_Q33gxdY}d zm^)$ahWP{LFPJ}J{)UAEEL_mTtx&$CP$4OTiGhJ3rJ_(VgCT>_o1wK(3AANIIjB$t zv}J?|yk$f+L$t$*A%%fqUd`2_bEi++!{{9f+uj-7K4tHGGvV0$s+0Eln&}+^Pfyw} zc;t9h!}hrSp<9doE5ABnFWsrY^jPGm{lO%T+=SQ1?N8WnYx^)Cw!fXpxboY}WA+@2 zw_DxgJZP`eCl(g-{HVR>l+ZMekbU-(rwRNBe0s!wYTn|<5~p_A3pri!dH(dUeWDJJ zq(Q}Ydn=WbUus?-vR{*P?i63b7JH7AcOi#<9JIe{u6g=I(?<KAll4>92^_T74BGKY z=iplV|0g*&6`LKfue_nn<@0fcy+ls^m#eq-*_&AKl|Puc#6EG_?5wYad+qoBaOBaR zv%r4#iV~047JKa9uTo-B<eg>DVD)wVYt!BKUw1dnSM8i^-)`Y{CM$23y^5gFvOCh< z_76I8W?gu;)BfJPgo95X)Z6#3+Az&*(N6mw&iE<X9I0@=r`JBV=^Uy1Vd`HWUUpFW zVf}uXdAEKac=cVndq2!R)9!tR#+{S*!`ySd&rCp{ch-J3-Bn=sMwBczGMcksKg>UC zC7bp=pSWZ{%)dTLrwfBWuGkL?k0Zq;<t7K$?uUif)a$oqo@v^+9~Pe4qSA|n61MDz zh4;hl&%PU1Y~K%ykGt;kj>??cwI3E=Z%=-A=Ly-j9~Pe*vUX`c;5@h=7T-s;K8yJ= zAKrh!Yce=JaLo35`9$RC{tXAp!Rh5=cI!4R+vEFT>B%-^@`8m`C-=kBS%jmc>xI*& z_ruZ|Or4(ium3ZfPuRm~)~6RtFSj1IulwM}736TtzG7|1?N5!z?2XwbpWqNaVt<A^ z>%*;zqxM&BWu2?|d(b|cL#tzB_7VG6oa%asj0fz!Lgf-)MIW{|<y)cj#bJ;Ax`TJ_ z`FkI-{}duU^Yo4#_6gD<HC#3a?cI+3)wIdoYR^|{(YDm!fPLDr-lkm%o9uTVHV>7L z*k{kk?^VgKwcg&i=X0NS@*evu$$kgkimtMMY&GZbn$BJJ^?Q;R#J4Q9pUP<QJ>v2X z`~40JC%!0LXm2@N+Fr?HyM6V^dY3KhXWLJxto`%+{Z{)b=l%$(h^h7h-*p#AUfN>+ z^8cLGO3nTDeiJf-1HW#z-xlAQ)IGV?{yKL^@B6sT_O%sHn-mo*?f+gqY=1~-v;FS< zLi~s865#yRH=KgR>l60F)Z5?A;`LFi+z&I4wXA|?;pEo+F#9?sbsk$Z_wSzqbx-T} z<EK<3rtUv*K@{xX0_QVg2iDKt5A)B8OIOysD_poA=HI94)83@EEZq+a4|6Abt1qIf z_QS#}T;YkdnAZCJv!>{S!!u^(?k9T_HtmmvhWFdUVKeP>x9*3<hiG7p;Q1Xp_QT>U z^}6kX?+$zR!{U>nKj&={<AMFK_@0}&p{MHa!TqrG;AZ}#nM?S{{_K*2;Peu`IpScj z!?FFabX2&oedFxr6Y%urd$)Y1$?a47VdWaEoP(t^Sh_Rb^s_J3`lvmO{ulc(>aE@p z`<3%-%Ga<Twijh!w%aIu*q*!dF#nP_2kh<Sl&&0NI%J<dclxvSH}=|FZ26(k_4I%} zM_BHb^@6+Xzu#f?T++SYe)@I0bN(IM?A5NnI=sSruYFkO*~j4yo9(xA|2^L<u-m@w zc=LWK+YR<c4=n!mKiFY^`%;~F9@A?3{2r?bGKaU>7w)~rI#*_y{nWGnu7n@hVt;C~ zLw)Yah4!X@&bTT+++^>)^YhguJag<N+2s<>>20)cPvaJTIdiK0oT5w1^HbN`?^%Cs zS*z*<`)?Il$2PoPW8bx1=<sW$cKc+O^{eX`R@>k4<`jH?yUzY-)1G$8jVtY)Ew=vo zJh#wZN#!cbC$^RLm*uuR?Rgjm=SMt$eT4B**nXJ$$43{O6r5MMA7<VWRkc67ck1@T z?3-w%#`H(2eg6`udsgk>aqLr_u)ios3GCjs$LCEy&78XbyWR<~e{!yRwXfruvmfSP zz7t0CicT)v4-1b=>CEYiWR~rRh1YuL(k;bItM@<p_ZA$UcX|7k%i3<(4-0S0b$1>_ zIc(k!i;s6(18)a)Y}*fuuQ^l8&TJCgwI3Fr>vxH>ZoaX1KP<jK_^K>e`sToX>ly#R z>7lFh4BI;P!~0?BWVz9E>2&L(`(f!R^T_&rS6Yw5)0@Xf79OFeC-=kZQCPhSE7xG@ z8kWvr>F(hP69vQ12kl`rW9_Fd%WDVh6@In3{9m}=zW>Yr{@mpK_8Vt;Mdwf5ZGSoD zLE{JBz4kiW*Uq1Mbi2Kv1V=#S`d#+nbNw`<&9~SuUz3v2>a)W>D<$t;)B6qf+?tcv zH@w<vKSjRK?%tEt_R%NjuRlIvvwh&wkEhqRFSE~b*q(URZKHkBl?AM~&MdV5*}pzp zR%e}kyVd;-FLGzw2hDri94@%Z-nlEM=V<2?`*jYN|Ck0Zvwx6e{nznGul;p*)^`dV zi|va`Ro&PRx7e2i&7Aktb%8zCOh%#Dxz+Zobh}^gb(mxSMXx!uNWRd%q3Xy_seqaG z4~kbg?#oEB-@RX=V%o*&_G?`l{!Ni^hx64!y*_hEy6=anm%5@BcQ+$xKg_&253`Dz z<qP-2?0e{>7gsR1djI*!3}E+&^K&~09d6kVbMM{rSp|MadiUFWcn0>5d9!TSna(Nu zVg6MwQ!Mz9J9|GYJf_d@aJYA7;eJ?n?MO<F+|<5oKP)_5%6AGsezJQ1LTGqbhkPh& zeZOHpEIuaL>^K=~zGXiwzHa`y9600X_WiK<{G+!*vS{k={jm7X`TmoWY0>`uvg<U! z>EZP&^%kSg2lvC$OWZTLeed*+?1!Z%!3icdUmB0W(_3v_XpHKs6Z>K9Hdwt5t4Cqw zEUa9ErF&R9gSjKsN;TtW=w5pm9d`U=%&u>{>{rb;Kbvg1(>~KCOtQaihdt+xtX<hL zTkYLByEgCS*k-?_Oeplnk&X7Y_wR3*-?-VH^TVdf#~;_&^SMOpe*d<?eoxarTjoE@ z>^E%gy>aXQ8vB<qon1AX7TTwVUr4KawcOtJ8}s(xR<rEyIP!HBDJ{02bLI5poR&%U z5A1JyM%c}@pYlWQ>8-X7`&Q>Ho1I^#*tZ;H%*feSYp=nyp?kG;xBapScYa7N$hQ~0 z_f%79Ppy55f@o#I?I?TCn>}H3&ZgL3z0#q!<+_!9N=sv0jgb3(8>3T!TkZ(%pChPQ zW9*x}|6xT>eCiuvIR8!25!uskgyHHVe`-`-x+4fTPkRFAvvb$2_QULZJug_P`F7NP zcc^<#PI%I(wIF~0kqP-=_s)L(cwNQ5+Wp}hE`$B!e&Dvm<F=0dF#l#+U)WICGHE|7 zJmw!xJkDr6Yd<W!ShZ&{G;CV9f0tSUI6SY!R|&KKS+;-s%BSG)K6vix{O2Fn?1#t4 ztpgpuj%?fyi!TE)W%j(7t@~l|dHlPlPMYP;{jhlEN{GJvJ9O{<kBoD{<xaizuftZ? z4(#WE)`O4ZB`P*aAKnj3PyB1wYu>3i3Qup67AF=bzCOMm)=z-7GhyvESiKLcM`7hK ztXzZT3s^dXx$}9Z@3aX=x7ow!IKGCRC#<*Ft8;7i_g&j$?{e<U|EKdd+Ka|C&CEz# zXa9LZ>|>s1Ywe9gd%}yiEw}e^&AL0YdWHQopJ!`l?OJGW{<oXKj(4%WyUg>`cfQTA zSNLW-DJpuly-=>>f(*_H_8Lx`VsoBPun(C2mTPfMqdf!D>ql%CtL>j#O74mJkYOL8 zR&q3xFKK^>d-PH7RPX()B@XTnKDX?bQ(fgXg(+))YH(a@X3@m`w(IVxCEO_3zvRCO zlfm0*`)kG3_QuW6-(NbtjM;eWto^tCuA1Z{mA2oh_$XWJ;yL?8lwVpbp6L(giyeD3 z(PfsuJxu-6>W5bgrPAzS=1t0)q*OOQ-yUY)<k@?VW#1^Ue^v}}4}(SAk%dfI_A9#c z!S2<1KK)r(s<-__w@+aIEa7Q>mh>TGKg_>9Gkf<fuW8&53lEMzZ00$f6ZXTxt9r`s zx`*Fp><@;9XQtlVrn$Qo?tcpn?^Rq!CYNnnz8@AJL3QUea}w9>hs778pMAxIYn%4N z;_>MWp&65oZrcxw@6Hv<1$(~j+CQZ|4qVS^%Vn?2Pu{;j3EJ-3rsDMMFw>#^u=Mor zD8uYc*+<~%?e&~3yS~3Xwjb6nh4p)2{RCJ$6xMEowF_YND6E`@m20qk150Nxe=L0I za7Nv9gFTFntnw^<_+^#7#;IBCTfVNaXSA)=U<_Gm&vW?eZyD`H_IJFrIVP#hw_nS; zI_Z(c4Er}+4R_zwOt!zcc3pi*Znyp6z}s<G6YA~de_mAIa<JI`R0sPD&iy(2kIV_O z+VmlQ{}s*)ZTY&r`~4dQtd2?4?Ek^fe00^_8T$iU`m8TIwd{Xj`SprD|HAzptXaCh z_O<O_b@%1GYs^dcuh{umQQ&mz{)2BfJu_Ebv0torR_xM8jr(QT)V}ZruG;_QOG;~2 zO8I`dw6w2ha#!z<i`c&AxkCb+FLP`CGAqXfdzktQxA>k;OewdAnYU(YBKP)3jrK76 z8XXI{<xjWT2YvVec8}0IwsWle+w6Ciaf01DP4A)KZKoFd=tI}R{<*bS!|1$34cxz< zdJcB)`VbEfk5hW@Gj<#--VY0}os8o7mAT#fVc~gii#Go=iy8ZW+d#tG<T+EYy!N8~ zu=rqpBX)e-*A@F=@znH#IaJeh!+uzN?wwVXc*=Uqepq~`SG~VJp>4<h#TRCS+f_bd z@BKgV?%mG}ZRdY%Si;l&^uT^tda8ON=KnhSFg(49bXdB?K0mr2HVy*ohr{}%uznG& zp8#v8!rE=Hb_1*)g_Yy5at)TRVCf9z&$(x}i`RKBwujO8rre6^IWo&W_U?*n;R|Ni zySv5b2&wkkYXzB`_*8e;n|;>~zHuzU{%*yN`x*_A_I_UfJ=vLB_g`w2Oq7{av)^H} zO^*`W^!*Cv+TEWQbnjmtTj5tSbK!p0kiwlSeJAaASvxafPSo=KOS_Jr((IkG-#~J@ zMELYI``;_Sxb=MMl>HHzCdR3K8}>8S+D-V?HEBP)!mb4-F`M>JajSUfn$WxdrpMk( zUGp~Yue<%wXv)>5{Vk{2e1l!K>_7g-$o0&JOgMl29HZ3%A2aP?>h~^|xwiRglReBl zL9zTpNeR97F#FD#pO$0qo@B29b&tcty(XWgPO<+veKOd+KJQgd8}?4I|FozW?4QjR zr;^wEPO^vjH!!iK{MUkRczEdF5q4CXR09vMv`d_px($*0Kg;`q!}DzR#!dd!9sAF} zxCjpK+1{;ZqZZ8A4~vKBo2>#3o{RUx;;UeD_Jzk^R_%wy=gP`sLC@!H+z*TIr+lJH zyE(S)7oQRbZpT-1{J2=Pe%JmB!GYj@<=R)B*Ozwh-w#Vq1`5ko26!KWr?*rw*CRns zkL-tyyTQhVVB;XLaS>R*6xMHo^%G$I23Wfd){cYKqp*4vR<6OyIaoS_r8}k@U3)hz z>ad5=S?kUvTgf@@KcF6aFWAX`|J|e$&db(x?+=i-?d`bJw*Q>Zo>fwXv-j`$nz*&h zW6J)%hRL_8CN9~3?$}yMrYp1d2O9kLRR~<Y|FOaXWv7Sp_D2YmDj7W4us`puM3A<| zg8hdLZe%&!-@M;sPFuzTsrmaiyYudPS-EZh9Q}sbyJBYVH;GbO9zS`<{)$EF8av{q z?zdie>0T}GuKiruA+?G6-TP-{+@1R2!mj-~fBTo_EH8rdFI?Gr`SkK4dzkvKsy>|k z`rY<0^RA}dVmTW()gESFeYDooGcmL61Eem4-Lpq^#jf>I^X*mIE5Yt9H_Xa2)L392 z0rgMb^Zsz}hx6=V{(Yo$qnYE%EPGgZOrI;eyvk#WJuJM${u`P0+-b9)r+yL~p6|jA zvP3x9+rPOE3Gc^j=~Y`7b?k@5NBSGJyM0Gy?T5wJ(KwUaOd(76!{ReP%k>T4v$gwS z@omy#e|16Q=KZkrkgHU-xXowBepq@*mhtCX<-K=5EIpk*WUa<+a}b{1iY+yMK7V?6 zKWv;AHckl}cY}??!Nx&g<07zrDXd=w>nFha4X}0_tQ`ldM`85<tXzYYbFg#<OLrz* z4r!Afbnl1Jz8_s4=Wxy0-{JSS_i)_I{THme>Qh6P?(e_9eU*ai!u@rtrFe?Zt=Z52 z?1XFS<R$x`%go7rsIqB)gx2fU(`C!{8@FYbS)bgx-{OK3%h^lI_O~>=ThwN>bN{hU z@mCBpm+pUC@p9KKi#_|Ry3gG0zOZP&)w6wUmK*o&7s!@5#cVxyzpZ=beSgaX`|X~U z{+-e_X@8%q+e_1R2lg+B&t%+_RSD;p1$!;A&91bEsee2Bz~{|flk8#UJ^Ujqq+mVQ z9%kRas!;107Z%w!K;5$|!_b|3=2Cl)<U+7}1!L}qUA?r-KFW9@*gq@h9D8xGY?(dG zzb=2WUsX<CVh;-sQ<Y0+m0cIw!@_Ig)XRRy<7V1lUu*^r&tI1>otpZf+y0JB5IA2? zUa~~aM$Qpl&jiYqKNeB#-4Bbe$HF3xQ&i^fhs9?(M?vK5>J|H8@$KTgnD5884f|o~ zp&(C~Y2&M{`<Lx10gszyBwW4SBCvZuEIsk;-Y&Sz-~c?m9j(r9u6=!I|6OxU@Vo+S z+#EK}3meCTjl03dm0;r_uzo(QUkdAY!TJfXb}+2n25XnV>QPv|3M<!O`4*PWVBt6^ z-hAg*i6#4Cw9Lt{e7_f~_uDOYndbUq#r~&#Qt2VgoAxvAxLd<jxMsg;^35q6>f84J z-k*6nGGN{Q%|)+n?{eC;zmJ_Q%P3>r{+uI>HQf69_SeXyoerI`X8(o9wvQ$69@y_@ z_-YRGtrhzxoc+CYL+GLXfAshtZ}3^XpDB@T`t6j%`?FcA-v`#u*dNT(*QWIO@c!cw z&pu8*)384T%0JYyed^0I4fZhgPRr)3tg4@34>Rv)jB)s3pT+ht`%LfkPZz$m!XD-x zmUY4W$urj2Z-Bb@$}6}3HW};et54ns`=_Vj%g4lkb@nj-vIVu43KXufhlR(gc{dMw z{#aoT3nz<qzuX@ZOYEz?zJT-X!EG15=5x)lw_fK8u16Xbst&DK(+zLe@b+z)a`RZi zepr0nIpQz&p=R=aSbQ>l*uiAayLdk=zL|6LRPWqhvp*C%E++8gRKbY}oA=MYI1@Zh zn<n`E;)Dk~_QTRs+**B8g@}Fd^tQre+KEFy4(^A|W5MQeVDk#Fadp@@FKnC>Htq%+ zcY=+Bz{UY!{Zd%J4AxJ8^&4RAHds3jR*%B!5m>nf%h#}U1`Fr8>h%GS=WN&yqgBn> zCi4B=vOnm;f8BR2oAw`wdCw)^vvYr>oc(0y&dvM7CX_!3)!n!M(aF+dd)hYde_<YF zX&ZHL|B*`JWgYRG_D5ag;4>;Zynj;AHpw;28}@sqJ->R>;^_Wu3D2^+zpvP@wPxj2 zrnY1I+bo_*tDTv@-}VIev}<O^_g`I?cbe;F>;C7xd|-aZ<%JSeH(Tvt>SfjImPehL zZx1uibBfl>Mc-H0!|dZ0>X+Qcyup6^gK)5W4lO=@W<vZX`yi-$xBPpoeYkD2{TBA` zVE;Urrk&%_x!E4(-_Y{BT_0OE*~7x2Ao<Dpr*k&g!@}!<663MJ7pv_(LpFlzi3w9~ z_NIp}wLb=Jw=5A#x+GIL+a4Al#^-h|<zQ-s_v>bCT)=ldp?*IsKGQTlpN)>5y&o3e z`W?0b^{<xizb>%~JPzx6F5}}>w~hOAq4NUP<!g;(4{zHKOHT_k{@Cgy?}4W`p_}jD zuM;@9A2x3an^%O*W5MQmVDk#Fad+4_FKiqXHtq%+mx7Ihz{Uw+{Zd#z71mFHwWDF} zHdwm`R*%B!8Cban%lEK!28#!;D|U0eb+_+_(USvKT)U~ZXTQ5{H|P7*9s332J_l;A zIj}!>p(@AE;vM^YCTfYF@;tns=UKm8sOt9pwjcYue;+@(zvU|XBCC)s`|lO(U&ZzH z`2On>-RYqp*Y96iBD7y&`pNw#mN~^-t5~wXJg2X3Va2KaE0c~`E7<hxzvk2p=BuT} zy<2b7V-HjR%KFE%?-fhzVdhO-bK!K+$MyCw`{L4WU2+TAVlVz!9qgXOtGXxSRJYp; zK;6rC=xH`X$qxHzQ2!jAE!^@Yb%#C7pV@w{6a00z+rz?R?^>r2!M|JVVd2HC+2|a` zyvd#m+AeWPnUqy}ZjJqQXgk^I!3Fc0iA(HZ@p1U8*J)+8>GrVrdOUBf0r&nK__$W- z;r0J=o=?~hi*N5IzRJZ)i}&|I$Nd*(&9RWvS-1ZfbRO!qSBh`Mfi3%C>1masl*ii6 zUGVhwd9vJ(VzUGLVe|B`d2ZOeDQw;mHjf3H*MiL}z~&iX<GiqOU)Z=CY#a_Y4gwn& zf%Qva{Zv>#0oHGTwcB9rI9NRjt5;#=8mydyr88K%yVG%V!{)tv_rvJWuFLH$2M+GH z_#WrIoonBI`?rs~d=4MkUm3w1zi|Ja{mQRyB<k!xzJIa)J<H>hcJ7y0_(RjI<kbGB zS*%7ekz4k!H#ol0r0>lBW%1w6pNU?zze+hO*lGLO{kKZxi%mTz!TI}_Ur3nbImsTT zz9n;U$@l10_Av8a&YvHc7rDhAW}mmI>XxgMcG^#Q^#JT1`4ty4*6iP7|Ihy#*uBw9 zb@KLe?X$lJ^-q=I>CT;d_u9k!`)-B!m%D0v>|x=t-DKGV#h#t^u<&}`l3C5AzRlh; z{42N}To?Ip*&~%r_OGD*GBeenV~T;R?P2j@J*iKqbml^PSbT+E|39HXw--K+dbR1( z7WRwP`(g3DO{bnI(r)g4Pv|_;M(#6w5rV7syF%wZo7j?PD?i+{AC{iNQg((%Ufuyu zZ?eKaqOaZBw;#4H1h!5DHct<mw};J}!sbz7^H{KXBG|kFY@Pu&&I=p&g^jzx#^GS& zAh2-}Sicn3PlfdpVEqPIyA9TkgVm$3dKFf#!OA&UI)kOVO-j7mc-0T>htbb2aK3e$ ze00Cf`CJ9A<p=lQc{yd%bCr|(_v@>=7&9E$FD$u|$Fl#-{#%!B{b~EQbH7&F&j%I} z=l9<zs(f==V&neQ-^x<*6&Lm|{^`HW_seuRUo-fl<nb@l?P2Pd7Rquxm)K|zGcR{p z`=2@AcG|=26Fu)b)tcdey-pn)*gaEiKPd_<KWM)Z>R#`*t8WF>583xMrGWFz^uEHv z_5%m)VgBWGTHvR<=72pcJn|O_PKwaoXAcXnmnuIO?s3{>ZvyRiop(Ms%jV=(`y%K# zfWLO~Rf8uR>|yc2x%9r;+^FUDu=uk6Ijdmx-5K_<_>_77R6C0=2|kbHeLrUY&M#B; zzk<$d9?^*lH49$0zX>|e-dX>6<~hBM`(f#cf8x(ZC6Dd!^mguP!sf4qd-ucEA;H!G z!PbSq)``I8>0$Hsuz6G1JSuD+3pTF>n^%C%Gr-1qVdK8AaW~jF9Bdo}HZB6|m%{p~ zuzmupoepca!P;@KdK6Z#!pb#RIR{H;uyptB2CtW6{?Yv~`hBBbmGjb*`}Z7HX!HAY zcz<$I%h5}J&hFnRU{d{X#ew~KW~(k=es*F11Le!7IXk!QXN<J_$IE|dKV#Y^A4}c^ zaQ<?si&0&?3+!R)qgH>4(&^l04>NDkYlZrED-PJh?7PZ3eP;Nl!}edH?zt_Un&6&) z)E?&EFgcG;UXzd7S3%3^aCvb9|HDV@Vg8-^dHR*po`>yW;qkEKc&uI2L3>zuxr8wP z<k8<}e;C@&XZtgATD#Fsdk5&a)6TVj8659#wui;XhnT=4;nUaH!{X~o&3ttS{)P6i z`1FhMob}*y3w++P@R{h-)z;nnYoPP!vjSa}6gd{}|E%(C8E9R{mn};^<fpFR4@*z) zmfui+@qQ~jz5NMwtaxp)Xa85|x;NOmG}t;M*g7HDx)9j94cI(AY#tppZwi~&gw12Y z=9ysg3b1j1*f=k292Pe21{)WHjf23(DPa9lSice0Pk^=KVeK|py9!p1!s<C#xdtnD zVCf8&4!U)MLrb=w*bk$vX54(kp?hw>>&~@WcYBZSKOnbuLd407`z!t|&|OixZ@)m; z+}ZInFYiB_x<v0y>M}Tg0+&TfOxiMgnED=lf1xY2`|M%n&5SoV^t$(`J<J|)kqyzc z+fUf<YSI9w`%{y5fA?5=(%yZq2Dluyns!RywaQ8R({rDI+Xd%~b*uIEAGe43x2`zp z{-5JV?P1|@R!R7|QPE+0Sa?}-3FJxLJz&orvkyG3bZo-OtlJiQ>=#1EIYTol)LvC? zvxmjUgKD#))V>Y&u=wiN<GFv0d8s`tKK<plZOAN|2%jexOF5mUw5N8z1a#fThOU%d zPhA)6Z-K5OlA6l!VZ-Y+`(f$nh~<i%l9#r?(_4mU{+ZXNyZ6J^S;E#y!q&aP)~(%v zu0w*YBZ939fvxMf3Z18i&D+D~O=0t%uz4)lJQQqR0X8oH8|O8Fj{CyK-C*NpuyGLB zI0md=3hP(G`U$XhKCImaYj?rwQCK|)E7xG<5-gp;(n-mV6UNK_pWY9n^XFWNNV|D) z|6=Z=61Le#_cwUd1q(A=-M{^HrSIl58{qsK^YkCGp50&%Q}287@VoHrqxLX!Y?BW* zZ~1@P9%f&j$$>Zl-E;OEpyh1b{*NcG{W)uY6I$;dT2&Eb-G9dZ8MGaGXs1@AdC4hz zn12nDk5+SkJ#G&RkAmWF;kPV~+QY(Y!M>B<H-sLt=Yfu6ruatXT5jBDUkx2sfBW0< zS^VT3_OSSPP`T}>Nz5jDSbQ;RHt|}huCRy2=lM`p1%tQK;OjE>3Es$=b2erF)OYEj zpmjpsu^U<TIn3FA8Jb?s{(rE)j$!qFSbCE5kx2;rx*49{-sjA_m6f+^KWtquY~3zw zoh57?C2ZXrY+V~{9TIGv5o}!uY@G;fo*p(251TiI&5Oe3v0(F5uz3a8ya8;S7dDOy z8+U_^tHH)WVB;LHekrWq3F{}o`T?+Z8?0Rht4Cq=B&=M6m0PfM21`e)Bu{o<aX-Hw zMlWzJGgrKPb^kuk#hQ+id-ucm`<MRXJTAG{9;U8wnXA!z_w)8J^VaPz(#^hk(H>@B zhsWwEQ713jd#W7**XySi7oEKF?1KHZ4;|ok=8j!jk1Zq4+dHSd0QY<R-m5#C_MNeZ z`FEGvcM-+uC+%V3;o-;2#oTtx9u{8f)Z4e;PC0CE4IMXMmbL7&zvTgYS3f22JPupR zqVKi5yX;}{VfZ=nVb{FP_OSS}zp|T|FL0GTEIxTz3<Zq0&a#Kax527O=`}*`_OGGq z(jJ5}GD!!_+;0nAC)t1XXrJW9mHT1o$>H?XiSOe!!_(U$lg9ZMp6%R^zK(o1+PdD$ zXzMK3qOE&FUx$RgE(Cp^9(~>veI5&aUI8}F8w?%ig^jzx#^GS&Ah2-}^nNLNKLNen zhF*`Nmuu+hZ1*~A>96h=_oLHiZmx7w^tfmbQy;`$^<L%jReP9u={oP#L>R8x!|Yq+ z^5RMS%**!Q_V$C@k%t|;<@oq7+3!+02kzJOSz2hxS6r}HdSngmXV+Z&o8q+ntUb)X zGlWjgSXgn&9u^*Mg|>69nH{%>g;#@cOqbH<!}im^zhnT-)9fi;{MvNg0ed~@yy2t+ zhi81au*)76AFpoj^9y#_Vh@Y2uM=&*pUGWq4~x&!54<>97tgVW#rM;4##&?FWc$6) z^q>>`4YV$D{~2g{(bAm4`iX7jepq@+KIHCjNN6)Wz42<9Pd8h%b3bg~0&HKyYUnz0 z*gA9Ax?b42U)VZJ*g8zux;NOmIM_NQ*g7THx)9j95!gIEY#tvrZwi}Nh0SBZ=DA?= z3b1(w*f=k2+!r?P1{;Tijf23(Nnrg_Sicq4Pk{9!VC^<oyAD>5HbCoHSh)r(=V0j! zmj0^Nt_P<(7!6YoOMfu)VCfHLA1wWCI`kRbZ-J#hHE2H?mi{WNw}8iWVCfI$-<i<# z2MZ5a`h$g6BQ*Ur>^=#e7lft1%g}jhSo(v-2Q2-;;tQ7kVDSk{f3Wz5r9V^Xx?5QK z3xuu%howJQdV-}tSbBq{KiEDQSp5szw*aeuVe80Y^)GB)FRcEJgs!uM)xWTHZ?O6o zwhjqa|H9UV!0KPvJUy)bh0UA7>R;GA7Oeh-%`3p_U)VS=tp0_KyTR&T*f<ER{)P2R zVf8Ppp8%_WVeK|p{R^u{Vf893U7?qAuyhAg52In`!RlX_eX#oX(+cKjP(K+~|1N=! z+raAIR_Hh)tp0`h7gqnm!UI<S!omwy|GtLKL&NId(@u5ZbrG=o7Zx9|`WF^ou=*Dk zpRoED7T>V?R}H$(8&?0`gYG+k)xWUx1gn2x=?zx@!uDms`X8`;GO+#!Y~KQ`{{dS^ z4(orw*7d^rAFy?nu>J>Z-5ad`0b7Rz>wm!3g~0kBuz7k|{{uE}IsrN_3!BG+^*><q z3b6hMY@8R?|A39V!TKMtaS&Ml1J*Bv^*><!1X#ZTmTq9}I9NG`Ua!K+IhcA_x`UYq zqha>J`XAYnn&5F6SpOpjI!*}df4qW@Tf_PvF#p2(AF%L%^*><Y1?zv@g|2&m^*?Sv z*WtkWAF%j<^*><o1?zvn;uF^YfW<eg|FI0ZuK?En*bUw10_%Uk(i5!z0ZVVN{s(Ly zAZ+{vwl51d{xT7|PX;#r0^7F$8-IbVBZrN@z}EG`#$RCTEMem>uyt>+@fX-SB-r>1 zY+VR!`~@~o4;z1h&6~o;UtsfCu<;kzyaH_e1vbtL8-Ia~yTQg^VB;XLaS>R$fb~;h z<rJ*n0IN6A+i|dZ6{a3m&cV!sr8}5?F#4)nKX@DvHvZD>=mZ{@hK;}YK<7zd<1aA( z!p2`<;Q<?efrS@r{3Ul)GI*U0Z2YAKy6y`${sM~+*!T-9zF^}ou=s?Hzrf-fHvaO? z39^p`HvZBLO)s$V7g&0NjlaOs8*Kaqw(k@+{tVj(2-_D3+n2Qjx=#zXPX@Mc2DWbj zwvPd}jvThG9Ja0(w$2x}&Jwon61MIw0J;tiwhjrlE(x|S1UCK*o2Q4(+r#EfVe_c+ z(0MG_ycTR;0XEM78|Q_M`@-TKHVy|XhhXC(uzC^JPldJ9(EAOrb{tGStX_qg2P@}b z_QBHKE$Fy2jK2Bm9e5r@pwJsU{(NT9Y4E&U7<Bv@=3mZ5v%ureu<#H(a}qrM3=6L# zUfkgEXIOZ`#-G2Ro(x`R=rsKYc>Ea_A1|Tf&#?GPI<^=*{tSyxedzczEWRf|$Die) z`=%DLdMSa%pO-=R@%@30Kf}_~C+PSyEWN?TpJDr)Vf&t8`%Yo|P+|K3VfzAM`?6s3 zZ?JtbuzfSIeG9OC46t?Nuyy6Ib-l24zOZ$cuyvQPb#JhBaIkequyskWbs?~IBCvUS z*t|V#-V`>E3Y*7*&1=Et6=3rWuyO`A?hC8;VB>JGb{K421lBJ=@2A504KVevb{xz+ zSiK6f4_3~>(itq>c|qshV00XG9@sIY1U&x+^DhT<{tXr$InenxSa@A~@ew@##tU6X zxH)Y%c>YaBz8Ji&6*m6{ix1fR8!Wzbp!09A_&n433OxS?i|=UY{M&QrK0UFM$&mRs z=)O^Nll$QLH&}Xt&A-9Y+coI?8*E=aZ2iFs=ssuI`UBX$Q`kOK*ginmzChT%EZ9CR z*ghH9z8ToQ1=v0Y*gA69y7D~ex?b2iU)VZJ*t$#Dx;NN5IM_NQ*t#Uxx)9hp5!gIE zY~CI=Zwi}7h0SBZ=Cxq;8f=~c*3N>B`@-7&uyHt8KNNjj1lCW5sfYC&VCKQvaWMN} z^(w4fgOzizbOuXzwa|43FdF9HO<fA$^#`!<INUi8y#4?dUJK=t!Rrr}!q%-q*B^*} zYyhvL&RY<947C0L79UQ~^#`!{f~`M*#ixtiYw-F5SbQt)X#=l6u!ZgmT?Ji#kPh9a zTnAl$083A>^#`!@23vmsJI4XG{tC9Q9=84pw$B;1{tCA56t?~fwhs`t{tC7)3%33W zwoeAO{tC8l0k-}MwvHUO{tC9P7q<Qiw$2i^{tCA44YvLYwhjrlE(x|S1h!5DHct<m zw};J}!sbz7?I_s17OY(ln`eOaGhyStuyHroI2>#o1U4=LQxEH>!pwv98({Xq+HtUY z6jrao$~9Oy2TNzLbO-Y<jE02=Z2c81ykP6E?4j$VVe7A6L)Y!Y)?dNm1GfGO7GJRS zSFremt-pfBH*Ect9&{fvZ2c8CbYDDd{S_=d!PZ~F(i?3373|y)*!ol0IS$6qb1q=} z>S60or$hHS!`7d|_MO7kpThP5!q%U{_GQ7=pThRZz}BC__AS8HpTgFW!`7d|*7d^H zpTgEz!q#2F!X35_4z>;nwk`>_E(Eqt1U63(o41Fx+hOyluzn<LUJEv^0h?!ljkBST z`@+WIVCrGxA~5q{{ZyEJuzmxq-3DvN!Rk?1y$UPWVC5Xlzp!)%3lA6#3oqFE)9Up% z!0Ym1>rWR#_ldyPpTgn;w*C|rU$FJ3u=s?nKZV6NZ2hTYbUk?AJ8b>wIp{eIu=S^~ z^aNXf3QKRW^{23NYGCW{VdsXx*5AX<ae%GAhwZD6hMr3R+vf~he_svVcM4m758DR_ zTYsN+39>H>w*DTrPX@OB9=2}*w*DTrjvThG92Tyyb-u8C4qJB#TlWTA2M1e+1Y4H` zTNeUbCj#py!shK^;})=aRM<EgY+egCuK=58fQ|FQ#(iPxVdHQx^I+p5F#BNrR9HU& z)^C8d+hFZDSUn1>S7H8zm2<H0fTcTFc){p4=spnG`g`pHE%3f1*!p`|e8AS<!{Q6J z{vH;eOwjfBu=s|pzuyBrX92eU9+qBU>+fOd3AX+omfm3N?_uX6!S)})&Z(INJ+}sS zZV2og5!g8nu>D7{ef6+?_ON}<uzk<4eW$SfN3eZ>uzi8BeOa)5TCjaGuzfSIaD(k* zfaPP@x^h^(0$b+`TW1MdcL`hf23rRQ>j%QtCBeoeVCzI+<7BXTd)T-wY#tRhkAOa} z1)FDpsfUgG!pwt>!@=x>jf=qC1M8>4`U$Xp1I$0Lb{x#VuzD319<Xu_7GAJ)=U{aQ zyzdA`AFy5x-UsI^GatPF2o@i({YS9)`dlpw-hTv(&pPP-BUpU9L-!xS(!<P)^5FeP zoQ(D0b7p1-9|P|{f~6-O=>8*Edb_);61@Kib`BS8|K(fgxk#{cl3?f5!1iCl&JBT` zBLX|e0d}qfY+pTWpFM1!Gi={8Y~LwtA1Z7gAZ%YCY+u$L8OS~@SU!dAn}O9Euzd`$ zb_{G?Ic!}oY@IKxp9foa30wCDTL%Xl2ZOCkf{n|<)``I831IW~uz4%=c~scE7EC>C zo&jbaY}^-SA8Z^BHVy(C7lFAK)=!1?6JY%Yn15mII9PbV>Qz{H!OA&UI)kOVv(SB) zFd7yg(a`;uu=v^!-G2#-&#q=a@cv6!eDB(63EqF%5CS<z<|$JKc>g6Vy}<Tg!qU^{ zd{^-POIUh??Z1Sb`v%*84m*bnw*MS<E)s12IqaMo`2KU~xgoIq=dg1eVEfNu`|4r) z&tdzVVf&t8`%Yo|P+|EHwoemQFTnQA!P+UXeKW9jGi)CNY#ljlT{&!A0k+N;w$2i^ z?h-cc3SS2YU55l)mjs)af~^yQ&6A?f+r#EjVd`P?TCjNq_&fu2ogHl47dGw&8;66r z7d9>e>zBg%sWAV-`VFw~fVJab;RUN#Vc`iY=V0j!p6;Og&SCKai!a#zb69-B_n$-K z8@B)aw(AP;xka%3=UXl9!RKhf_MgMj6KwxEEWN?@pTo`>gq^<uJNFHC{s!zEE(_>6 zYOr&WVCQeZ&Z&W&zX3Zp1a|%g>>LN!xehRY!S>n1_Bq4$J;U~$!uFxU+99xgfv|Ql zY+oI$p9b4E0~<Gh?PGwgBZsXkhpp>{t@G`NuCs)#yM)a<!PddS=0RcWl3?@l=<7sa z^Y$?Huz6INd9Zmcn0>H$2G}?+Y}^+%?gkr&gZT$GE&}r}te*-C4_Lne7GAJ+9IPIN z)vK^_4OY&<;sciMVDSZ`Vett&e*+fZu=6)y=>c~B#&qbpZLsq<VCe~V{st_)!Oq`+ zoofj@e+YKYAng1h*tu`8^M_#PaKX+Wf}M*5JAVjvP7Unb8koCb=ZL`a5$s$C*uHw$ zK6}_cXV|`HSi2Xt4;9uAgY651jSIl`-NE+B!1h7G_AS8nF~HW5!`79<=9OUUd|~sP zuyvQb&~<OHb#SnCAn5CoVCzI+>YqZ-AA*?&n@5G&2b<S|%`3p>8DQhQuyJ45xEpL7 z4(4CjxCksfVEt5Bc)|J&uyz})9S5sNVf893K49e>EWTjr4i=v<8W!KM^M_QNrhw07 zgPlLL3VKc;?EE2EdV-xl1WRwQ^M_#PsKU-)gPm&$JAVy!&LHgkHQ2dtu=Ceo=WxN! z<$}2jc1{v3pKwCYUxSs~uyaIU=QzO5b%5=whwZb6^|N66o?-od*gjO)I23GOAZ%Y2 zY@ZfvpA2lD6l~uDY~KoO9uu~Y05-1<Th|R+=K))H30rrBz77tyE(xX{woU|Q9&Fwo zW*=-G6*i9ro7aNPE5PO%VB@^7abK8!VdHSH@PLhrz`_gGPlfdpVEu+&&~_WF9S4sO zXuS%HFIYJTi%(d(gT*(DE`pu|2s?l69rRpE*!gR)^aMM94VK>E=dVG}?S-A=3p+;@ zcK#&nTuazFm#}jNVdoCQ?1h~J2g?VrbGcyUGVGirSUUlBZVl|*5ZE~)uyY(>=Q_ap z^{{>RuyH2XzGv9DCu|=oY#$(OUm$E>7HppuY@ZBl-wbTt9JY@EwvGd~t{k?m2)51_ zw$2TG-6d=t985iIT@uVZ*g6rIeXx0Z*t{uh9u+o^1)JA`%`3p>8DRc}jr+pF12zr^ z3oqEX2&`WU>!-r{39x<xEIwfEI9Pna>Qz{L!pb>Ve8bY6)&*(sxsx!uy%TZ{D(w78 zSbCcE(FuJ1B<!5%dcM>c(D{?FbE0AAM#IkSg`MLIJ4Y3Et}4u2*g2Q5bj|@icMw+Y z!p?z%wF6-1a>3fAuyc}N{Uq4AHL!C-VCRUy#*tv>I>5#?Vf*Z1<Lt0~&#--`uzjeo zeSolifv|m9uzgyvbsDgJGq80-uzd`$b!@P8<*;>?=<9r8>n>sHVe8;v=E2q_!R&*r z6M@as!{+T_^QKMEc~sau7HnP%=3m%611vmX<G!%)f{nw$#zA1?BCzm=^;2Q-0qZxw z;tSS}gT*JTUWLUstek_T2Uxm;r56|tOHXM>Am{JG(i`mjUD&zgu=B@Z=S0KKjfS1u z3p>Xbc8)6STvb@P3OnZ#R?oxE9fY+zO`+$H!}>w6bGcyWBEil{f{hcw&aHusTf)u} zfsLcX&UJvz>%jKe!{!-b`<`L*-mrbBuzi5AeSxrbL9l&Vuyty%eKW9ildydZuywrX z>&jv4d|~Qg>n>sD!PddS?1QaKf~^aItrLN{7dCGX^Un(Cygbamuz4+5c);cvVBrND z_l1QgY#a_Y4gwn&fyD=`p9+gFSib=lpRjfuEWTm&Dy&?Cm2<Fk21|Fa^aP_}=?!-N zIP4sJ*tz&HacAiL3$SzwJEt6WZZGWIXjnZCJI5E+u7sVd3hU><&bfs3yJ6=J!p4DM z=fJ_vMS`8n1sf-aos$Hcw}G8o1Di*Log)I9*M^<z0Ndvb+h-43_XFGa3|of=+lLBU z7YW-J2wSHM+ouIvHxAo31LMQimBZA-*7?HBgRQ%S*#}z(2XhZ>T@q|v2yC4Q%s;Ss zdzgP=^Qf@!fX!>c!V5Of02}9pjr+pF8#WFHix1ei2rRx}{Zv?d!uk!c_=dIPVD%`h hUWKI>SUCqvPq1_cOK&i`P_3j;J*hOYq);PC4*<k)=wJW< literal 0 HcmV?d00001 diff --git a/data/sdf_circle01.mat b/data/sdf_circle01.mat new file mode 100644 index 0000000000000000000000000000000000000000..35e6cd750be67e1a5df63f863b62013e190d9e72 GIT binary patch literal 27305 zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2cQgHY2i*PhE(NS<NN=+<DO;O0t zvr=#?%2e=8EK)EuP%t#NGBU9;vQ#iKFfvvk5-`93qo*%F0|UcVAqIwuIVCd<`=vug z+<1iqMJ}{-b@g#|GOZHfj{V0Fx*|Yo)ms7A6+ufEa|JC42;8f8OX-EE;HAkLYEDvy zhq&CtjvY(jIi|hhQeys(eHY7xWY(*t=bpdw``e3q<;U0is7>C!&R=ah-*Itu|GwkL zx&MFs+&;gYe}B!p5AP=*{vU6(ucG$T{{HLxYTkd?Q#>i=%=);U<v&mDHQp0<c7E)y z7a#sez0ZGk-e3N{UF{!UyZbZuSCqVdw>bL#qnf=xcVCaUvn%{}@Xx=akE>(yv;IA- z__e#g;P1?z&FrtU`LEe4&D;5N`TqMhc7pd0)$~03G<Ck<`%8b$A3Z+3qdv9n`Kz<r z*WdfY`g6@$`T8Hn9{$w+B>i)>H2?ei6+fgu-JQ|@eqZe`>reM*=^xtphy7ae6aW6@ z`~Rz7OHbcd@q6q34RO_<Cm;WLa<lF3;!kgq4Nvc{**oiL{J#A&pT^hiGd{iFCO&^! zedWJ{I{!aa1#kXy`t|uUpY$Kkoc%QZ*mi@{`+LeWr`1dTKc(~EX#Up7f2noKr}v*J zPo7pk^ZzxS|3$WL)9QEHS#SPx`s=O3=0C$9#5=}3U77rA`Tg+XU*Qkp8|#17r2GBn z`lJ0}{ek^VV21zo`7<v=*>FL0M)1FAr~$11M1O=oh~HTMX7QTrFFt1${}DdDzb4-J z^!`6J>C@^f_ubR^|LIS2<iDq1m!<ujetrJTr}4*(%}?*|$xobCFZuqG&VQqSQzQSS z+GTJ4bNY|lwECIvkLmn3`oA{v-_xodo&QC(WobXBzuwxO|NNWU>-#s{Bd@=R+-E1Z z`QOXGXXkuff6Kn`+WxLz;er2G{S{9CyLv^vD6*i|t9VyDg8o1A&qn@7Qnct_=<nl} zjla#l1zt0)QhOHouM-v({#1*K?Zwp-|CY}Y-jSaypBlS3{pY`rGyZPBX<zv5e%k*d zH|uBr6;A)3`(i&YRIu&e?F-lirGMo+W2l4)p8faHj^laIInF!VH}-E@chT6e|C7JY z|4%i0BmX`9$qq_RK4(AaueUcoy}u_tZ(6<N{X;tcjqWXt{FnMid-I>uHP_O9PXF1f z^WW%xZ{)w!f6<_nBoX=V>903uKIvc2=D&aLZ~6JZ_bcR8b}sOn`|#_p5Qu+%ew^`p z|G9si#s9D0v@d-3-|Dyexxex0|Bvka|9S%$y_s+Ey<YYk2+jQ~{QTeg!(g<z{?VQP zr8huG^Z!~)m=cIW*>Cl8|30pbyI20WR^Z%HMEE7iJ)Iu`^O7&zOUdMV>Fm$sJ?Ce~ zY~Z&({QB{nQ#%`fTYL+>W?6MS`JDXTzrxS|T)b^xh=rN_?}cIgmv8=aKD+n-Tznj= z%=THo+-LXv|3@tAzs#_&Is-K#`QMAT^81fJ`p&iI{OR9$d-xCSoKnp6Jr_A6?YVbF z=l`c)#+(0~u6UOAbNWYJQ1*|R^)&w2exuX-&!p#0tDkxQ#3%jFcV|3}f41NB^!}Ls ziBIF#oj<4Zzi98b4{z_(Z=X|t{Qbrsdzs#U{<Kf}eRF}{(>VufrGIO@ma;wlx&Qc_ zdiLMy=YE!NsDFf$d1gQJv-=O9#4|83{QsXabASJx5D~ZSixzW<ipq-0hD5J<@rgg6 zwW}}7dR4$uFZZb)%B@|eVy7sba5b5<R8yumWi4w|Q&{9N3rV)NjSIYvzvFpc{cB%j zZ;I^CpFi)t-dtz@=f|JBi$9m=pR=>5t$O|?v;N!TzthjZ-(REs@9D$Ir`!FnJ-U6Z z|J;9p`X7g%oA2+fsNeB<`H$1$kN1ivKg`eHUH@<TkKd<nOXnQ@{p|ZcvHky^Po90> zuV&N#+2{7%`2Y8d{roux+AGr1&rSVt;rwQc`0oAH?_bZpZ*$B1(ed9^x{rQ^)tpF| z{9F2{yZ_$)zpnrHTJ*R7{;OkW(|>x$f9;P8@9lT934i2l8}D9o{n6jY)_-*$#diL` z8^5pX2m8nLa`~F7?+4#cJg)KY+3NoMx_qVmhkiW&x$fwX=Rfn0{CNIVqo+P@{#uFu zd*XI~$h}m2&|c#I9BspYKYvE`)Zg1D_;`N&{jQJaf4@HR<N4q3?(=^BJtg@6&%gCY zemwse)A{lI<K1G9=O5oK^mzVp`;L$2KT96{@w`&>nEm;3-DCFW>y?k$uXj)U^Z8fK z(I3yh-fGRCJdgR`!5_~-TF+lk`t$jlug!PY7-TEnZ~ek^P8p;Krv*Pj78K*K;8kC} z!9UCJq(7fOZ<YAJr&jp=OjHXVXY9Lv$#{Rg<)it>i&c-=_rH(osh4~2)l)C`U*VYj z`E-%T^Ph)H{I|KkRpP(R`(BCvd*oU_o`1br^zr=b+2Z?rR~i4Syqo?}MhhNHorSrA z*S9DA`TWny@L%QszZduD9=}+#erbAv+(&*wHv9+Kz=vvs_P?M1f9m{ukBEuR$}39i z@{j&_{<jol&+o;hPUmhH$(;QjRrT!S6y<-FHK5?E0Qvg|yWzh|yTua!ZLq}kbXZ&m z$$mV4d}?Hm@l=Wbe`<7(+2`LsvZq$?`u2}8(=W!WAG3eI-#%*kGb_dCxt4v;s&#Ff zKX=p{{HxrfZTPS9PqV~-oB4q~_4DEr|9l2p{l`4%&*v(UWA^LMhxOFQEsy<SwleR> z^UpIP&txNoROut7b@`8NH}$}R%Fbh*ck$c_pSQ`lf8IN7US-YG4&1Q;4lA_S5dQab zPo>!B7hXN}b^H2uonKrl{9CAd$G+aeY)M#P)o#DTezG)kop8~;X^{B0=>NfvoM^#m z3Rj|i)l|RFZ}!n$YlR-qpTGX`{~Ecsvv06RANiVlX7V&xVBLCt7So#k?J^(FKbAl8 z<N42?q(7hka6O*?9F%3A?-zVL|N3Hsf0eb{K;icG(00$#1NQgDXB@S2R6S;IAAfv@ zQU5pJH@o#7{W5#@^z@7M-OtMe|9*GBFIqI;{QKP^d6nx+clw>EjP6_fGtT?`1Di+k zJ@xa>YaFvbzuw?q<-Tmgf0h5F4gXczay_1ZJze{l{rc;AHM6If|NCioMd#a6aGY2@ zGyI%!|A^gxrFHqqvO9YI?W*vParyuEV)>ds>n8l(v-AD5-wHeR&LSta>B{ib_3v-_ zg!&Ej|7D-Z_teMrA42fw|2EU~|MR)_eB_(vu)yi*h5Py{Kik|g{}`JoczyfE_Zxct z?UI<Tx9I=B7rQg-j-PtZwe$S!-+4RvC%q2(VGk?K^JPLMT!WW_OY^DLD>QU={b3GS zxXAlypH|4qOFOwPc`aNtOK%a!EF~wY%<fK|-c7fz8eKben1?MnqD*8}$&I*kvUhi1 z^IH74^7Fm#`)q&LoqhTD`10@a<;UmS*8F;Tc=>+W{WTwco_!wv@pixcj-Ov1)%WY& zla2NH^Dy`0^k$3o@%R45cJ2T9Yx3;#?;pMQyJv4-Y5&pR|K6UePtW*2>YtBWGk@<6 zq5IS7>OL+1u~>fJj-N5z?+;aenEqq+_38a_@qT~a=6+rOPpp3bzT5=uecJzC$bZbu zzc2gt$oF^u+`hiP_rLyw-MpI{ixuxZe|>yMe3$?Gef{$n*Qnlmo-KF&`lH(VbBp%5 z?fGwa-|a8=$A7zezJ)#d`+Z~Ix9E>{_n!-%ciVIRx$BYZppy9Q(I3x0^0UvYtY4=6 zukydzG5h!Lr%L?a^B+_&)cEw&&jXe7Hu1YZ<c8RP;6G;n+^MI&ZokOm`O`tgIk-6f zC1v=p@^AN#xM>&fdG*xatJ6Ma-=D60%)WoU@Z<T%Z+C+dGv6^#5v>g>qW|s{dOZL8 z|EQk&dHtOq&wmD$k<YWm_w{y3|4;Z+tas$c^Y7k<|0;ia_0<0Z6*BABgVO!}j=}`V zw`-65c>WetYF3yf{rUU@RCMmyEAijv-@^pCvxU00pa1!W|M~ncuJhyh-}Ru<{x~Qb z-R$~!{_|Gu`))@&K?S{d(x1=2mrDHKQzKeA^EkN3-{;m-ANODB+1%7E+l3y_-wvu> z9^dTvc>Xb{7W?=OWLbN2Y4o|pmg|<L7sP>z==c9Yg}uD+<N41)S>^n7y&7pzagY~d zIzFC%KRxNs=TD)3_?~}6_C)D9qg&CPAJ5;;2L*i&s7?eG&K1?567ujtJH_j{mi^Cm zOq*9(zZ?`L>p`{d42l0X_R)s_D)*H=yuHE#Eh3Wc%t!IWwX3Hy@<l7--G8r@_-`|> zSK`0TKEIxNxqSU&_Wku1H>N*B*mAJ)sJ+C08&pqR)kE{d8slduo~XTJbSwJk`x;s7 zQSm@V<@?f|{wGSI!FBKd%P#cD&SL++=WgVYKtxEb$wm#S;^x#X+q)`LaanRl_q=0K z{H)(CG25To%z7ME#r9nEf5M;7e?aYqpSN1`7r#1c|Gh@+9#VAwUQuZi-~I9Y>ugXv zM*q>AEXlWPlVmsdeFKHA4XEg=nJe+%=3X4AW-ll>H~V>BCEujqm+q)P*`@pZqRjP6 z#^;Zh9JN0Ws<7s*kNsh`O773+!cwuvIf}>Z_s4gAJbyZTV?JtxoO7Li@t)STv#&Gu z_wJeJj$zG*IK}gopK<x8=Ix>Ff#)C2KmU01jO*Nn|9<`r>Zy<0*Hyd-)h{x2AivDt z*7xo9583%xY?0Hxul>X_f8Oq=mLBiCi|0>xy={(rag6tQ8~e+m{~?*&^3Uhe$2!+0 zwo3f}Q>S^%UO)cm&L=kW)+@i;yaJjHzuSS5*K;P2Zz><#Ji^E&-<RE+|8j5GzVs8t z_ifrw7f-yi{MkOm>xGu<pIJ=5_f+)%o$Z1DkK5m0ekSRD5F-BPG@pO`)2dEAJuqJR zm_4*WPo1>apUF|+=-L-+`npb?>+)NrtylJkIV|F0WNdmXsKKtH78Q9>kc)%W&~W2~ z2?8A~1}*yjf=`sAE^gfZde7OO*`8+m@7#H}XS@EZX|mhzd_87g7#jWFX?=LX#rxua zYWZzW-2EEbn9TotRojEpZWFRt@*n)(U~SGcyO!si`??U{0}oC5tXa$b`a;hh*#E2i zbfE2jnK>S5=O^7PU#QLyzuZ5xmZNrB_RKWX%YQ_w#2;+f{h*%nulEMI<B^>9I`><S zhwpE-&)$&dpW6~&W`8MLy)0_?{>nVtbE+FYu2$N2{opfOIoJQfH-3E8`B%><ZL3=t zC!KR|v+}w3O`7$8_Q}oobo{-1lV*L*Ji!^Cjz6E>pjlsYPZGqO->g|5<Ik|{Phsse z&H9-4Og=v!|C*j~O8)y{qo@79H?sx*+Y$f2xSU~?_`DF%;MP;{;MSX&*(K+~#g5(I zVY+NpRAla|X=QHFT@RP)&E&I<nqJwhxTa-=g7<`sl38urUI|?c;LH@h%DGm-rNgUI zS-i0GGne9bHOGMFmW<?ChLYDqN_kHIo4e-g|4Zu(Kd-&E$xr6_o4+>a>*sI&ZkL(6 z`R1FwIorg_-OF$P{djNwvrqZkch|FT+hM-zrrnO;#iy(PSC;R6`{v5@N!w&X=N-Jf z{_Ot^PGA53T~{9NqaB{Ux_O#guU_yxQ}<&fzn6qm`p>xj@YnI9SqqX+Tw6AqYwG7a zAI?XIJkdVBe5>ZV<$W6U>o3dB);|%fSs(pB;+H>rXu*SzyWVCOhaP(RdHKGXB@RD# z6s6f1oml!z#x^&9Z`J-!pR$&E=e-H;dn{`Cck=gDzZP_QcTfJlqT!Lm{{_*zr`|u^ z{LTN?6$_*DR^JlUo!&fR-~WI?+~OcFf5k!JdoLRFUZh-oYdn9(V(*&MQ)_+Fj)#6T zntJ!`?#;igrmx!cQCh4tKWa_i-pRFZI<;9>o%;D=?(=-PK;!mxC*D3Oj5s^*e&xj5 zADr2i*S~r$%f7{Kfzx;U%PLt7VF!0#p7s9j4*#u@kH0&{zqft0GgEN?&Tm!QK6%+$ zs;8z-J5|5S{Y_f+GB2OF%j(;c+djRh%8Q-9{qof68_q9J+eo=Ti#)#XUlRN79|!d< z4zQjJfI0a>Y>A_K_MU~a{cA3ir|NozHhk`0v%R?M%9XW&pOfvi0!kU*&RfZ2vnVJ{ z%t}sf`tq!@r{BNKDGEq!xgFZQls%-Z?)*XCjZu|-{a?7Gcluj?-*av0-xGJ+?yEhV z`!VI^i(hKiQ$O!+*}D0+#nipgOWVtTglp~AveztqVJ15D`)$p?T-&t1d(W%es4ES2 zB+tglx@8x|k;66EmA%AK{k5fQdP#T3nO7qAOAj2W^f_^2%ad8V3%{M%aozga(^W^C z<xdA&>ZPuh+*5RBmrU`#1HOB<NSXY4u)=?CR_WxE;?K>!Lszf%U-rX2H}Ky3>ASbz zoxCJV{cnWOe#5|>X?Nm`b5fg%V(shRy_h|J^@WqJzD)gVjg&vHNtQ3$Agx~aft#<O zar47NwRb)=h}%5u<$KZOckOp?g=_h@<!xsC^EUjSz$Lr<x9#`1bvxsK%$p<azWm+0 zliWsEUOn3zbS}PZmc_?EPtI@W-!p5ClefIw)d@S}EEhWqrG4!WU$>Nfy{msm;JW?c zx;it&*D0NQQ>*mlW9F%6GPZXX|DL}k_m+Ly^IonTnbdiuZ&n}byTL19|AUo%hZ~1q z!C~gO4=vw0E-E^JA~Np59pU{o>~H4Ic`^A@j&=CsT6yb>`nebMu6n&&u2p|O>i+cY ztKS@M*Ye`eH>=wEbJwos-M_sKUw`p2=tfECewW+gH?F_kvB_!sP2u&Q+1FlToL>Fr zPT?!=-2HQZ6<(djVhl?8p<j39LsL?yw#Clyhm86ahphcdR*1%4I62FIr~OX$j`V%p z-)0}>*VX@J=lyd<eIDPuEtltrpA&q3)@+&D+y2vbA?uzdr&a6~UX|<UJ!{)r$>qu} z?dkE7kF5V6=iPNCtj@~TI8xm#rS$th71#NDjLfWx-c7eU8Yr^OEN8al%yik`8-%sW zJ|_R+;wxx0_JRb1xmAu(?1ht2?GxH}`0o_&$loXYZ%&cko_Qso&%D@pu_pgq_<>d2 zpC>gf>Ww=%HEaG}BTsAPV{xUfc~|d+_t&TJ-&t49ul?6>Mf$9~%i3Fie7$SE^|}14 ziuGrFugAsuwbVX1VE{@I^PS%$n1}p%aQ%Tm%!4i5c@+n^=Y8lClD)t>cg<JnQjW^b zJqK(LTB_c=uy&32>BNbO+GqA2-I-PK>rBHeck>w)FKg^J@3Ts1&r|(VD>wiC`!hFu z-c}@D-q!kXcS)|thwqEyZ>&vzwM%B7ZAyRK9pgWT&ItE{?3FLOb97zop9kVH515vF z7aTtPMD`(Peu;y4$);<2ued%b|0Gh`_`CPVbL*z%S}BRYW~`ZNb$kEz*{klZ5Po{6 zS3>a98ZZCLb%%-za+Y7@Uw`h}JNulx>f&Q@vAb<^Ur)DkkI%MUb$zD9+I(C2pVxOb z>iyl_Dn2P{sd?X-&3(0JkDPPZ*e(C3i8&^<>GHPEb`XENu`hp-^eV=CzeBC(ok{O~ z?Ns$o*dD5Uu)FE}>D-O;kDmz<PX8^lxb8)v>FFbZf1|I)TK{;Pe?|S-xjTE>n19`V zBs}rq_2Sa$f{VBL+^GmSE_SW-MeXYaeUm3vUfalTU792gPU~};gdsW4HURA5Z;@AC zUYUKx_euRHiSLHryDXP~ntSl{&r|kSzg}sZ|LxVUmEqn(HZi{}AKnf9qVe_B+u1(@ z=FQul_vf3mEZ^&^|7X>1zPt0VRoTx=8@{cVGlT^11!1sV2Qstk<;`t;<|)Q=zyG$# zdDV-gUs~SvUi+Qz1l9++&$M#+B=z0+JJ0jh_k4BdV((jCG+(Cv>{aFR<=<^4R(agF z*}k&0;76?Nw3v00SB&D>{rG;Y$w+9d(^|gkPfz;7ZBy^hnYwFB!0wYLY+m1JpSv^Z zx;Hpdi#X1SfD7i!AH-x|m|r}3>-uH$OX(N&KiK;GRr|{LbIP6rdy~(tJi@qZs>dDy z{}%@TE(WhFm*)EVfOYC!|6c8DZ~qnTo&ElVknEe;KVROPbPN=R`CT`DM!XT44k|HD zm_L!Mn6b;f;*fP<$pU>B_1D!LKimFIsq3g~{59ph>)r!974!@C7y5OYmv~FnU0hrj z-MH;+m9p}ZRVQ8;?|B!r^!StYUxQZ)Usu^5wCp-B@2<TXy9!h0mA*NAj%_2fn4Flt zqWJ;1Oc0cH5ii-aPg!vOm#<Hhf9d>UeHHr3JlNht?up&=xJNaD`5)#QYyOK6e;>5i z`(roX|M_tm*EPe}$F0)RzJA|C|IO|d*{$H*eox%G{Ls7GbG7O|WWRGi=(_KlI3zi| zRkvEwrS~GGw*IkgM>+51>)tk#*RS~P``7EQ+E<qE7SB7XyFR<!IZ$6*cKlFW@#(GU zlgpyF>#n{adQvx7ug`wf?LB3`UEe<f1=Y;jr1V`kYMsH6!StL192dg(zAkc3zsNdw z&HHI{@;~vtmsyrKec$97$3Ma6gQf?@Pf+i)ey4j*zOY5kZT&0PT<=qTFDufoEqR(R zXAVkL@e^-Is+F~Ofiuq6rp<l@hYx=cn)g9ZSoTHmg_FPTdG9@BE&bm0%JDDnzbyO` z{Dt+2^(&KCyr0y*+ZP_MH2J9Z>x{OQ*x#^ruf28G*G>DEr+x;UZe$;Au&q0yZ*kzN zf`vUabHu7yxpaU+L*4xS-cNCv_ZXkpJ(X>#o40Suo@0NV?t1TA!1wUZB>o-OJJcWX z?A7~KQ<@q1!qeC&|54z#M0LkE=M>g=#0R@`=@qzfg7ddp{37P}I~Ps<`kni~_TNW3 zzjoR#sJs7f`mddlU%G$Wzx@9CXYGd1zq9{UZ20pr`|qcOuiT{ywfSvJ>z;k*-ZMEk z-Ad<c@#gs}=P!?VmHsO3hxp$^zh;)csrG;KC;HX%1Mv~|m-DT3z8c%C*7<5&>ss&q zS0?f0{Y|=mXP>^bzcuQY`90m2-4*VSsxM!^JRg)71kC?4f=23NBdP<gTo;Y|Ah^cG zH)T?9Vh699RYrlK)OO8hGZvJ%YV|W7<v7LhNn_hhfy+{kM;0&&DQSoqE^<sL?Uat0 z_h}#Z32C!#{wYx`Q$j8Uc-J~kF?rkbeD|3T*Jbsd-}_v9|NWkKRXqPM<Zgd!z0ibZ z?uEU{f|2hYKK?yfJ?u^OhS(YPFDI@#=lAjci|^ixe%*Hc_Urvyw{u&&uILLbetSPx ztU9k{XN+--$oge}HZ3n)<Hz&AFkz#!!`cGh*cZmen`6#IroD&=4^-*-DE)r6(5KdK z?=CI+Y#Yycw6Fec))y|pWlvVD>5nnisD7}^&fula)IAyBBsc0UVzqu$|1}}j@@V|N zg|ntcde<$Rbt3()kI8eF1M2ye;X&Kp9qFn%F66cO-nZH}A<V3^)2n)?*7+{0|Mot+ zw^P~m$3e~F*h?+jX0HmpyVI$3bxDfCk3Ad}c4evWj|CMTP<{99T=vu?l`nPa)8=}v zC|9@VI=xOj?Q7BhOX*UsQlGi+iq&kXEnQtvU)KNf@kX<U`?lZSnI7D};s52ss{f>| zTvc8-i~inlc&ETV_Fa)T`2R4m7yLBzGkwg)|AEQ>!?_gMO%dWRSnh6mckk8D+t;kR zR%px#n;Pu9ZrAJGUtgxCync51LFo~XQ_lP4*H?$T7w?lU*s#%4{j&JJiLdy#D$Uy; z{w-qW=I@GoKYZz#<$ttwQPktFmoNQJU9>d4@chg9tB*XsF7$7g`Vp0TNp^-Z;)`o5 zH!D4w)$;SfUB~T`x<T#Jv)Fc3P0e+>X!>iS;zr+yk6pDl9?HCP+~a8egO%^WB?CWO zh*M9_E)CeMcfs?OSW)SY2Pw9h+3_3wqEt^r{G2!WZyVR9=sV?EN7;{t+`B(1yTIN( z>Yp6f{ey?@*ZtFEmESEpS$dw%vAB1!rWU`qoVd2deSVLfdXA;O$@=GWZ;O{ZK601Q zc~BGe=4Gjw#&XAFsb^N6W$l-l<?^h)toPjO>`5P6Kw)Tkr`e`e9^@=0zt+FOcOEdt zJ>b0eZB27yX`|l--~Zp&&)5F9vyqASxNP{@RmxAFv6Oyk>$EsCXPZy`>va=P&DM)E zGM}sXc3au&3EIas_HBQ~`u4N6$JX%DSAX^%USw0SQ<V29$JVQj47=yDb&B7QeP()n zVfgpE=eFvadEWc~V#hRr?;GMQHp;EcJ<NAn-n_#kZX>_-oWr&&@;@-|KfoB{c8nA1 zsMV3|FPP>o$$KT;cln!{m+k53m}@;txwzhWaPrj0IL_mqn_E)touP4|{>1;cOrQRF zPd|I-RLsotT#L_FSwFwe>3wBi()~>(o&5JDe>UyfJnQxMEA!8tk?mb0<jmUo^L~ev zq0#B=$owvW<T>R*@la>H;D$Kk0L&E^7<W0SzgT8zw8VVx6`xx-Tz{{<7QbYYg{ZIp zn#CJ$XJ4xIzj~rBS@^Hd$LD27ryThj@^#VFO(~gtGv}mLsJO_gxn2ym`LXT<{}stE zl3(^ZJbQWar>TtS$HnJgi#~d>A^7r-TV4~76hb2Bi|(E7eGK~#C>Czab34W>_keZY zgS8ShtC&j~{VtrnvHjnVFW$OaF8?j~>$F=<YR`p7i|*WWV83Q_-NRH(a=Y?MN5<AX z@e}?3YSogm1toXdt+;8wPx%vLDeG6urr$EpSl|B+w-0U=>3AylO#Pu=>!kgk%&*;h z`SqyN+Y;Sh!oHm`pt!jF!M2Dm(f)&Wp@H6mmOMxb@mvFTL+;kJ+<#N{n7uMT_wVV} zi)okM=UhH>=ga?^y0Xhw$B$LrTy7e!e|5*JRT)`r)#1PUQ!mLslR0AlisO~xFR5qD zA7=eO7xn*!`tnxSi=~TqgsA`i@wHldvhcSW<u7x}BS5wvwk@hVEW@<_1$Pg(g<w9! z_6v&rpyYYB+PmlQVSVctr`n%}$L+Z8b@<@txR-pH%ciX0ZQXUSd&ymuFufnewz2%5 zE}x%ozEFAJWZmP{_SboDrtj!Kxs&~~^7-Xm=8F$!JXvsNwdGXFMMvMx;7F4VdEBP| z`+V9-Q2yBYFV+A3#A5y5jqx}ACB<J!e-g1dz-@ETy0=s!DRzNyNid(;Z#~NtwbRco z|DL<NBp{<MDCbrBTBH5Ew=xBzv(!tDlsnce_RH9D?^NB{%0^2rn;x6%vxTJPRPsG1 z?o!$xFn`Iph)Ve~i|hT>KjO~p{lG1ApmiSuG|kUj<81JnY3_pCXB$5pQIY4Ks=u>% zbJfI|PQCF{YXdI(edJ8#Hr}@V)aM_39qCtCuY`Zm{8D|!cZU6D%lHp+Tn}udZvI+g ztgJp;XX0aZn@gKdeEN`e|Eayq@`^K_dy?eKEZVCJ^zU(lEM=W@0FuCJ8%jCSKwc}Z zeSGM6dwI#k&G{=I8%}NgV|3_!#Io5-ZA;w$EaF=*Z;`&_&i#+8lkOe4(^tN5;{Azj z){*MvYYi?fQoUPtukijUE4lkWUcc3j{_Od&*+`-E*paT<8z+lv`0p*RNI$4wz##ri z-~pt}ZsWVaYUOC2>T;v`hQI#B!2MHBv;W=garA%txsLEBwy)HGUC?{?yh!iS&!BTJ z|J*ycTd<P3PFG*w*LJbi)L)P8HvfWTm@_|jh}S(3w`lO*w_U=#qCvc($-JZ??*iv8 z&AhgKGpoOC`Xjwd$!_BF71Li>zqr3*@(gvu>CEc)I#=lL{#nkw^;TGk-Ku#n+#fo2 z%I@--8PAx%`)c0V+0c|Yxt@RKAG^7K;%3<F_kFi~^6T4ADw^~DGjBhxzy})a3U>jm z<g$rmx48E9ySbN}W|j1*=r~vB*Hhz+tb#6GQ1J3rF%s6*4(X87R<2G>Ti>HDV`Jv2 zlq$mM+8ngy|KIuZKbOR=yS?`JuJ^@u^Ur_3^L$S9jPCC2!<l+|`|OqNum3;x|G#-{ z*16Bm-M^lf-}(Pf*{!~A?+5??e?I>!V|m%_XUE(YZ1h<4U~f)1b53<yfbio0{%fnx zRh*Mt(SBok^zHCfSDaFx&%J-T?RCP&=`8;L<F5u3cHKXfute<Dl-3lJ%E$|;$7;Wy zo-OsDUq7TqFMu(-=DEC;T?p^)rV`&DrKc1=e4FFBwo!}2PWt}-4UL?OvT|oRGwT?{ z7zgJx$0};*t(PyCTG>)^`WahdbIk+$zj4cpHZfH_2{*kKdz$fm)$|2tEJZhX?SAoi zRu|8YSp5mr|MnMcsMk_nVL6S7TQ~ODZR^txevUti-mhgb*dv?x@8&WQFZJ!`rgZkS zc1+>6b{Fw_;JoU|min|-=IM`~ubJrk{o&`o*@i)JjRnq$?vsB;s_BRe`o4-g(-7x> z`opw?NfmOr0frY|Gu5rUx&8m{V4izS>dh%XdM!W4?mqp)=|}IxYyBH0X^DN`sqpst z&%HGjNrucX_xSc0&F+nB^!Qft<DYr@-!%nKqGzAl7XIF8&zhVRok{kRN!L%c<*6$C zT&;OCqqsnam9u}<XK}sLp(zRvMekgcyjG&i&gnn1viq?s@280$-f7Fsxsegxam2mf ztz>HjvvBA7@~^^D#-*nw{5WHlzV4Oi3WrC#tG+J(`MV-AVxRSdwe0^oYV+g%hW6fk zI{jvq<mLkXh_yl8OVmZ@9=KuAUZpbkWv*y&erAr#CCRnX+aKJ1zu|T1oWfrYr<Wan zxAlR|F6V>!!_MY<Anlt!?bBX)OIyyn2D>iJC2!*Gb*7Iy?JqU04I_|E%*e-C7c2 z5gX4<*ez4`TF+3%?pTxWA<apTSE){jy!f}~|2(fBD>hrK*|>s(N9?r)r&joZzFG5E zPg(n}@q;UC{}JA&2|<sGdow3+e~<jM%lqr8hU#y5zi+PIoUr)2dc5Rfp$y|&;s5gD zw7zJ&ZMw8~#)F_Fk?Kpjf=xR`)xOAyGp;&zw<=edeYUb#^OJK<1#8N0&YLgF`0rbI z9Z&5JcD~zxzKFeKJv^s#>&B$DKR-MxPnq?k)c(fv{IY=8#d0EZzE~aQC~4<Yl%3E2 ztFqu|-h*5A|8E>VUnUmP_3U7+$lS9<W)IF>p0_GfzrfHTia++*W|RLL-EZapQr$hX zu<@P8^UBV(FT;1dYj+I(w!{C5na4}k<jT`F+x9-+W_nite3!!~^`%>N|L*mbnY5o- zJ8!>fnYYU!=Qv)G<Ik5WJuHoQboJ5C)te1$GHj*KKQ4d2U(58wQQ_c<{mHUi%&V(C zYR+6P5~_>{X*n{zu-)rTOBScl{Fu6F@BIp=D}R(to%Z(ahs{i${MJv@*K?m+>(tFH znzL8t^ySqSE25g0D9_#@Z7F#DwMEtK+r^eX>sB*rt_z%_Ja^(dvuAZ*O~Y0_>)G$u zBU`(D9oIWI?RdZW2b(Hm_CyznZFV@(ZEF26*D@zk>1No;V9xYKPHda!O>1)cx8%s# z#-a!IDO2aG9%|`)_j0$81c#z*^sMjIVhtt74YX7JHDx|kU8s*R<>zXxv)g0$+|J^9 z+;LZ}n+msAe>lB<JIAi5x#Dx4dNOC{yVe(mG2WPBH&ae^9?L=gx~;N$=iA#HLd94A zvbgbOv1e^$yxtqle^=%i7z-Uq(=p75o6c<R)ADQ|CwH62?Q1F7k99OoR!w@M_mn^W z&6(^??w1~>E!4<%`NSsSd)lk`nOE~WTP7v>U2_)xiaN==an9}yz7Oy9`v;mS?A%j( z%YEPO4;7yO8p`i0o{P!PSkm?Tsp68af3I%*s&vuJ<9Dd1=l>a3-}aQ;({;_9CqK#D zAo0|^D~o1E*Dn2YF)Vsg?t1tChkfM_Hx(wu%HEbg{q3%~QM#MbPwg0W>D*}f_7?Zo zlJ|_uk4&EM;nb}%!~Zt^Qv3dV?dwfHuXF#&{?JVmKir#h^KNyBDSv1G_wQSFZL#H^ z@S}@o-<2uahg=>RckfK7x_3E9WuN<Zy?)=C=svH7`)^e3czN~6bj#+}C6iNbE<Bn( zW!1zlDPdDXuRV)9eew3WjhV)-?5`%i3)~g1xo+Y*&-5j`_gPtX*4u{`xV_Vs&(Et1 z;uE@Wv+aMM^9SkwZ_`s-CJ68QJ)`tPZAHzV@2_6035p1wAt)RZ%b~>mlmAA#;lqL} z`$CUD`K)O9cgxD%X@WE6+x`A`t35{`O4LBP`iC5g?$(>1{)9@&tnNR#>g##twMRA{ z>yh?MSQnLXW1lj6&h4b-<z`0@hVR$u-hVyi5bsg(o9ovqIavR?`}MHW$_;OqEiK;~ zE|9Xy?%&3j2|{1y1^wE&^%T36$+G>^AN03>7GGL+CSm76*0&p8Fn)~?;CU=&EBKUS zWwxa9(p3`;_PyJ+?ro5CLFtNLyIwd*e3|EeH!7{_mH(aSjPl=s94R^Fe?_+hG#gL9 zcK5GmfDd!#-$iDbf8{O*d;E6&HuH8|d{%Qw-27Wse~SbAei%&X&OP6B(`nTr9dniH zYPnhSZ-+l^KDD?xr~b+M|FxS>YIB^tFFO67rFoWUV~e{-&+=EFmszOn>GG45J-K&c zOUv_br;ir@j$fL1N910g)#B-Q>{FCK3NDBVsduuI?_8gLUoN{fOkC;X&*a#y8Mz;7 zRrcws@8kM#$MKJKv1&BC@?q^eDfXcXFJ0>kJa7Ik4fUuJ`ZzsUZhzF_i?>D7?>-JX zT<TFNefM=qf6$zzJKT1u&z&fHtSe|)e5tqP<At|=JwBuF<Iy{zK4HsFi8E99JAYP+ zG&)SQy!_~g?FE4yGsDlsGyI(HQmVeX$l>YH?l&`5S1K6H%D!(ujqlLy{<U*|OyX6& zek!D*Ud6F&-ucVt#9w}Rw&S|*@;ep=c6X8<+05b5im%|_uifVSJ?!D^ZIeC<+`J}{ z@luRo?}e@2`Q9QdzIu``zI<T$p8M?7(Z`Mt!VczUMV7IKeM&LDG(}><ll!)sV%fZ( z=kc$Zq_dN0x$3i<^IuzW?VeXD-&Io}b?dgx($Cf+54YS7tyK@?4VjjE;naJvgB}w@ z{iMH3CC;Dw-ShGGV$09lw}sZ5@0)u%vi`|yuiw}51kcqO{bOoAE8FGqyw^p3%AR|r zBBsmUto&a8bgR`hsnvhh-u&j)^CkV()YyZEZ}#8(mocSW_s#bYE}yRZtW>>yOqYSP zUt(IK?)D%Hp}5PbhjzvkKc1kOd*ImSw@VJ+k+h$*;ln!y>GxgFzRg&;vE<2zi;?es zuGM1iT+jP?`dN$WKXO9zw+ipRHX}#R?NQ`~2R9zS-?h}@nC#JO&E1#Gr#+fn^kn`e z@Bb%T%XW50dQ{DtJ8|o(+v{d6HaBeDRdIM?)Do-4Pd=w^zqq%&ddee%qo%*~JknF1 zPq=yS_P?vFa{}XUY~A6fDIdnM_S5IuN!+(*&#vQld46Ee<Liffc#R)U{@?j;X0zmm z!0W$1{A^0O5MutGv)HHc_N->fhi5D%{6G2NmXWrc$&q%c-TPabj)+ffSbSW`;YwEa znw1UqYeG)jCR#q>XLq{CdCWSNRd3$J={5Z_i5Z(CPF@TQzA$-TK(1O1SE)%~=zGI^ z&6o98WPPbCI;<MBa(8pZnTtOb?Atx_t!~K48ok^PUu4w&OnZFu!#lIMNQKx{^2xp8 z&Uzp7_K5e*wGU=C5udr>*3T@v?>{{LWZ#r-{juTioS%PYP2aw){XKh>{ev6V^ERJg z$=J{EnZfh9MD~KuOv+x_I`if=AN>4(bGx_WeyzUaGC>m8a;n4mTY2P+?lEsKt4aAi z+2Hf92P*@ERqFR1t~oVn?mDH12md$w-wQe*ZW!TWQy@lGEfh|EiSy=a}Q;o}c? zMeeO^3T>RuHEH6H3uSq6shKl6kBEP_+wi%5QumR^4l6P;97-95?#F4KJ~y={G3f$t z?KZiTX%Q(=(u(tDUGhFywAE3}Uv<*Fy3k_RCs84VU(AGNo}W_@+#E7J=&o}q)6<1h zwg<*JeAhUd*_Hj}?u2tE`gfSG^3^s=Ek35qw_ER9@A5;Z+CE7BoFCi}{q&Sq-dvX% zCx88kjd{jb@_yPw{s=i{$;+4h_J(k@Zf{y=`+hFd!$7xXr&+bxHf-7T<u-3j!prOK z+Ml$~IX6G=xz3c)89MP+*UrW28+{iXJzTYx|J{Bo_l;=_h5gPwTkH3h(^c0ze(CY7 zyc@fJ7E7P$b3Xs_S@!;^l?hqrXTMP2(>WvP&j0gs+{)KfOt;q1pW5`f_K?5IzCB%! zYA!wAyd|ix=FgVhKVP2SapcI8;>n>>uU?BT<6kLN@a)mtNYAcl-iB-2OSt&&Sk+2q z?EjMByztQU@WykycU^VM$Xsll(cib}tzhl>>o&5s3w}&uxP9;T>tEkWSMEHx@f6=S z`?(*rQ|0wa9{u=zad*=}8QYTw7bo6jzr3jHvqyz!vHpf_+mjR?{<Mit`?kPIqb~fP z{JkP)8}{j^?l?Yb{{HsQt)}yu9yYPFpInXIswv%h{!&5oVYa0;D*Mh`)V%ocUE|P? zX&l?+OrF=@_+J>N`N{8P>KW5Swd4gdp8kC6pS_p2Jof0tTmAEAZ~OCvDt~oc!MA+Z zlvhftgf;w^?&<x+w#!|!FF0oL`>qnVIyYO5m<M~b=Na4azp{B1E6S%{6*Nz#kC(Mq zR$bS1X-&e`ZF}Y{Tde-aT(S5=W-sHH8=`w7URJSud;4#1He0{-wq4zlUF;SMv*vBu zt`g7Tb>g&v-BBZhpDlOyaVHj9#2#Pv<J~=mJM#q71S(A=h4ZhOt@*@xh&M@*o4<X@ z<M!0`*6mHr0q1@47H`b|Hh)iA_y6s)87-WfB(=+J#GA`*d|fI1<H6>+OV^q&{(NEQ zYJam8f2DSretK)W!&W}~{;5^wy|W6RMNZt5$$soi!?o$hy6sP9?`PXG|3uB}Q^9O@ z4$n<bo9oRy$k8x0_}hwaQ%}|M@HMShY|N=Pwb#l2b6uUo?w;lOZ|eRF4{BuRGfdbY z)s{DX|4}20iZ$=&?G80M(AhlQbM=8|Ti3aNxO{1Ej*Wc8lPB-WzZ_{_FTe3)-;Kxr z9%STb2=O1ExAT1AUD5s{*LT`<Py6>VM8)p?t?L>5uFn-d+A7%if3mA*)Z@66{`%m= z`cD&`#e}2zvNnF5WccRl#1jrGwz(5Lc2y_%tGeeO`+Rlv7WP#We&t;eUu4Je)pdp5 z(p{aeSXN1E<S(}AgS2}f4cv-(_jJzgeSA{#tnut#k?A={KZ&Z(lMi8jV;Q<#S95M7 z`?OU*D$1%8Hm7;MQku=7^E0R7{iSmXUrKr3l<b(=aMn~Z$tL&VimeZS$8oMTc{9Bu zcG)%Y-l}3fue^$>z3KOLuS|=S>|Vp`e>OD!$-<Z=?U`>DNu0Sn<7iECe8a@ovC1=! zD+jz=k@D%}0i`XwHz%b|tTWI$y47yKQFmvPm{mB>bz_4&ThuGIw036K{ZfrmyIu55 zH}|`ox!P@}wf}XNYChOLk8w-<ZqMgi*nb^#c)n96+V(<;4HKVE)tPO}U9?(fZq;C$ zc7H?q?dqTLpEv(s*TcG3v;G|4PuY8y#d)8+_ddD&`p?^k8b5G!{<6NAm@Xg?_qqOb z-luop8eWKg&9S*JziP&^<!m3eJ+8O-{&{ZUj`HKeo&Co7ORH;EADZ}M)wSubA6@ix zd35<!%-^*i64-?97v)`$t2X8L<PV!P{gJIn_k=G=D*NIsO8qq_&6Coex}<fAwndfo z{>rD<?W>$$eGZyG=~v#B_QkMh>Ozhtb`5Ah<OxftFH@D|oWl1vzSVgj^!Re;`{nkf z@@=|a=CM7c!WZAWCa!9!lCfH@JF9Y^rEEI?wUzNaJ_i&(daT~}bmxlgRY&%I3=+!w zvSG)Pka~X6s4L5CGkl%QuYWJMN$K8rEc>8Ac!biv&DDqBglIUjy9M(ve#)oCH|=BY zInBeCVvl(Dcm=Lll=v;_|FfbyMin2^FU<e(V)eeqXP2ER>aoo{7MqlJerCMY4r|Aw zhvYK;?(X>HUun+V@a^?)!|)I%O|BPv;`w%YEz~VPP<l#k*~|S>PKlixE}v-2U4P=E zc>AsB$q!x}v)E=|uACibziSH1llR}`#SVXK<aDUWiMsbIVP=(#K-~5w`S+`@@;CSW z_<zHquKLa`fjg%g@(Od!?lp39dQX|2^Yp{Qtxt;n#|PPY)O<bvH#_^nCkYih{m=8Q zg&xjmI`aKn?#Dxu<2OzGQt<9g%)B`dt7jb9bKqL{O2_Kyx99x!sG3!|KR_kZ+0s4# zYTvc%>w7|!_G$;*pYv<amE%kHNq?QXg5Rms?&{Kz^1!(cl?G39r%b<K>-C-GsOhgf z&-9k$1D-!0o#z&|SX}W}@pbL`%LU$lzTJK`S<~Xd?iW3L`Fcg-UoALqEqScAINkA& zHrEI1S+;gX3c{ar+IP>$^7%6NV7gxRtj=kncl?}QpT6LEXGx&BiEQ+*&|gQZP8A;C ze|X}axtcHjKd4UHmZ$nJ!hE-Vw6R<z&rhdaew!~O?Rw0vJKg^0AD?*slm0quu1d>r z=M=flU0Jy_I`PA7weSUBvpgr3%f>i0e*0Zj^Eh!@n$3mE+4Fw4o_XuTxG6pA*}<D` z3tvvKP0hdkq3`~wWiuZ{{BL^Ozy0-(84nwM6zYF&xvPI}qbK_zWz%1gN7tPcQJ5o= z_h$2~m>CRxyN}<j2ou>-vrVjrWvNf8*}ZU88Li+;k}`H_KcBJNt~PXdw0d{e|387t z9aQWp-MT9h!fb@@?+Hs^x+{Lk;&)Tzp6$u92!B4mU9kVqjVeRInH5j&N4G3d|J(C* zw#t`N7FFI8cgKEr&-ZvA@ukoCJ=c}+<#zgCEmyvGt5tv193bEDDrnW|6@CkBW_~g| z75zePX?4pZsoLvHj1O66wVrPLcOvJUhU)op?;q)QS=;}_{l2Fwe!9P2Am&wUkNY%- z2WIPSY@XErocivTQB7q1#$`v8`QuN^-(6F*g=zJh=W3Pn&KzCr|D!cEH&S9jtirT& zvwL^@wB4Go8Is(1aDipA5^MPK`?=3Qu5jdUJ$Uqzui=Uy^^)Hw>$l#RSGg_l!%S<w zpB$T?&+C6IW@S^%9`*k_*YhWfeGMIczvYVh`|O+x*M&-_MRiH6d+Og!y>~dLTI+H4 z>_7g`n!Zih`*Z)V)CjrC3HHfNZ~forN}iq@SLyKYRh->vnc~x|haS$kxVoKr;!TH& zuU<vhGxqv5$68tD*!*pt_jArg-<DaLdxTy5w|3ojDfn^Y#u=X{&kLV_5S=BgeJAlj z(~;#TPx<rr$#y$Ds@}W&6yx-VDQ_ZXoVwOrwfW@@dD{h7r^HRsOF!}Ntd_&2;5je( zr<%H{%Ixj5`~Bfz|Gf$Kd|$>+ytnbi-1d8-E8hFn#=U9|y5IXN=nDHnJI1R^L$-r* zh{@C3Dc3<cgy*Q~uXmp5M-IQ+rW;s&O|0nH0cK7Ad7jf39E)l@!2VcC<l_G&s*~gT z!U}`_#Ge;m_<ikFZmF`31<w}T7m8i!6e9UKxn})Zx9#5~%SB2*ef(+ZBl;)iH+%Ro z$-5hMeD)}G2`PS!n7*X!PiHJs+_o=YdpBnv*u8w!6qPA^KXytr+fFyTSU-R18{RE- z5t&5=hx8aT_D8?e?cJ6f>G1sJ?2>KoE}U#<{i}WBo5Z&DYp+k&YOxId#ud&NQ4p)q zbUXN8i2IpRX>&cMPnVPC@B065YpX-WA9?70gldz}*Vm7GiUhrwBm9cPb!kcH#5tO= zOLi|&VOO6X9Hb=b8KiMib^4Q_ohg%B1O%7G^lw{qgG;kWKx^ufdm-<Ctt~(Mx9{Tn z^`_q>ve%f-fB${+{GHFfYwpY3;rlLHuj5{(wcgW7eU~=Zgyo5!W?T8LK1`l(f7|Le zhh+HeecIJS%|+c!k1a^7{(I@M_`1eZmu@yTe{77f+tO$6tHh=CXq(>Yy7b<;;d$q# zyQ$w0`s$|jatFVR-uXZ8t_XJ(_WmhQzi$xn|NXJd_A8SNDsFbK3ICm}c+k+#z0!{@ ztW~<={d&)aYo%q<XT6Ir?ml65@0PjbrI4xD?(t~2^e^AM)NGG`>%oGj=i8@ES+wok zkHv3&WV!tA$+POJHn4Gu$L&;cShC{WhRVaT5%C5ZJIwU>#QLwSn|Ey@ul9M%x$Aj4 z8rFv2n{!WbI-B79sk@{8oOiO`A~pMep=4EwXt&P2w~P}5|LUvCRUS%{znL1!b$nys zrhE4zwK<M%-zgGvU!qO3eph}`x8DK3sXsE?WWz6{{JY+(rDlF1X7m02`j&Rv;78xe zea;_Szhisi^JzRWIk#mD%iQGNg|vN`%4ivpWvO8LwT4kDlhupw#{GjIbXWSGm2&;@ zCZ{;euW;jstg>|vd-7k}UXk1uHD5op>BF+mciKL2eK8Q<>~}zE&b9Af)>}sOdVZew zI;>rEfzzSqk=}RSB~M@c!|C+goii-=|C&}-p!;c!{ne}QybdQ{*74!}ET3@t;mUR~ z>2}WgWUsxxGYTax{M}bQ`A%E;QbybL-(Q7JX0|SFytln>W7Rr^H4I9>Zr!X}w#@F< z<UZFAwZ+2G@7F!CWDi)yGOzfx#kL0dfPfgDE%&ysoNTu(Q-ty5#Pw&h!)@akmW21e zkMmzCS-D_$a>v!rN6zuxJN@$Jh2t;xy3}S*U3HRI!T!skrl<Vw>KuQKsy``zSr9hi zg8#vl{5zB7DwTfDH9uZ`<($3VV#DKybzk24_^c>%$E?Rsr?lki9Of4Czx|f0>`l%N z<tJ~ycg`<2-8v=V<8#r+X~&Nqju%gsZI9_XxcOLhXkMDmi(`IlTeto!IObTDlwYEv z(HFR_(&D?=miO0wZ(F3qVZ<ZlCN}%if$#HP$SoJr{}BCh$DRktEIG!vLU-P)vz?W_ zBlMl@ereyD=C4vbwI4sTY&^APoA->Gb-(_tS-=rfYJ5xWoDwVlp1>Pc)oC`fPRoVO z{+}=Q>WXw<-G>?N-G9EzUR=Gi_h7P$jQG2g2V}1AyXJJ*sV1^`$Ep+d?AzPr{$BT4 z&GKD+MgQNT7ksB8d1syR(-$pryOX;5<dc(UG=$@JzgE4!hI9Rqr=`C`mju0gxcOek zn|UuyJ3e0t*A>z{+Szse_{+8{Nv9uA6a4Sby5i-BrLo8N@Xg;<wda}k9jD*xFHc#v z>w5B~>1X%!$}eu&pFVS6_?mwbtR?lEPQMSb<=WmbZ~Of6G?jP3j2C{_t&NslsUE~t z?B}uJtIMnUfT_Ma*|+boy8fH{-ifms%@2ah-c&7Wij6rp<sNg?`Xz~f)Yg6pQm~yP z{PfV7tzHdZUd{cdy70sXfi0G{ufy(jN5&OP$5`$)Kfg`Xpo;mY;><rzk2kq;O^C0G zPvX0mv0LR!@hOMvImbWlZmnBipIABdqx5eJIjv_AGj$E0KGf`+&hY+RMw__GJhhk~ zj}CuOdGa=Kx-Z{hX}%>Zd&8!7Hr@!?zV)w_Ye7zG$IeNWDlJtNy<d3lv@JTLHgA5F z#*v+Neudl4w>Y<7z22kwy;a@qQGVz#HAej}H`9++q|dO_(*DpEr=YQV`Py!YryY;- z*xcSW@2d#A6ZpE`>(9p6@)d$=Di0293D>oqC}<IH=#}UmVwH9L=ljm^U0qK<&$p7= zbK>Bt2G3fhizi=Hf686-{}7Xn-O<MUGELdJS}y(G=d}-V^>^;seC_6std*_n_8qmA zxh}gXXW6{|=lU6byYgf7?p)89uY9g!x4(btwY3HR{=aTGe<xXca`uK5F-Pk!w-;@> z`!@P##S_MRldgVy<Dq+$+53`|)UMYXYonOsjTHD@%m00hO*@n<QhMUv0&QXaKXVs- zIq_(YkIjX@V&--cXE;wVSH&gm@A`iGO3;GuoJSTfDzDLKS6tV3?c_Hzw{ue+zUGC; z{Z5RnP?dRnZFBOgB|m4M)^4~aZ7Ff}gt;F3rNxKW{yJ~3^N#W5)~wk#rWbzpS}^Bw zaFWWeM>AVl6ArwLE%-e9M3=awL*1UuCMy&EFnn75@yAt3HesQby1Y#$eTjcWe!ac) zH(UI%x!m)L4#UW`%dE^e*_W4Z-cxgC@3+2(tewlZYx<|}k=mjAB=<;8tMy@NzPq-i zp|%P#F4K5Zzv=Je-fRCTZ^c|b#+!+2D~=xDS)docZ%_PNo|>ZaHK$(rv#-nu{AsUW z-E{lf3+>p|j}Q7i4exJWEB4{_=XH@;u5k<Ie~dZ&?4{!t^R1t5zpi2o=C%q={aH8T z+O_JXo;s4Qd;f9#onM@ou`j%OVr_xgv9s2a(HWB6mTpF^d*v!CbT(|$-C!qpf5&&R zZ3bW0MqUtCzc#7s>)c;+yK<LGY%}#YjMd@Ropr1){zX*Zz9z%Y)7zdd`uBUnf%R9X z9<uYQ?fYhWtGvVbkwQV=%9+vs><_)(<b8Ji@hsbvUnhLeZd>4`6S_|_EABJLUH3;@ zgZ6p_l>L4h@%jZ@O+tO{yX@+DAIeuu>6^Qhng7MV&!5dsn^#R`u+_g+RTo(wpn9QN z?eC<mr?vmfYuK1hUi-Z6$K8ipPnsr#*i2qqe)H(>O@A+EH`M=(FS9Le-(vGT$@D^% zzs#+_NlSJKdGxbq30F(lHhj4>b&~CZD`8xhrm3!(7rM~2Z|#|XlMcMh-F7^7@p)zb z6a8i9RoC^%JJhwFk^SEE(fh^oXDcfsd*(asuiyOG+MsSR({bl(mb;zW*uOlj2<#H6 zJpQ5SNN{Uj7USW`6Smv#x4Le$=V6Iu<-N@(8w=dheb+v2bULw_?}n_zJNDU5-)GI@ zd3s>;*GWapYt;{wl|1VeRx*F^_|E0Db-}R;zS<l1-%#LV=}k8+ny<tB_}S;%kL7MO zURUW&TD|&Q#oP98D}LtdJ<4}9tXJFr<zRLGrrpz5El)Zk^}X+^nWKf<QiHhLM$&hM z)6eWl7VW$xmZf>^VR=B={VBhbs-_o4mUp(iTPq|p(|VzCwzE*1-H%ywQp4IR&F=O! z9$uO8=7Ms@wUtNIxf<(d?>&FyIlG~6){3YfciryCUU6T^e`%lFZ=piBIXmwDOc#95 z9$!(~mu*s>QP<*mX)4#R-&!A^+!9_U?8-9Z{`-xF%gwxh{C+vLQuY!5riOR-BZB{x z9pC4*;Q5Q&_pImsZsz~8kNMWpuuW0?iruS^bKSpg(reR@w(j00lkaP*6t`YDyLR~+ zhs)_JZ`UsPdT#sLFH_`Jv0ZxnV4=TtREbvOmoEL4l4)O>3g^r#lh*xaT-Bt${QK<n z+Pft8HvY-@eqwdXABIofCD$}1Wxs|joL`#uNDCBw-cROP8P+8~xxGJ9j_Lm%(_Kdt z_@pi3<o;y1g9=;`bDxcg6+)lP9+WI?eaJgMe&f#TSKslzO`n_eYB%?R(|cb}S+$DY zX5u|DlZ2_v+NF8@JaJ`>r-EJ6Me<|Yzd2`WdwqM;Hu<>LK8@VYwx`_hXFsu<>}(w` zzxwefw;P{yb#t00G0XK?DQ$aE##8qv<N5Q9a*^d888uT3Jr4`7mE~4uS-zn%_v7Q@ z@S=dM`)2F+KVw}Kcj)o+IHT4d?%yBg3!eEOct7KP_t{-bclh5aK5*o>P~R_arBzqg z$^0$drJj8y&h_k<*G+-{LVxg|zp&}#Un#ZSCl`B#^cMInux!|O`^D~_Um9OMZ(Yy0 zuXL_svB}X{2iz80Jv6QRI&tYLUM3wwzgb(IHb(rAp1j{XOk0fk?(@UZdkYS%;CoRi z`C`@9^*;qK{64pCuV7-}XO14zQ#Y#KsqegIQg+?>^6t$!d&{{@txM0oPFQ+ayUcbr zTiTi!nHcUZ_ndkk#)mA<YWNbd+tbbK%ziG9^*xg5pCtD>^k1I*{DgGMmzJxuuSD!s zeIwKkioOVg9h;j!GM~`T*rRQBPlfk{|ICZrDU$ZBb?5!Lf6x9X{mWABR(VI?tUany zUoGRDSf(E2(f6>C{J<(MF|TEhYS!<V)QZEOCtVZVeEgQPN7S{Ws-3I_;!oA5%-_(M z-&CYA`FtbW(>v=v%iea(-Ief#z4LWJ@zs;>1-abL)N;KOmv?n+|6Q5=`?|29U*5Hk zuF|<5r)Tc9^1L|v;)j5xE%8AzQ8S+}|1NQ3!KUWJ0dN0FNBpgN`07@9^~LG>8!wiG z>Ye*3HtEfjM_q>1x0LQ*6;ImpHSVI4J3mj{-eu*5Zkl1U%lGrIx9wc^xasdzw$p(V z3LEsVzG?sU-y=RJDya8*dH9Q!I`;kUbDNrl@4WwiJ?7WvIUhck-#%aa<$3Mf^mFgh zUiiuxReg!Se(Txm?2zkrOuOn|&q!vzEY0F@{&)4chGY6lObfoNuB-a0|Jm+Gh|HmD z!Fe0vD*v<`53u;F7@faiU7uR8{))X<lCKx<c{y8lv)qEOag)~XJhDdIqu*5TS7DEw zL)F|l8I^q(xfiunyZ&38esSTFIye6rr(ep3Mb&j)ls?hF@7g3^P!Y)+dqF<!57Vc~ z_xI)LJYIg`{kp}A;(G1}ocU;Q{LAI~iUU!7;%6fTAAzFi<}^u30~^lD$rod!j(%u6 z_S#Ej(Xm^~AIl%Fe0a_Co|(s=ncu&JckH{<cVvFs&0UU~uF5L^;&zAqIahBh*s=fq zs<^la(an+_`!#iE?n%|oY4~xvC%Pi}^zXe6f3BL}%X_eV2h+#$BPDlD|JKSm{F(am znfteE^ICFc>eHo#MEXpBJNjR!m}bw*zIyx2)90NhCn}tlf2wdq{_N5t@@F}Z$k%Qc ztiR&kX@Bc_hyA0T9oA*WlmGgL@0$N(@pHlYFLsagcjru4+y3M7gIB8m?tCq}`Bk8P z%e{;~`%BF~l=s$H{PVlFHu0;<KPij7o%Wxae<<(ou>Umqhw}c;{YPc^kId&h&i&Ed zMzH>f^vCP*=XRU#JS@HPztif)XV^mdm;AdODH*R^ew20Pe=XCjh^=pL_FkzkubcVu z<Byp#0sr5AdLyty!;<0a{a@2Q#)`k`WBPjk(GpEx%cfw5zvUja%_%4LI4t=0TWwp3 z!Sfi2fd8gXj>g(QtL0qr|DWAsmoC@0p$q=qUjBV&^^Lz$0soiY>HVeuRN+c}rry$j z-o;E`<A3EX{^wmg<(F3B*UxtDcP{=dS^7Ni|I+V)|Nl<V+V<N1Yry}#$F6$+z1|<q z|1~~J|7-lE=3n|RUH@*Me|hpR{rewSR-1o~|J3|T|EKHU>3^2|^R5g2fAvAdf`7l~ zJZ}A^|Iy{|^oqs*ZhzXkCBLFb;nE+)x+AqQPCMi;{Rv#S&QR&|8<tCd1g}4QQB}2y z`_dmn>yS;qA5Kp7`0PF-dB&~I)#n^OyH`Dlzw}eTw&5B7&iivtm)&MzJS)%CKJ6m^ z8l}c({e2TVKWSA>Z+zCTKdrlLO1=lvEc+g1lWO7i;KpbD+%Lb@T>n1P;WPWwJtp-5 z@7zAKKUFZPzYyc~nf)xMiMfFJ7g_h;H_q2i{bspq{*J}ZFHXJZ{`t7w(a0QIlX{E7 zub2K@JifO58Gn@iS@|XGXXTfOpOw#@GtZ>n=INJZKb!w7_}Tnt(a+{TOMW*0S@yH} z_(zA&$7B1%XW4hC&$90@pJgw`K5d5%7l-?)cqWnOC9mf2Ii8Bw)442Rx~Pmzv;L*Q zob1vMo798;U8+AC^0s6(YtTO-x#*a!?_WAD_;WP9wfjrRKKX!up9(bND}53q0{%VP z>%DpPH6t;N`k!@YeBKAuDh2%8^rO<_%8_Xy3;wu9*|`|{>u_t-Z!F5F+>;iq7x2&N zIoFi^0q4AbazD+RvOj>&X(vaaW3krTU=#Cz`+t5~_x;kk>;GHEwtvUs?-xTisI{zZ z|1|x8Xj{(SDf@FCTyy@(y>|L1@k`C0#4mM!n*KL)+nRY(_WvlV3jXI*7x2%iF7TgI zUC=+Ly5N7Graf-_G~J<3{8W60`l<L1^HcHjSo9~(RTplaV*i^r{C>E3{}B$odcWPq zd9E9(G~ynAdhRK_g!jqm+KqZKEobLvJN@Uq&{OWFtXD5~apNiviKiTQ+yCj;PpjAd z7$EaE@%8TVoPRkN0@wUiG-66C`dK`&Dl2i>rO+FTt{ixFzO`CR!qZ*-u%}+V+}uQ~ zgWrOz0{+EJ3yGD^n#=q1t@VvRKlK~KdZ$*|`MjEQ@Y(s2KG)|_CF>lX^Gx|1^v-)< zFrWKxfu~AW#54VtzTYzcIoq$k#lKgszw|j{{^IX~<+?GpCAsgaPhZo1y=mUr`Psrh z-&AcY{%O26+{m!)?0irDowuaS>ecpYUGv=h)A*(5&(8+-bFSVtcx61>ezN@Dn14QX zZ`XMxeNO*WEH~YL^8JM||7@l{Rt&6>oo)ZwBF=5cpP3cE4m>@7MEdFZXL+97;>V6I zZV8jB{#RCz`Zj;*%#9II|Fx_CJ-Z*d{fBdj%C5Z7<n#A$y^`V2TgBmIRp0&g{h0#& zkB7Pf(!c%j-K@T<diUQd58jfr^Y``6y}fZI>jgKfyUNYa+j*HU|Kx5`df**mQT?y1 zr15bhtMTK+Uy}Ua{`j8zbhPYXvfP@jrHSY7OZ~p~Ct{yE^QIl2x8Hd`ck}7L#R=Nw z|IGHS$)4YFTXw?JbN64|nQ|`hrA+s^tuJC4&q-#?UwCKwFS(}$Uu4~Oa-3zgKfm%n zHv_}}|1Oc`Hdo*N=krh9@?nXp=fX=Nwq^NKs)96KwI)qSX)zIrRGi)OCg-G_y7r#J z9+#zy5=A&&J2j&G<M)5IdjI+N>$<z^uRp)@c-Okma@TcBGfDyvSFlW(vi4=D#V3iN zH`6}x>zrwy#nayG<#;w>YT4iLH`{{(_WVq^cj1)Vhu~oT@L6H3pWl3|xyZA&IlTGs z>F0OXvH7mw=HGP5bk0SC(zlUM*0Y$-{aki^^{PZa-}z6sd{#ZMdftMqdaJo4{(n7X zu$SwgmZ!m!6T7@JWbgcRziGohZ@;{+PPl&aO}2c#?Y6$nPMfPQp6V@kESdTDrtfN@ zg#A69w~v)qZutKB)G}Yj!~@fRxoN$-k-KB(KKIqmcbTu%^=~`-rmMx2S@3VW|J%m2 z`Hqs^md6B)ugrcNP$~NLzQmE~eeDkqBy&fvy~QAsX>#%J^8l$`-|g=GQq_`qx@}Kq z)fPMX-0kMwAy1DzxFdc4uF5~%W6IAbr=}(gf2f=v_%M7o*O|B{ah3(jtm08Sd$w;4 zclvet+M2}OE<DdC&A0fphGoqyll9NfbRF&d`rS0_i%7$rtxNZvwM$%KbE)3f-+;L+ zv{dKK*H3PH_oVOc&*^TS5_$bsult#g$Ij{Jr`r^^e0vhQsbAcxvF<Fl{_Il52#x;- zRLbKyYCJbxx00XQk<Q-q<=dA#P0G^s?}V?PJ7T@@>-+gJKF!Iqh1XVxi7nXw`p1it z%$5&s{7~x&-^H;z|5U`^*TEb;Gll=0{^Pw-&+SI?pAgPFHWU8+nHpjlck{}V;HW^x z^?T!=iiU+GxGa$E-tJI-+4kt4sVlbU96a{-Ow^8|mgg=?(;g@PF7120k#W|edzvOK zQ<8jNN}2Ec%)j7%(ylcJ1KF}_r`_9Ux<PW=?7j7O*3578J6}J0mAE}i=XS4e(dF$0 z&-lD}7EidV=y1P!o~NDt!ugC_?j_Xh`F3+L&w{_-@}{P2X`ON?@muAr$+OOE-Y{pe z@uj<=-)b|K{1*Gs)nD#8mHFYpFFUX9y|Z70W4()^uzTH0C6Cm{FKf1_zkPd(X@>o) zQ-^-%9J=`H>8+U#H7nX#xMDlJo_tvP`;q^Kx4Rhl=38&$HB7TBba_-=xl{c6#oWRb zhpsudIeea@+$32sYpwWmNj1aZD#uHg6E6qtU3}ilqAFU;X664qdxgGxncS3>@;z?f zTf_S8e~9D!FK7Fg+XQV4H=H)}X)s&yyVLiqu3IcU-u{2<t58FSUyq7z3%-{4BR3^` zXGt#KoWQJ;n?G|L^xJB=@43pPgZ5erK5d&Uy)H=KBF|KSC;rdszo*>%4luu;`Tpkx zuC4x6F~250U|-dG>YtFY(ucK2PG|2~dr{zp)UCSIS25Ogk$iD{zKqY$9NqjzKJ{Tw z^E~lqH&Y$sA3VBjWy89n`en||74D24r_Sy$d~VLHyVa!W`+9A*v)6J9tFMVZF8IYH zroa2z&yTfsp-nY!EE4Bxocw)vAIst6N>Q`Rj1$A^C(KxvVrX?+djHR>0qk<Q#*;-# zC3aj~axU;eHuIi{7v)jqCgP`;sQrtRyU6kV@)h2B)=LtW1SF>{I?TKyIcBxWx;OoS zZ+1?Z`LJN!Ci@HDo4qFg73R})czyRoR;HEtEbcAy7e@ZRJlBfnfp=Hh+~+B`cey#A zO!A9aY~gODeD`8)rizW$rmLka^~Mb!gd*N<`)lCO$RK@QcA27>xLoQ%0qOj&@5C&Z zOKUpR80^n7TRTJW;KSA5!@jvOZh7>fvT^3``y%r_YTj3F&u`u9Vze&(wc<sku(+ID zKCTU4FD=a6_1;62+0$S0wD|woZ-b9rydAOP#*-=cHeCDnXQ}Shn=kIa5SnJZUh2se z$HJOAw!Clh+*1rS+m|n|^i%I#wy<}9(guyq97>P&hQ^z#FKW&Dcyim^p91G?ZXQ@F z(!2ai)R{T#`ONRvOBRLI9Mt&dtNLU{9gj`K|C@W3+>)DlOSUc7O>e`(;0Rx(Bb-}I z{=b=b@@=PJp?7}#Pjeo*x4P#ea>SdX!|iH!I4KK8O{t0Xd+9g#_qElp?oNKGaI4~m z+v{4jTeJ35-w~h9UU*55y(;<hgaZ2;pVIGcaAy?D_e?(8xHw5**WKj;xi>0b=00xy zRkYsejphA<n%f(*R_xW8^Q?9egS;09|A~3YPA3j7a*ed9%*g)lyFPXwUsumB-9KuV zKE|)T`Nt;i>jTyM`(0L9%9Y<*dBpDja>iRl=gdCt|McX{-ubF}dp?~Ke(`(hrE}Z+ zO^-9L-EN(8wtY!w$%4NuElpygAD=WP6f3{BF+It9_MkeKMdg%FT1QR0ic8`z$v$*c zet1>)VCee4cEU2n@tZ#I@KuWzM45{U-K<#7A!#qVwfo_RjXwh;we+{Ex353(r^Qmg zCuYJAlW5a#0rB>1o$q<?X|cz!<u9JKkK3DHs{f4p@-GejH>v`6P0`&wC9-A7a+lBX zTG_s*EUSd8zW)3?<(|PG*XoeX)pnm|@x8v{6V^QM$8J@n?4ply!G7g8CwiT@Ti2L) zes1TcyPy8EmQHu*+<5Y0)9cI?_a)}v7Kq`voVsZFPBE3FxWogeznMoqkob0a=ksWl zQ-8nB-j-WtG^I#)!_J(y?eafkcel-TI`~F&{voUDi+`MNIkRy3S2u-E?q#d%my4H3 z-oHO@+S7LPL-R9DCEO0~yC}Oi?9WDx`q<M|MTuN3b0%p2-DSYac|G@}bzYXwGUdPL zcY5ub^W*KNf&-6vwBP?ZVtl?&?!-ZU?fsiyEjTdKSSRIdaX~xZ`R|d^S_eI^tnoBA z;uih!d~WQWP*!Vm-eaG(EKhbyee5zXkM~cM{dV*Ge4}R`@%rx!Qk$NZ96x_cZ{iEB zoCn&5$6DHy|1Ms5cHN6(^R8EjwAW{pvHps))%|u(bh^iauluBPt|h<oWxJLC`t<^i zhu=z8Top-}`+WTz=VD)_%@(E7R|}`OYp%akt9edPZi)4wZa00ic~@O^<8yx`J02_J z<<2+%sBQW(vQEY>jos)_+*N7Ck5?mi$EO_86zTMTH&G~W#@wDG=KWtIPo6Um7G;j+ z|9sDO!(REAjwRy82Rqdd-c|LeGCFf_L+Bq5Rrmat(RuggUy?p`CVrAm$`Oa7Q+pr# zIh!UFJioV_YcW#~-|YuecWva_mBxNAZq=mr{Qtkd{J6EJdc)(Ol?Cs2pV7SMoc4_O zUh6gC@Rsy{Ub9aHyS-N1_w?Y*uZ8alQ#&g>4lmL^WhfhGobv43Y4@b%%XsonoT{kE zm>Iv}{zd8I)81|tn#9xa@bJ}o*+BLks&_LMX)J5ApJ!VuaPKX<{-png+%La5A37UY zY}1e_(fY;y-rC9Y71vGdm=gZu_xtP^4qwMr#Ufw!Rhh;596o+~b!6%HX<}dZ%h~xJ zmfXK6FK)W-+%9cyw&ku%lWGeVySCR@a!7=f{oJau@a4YbYg|948k`RL!q^`qvn+B+ z@7wSdOiTZ9?>}7Qyy1>rWk$B2sEbd{jIg}Ina1UNug_SNR4Mu5_re*gHaguGSnzdM z9Vd6E_EE;S>2nhoiQJsTQZmozv+r(}MM0LeOka6R<X1g@v2RK^Cx2IQi}*{w(2bK@ z)q6v|X6o$u-g&HBvwLOXj|*bEul(Mg!XOl%+~<AuP}ni$k1x+;SDPF&l<j;UdwcS8 zpSsh|mz=LmPQCN9xqM>(>DHywXE=YAmHYfo@nz<<y-5eR^oGjpEu9*>udH+hTXolO zbu(*`#S^Mt-)po}t4WgUV|u;2qNcM?cz=<p#*>mV`Cwg_mEz53zSxRBSYcdr@pG}9 zTH4H6)$ep0U%BUgWL_K@QL|iM|K3^wqaWGM`eJvF`5!!9o&Gib@HEDE&Zko*eOUhK zYG{nrN`cMCKR?{q;!*SG-$}#&w=!PLE&O%ru1fx&ABVOtQaCB?6dV+>qohl&rq)At z!oTHSmQsxO=EiwcF8=rW<Zki*yw0nV<F2IZKg#PpR<FXp>D4R8k1e9FY`0JNwWoN# znB0`esSQhS?q9xguO`dY2>s73M?bE}ey3nsv6Dl{e8;gVl1CmapP%(t-O2lap=VWS zXQIV91K(LYtakI5Y*AS^@6LmE#n3Y+CQdVo_;0Oo_0>1abAKXV-%a=3ex9-R{*rfR z#J3#e&)7RX?QOlS2@m_N^Un3Y@tYlaR=E3YI%j-QBx>4ji+Pu)2hH~a$7KAf#~l2p zR=0>BwY((x@|~h|zEmgw-~OQ3#;q4UY8F16#+NsD<*)X$%io)>?JxeddEyT}?a&PG z0&SP@ZF0YMoSnS)+cbrjM_o$(&*`xh^{A5Cd~t$Mz`AEI&c-EP{jRw;@%p;x(nb8I z48A<u<N0=9)Qt6SRetPP!cx4-I{y9rVCKW+Tdn!3#2>mm|8iHejzjU3YVYx<5*j-u zd_VB?j#0vv@XG%y&K{4g+1I*sk88n&pBJ-MPY$kp7x;PpzUZF^J_jyY_#)d%Zttew zPUkDwUhMxg=lL#c$;WbH`UUsq@+`mk-R|Dhe=+)sLNDrfJX>&Q+TOl<4wVZ}FMM)V z>-yq3R!>(xUo?4<qMt|U<{fi-?ngAqYi<0!<E2&E((-Su3O~74?&IYTV{dub(O9-N z&PA*y?&H2%1sk5-vi~j1wG8-otjpv6XDE7qm%jMbTtkh94RiTyrQ8^IWWLGhf4|nL z;N6RfyxGOa5B7Pjb+(qc@F6g4ape=9^m)?Svnvz6@Lc#VdUDZu?%Q`7uL+B6Txn>= z#c=EW#F)u9eD(|cIAGVG<$rXV{=@(8n<wZy)oOlSSrY$BnuCAUV-EdOt6S7T5$0Z= zW0#?_FNa$<y6U*`gdcl~cm8{DU~P@czTeNzOy6-lx<yv%rH-w|CqKKLlcL+oMH`k* zf8#h`d}Uq7lI6XRGRyhitbF+5?UW9;w+Ht{{kdYXF?*Vy>N;K#%N&CM_T_OaRD4!U zTfA)DjgyYGYDe1QbMypsoGXtl*^wVFy64wNAL~u)B<4$&?Cw3YuAx^+J8honj_&j+ z{Q6v<v>5et3pncH8rq6uyx(gZ9;lt<UjFsk%8H+kTO+Gphs-^*@A@kCi^h`imlyW^ zwV!?T&#&MRHO{Km|7l16zIi$GSpHMx##a;6Z(8a-U;DS=)uYt=^F54Ly|t{3%6;Q_ zv}vNS$G(Gg+bg@ia?cIue4CtE&7x%+{QW@K@`}u>J3AlP-c7pO@1@TZW&iqY@w^ZZ z{kO4p>rV>r4BLHv+H@cGi489vugt$I5cU4!nTzKdJ{AW1b6@uU@wM%g%jQxw2|3Mk zQURuq&bZz#zASRJr(s`y=KRa$hrbA2sET~3&f0Z<7vtOK%hJ?0tT`#%^LldRb3v(j zA2Y0~3Tw<iJ$zyJr29hcN`F<G75NH#R`e~l@v5BqNb0ENqvjI#nu{JaVJ`Ene>!wF z9{K)Kdwx7i^;CtIiR({1{~Mncp1=~Vzj&SI=0CqyP54ssbnc0NinVI)`SsP>Z!GKr z+AGCUSMUD0YU5*`>vQwvUnNFOxT>DKPHW4P*vE78Be{2`axW9#x8^`@lF(YyJO9I{ zFW0vcuj6k%_Matkwt|Vzo4D*{0bT9VQoN`3O?8Z45|H8%E66wJX1Vc>bsR75Tw+?? z%rsw|;cx#l=^gGf=bvGe=H?BHxxeQ1`t<k`X}yiFduCVv-n{s<P~)q6v5S~}=iXOp z{AHk+y_KbSUxGsX1BWB(|5DV=Gm8Ho`m~ZwL~N?Ov?2S&^R+tnL(-?$-JKU)9dJ6K zyyW`~_72aqjJCq)2mexCE7Q;2WtckWdXSrb?#GSqN<Q|^&pULLHAuAWqnX9;^}e6Y z?)SERyDzC&Ej~Ljc)QyA^dg-<dmlB|6kVx17qjHw+g{m=zZ>_xFMgA`IXZ`JOI(L* zM@#V@-V1dL4Mo3B)6Y$0{%v#l32(OWPF)K-{h#}u?4HQ){#`Y`<--EIsUH^YYkXz4 z>aj+DP@Jlr#yP?GCDxB7yFI$CUfY`cbGuQ@O1Wc~ntNuX_Vscmu|?~<YG3_5^Ny)U z)vxt=Q9g~Sy6*YT8+$f&@_&ukp627y-QIci*q&={&b$KUYUlr-dY_dewm~AZd`0ub zT5aJV`>=<9GHq^unr?6~>yfqm_mIj1rC+<(ANs8I!2DAEx8Ute{+eF1>aVY5unKqU zx#KL{8noa)XMCsY%%ulo7wkB@Lt}|`<ZS+En}$&L54|EW?_|$5$-ddO?atny*5~0; z0rS5eJaGK_N(I3c@56r{3;q*qD0QN}D>s*G>nD!JMO*cB;!Sdb)~vsKTKe0D+N%C6 zz2eoURwrDy->)BhlK*m(wYs?KjfRfQ$unJkhMI`ZufMtNu|x>hzAas~8C|(+?*;XU zcWp3KezIlq1-*c2wF}O^5IM5XJ0{%xy>X))`+wo(c3j2J<j-H#%l)Ns=I6f`-$SR* zFXJk@HktF;>-sNFf|>glMz^c1E0sB%@U2p+`Tox}o*N$S(Y|2K{BzHf-<N7nwkzA5 z5Leu@!hfk<<7M?tugz9Hwy4?gS?i-~Vc!pxyMp&k<|V&A==zd9_+zCg-{!BEs*@Hv z9~4rXH1Ap3^L6h^Ry=-@JMquoCzG9dPe?ER`ta`k(uk#!ztWFBsJ)?SUfCPJLriwl zUh923t~Rgt|J*ZYrcX%MpZ7PGf8A<eqVZPII;_{bIIJ^YYBk%f@+sPR$LG7S{E}~p zHNUc*b3ePt^vTglj}E`?YJAngn%K7b>r;sY-NK2X+oNX8lP&viGj;xXuFbzzY)CD; z9k>59%d6L=!S1s+tzXXg(??${%_l$KkL8msv*`UNd$;#B)@<Y7fBX8W{pZwqvaU;Q zy4T9qJGE~0RW^q$+uDv4-AgIirLePD?K1PH8sF2uOoR#&MW$c)XyN|Bdfxx=rCB;{ z{PVBbZa>4f>+6mk(`7miDgNf!>Q$T^QT<}};j`+35gQMG-M^uAx5ekRyJjW-=`K^O z)xWBhdLuQjXG>gE|IwY<r{7HYCLI=^pCSFc!}GFf<nu_gaEH*2mj^44FDloq3SUsp zwoCu#z83bM1ufz~m9D@V3`f|n79J6Q8v1DRBbkrByP_{|YFWxIwSM#0q#OESmV5Sz zdHP>m{OZ)>yst8Qcb>6KRj{m(ch|M;zP7?7A?~W``^wth1@RS2eoe@-jheUq-mf#8 z8S_r_=~%7}X`ju#k9}UyD(0x~=X9TXM;gy`uzLPA@tp6^2v(8u-p5M6isQo=t}@!R zua(@{l<H8Ld0O$U`|rKo4ZoHgw&Sn<8p3_>&Whfv+m}zBd;F+YrGs72naL|y({~y- z*7&~edVcw_%jJfeZ@&UU#a|zDDM;V{J7Y@U)}?2<Smme3U9(=YduGFj_7w#Z=}UTU z%sF**<J7;I^ClhKU{Q0+<6eMRl?b<CyUM!5hwbmau+bHs|MsEqpF0nE9CJ$>yFa~O z(DmnTd~#vghP2WfU-XUNC2W>he189PgY<Q_1|g@velxiAef?`8Dc1a(ifhaLgI-;K z<CnhV?;nw+JL)}Sj>o_GykzQ`AnqwV&-To|WUaL~`QmP|yVLJUa`1nu<<$T5@C3VJ ze9MPrc2j>W*w^@DkyYQ1C40PTg6AmPPI%{`Z!&4#nXF>%Xj^HEs^8^b)<-(k_O4hY zJZnP5QbF${-EFDXat*PMzW$k!{k`ySlYIT6b$8p1-`uzCfAhqCMHGM2!RQF9dp62n z?A90P2E^p%Zm)RqVYigVdHzqH>+Mg5vWArB&7G`gxGswC3VZFvNTI7CK}^1;UTgN% zr2l<?_0+ZU#u}x~zZkw7EoM@<KPh|Nv!^?=iyCXbgbA-*<^S2Z;a=%KcF@qXZBVt( z<#K`Y%et=I9<6y>^A@J~EL0KF)p!2<RFQf6IbW8yO_Q}&7oJFwO`D=+wq?@VgB~8P z%&8ZmQbpftZHvlScs*YJi|dR+rH-V6072^<0nhAba-APf$4}oK^UhBHqeJVaKUY+@ zJfHn<Z~Whj>=PDeRVoE6d@*f0FQ==U#XG@#^A{p7ZmxeUJ)b4NL1X{S-ScjD$wbbo zyOqm)=I`}FRmWRB8FP5|#Clilh}5&bp_1=nx$ok%2X(erxc>^VE#?co!Bw0iVmqTz zDA)dE^t<&v&qepEKV!V-{#{$yvd8Mc_j^ABJ;F{(Gyb1^F=(0b(z74_@Xk9n;q>ei zxvV|kfB#8b6}p_^jC=9#b63-uM86nZ6VFUl`;c$*UHtzQ)d1PLqkbzT8hUKDymgJ? ze{((0ulspaLjE6~xwF<VZ@MXQd2Uu~?ew$8&3j&Nh=~30aaIn~v$xf^E|i-qnC|=i zY5M$|3VGK%XWak&%fZC?uTnr7`#*_j$y1UI^>ejFuJ>CcxE5dib8!0_fwfmvZ3V0z zPS0OISM0B|!*8d!Z{2N*KU9=;V)ve7Y?rY*J=aZq%FCTc-px@ISW)*X)Z^~ysj*7o zt6b!M)NI<;TY8`SPE6E6-TC}h{yqN;9(`7{Iluq7^edstfQyH;jE=r%d^&&U;@LO! zLh3h_#6Jo)zmWPd-~YIkpzh_hb^-;B*44a|H(%;wl4T3d$m8(Hdc3V&kfrOdW%lpu zoBeJSZu=sz@<`K<D@^Vi?=!A-zk2z7weJeG?8~Pk_sjAvnUJ~CNx4XF{j_WBIR#s@ zjz}LhR^WX-=eC$yf|2sc=hMC@o!Xo5Vf$7V718gu{)$^N0{<oMP`k@j`uNP%Ijp?D z5@L6}dDp@kF*9&$u&(&Uj`_>ug?J7nSWG|jCSXzc(w&DKo-wdg3xydf-M;v6$E$}9 zvm~dTnY7MjrVQ(9D|M~@Sq~OVcztKC65e$n>&33GJLKicK5_h?BH^l*KF5Fm1^(sA zr_LLB{&`w=eEGEU^OJuxd}sca*YwYW=f}Mx_pdEitP*}upXn=hFl?rb!Mu&eFSE_v zw;jG~Ua|A|mFR6H8jPPWUQ(HMb6R)9&my6krn}ry2dXTx`mS!`&`sDIbzJRMo9bSJ z%vIHE)s6esE&IUJ@U<?bM&_253DdXaO&c!nJebY&=K3D>5D&eL4BtLSE!E7_YGQb^ x-MDh%uHvP#M*{y`uT;=1J9qYo731^NhaQj5rfkc$Ot&bv$lvjgpS5}Ae*oo?D5L-Y literal 0 HcmV?d00001 diff --git a/doc/rcfs.tex b/doc/rcfs.tex index 7b1b9a1..ab65444 100644 --- a/doc/rcfs.tex +++ b/doc/rcfs.tex @@ -3106,7 +3106,7 @@ For proper scaling, the heat source is further normalized over the domain using The heat source is then diffused over the domain to propagate the information about the unexplored regions to the agent(s). For that purpose, we use the heat equation: a second-order partial differential equation (PDE) that models heat conduction by relating spatial and time derivatives of a scalar field with \begin{equation} - \frac{d u(\bm{x},t)}{dt}=\alpha\Delta u(\bm{x},t) + s(\bm{x},t), %\beta + \frac{\partial u(\bm{x},t)}{\partial t} = \alpha \Delta u(\bm{x},t) + s(\bm{x},t), %\beta \label{eq:heat} \end{equation} where $u(\bm{x},t)$ corresponds to the temperature field, $\alpha$ is a thermal diffusivity parameter, and $\Delta u= \frac{\partial^2 u}{\partial x_1^2} + \frac{\partial^2 u}{\partial x_2^2} + \ldots$ is the Laplacian of the function $u$. diff --git a/matlab/spline2D_eikonal.m b/matlab/spline2D_eikonal.m new file mode 100644 index 0000000..5330fc4 --- /dev/null +++ b/matlab/spline2D_eikonal.m @@ -0,0 +1,378 @@ +function spline2D_eikonal +%% Gauss-Newton optimization of an SDF encoded with concatenated cubic polysplines, +%% by considering unit norm derivatives in the cost function. +%% In this example, the points used for training are only on the zero-level set (contours of the shape). +%% +%% Copyright (c) 2023 Idiap Research Institute <https://www.idiap.ch/> +%% Written by Sylvain Calinon <https://calinon.ch> +%% +%% This file is part of RCFS <https://robotics-codes-from-scratch.github.io/> +%% License: MIT + + +%% Parameters +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +param.nbFct = 4; %Number of basis functions for each dimension +param.nbSeg = 4; %Number of segments for each dimension +param.nbIn = 2; %Dimension of input data (here: 2D surface embedded in 3D) +param.nbOut = 1; %Dimension of output data (here: height data) +param.nbDim = 40; %Grid size for each dimension (same as sdf01.mat) +param.nbIter = 100; %Maximum number of iterations for Gauss-Newton optimization +param.qt = 1E0; %Tracking weight for Gauss-Newton optimization +param.qn = 1E-4; %Norm weight for Gauss-Newton optimization + +%Input array +[T1, T2] = meshgrid(linspace(0,1,param.nbDim)); + +%Points to be used for the eikonal cost +%id_eikonal = [[1:param.nbDim], param.nbDim^2-[0:param.nbDim-1], [0:param.nbDim-1]*param.nbDim+1, [1:param.nbDim]*param.nbDim]; %Points on the borders + +%id_eikonal = ceil(rand(1,400) * (param.nbDim^2-1)); %Points randomly distributed + +st = 3; +idTmp = 1:st:param.nbDim; +id_eikonal = idTmp' + (idTmp - 1) * param.nbDim; %Points uniformly distributed +id_eikonal = id_eikonal(:)'; +T_eikonal = [T1(id_eikonal); T2(id_eikonal)]; + +%Prior surface (SDF generated by demo_SDF01.m) +load('../data/sdf_circle01.mat'); +x00 = y'; + +%Reference surface (SDF generated by demo_SDF01.m) +load('../data/sdf01.mat'); +ox0 = y'; +dx0 = dx; + + +%% Precomputation of Bézier curve in matrix form (as in https://youtu.be/jvPPXbo87ds?t=434 or https://dergipark.org.tr/en/download/article-file/1633884) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%Bézier curve (as in https://youtu.be/jvPPXbo87ds?t=434 or https://dergipark.org.tr/en/download/article-file/1633884) +param.B0 = zeros(param.nbFct); %Polynomial coefficients matrix +for n=1:param.nbFct + for i=1:param.nbFct + param.B0(param.nbFct-i+1,n) = (-1)^(param.nbFct-i-n) * -binomial(param.nbFct-1, i-1) * binomial((param.nbFct-1)-(i-1), (param.nbFct-1)-(n-1)-(i-1)); + end +end +%Matrices for a concatenation of curves +param.B = kron(eye(param.nbSeg), param.B0); + +%Bézier curve constraint: Last control point and first control point of next segment should be the same (w4-w5=0, ...), +%and the two control points around should be symmetric (w3-2*w5+w6=0, ...),see https://youtu.be/jvPPXbo87ds?t=1188 +param.C0 = blkdiag(eye(param.nbFct-4), [[1; 0; 0; -1], [0; 1; 1; 2]]); %w4=w5 and w6=-w3+2*w5 +param.C = eye(2); +for n=1:param.nbSeg-1 + param.C = blkdiag(param.C, param.C0); +end +param.C = blkdiag(param.C, eye(param.nbFct-2)); + +param.BC = param.B * param.C; +param.M = kron(param.B * param.C, param.B * param.C); %Resulting transformation matrix + + +%% Encoding with basis functions +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +nbT = param.nbDim / param.nbSeg; %Number of elements for each dimension of the grid +t = linspace(0, 1-1/nbT, nbT); %Grid range for each dimension +[param.Psi_all, param.dPsi_all] = computePsiGridFast(t, param); + +%Keep rows corresponding to the points used for the eikonal cost +%param.Psi_eikonal = param.Psi_all(id_eikonal,:,:); +param.dPsi_eikonal = param.dPsi_all(id_eikonal,:,:); + + +%% Generate training data +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +w0 = param.Psi_all \ x00; %Weights for prior SDF + +w_all = param.Psi_all \ x0; +x_all = param.Psi_all * w_all; %SDF reference + +dx(1,:) = param.dPsi_all(:,:,1) * w_all; +dx(2,:) = param.dPsi_all(:,:,2) * w_all; + +h = figure('position',[10 10 2600 1000]); +msh = contour(T1, T2, reshape(x_all, [param.nbDim, param.nbDim]), [0.0, 0.0]); +T_contour = msh(:,2:end); +%T = size(T_contour,2); + +%Compute distances and derivatives at the desired points +[param.Psi_contour, param.dPsi_contour] = computePsiList(T_contour, param); +x_contour = param.Psi_contour * w_all; +dxmat(1,:) = param.dPsi_contour(:,:,1) * w_all; +dxmat(2,:) = param.dPsi_contour(:,:,2) * w_all; + +%%Compute distances and derivatives at the desired points using the raw discrete data +%t12 = [T1(:)'; T2(:)']; +%for t=1:T +% etmp = t12 - repmat(T_contour(:,t), 1, size(t12,2)); +% [~,id] = min(etmp(1,:).^2+etmp(2,:).^2); %Find closest point on precomputed grid +% x_contour(1,t) = x(id); +% dx_contour(:,t) = dx(:,id); +%end + + +%% Batch estimation using distances +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +wb_contour = param.Psi_contour \ x_contour; +xb = param.Psi_all * wb_contour; +param.Mu = param.Psi_contour * wb_contour; %SDF reference + +%Reconstruction of derivatives +dxb(1,:) = param.dPsi_all(:,:,1) * wb_contour; +dxb(2,:) = param.dPsi_all(:,:,2) * wb_contour; +Jbn = sum(dxb.^2); %Squared norm of gradients + + +%% Gauss-Newton estimation using distances and unit norm derivatives +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +w = w0; +%w = zeros(size(w_all)); +%w = wb_contour + randn(size(wb_contour)) * 1E-10; %Initialization +for n=1:param.nbIter + [ft, Jt] = f_track(w, param); %Residual and Jacobian for SDF match objective + [fn, Jn] = f_norm(w, param); %Residual and Jacobian for unit norm objective + dw = (Jt'*Jt * param.qt + Jn'*Jn * param.qn + eye(size(Jn,2))*1E-12) \ (-Jt' * ft * param.qt - Jn' * fn * param.qn); %Gauss-Newton update + + %Estimate step size with backtracking line search method + alpha = 1; + cost0 = ft' * ft * param.qt + fn' * fn * param.qn; %Cost + while 1 + w_tmp = w + dw * alpha; + ft_tmp = f_track(w_tmp, param); %Residual for SDF match objective + fn_tmp = f_norm(w_tmp, param); %Residual for unit norm objective + cost = ft_tmp' * ft_tmp * param.qt + fn_tmp' * fn_tmp * param.qn; %Cost + if cost < cost0 || alpha < 1E-3 + break; + end + alpha = alpha * 0.5; + end + w = w + dw * alpha; %Gauss-Newton update step + + if norm(dw * alpha) < 1E-3 + break; %Stop optimization when solution is reached + end +end +disp(['iLQR converged in ' num2str(n) ' iterations. Cost: ' num2str(cost)]); + +x = param.Psi_all * w; %Reconstruction of SDF + +%Reconstruction of derivatives +dx(1,:) = param.dPsi_all(:,:,1) * w; +dx(2,:) = param.dPsi_all(:,:,2) * w; +Jn = sum(dx.^2); %Squared norm of gradients + + + +%% Plot evolution +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +colormap(repmat(linspace(1,.4,64),3,1)'); + +%Solution with cost on matching SDF reference +subplot(1,2,1); hold on; axis off; + +title('Without eikonal cost','fontsize',30); +%title(['Cost on matching SDF reference (cost on unit norm: ' num2str(sum((Jbn-1).^2)) ')'],'fontsize',30); +%Reconstructed surface +%surface(T1, T2, reshape(xb, [param.nbDim, param.nbDim])-max(xb), 'FaceColor','interp','EdgeColor','interp'); %SDF +surface(T1, T2, reshape(Jbn, [param.nbDim, param.nbDim])-max(Jbn), 'FaceColor','interp','EdgeColor','interp'); +contour(T1, T2, reshape(xb, [param.nbDim, param.nbDim]), [-1:.02:1], 'linewidth',2); %'color',[0 0 0] +%Zero contours (boundaries of object) +msh = contour(T1, T2, reshape(xb, [param.nbDim, param.nbDim]), [0, 0]); +%plot(msh(1,2:end), msh(2,2:end), '-','linewidth',4,'color',[0 0 .8]); +id = 1; +while id<size(msh,2) + n = msh(2,id); + plot(msh(1,id+1:id+n), msh(2,id+1:id+n), '-','linewidth',4,'color',[0 0 .8]); + id = id + n + 1; +end +%Observed points +plot(T_contour(1,:), T_contour(2,:), '.','markersize',22,'color',[.8,0,0]); +axis tight; axis equal; axis([0 1 0 1]); axis ij; +set(gca,'position',[0 0 .5 .9]); +%sub_pos = get(gca,'position'); % get subplot axis position +%set(gca,'position',sub_pos.*[1 1 1.2 1.2]) % stretch its width and height + + +%Solution with cost on matching SDF reference and unit norm derivatives +subplot(1,2,2); hold on; axis off; +title('With eikonal cost','fontsize',30); +%title(['Cost on matching SDF reference and unit norm derivatives (cost on unit norm: ' num2str(sum((Jn-1).^2)) ')'],'fontsize',30); +%Reconstructed surface +%surface(T1, T2, reshape(x, [param.nbDim, param.nbDim])-max(x), 'FaceColor','interp','EdgeColor','interp'); %SDF +surface(T1, T2, reshape(Jn, [param.nbDim, param.nbDim])-max(Jn), 'FaceColor','interp','EdgeColor','interp'); +contour(T1, T2, reshape(x, [param.nbDim, param.nbDim]), [-1:.02:1], 'linewidth',2); %'color',[0 0 0] +%Zero contours (boundaries of object) +msh = contour(T1, T2, reshape(x, [param.nbDim, param.nbDim]), [0, 0]); +%plot(msh(1,2:end), msh(2,2:end), '-','linewidth',4,'color',[0 0 .8]); +id = 1; +while id<size(msh,2) + n = msh(2,id); + plot(msh(1,id+1:id+n), msh(2,id+1:id+n), '-','linewidth',4,'color',[0 0 .8]); + id = id + n + 1; +end +%plot(T1, T2, '.','markersize',8,'color',[.8,.4,0]); +%Eikonal points +plot(T_eikonal(1,:), T_eikonal(2,:), '.','markersize',25,'color',[0,.6,0]); + +%Observed points +plot(T_contour(1,:), T_contour(2,:), '.','markersize',25,'color',[.8,0,0]); +axis tight; axis equal; axis([0 1 0 1]); axis ij; +set(gca,'position',[.4 0 .5 .9]); + +%size(T_eikonal) +%size(T_contour) + +%print('-dpng','graphs/spline2D_unitGradientNorm01.png'); +waitfor(h); +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%% +%Residual and Jacobian for tracking objective +function [f, J] = f_track(w, param) + f = param.Psi_contour * w - param.Mu; %residual + J = param.Psi_contour; %Jacobian +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%Residual and Jacobian for unit norm objective (for points on the whole workspace) +%function [f, J] = f_norm(w, param) +% dx(1,:) = param.dPsi_all(:,:,1) * w; +% dx(2,:) = param.dPsi_all(:,:,2) * w; +% +% %Fast version +% dxn = sum(dx.^2)'; %norm +% f = dxn - 1; %residual +% M = kron(ones(1,size(dx,2)), dx') .* kron(eye(size(dx,2)), ones(1,2)); +% J = 2 * M * reshape([param.dPsi_all(:,:,1), param.dPsi_all(:,:,2)]', [size(w,1),size(dx,2)*2])'; %Jacobian +% +%% %Slow version +%% f = []; +%% J = []; +%% for t=1:size(dx,2) +%% ftmp = dx(:,t)' * dx(:,t) - 1; +%% Jtmp = 2 * dx(:,t)' * squeeze(param.dPsi_all(t,:,:))'; +%% f = [f; ftmp]; +%% J = [J; Jtmp]; +%% end +%end + +%%%%%%%%%%%%%%%%%%%%%%%%%%% +%Residual and Jacobian for unit norm objective (for a selection of points used for the eikonal cost) +function [f, J] = f_norm(w, param) + dx(1,:) = param.dPsi_eikonal(:,:,1) * w; + dx(2,:) = param.dPsi_eikonal(:,:,2) * w; + + %Fast version + dxn = sum(dx.^2)'; %norm + f = dxn - 1; %residual + M = kron(ones(1,size(dx,2)), dx') .* kron(eye(size(dx,2)), ones(1,2)); + J = 2 * M * reshape([param.dPsi_eikonal(:,:,1), param.dPsi_eikonal(:,:,2)]', [size(w,1),size(dx,2)*2])'; %Jacobian + +% %Slow version +% f = []; +% J = []; +% for t=1:size(dx,2) +% ftmp = dx(:,t)' * dx(:,t) - 1; +% Jtmp = 2 * dx(:,t)' * squeeze(param.dPsi_eikonal(t,:,:))'; +% f = [f; ftmp]; +% J = [J; Jtmp]; +% end +end + +%%%%%%%%%%%%%%%% +function b = binomial(n, i) + if n>=0 && i>=0 + b = factorial(n) / (factorial(i) * factorial(n-i)); + else + b = 0; + end +end + +%%%%%%%%%%%%%%%% +function [Psi, dPsi] = computePsiGridFast(t, param) + %Time parameters matrix to compute positions + T0 = zeros(size(t,2), param.nbFct); + for n=1:param.nbFct + T0(:,n) = t.^(n-1); + end + %Matrices for a concatenation of curves + T = kron(eye(param.nbSeg), T0); + + %Derivatives + %Time parameters matrix to compute derivatives + dT0 = zeros(size(t,2), param.nbFct); + for n=2:param.nbFct + dT0(:,n) = (n-1) * t.^(n-2) * param.nbSeg; %Warning: n=1 creates 0/0 + end +% dT0 = zeros(size(t,2), param.nbFct); +% for n=1:param.nbFct-1 +% dT0(:,n+1) = n * t.^(n-1); +% end + %Matrices for a concatenation of curves + dT = kron(eye(param.nbSeg), dT0); + + %Second derivatives (to compute Hessians analytically) + ddT0 = zeros(size(t,2), param.nbFct); + for n=3:param.nbFct + ddT0(:,n) = (n-1) * (n-2) * t.^(n-3) * param.nbSeg^2; + end +% ddT0 = zeros(size(t,2), param.nbFct); +% for n=1:param.nbFct-2 +% ddT0(:,n+2) = n * (n-1) * t.^(n-2); +% end + %Matrices for a concatenation of curves + ddT = kron(eye(param.nbSeg), ddT0); + + %Transform to multidimensional basis functions + %Psi = kron(s(1).phi, s(2).phi); %Surface primitive + Psi = kron(T, T) * param.M; %Surface primitive + + %Transform to multidimensional basis functions + %dPsi1 = kron(s(1).dphi, s(2).phi); + %dPsi2 = kron(s(1).phi, s(2).dphi); + dPsi(:,:,1) = kron(dT, T) * param.M; %Computation leveraging Kronecker properties + dPsi(:,:,2) = kron(T, dT) * param.M; %Computation leveraging Kronecker properties +end + + +%%%%%%%%%%%%%%%% +%Fast version +function [Psi, dPsi] = computePsiList(Tmat, param) + Psi = []; + dPsi1 = []; + dPsi2 = []; + %Compute Psi for each point + for k=1:size(Tmat,2) + T = zeros(size(Tmat,1), param.nbFct); + dT = zeros(size(Tmat,1), param.nbFct); + %Compute Psi for each dimension + for d=1:size(Tmat,1) + tt = mod(Tmat(d,k), 1/param.nbSeg) * param.nbSeg; %Compute residual within the segment in which the point falls + id = round(Tmat(d,k) * param.nbSeg - tt); %Determine in which segment the point falls in order to evaluate the basis function accordingly + %Handle inputs beyond lower bound + if id < 0 + tt = tt + id; + id = 0; + end + %Handle inputs beyond upper bound + if id > (param.nbSeg-1) + tt = tt + id - (param.nbSeg-1); + id = param.nbSeg-1; + end + %Evaluate polynomials + T(d,:) = tt.^[0:param.nbFct-1]; + dT(d,2:end) = [1:param.nbFct-1] .* tt.^[0:param.nbFct-2] * param.nbSeg; + idl(:,d) = id*param.nbFct + [1:param.nbFct]; + end + %Reconstruct Psi for all dimensions + Mtmp = kron(param.BC(idl(:,1),:), param.BC(idl(:,2),:)); %Resulting transformation matrix + Psi = [Psi; kron(T(1,:), T(2,:)) * Mtmp]; %Surface primitive + dPsi1 = [dPsi1; kron(dT(1,:), T(2,:)) * Mtmp]; %Computation leveraging Kronecker properties + dPsi2 = [dPsi2; kron(T(1,:), dT(2,:)) * Mtmp]; %Computation leveraging Kronecker properties + end + dPsi(:,:,1) = dPsi1; + dPsi(:,:,2) = dPsi2; +end + diff --git a/python/spline2D_eikonal.py b/python/spline2D_eikonal.py new file mode 100644 index 0000000..bc8acf7 --- /dev/null +++ b/python/spline2D_eikonal.py @@ -0,0 +1,320 @@ +""" +Gauss-Newton optimization of an SDF encoded with concatenated cubic polysplines, +by considering unit norm derivatives in the cost function. +In this example, the points used for training are only on the zero-level set (contours of the shape). + +Copyright (c) 2023 Idiap Research Institute <https://www.idiap.ch/> +Written by Guillaume Clivaz <guillaume.clivaz@idiap.ch> and Sylvain Calinon <https://calinon.ch> + +This file is part of RCFS <https://robotics-codes-from-scratch.github.io/> +License: MIT +""" + +import numpy as np +from matplotlib import pyplot as plt + + +def binomial(n, i): + if n >= 0 and i >= 0: + b = np.math.factorial(n) / (np.math.factorial(i) * np.math.factorial(n - i)) + else: + b = 0 + return b + + +def block_diag(A, B): + out = np.zeros((A.shape[0] + B.shape[0], A.shape[1] + B.shape[1])) + out[: A.shape[0], : A.shape[1]] = A + out[A.shape[0] :, A.shape[1] :] = B + return out + + +def computePsiList(Tmat, param): + """From a matrix of values (nbIn, N), compute the concatenated basis functions + (Psi = surface primitive, dPsi = Psi derivatives)to encode the SDF. + Tmat should have values [t1, t2, ..., tnbIn] in the range [0, 1]. + + By example, for 2D: + 1 Tmat is a vector (2, N), [t1, t2] to compute distances. + 1. It determines to which segment t1 and t2 correspond, to select the + corresponding part of the polynomial matrix. + 3. t residuals are used to fill T = [1, t, t**2, t**3] and dT=[0, 1, 2t, 3t**2] + 4. Psi and dPsi are computed with kronecker product (see 6.3 and 6.4 of RFCS) + """ + Psi = [] + dPsi1 = [] + dPsi2 = [] + # Compute Psi for each point + for k in range(0, Tmat.shape[1]): + T = np.zeros((Tmat.shape[0], param.nbFct)) + dT = np.zeros((Tmat.shape[0], param.nbFct)) + idl = np.zeros((4, 2), dtype="int") + + # Compute Psi for each dimension + for d in range(0, Tmat.shape[0]): + # Compute residual within the segment in which the point falls + tt = np.mod(Tmat[d, k], 1 / param.nbSeg) * param.nbSeg + # Determine in which segment the point falls in order to evaluate the basis function accordingly + id = np.round(Tmat[d, k] * param.nbSeg - tt) + # Handle inputs beyond lower bound + if id < 0: + tt = tt + id + id = 0 + # Handle inputs beyond upper bound + if id > (param.nbSeg - 1): + tt = tt + id - (param.nbSeg - 1) + id = param.nbSeg - 1 + + # Evaluate polynomials + p1 = np.linspace(0, param.nbFct - 1, param.nbFct) + p2 = np.linspace(0, param.nbFct - 2, param.nbFct - 1) + T[d, :] = tt**p1 + dT[d, 1:] = p1[1:] @ tt**p2 * param.nbSeg + idl[:, d] = id.astype("int") * param.nbFct + p1 + + # Reconstruct Psi for all dimensions + Mtmp = np.kron(param.BC[idl[:, 0], :], param.BC[idl[:, 1], :]) + Psi.append(np.kron(T[0, :], T[1, :]) @ Mtmp) + dPsi1.append(np.kron(dT[0, :], T[1, :]) @ Mtmp) + dPsi2.append(np.kron(T[0, :], dT[1, :]) @ Mtmp) + + dPsi = np.dstack((np.array(dPsi1), np.array(dPsi2))) + Psi = np.array(Psi) + return Psi, dPsi + + +def computePsiGridFast(t, param): + """Compute the concatenated basis functions (Psi = surface primitive, + dPsi = Psi derivatives) to encode the SDF. + """ + + # Time parameters matrix to compute positions + T0 = np.zeros((t.shape[0], param.nbFct)) + for n in range(param.nbFct): + T0[:, n] = t**n + + # Matrices for a concatenation of curves + T = np.kron(np.eye(param.nbSeg), T0) + + # Derivatives + # Time parameters matrix to compute derivatives + dT0 = np.zeros((t.shape[0], param.nbFct)) + for n in range(1, param.nbFct): + dT0[:, n] = n * t ** (n - 1) * param.nbSeg + + # Matrices for a concatenation of curves + dT = np.kron(np.eye(param.nbSeg), dT0) + + # Transform to multidimensional basis functions + Psi = np.kron(T, T) @ param.M + + # Transform to multidimensional basis functions + dPsi = np.kron(dT, T) @ param.M + dPsi2 = np.kron(T, dT) @ param.M + dPsi = np.dstack((dPsi, dPsi2)) + return Psi, dPsi + + +# Residual and Jacobian for tracking objective +def f_track(w, Psi_contour, Mu): + """Residual and Jacobian for distance tracking objective""" + f = Psi_contour @ w - Mu + J = Psi_contour + return f, J + + +def f_norm(w, dPsi_eikonal, param, compute_dPsi=True): + """Residual and Jacobian for unit norm (eikonal) objective""" + dx = np.hstack((dPsi_eikonal[:, :, 0] @ w, dPsi_eikonal[:, :, 1] @ w)) + f = np.sum(dx**2, axis=1) - 1 + f = f.reshape(-1, 1) + if not compute_dPsi: + return f, None + + # Filling dxn (of size N, 2*N) wht diag vector [dx, dy] of size (N, 2) + N = param.nbDim**param.nbIn + dxn = np.zeros((param.nbDim**param.nbIn, 2 * param.nbDim**param.nbIn)) + dxn[np.arange(N), np.arange(0, 2 * N, 2)] = dx[:, 0] + dxn[np.arange(N), np.arange(1, 2 * N, 2)] = dx[:, 1] + + # # Other approach + # I = np.eye(N) + # I = np.kron(I, np.ones((1,param.nbIn))) + # O = np.ones((1,N)) + # dxn = np.multiply(np.kron(O, dx), I) + + J = 2 * dxn @ param.dPsi_reshaped + return f, J + + +# General parameters +# =============================== + +param = lambda: None # Lazy way to define an empty class in python + +# Concatenated basis functions +param.nbFct = 4 # Number of Bernstein basis functions for each dimension +param.nbSeg = 4 # Number of segments for each dimension +param.nbIn = 2 # Dimension of input data (here: 2D surface embedded in 3D) +param.nbOut = 1 # Dimension of output data (here: height data) +param.nbDim = 40 # Grid size for each dimension (same as sdf01.npy) + +# Gauss-Newton optimization +param.nbIter = 100 # Maximum number of iterations +param.qt = 1e0 # Weighting factor for distance tracking +param.qn = 1e-4 # Weighting factor for eikonal objective + +# Bézier constraint matrix: B is the polynomial matrix, C the continuity constraints +B0 = np.zeros((param.nbFct, param.nbFct)) +for n in range(1, param.nbFct + 1): + for i in range(1, param.nbFct + 1): + B0[param.nbFct - i, n - 1] = ( + (-1) ** (param.nbFct - i - n) + * (-binomial(param.nbFct - 1, i - 1)) + * binomial(param.nbFct - 1 - (i - 1), param.nbFct - 1 - (n - 1) - (i - 1)) + ) +B = np.kron(np.eye(param.nbSeg), B0) +C0 = np.array([[1, 0, 0, -1], [0, 1, 1, 2]]).T +C0 = block_diag(np.eye(param.nbFct - 4), C0) +C = np.eye(2) +for n in range(param.nbSeg - 1): + C = block_diag(C, C0) +C = block_diag(C, np.eye(param.nbFct - 2)) +param.BC = B @ C +param.M = np.kron(param.BC, param.BC) + + +# Load and generate data +# =============================== + +# Load distance measurement Xm. Further, only points on contours (SDF level=0) +# will be used. +data = np.load("../data/sdf01.npy", allow_pickle="True").item() +x0 = data["y"].T +dx = data["dx"] # dx = np.zeros((2, x0.shape[0])) + +# Compute Psi and dPsi with basis functions encoding +nb_t = param.nbDim / param.nbSeg # Number of elements for each dimension of the grid +t = np.linspace(0, 1 - 1 / nb_t, int(nb_t)) # Grid range for each dimension +Psi, dPsi = computePsiGridFast(t, param) + +# In this example, all the points sampled are used for the Eikonal objective, but +# a subsample can also be used +Psi_eikonal, dPsi_eikonal = Psi, dPsi + +# Superposition weights estimated from least-squares with SDF reference +w_all = np.linalg.pinv(Psi_eikonal) @ x0 +x_all = Psi_eikonal @ w_all + +dx[0:1, :] = (dPsi_eikonal[:, :, 0] @ w_all).T +dx[1:2, :] = (dPsi_eikonal[:, :, 1] @ w_all).T + +# Sample points on the contour of the shape (using pyplot contour) +# to compute distances and derivatives at the desired points +T1, T2 = np.meshgrid(np.linspace(0, 1, param.nbDim), np.linspace(0, 1, param.nbDim)) +T1 *= param.nbDim +T2 *= param.nbDim +fig, ax = plt.subplots() +msh = ax.contour(T1, T2, x_all.reshape([param.nbDim, param.nbDim]).T, levels=[0.0]) +T_contour = msh.collections[0].get_paths()[0].vertices / param.nbDim +ax.imshow(x0.reshape([param.nbDim, param.nbDim]).T, cmap="gray") +ax.scatter(T_contour[:, 0] * param.nbDim, T_contour[:, 1] * param.nbDim, marker="x") +ax.set_title("Points for the tracking objective sampled on this contour") +plt.show() +Psi_contour, dPsi_contour = computePsiList(T_contour.T, param) +Mu = Psi_contour @ w_all + +## Or batch estimation using distances +# x_contour = Psi_contour @ w_all +# wb_contour = np.linalg.pinv(Psi_contour) @ x_contour +# Mu = Psi_contour @ wb_contour + + +# Solving Gauss-Newton estimation using distances and unit norm derivatives +# ========================================================================= + +# Initialize weights by loading the SDF of a circle to initialize +data_init = np.load("../data/sdf_circle.npy", allow_pickle="True").item() +x00 = data_init["y"].T +w = np.linalg.pinv(Psi_eikonal) @ x00 + +# Reformat dPsi for Jacobian of eikonal objective by spliting vertically in two, +# and filling even rows with the first matrix and odd rows with the seconds. +param.dPsi_reshaped = np.empty([2 * param.nbDim**param.nbIn, w.shape[0]]) +param.dPsi_reshaped[::2, :] = dPsi_eikonal[:, :, 0] +param.dPsi_reshaped[1::2, :] = dPsi_eikonal[:, :, 1] + +for n in range(param.nbIter): + # Residual and Jacobian for SDF match objective + ft, Jt = f_track(w, Psi_contour, Mu) + + # Residual and Jacobian for unit norm objective + fn, Jn = f_norm(w, dPsi_eikonal, param, compute_dPsi=True) + + # Gauss-Newton update (see chapter 6.8 of RCFS.pdf) + J_inv = np.linalg.pinv(Jt.T @ Jt * param.qt + Jn.T @ Jn * param.qn) + dw = J_inv @ (-Jt.T @ ft * param.qt - Jn.T @ fn * param.qn) + + # Estimate step size with backtracking line search method + alpha = 1 + cost0 = ft.T @ ft * param.qt + fn.T @ fn * param.qn + while True: + w_tmp = w + dw * alpha + ft_tmp, _ = f_track(w_tmp, Psi_contour, Mu) + fn_tmp, _ = f_norm(w_tmp, dPsi_eikonal, param, compute_dPsi=False) + cost = ft_tmp.T @ ft_tmp * param.qt + fn_tmp.T @ fn_tmp * param.qn + if cost < cost0 or alpha < 1e-3: + break + alpha = alpha * 0.5 + + # Gauss-Newton update step + w = w + dw * alpha + res = np.linalg.norm(dw * alpha) + print("Iteration : {} , residual : {}, cost : {}".format(n, res, cost)) + + # Stop optimization when solution is reached + if res < 1e-3: + break + +# Estimated SDF and derivatives +# =============================== + +x = Psi @ w +dx[0:1, :] = (dPsi[:, :, 0] @ w).T +dx[1:2, :] = (dPsi[:, :, 1] @ w).T + + +# Plotting results +# =============================== + +x0 *= param.nbDim +x *= param.nbDim +x_all *= param.nbDim + +x0 = x0.reshape([param.nbDim, param.nbDim]).T +x = x.reshape((param.nbDim, param.nbDim)).T +dx = dx.reshape((2, param.nbDim, param.nbDim)).T +x_all = x_all.reshape((param.nbDim, param.nbDim)).T + + +fig = plt.figure(figsize=(20, 10)) +ax0 = fig.add_subplot(131) +ax1 = fig.add_subplot(132) +ax2 = fig.add_subplot(133) + +ax0.imshow(x0, cmap="gray") +ax0.contour(T1, T2, x0, levels=np.linspace(x0.min(), x0.max(), 25)) +ax0.contour(T1, T2, x0, levels=[0.0], linewidths=3.0) +ax0.set_title("Exact SDF") + +ax1.imshow(x_all, cmap="gray") +ax1.contour(T1, T2, x_all, levels=np.linspace(x_all.min(), x_all.max(), 25)) +ax1.contour(T1, T2, x_all, levels=[0.0], linewidths=3.0) +ax1.set_title("SDF estimated from the reference points (least-square)") + +ax2.imshow(x, cmap="gray") +ax2.contour(T1, T2, x, levels=np.linspace(x.min(), x.max(), 25)) +ax2.contour(T1, T2, x, levels=[0.0], linewidths=3.0) +ax2.set_title("SDF with Distance + Eikonal objective") + +plt.show() -- GitLab