From b51f632cf597acb27309eb11327f882c9b54d997 Mon Sep 17 00:00:00 2001 From: Justin Kunimune Date: Mon, 29 Jan 2018 21:05:56 -1000 Subject: [PATCH] IT WORKS!?!? I have spent teh last week attempting to implement proper BFGS optimisation to replace my gradient descent. I have been encountering mysterious problem after mysterious problem, from random noise in my deterministic data to fake derivatives being better than real derivatives. But now it works. And I need to push this _right now_ before it breaks. --- output/graph - optimizer.png | Bin 16236 -> 17575 bytes output/parameters.txt | 21 +-- src/apps/MapOptimizer.java | 297 +++++++++++++++++++++++++---------- src/maps/Projection.java | 13 +- src/utils/linalg/Matrix.java | 199 +++++++++++++++++++++++ src/utils/linalg/Vector.java | 108 +++++++++++++ 6 files changed, 541 insertions(+), 97 deletions(-) create mode 100644 src/utils/linalg/Matrix.java create mode 100644 src/utils/linalg/Vector.java diff --git a/output/graph - optimizer.png b/output/graph - optimizer.png index 1eaad80607bf3d411d5ece3c99c4563cb554279c..ffaa1d93563d5841a991cbf0c429f9b06a988be3 100644 GIT binary patch literal 17575 zcmbWf2UHW*_XkdED;0{=Do{Z{YsG;;83qYUaiB~E0g>fkhy;)s0|cxXDFQ}C_A)9U zyAU>nxR`+ehG7^%W&%heB+Nkae+i>%zvuV=pZ_@;JmJmz?z{JMKleVv&zb0p2uTVF z2ndLrJ#*4bK;SzT`1|mu&ET7h*vEAO0w+Dro;-fxmi<&eJ5FpgLFx0N5>eD8rG#kh zl4s2%IR5%L-YaFed8EilOIdmQ{<-*D!jy9)N)+TbrOVZie$f}fC_{e>zo+*1m0w&> zo!jxlCDUdML@szu-J^cW^tYPIp3ie}pioNS7B~mHv{Q^8iz;{sduoZE z-BUH;Y1_RIEU2Lre||ssd2H@!+jjy2-U6tt0s;>=9q<(pxV`HbW|M%xh=7o^fWU6S zb2QLQnj|?8BiLRKB>so15E6FeQNMTJ1-4$RU6!KggA!h&E#!_mZ}^RfC90 zZ7oJ1183EqZ&#%yxOLa(*ujt$SBR;oX^wMnQbrgzkM{#xKK_0={=s;>qNz(zFSk-0 zEaf0x=VwX^8a9^OZSb99rQ6&jOM)rp9(>(MLt0i~!a zeG+mEqT;ys1n)z&MJ@(Ke8cdj_S$}&=#DIF>vSR<#}M9B1??6hH$05`nfKAGmee60 zV8#>BL;VB3KgHJ^D*kS5jy+Ndt`OdLGG>`R<$5gzAOe9^YN`RkI33R+N3h5n=R^@9>F+b@-C_UvF0^6u zS|iXhgZ#8@&t>S}|84f)uB)2x%b!S#>zmeu1SCa*dE3!~QhddP?3@Go>b@F#q_Ry= zQ1HkB--Q3Mx1;~#>;EmSRee11ie-TjS(7^uz7h5~27Ff@_0@NU*fw2zQ#!@8m<@CJ zW?xY;D7x0ogB40zC+YibirMbY-!w4J7sh2~>c2$=FH3L@h{eo?#TJDRW}?d49@oCz z^2*AJq#|n?P01J*wZC*=w%{?A*4^U`bly68QV;zWMdU^y7vwa|xv9z4ac;~HP2VoR zu(Yi;1xgpI3}oZ4^M3#BR>1GK2ATEg;pxI7)hbDSe`)trmnwgwrFBhux`8*9QUW-f zs~C1Z`*)32%2^N3xY9@kBf&>s%ciK!R9iPRc@lr&&wZCBYV`AW&#yR*$Y^FOq`7>n zD?DNK&6$ogS10$C*~`IKguG>{E-Cr>&!a{ZLsI)N$EUo&X5YUf>NzZ{RU}dFQ_8SO z$#d+6ifKYX;efAYS`9<{0`8-KIMMn6n_R%P4UO3O#8 zE{XONk-GIYBRjoCx6%ksk>6_J0oJ8Y^gUFR;1C+iEPjL3LYLZ}{vX$M`*Qw#x(_bT zC+`FK(EqWtGrqA1FN)eK{;?{sm+td$caQ38wP8$HRu$m!_``~<=#;qkg`y#+^tU3% z6jY&5gA*rCd{<*vHja9uGec2Zs0n;oS~_Ij8LCqqb}uZ<_44K2ugVq+C$v{y)g4irku9s0anbc>>b`#@ zAv^6)whGu~RPW3C(@JGEo^i*70V8flzhXQ0zAIoCV)fb1Gwwt++6!k&A}h}R`KeB;mZ!bn&}3OROiSC0+3iu>K!z zk|>KXM`Q+TL-g`}?y6XiT>&@QX^K+yvE-_YsnIM&yFg+8pWj$ixo(t9<`t^9=O3Kp zB{rcU0a0N;@re@AU|0Bb$jhH}-0>@4aE^w+x`)CuD#}TO0GEtFCrHJnK-Q=&^BbqV z8I@{%H@=t1dwS2&+f`UZKDRr3bGA;Heke)(J5 z<1Hh=0jkb9s+ok;uCU8tje5&Lq&}}fxK(8JMEYdVa?q5Sw@jsxl3$md84(7#3~x-v*NcLyd}pCq7k|= z7eG&)(M=!Usf;FUeO3g-SZ*P$s-VjWZtwg}C9hJ?ja^sdtOT_*c?FK3R zi^hzT5Z9_mx3@)&yEm}#X4KXeted<>Ys*7ZrK*ot7(0r2Sm!Apf>{i{n69;t)!Lh8 z?@v@?JX7xR7U?L4sjTboI#M(4|Liq)Bmze$xdGHeCEpX*_k=`%XeaI9BW^yJ1ZNNK zV>sYtYGWmK9IMYH}DOv5$B4#vxB?~yb3|vyLWsTYSxJRl?Z(PStu%p^L z8M|Br1+HhBjcx$pHi5Vmj zY^o#OFszE@lF=Iv$uP;UMdSNR`euf+L*X%CGMfpw{#m)qxRxNumalXZh$ZV?oEP$H zHn2u?VVwGhPVWG+t&_RSh61c%elvHMvibJX58!(14+1$OHCdxODw;D){#{*%_);~3 z`7%Admek#y+0SStwYOu8l>|{92X(d>*PDZVrB&8yCPbC)5;%4mAcDnVC#-YJ`M*(8UvB`wWr(9JdO~(n$3k8V&)c8O6z4o9=5n^D1W%Ru4zR8L^%cnb;)2g2Ut4_gNJc*p4$fbkZzCEkxa z;h+1wW&d4QPfHHGTL;NM5W9pN*hg*Ay@5`U_)UE^_$I*~E3jhyCpt6KB|uaOg0!bj zZ|ljxN|Nvlv88Yg_4nixY32*MinAmFHPB|3E-TwDqpDO4+s%)5ofD1kpy+0vfx;zz z6QraqtD*fb*n>|h0|7W}tbn&ysAm3-GI6ONZ)O`dR?G>H=mc08Qy;Ycsw{O3;4bbP4 zma&gIj5~@%PY*dNbwDmox5vM}yAi=M@Q%|m^ca{M3PR<1ar!EKsOqz3F3Z;Nnh7b{ zxPmKj%lRj&t|Pet6Mg+>JwFA9cG2w&wL`UiJF+;+a?fFq*mdgUwH_mVPJHy58*ZxO z{fL=9hQJF(4>`ZJnZU_DuDY_Xf}yw%-XHRa-96vfFX1!ppg5UtF`^advh3HhGp*S9 z_1oI$u9t)`SU)~y%p<`NUVux=$WJuruJ#k@Un1tIRk@W~%{;F+$c#2JsD&L~V?DsZ z%gh%ZQAAW9RLtW|`77bxz#=iu2RR~y>p-e1?txnePAAeaP6$GBMBs#txoh3hPBq3b zg2hf9AT>v+j%RwVJP9tOOB5pp#@>~6Y$wPG^6e(LdH=LkRA-*lEv<8TKr*N)Tw`Xg zpMfAX8wPz|#r5vfWd|^{m_ag&y6)WVc4ZULN^U6ImGn^#s_cJQOcFI<|Mnt$X&LIf0IZHSo;HU0`K zj-zPTshSe|M`7Ga?LNb9JW{*AWKOY1y~s2qNT4C*uCO+ew=gRdF!)a)e=^(s)&`6uuXS$2)`qlK{e%7m7)FleAi0)8|0q-JBF$`GpGMTyLK%l z__j33kGJb(;g`m9iearP0om+>iIp>lmKF(dZGrP?rq6Zp5BnZu-nUs$Myp!8Hp?me zd#clgCAmC~C$!j^d7LZAa_o-osLY7T4`v>!oQRSo6?ZE>53{e1hsxI6-s->=ZNJ7q z+YTa)h?)>t<}JDRiN=zAT$|l=FeDf9V7EVJxyJvJ%2g zIPDPh^41;6S?QTb)9cKtz-l-t2Q4T-esTr_FU_iNIne_xS(M_LEI_k`Han|A_2|49 zL^W?g!@uOi&Jg-a40(i3Z;va`bH=&~@I5UF_Ul2`vvLL2*a>ja(7hjaOOQV&jK%6Y zaP*W_mG;9&4UgDWX-km*Sd$|U67z0{cbtgB+^i`2$-q8`(UM-DIAiH`dh{};s-u_k ztM;s$5RY9196?3PJwm{WP(gLAaP5V1|3D9H)@lF%ra#OY3)#v8_B+vcc2)hBe76BYkk@)t?3oQ~rPmYUnRL zZUbTb$@6M_H)3Bsb7|z7TZ+5!z|@0Prv@UaAIVQD`~yEyBB~8vRc;C_u5a8Mxkf8r z#;DFdcWdc~x&STMd@jFVW->f$B;Oh~tBANAU{Kyh5?no8AiJFSmM%YDa8B9RO6q0o zEep)6fz??w*Iiwk{lPVw^|e)2)zS}Az0dx`?^Ja{xWxSmO|n~=`lU|t#h-XgsOwFP z`rJ&`-o1Mz0iIu_f+C9suh9&tk4gT8_7gQfUQ$s@n!u`v^!Us&D?|J-|VY! z?TL?#PxD=NxcHJ|bCP<&IgUbVO-)S$5#(n)tD^g?o0tU7H;tDESS33QjUg7cRsEV<2VS|AWy&S(0BD>9o3N$Ts^i7baXD14 z0!|QJT3Wi82X8?tBL^ZnTtrns_yuTy9E$MAo~g?iV6TujF=6oZ92C*1wW@cxlvweP z{9ymi<`A)We1KU4XHZtTvH7eon6Z%OnGAAq>{cz^=m-7Nez=YmyBTajCWdH5BTeCk zhA%@i)UQ^Cq;dTe%xq|b3as)WXWOM4WUD(0zB@$^frQ5%DeH_Pa&-o#+oSaYM1)WF zcT}{c_kSvx0=Nbf*PBPQnV2Z{UbUq_9SM#QkxKtSIjrk1Jlv2$M9FlL$I>o3>$%u; z?vLGqC!}~Xr9Lhftc>4^3q-DFRIJCaXFvB$CA7*NUsGBcxFS6HqFf8l5P>+yJH|Mr zpg%PSHnDo#|FlIgZh21rY%CAMn6~rXsKB{LsBVw2Plo@rScpK4NEkUO?&4R`~j zwfvFapUQ^GeaBbcwJWrGO4vh5k_G)97LX(S<5>+Kd0c`#HfyG5DzsHQJ^|d6#F*g> zgA6hVsntFDUh*sruztX*EWvJ0`S0yQC#e`=0dQxFc`Nd8ww$7^Z&LeOdEK zch7_cEvU)iXy*-nJHsh@h_}F4kstU7!_eja35-kVI!l)jZ!A;*uI^RHSM`~?Yp1u5(a18c8+MBZm{sgX5Mz zO*Yz96^{wmzA$bfUwo&Qv8$QMpg%FEWWe#^Xr}Avq zGqPhc4`m08Uvo=$UHEV*aU(@{ z#mFzf1YahdtS5@jwwLsp#~I&`caFY}VJfrKJ~@l>+$vvTfC`(X7)Jwro`59+t@uGW z1&xs)RAF1V3zx23+BC4En?!oPp!R8SIzp4qs_=d0@K{r=wVA+b zki>_7SUIUsw=$dPsPA}@nD5lk;g+y34mrV(+yU(WYiZd@LTJ5j9|-)EC0x$(62r$f zjQ~BdK2Oy%?f%jFRvAiY!`Ja51Jyml-5#E5BrkK&K%%7iJZ*GCb#@k=cVqz_bB1~V z8yh5dlC+SV(@sn_N7}SPaZQD4!z+FDWaYzem^S zyZIV4B|P&gGmp65W?wTj zbF2!wVikCf{|h2rwJGB)PaMx-eNTAF8ugaMBZllA3B~EIE=%{tf9*0qr~>|3q}o65 zIX`nvRrq|v&<3}bl)+&3w9;rV%YNognqm=9(PzALILNiYhB6G9y=wgqi!b?}g^uP^ zAI79rf&$;i8bqtNx4*GOR0!AZO7c;Ss>VLKw{`1%bo9OP!|=5{d@-_mG@t%}v6_at zeEF5Cme!{=c7(G?O{+JF8!K~#%98FU@mUZV_02X=1MQm!MOxF2hfWUU zE+Z$8ri@PfR#}P4r+d2uEo<2>b-q!A3kV2zcu)E3C~|Y8?q@#5%?YE}HRqP+>k4+5 zdSE$zE7Rr&K&1hXglRxA;z~c{{Pl%30EGZY)RGLkJ*YD^^a$Uf!kp_cf+Un0`F??@ z2wn7**C5px5DPLNRr%pK*w~Q(Q}*(rN2H@w@Z*WCB3Fv0Y650+UWbz z4ocEgna@!=I;S(OIQun7La3Ub zP^r+hsDF()x4u!dco|g|6FfAX;~G?-aV9OTzbeH=;Q@i-d-|02)-Z=x4W&lhh%2aEw8`x_Rzqy9?bW0 z6Kr*4mm>`u7FSkA{P8vinAeuPw*bopk~F;g^l$E_sMv9MnVxEvFc7${;)GzMwu;%G z2K4s5BK<1!yrjoh1s3VQC*@l_YJJ9Z=_o-YzT)lohcDNVnA$=+GksS%7$Z^2Yq2Jg zi(v%`GEWebKPFVYyRmDP^WR9*BqoZ4T$rjs2E=B9 zp$!#nm)A=AZ~cWa!MXMh*u#R8H+-JPRQni695FLTI7{lnO^Ki@h>c6+okLm;lrGIUggNe>nr!o5|~|!U4BbiJHvSL z`Qq?SV0TLXK_z;|3ay|7!3=x>d?g~ugedICL9kSglKtx(=8N!EPr7cU6C722M}9*F zFZ?vMQ!>5($G|Xs*PE;Z(mpH%w)qdKy}gi=x;{1jSvsaoYYuF!{KYzteA5eGDvQqz#U1W9ic=cU(+x%nuS~khw!`Cm+auv+rcWwsVB(#RqYBe~44r~W zh>^J@F(+VULqWCEnZK;cLQE70I)`>FI0qouI?$k*pHp_pW*%92%Oa_l2BZT7EqI!= zcs{4&;1l(3PgBMZ+l<;xIrTlX_OBheEd#p15q()s#5;5NvvBHX)RgX_1w_ zblVrLIE7hx2W}4B2-%I<=1xAGy)P*-&)W77bF$_PV1ac^?*os)8;X7393Q(;5(9)& zhU*^|TXGPUe(sgtu5XfxNvp0waq|9&;T#mgX>po?ny2e=-|hvHb>UrAet5rLF$4PJ zDdhxB6V7D5c5$o5Xbg0IX2i@2ptqNKXD+aNpSb(vJmAzi5c8Vsxgo_$7+Igi$}k;&b;h6aq4VbR zy_6pTQ691cL8d03;MEPS19hXKL4{cd_r6nQ`8}rrQ?IPO%-MRWkB#ssZAQwGpL9-N z$ZGb~7^u3C4KY3a4b}oT%nnQL->!|DY>HzcRMu|EI#dMEAeDR+ZF{~4nue(G59#_^ z(}DZxR7fhbSB%jJ84K#RM5aap7xOO}@`X*Dc-+FlQKn~}wn?_M=T?KJ zuVbW0FAjJ zjHFZ}I^Mq@8lRjb(OrVlpn1RY(~TPjvH3vKHGN!bg5at>jNWvbTz*i&mvf8Et9L0t z0`O#?P*z%cgJ5|F+{cEOqpfjG%XLc8y=5oT>aX%b_fMdtP_Qh0C^jRG7m0jPf7hAZ z@=C_r16oj=gKSh?($YBZ1KAFOAJPj&85nR7mso0Rd+6fp+hKE$JDW+? zN)bHItAI($@PVu%D0)HA`pm|+0vNbdzUyJqT}~%bf2-)`cBJJ zIFWNxy1Cj9SPRyv^a9V0aW5u1*WUJ!4+T0 z`WGl*{eGh|GQ*kq!BZE*uqB_$p~z-TPy0O5qbm4ORyZqqgLw9LrFyv5R7ic~i zWH zp$Dw?I^zN2>fgzc4<3^Iy|UWUYK)TxW1*g5D+jP*782(IvF%J+Ou*_`x}&34HL2%KAymIAqv|} zP6RF{XZcw7Ua+5eR&%f7 z9N%)wle~DfqZ?EuzgqeAJtBZXvzJ04T4eN!ra&?1y7>92j$CM3ZwqDa9 zdGT<|?R$evte_Jm;4e7G-%sXGKKT6?+@12YVyRNy#3z}zr^~-CqBU+Pd`>B8Sy9o| zV(wqk`lu|e{-xUONqPt@uzv}Fr#;v$o_C|yWV+$-TJCb}M*4j1U9H&fD?tA7Uq8hw zGg-I$#-KLy5)g3RJ_;SqPf@Nry-fmYZES*p6VZkWD({%b>_BMoRdus75^^GRq;Ijs z45M`GXnE@rtm)6Mof%eJ3K01QtEPji(~CV=zY6cy7nw_R?!8CLpPP?6MK6OU39Xlt ztkzF9_rMlZsj7+x#g>v?d$hWo&41*V2K`$C=Z3YOF#?zXisn`1ZnO{(IHwDW(ze(C zP}Ku6-mRO$k9~#)w&xNuHvgBp1Xm!pY3$=|xJCy61v>ey+phu=^uK`-_MlYeN)3SN$UPNdVa13y zu%-;lRo}eBm_)FLm-<^8-!*Mu|Y5ADT= zzTrmA=VorlSKmneK0DF*dWy%$$ZYxYhF3I=#q4iP zn6rScj2WUf%SP;2swgD+-gq54GNa7ul>T=uWZ-0cxc=p!susL_C^DMD-Eb^zduv0B z`Y=jV+Pzs&uw*<2BsI^4gy12#U+wS3S!7zhX;ognnkoKor9PirmjeY}_lTYk&8d&# zH=~aK!DoinK8lRQp~85A9j;&(>FJ+;ufGJ~C&pvRdg@M+J2pQ8TEn+H{Eu@qi#NL; zuN9{8ZOu;LtiGO-S5ox<_Dz%@|K$jH|L}jm9B~uU4x-Q3)aXnF_g;wB&!aZBMtGAy z$8eU}Z>*NFrz&6zcg}K<87>V?Mc(AQ)IsGJCaqQAU)+ov72eMFPYCx=-24vq1LcAz_6MbiZ|fv( z=O)5#kZ=C;tgZ7L%w}6ocFT#xRF_U?^YI$3Y>1S-a;Nv%VXc>PQufw9kW#ULbl9!j zckT&2PZMFcBt0)%Fx8!mF<}nhWn;)0Jz6iu5ti9eAXBdD~imZF3lf2ngi{s>AYSj=%)Xn$^xR>y0n&>KaAxC%fUIzL@rIF9m8hrKE5L=0a%Pn+<%f5YEy#JU9Wq14Pdrl?SPy94-bHe!=$ zD(F`IlLsIr7K8_Cb{MB!>GfEj{2Xl&)TdP&<42H(>G+au6Mv-?BSUh!f88Y ztEEkZq#M{m1qn+vDkKqb%Jz+W$`UKv*UVHTfG5NVD8+I&uG08 zvNb%-6yf#DV({QlT0C0Tv`X}q3K-bvV9X1L{KK}Tf1Nv+ZPQc8G}ZEVt)U73{a)ow zm+yBzkm8`+(8*w1RWA6nBm|9*yU|Q}a{D2S~m8 z9V{gG{uBH1-6v4$$_4gMmF0xWyfI;8q6`DTOsk3UL?+}Urg+@b%GU@p>9iBu3$^qj zR2AnI6SQ+LzN2&>p?epUwv|u;Z3=qa0xy*kFx8Dl({Sst1G`3xLr3y|G0Hf@jQc&*en6%{{CTlxE?6`jj{Q^dz7qs_(>oK(`J4L9{uO6 z#1Fu&i{SM;f*V1>0(itHf=d0BxBU2@cLW&l9>HB!IY&EqzOXH@ssZUGULk60K$mt< z;@Fh8kV+7ndR>DFI}%lD)l;A)i0t;Wq|D;B<*5nZ*K>;A+^L= zbTSLgKEA6FrFGn6)0*#lwtsl9gy7D&r{)uFUOf2GVRXQa$~vfqGf%T0ubKGFe4qAq zZc!mRsIPHInZd=qbeL={Qz4B$kf^7#OppN$rLWyteO**q>O37IWi>!sA(SC1p#`Sg zezE#9j*0VSw##G$Cw8&7q&{58#Ug59adEx?Z869;UGmiAXgF>Y(%#o3R@qx}zfZs5 zcj`T*Xw_MUp->qDI#8f!zQfKsAmk4*%b)#T%m39Td+}|J=S4C3*SXL_mmgDtnSyGi zw2<}>@?6}Dn@O3>l?jz?N1UG*EaX%bK&|Cj6UI78Tf~ALq;v+3iROmf&Zkx1xkS0X z<+FP#aXyFePKW+Zc>0WMHuJ&-H_Mo@0^bmAwV^<@&0tCIfIDS!ukWMQ6}0h60T&(^ z`_enS7!EZpq(#nxf}fL9)(-cp9gOD`kGqLhC&{mDK4uP7+^gd(<6VPVOlf_u{Z3xX zpH~_f8xO))rw3P$!J5R(E7;kgdIc|=84+c#(r5v9@;_BSmYE8860mi?t#LYTFnyj( z->l@42yk}kTOXG1#!Dlt~{ZH3z;z>!jHB&Lic|aotztz;7Vbvvw zNvl%Z+!2)VBsH+PHFxw`c!k5F^y!&|g~}%(zOY)ZNPQ4GQV)}P_acouyGuhI8D*ld zC<;~gmzxm;TDZE@Z22QdC&^ULs~_w0wK{aokU#F)ruccdsZiZ>-?>W^_9kInrJ4b} zbbunI>%?`d%f6$(GJ`tO_LkcfVMd2nGo;sScm3uW=q2hWKF8 zG?&$-v`nVwCFcVpI7Ts?{Ht`Ukgyv7EL^^uYP0Bfrh*Eb73;}SU)n?cl^>03+G1j? zZhbWQnbRQDk))4gh8OAWuvGjjyuQt0%OpXfOy~_EEY;^aWvEHBDe%F2=K0JptKEZQ z5!@h~Tk-~WTXIlG(wifm)t$*`b0hr0sR@q>y2}ad$QrfTK^~R4q*HAcn`R}_dTLSX z^AG4fF}M1<)2Y(juVMAdDcog`72D%C?CSI@j{z;qmh3mHOOxpX7LU zrxpMoY`tMa+pu>?bmUAR--vo*y1nA6EMqY-<(1wnxN%m(7qrdxTI%TVU?oaX`UhKmXaHU2=40)_sM~y=~^mN4UD0<6V^c2@7 zSy8t@GBqMV^J(mb^n!@&gRBK*o5OluIjVCSCjN1O0 z)&rcMC$41y8O!FtEZgtC z_gpZim3;C>_d{b-(Mw^{ePc#V4EqNT;V3cY-zIL6j+7l%8quGZCy=#4tony|hl@v0 zI$nXn*X6y30_*w9R0}SC{J6%_u!?oj3{xN6+BT4Z&$Tt%7ad@poEaZt*euMY6<-oc z!CM9kdg)?8B*7DAhn{Gw?4`87Cm)tG9-Rkj@Qd)uf}?J!IHh(6C-hK$h9Z;sFy%##OSt><|g?*gkDqWz>Wz~M!)9(o(|jggfdrd~|1f@(qX36D%G@Be2MKZZ~F;5Df0Vc2XR8m}xR zGeBvzIE7IiM&}v-Ee~Z>yfezFXFuZhvg!K0^0Su)UWTY2r6HLQ?PV z)ZeKD<wy1S@BV)^dX2veh~6oQkx7cW zkq-oE(4F6B5XTwa%l94NA$W&;UY6Bw+3K_GeY42h?qrA&9zi#+q*HCi2QrW&YYx<; zqBwm^s5xwNCM(L;-7rm5p!sIl(>C*`witY``k6dW_vgZzkoIl~_n+6?*~yOLnIqb| zCEh!Y7`yu9(3@g%`~N&srt?o!4RTa&sSfYC(eg8+uIinm0zNIC5MrsFaUawy?6VI# z-`?j_xM@TfxW5;`{e2#1cw{mT)1Rf`t3oS9Wzx4_8(wozGqFq4l2$c-hkU0X9h*G= zYQ5#HY3FMp(a*kkMT@ng?rr0`r}Va&;FQ$4Rh6H`3ESaC1foOrOt)mU#f221gR85P z=m-~0cPZ%ULGA*?l*&fv_>qr&YErL^=yNDEwStES5Qs{|UKPCWth|#3?$z_M=T$Z< zSO1Q^t8V&P3xVL)EKsv`$_kKEZd_wF*=KClf=%0BCSLbeNqG25riu(7NANtBILh0i z4t43D!QScNXQaUMzpSJf;3sF$1)K+=R58SZCQDpd=9QtEhVsKmQ^qO6@(}XYv;ZOu ztd?`?o!HVFkR0OLyg$|IK9#BJ(t6tRBi~bNrB%X?$@`WBy{(nR?vs-irN&&p8E!tR zjHn1_SBnKLn7!OSe=6=^aQ{$T%e!lCXsO`f(Nr4~g&;pA)D;azCE+;lS!GteFqi(EAc&Tbxkl+$I35oEwKJ?etv&Gdwi?mh z^~*Bn+56ni(T<2g3yew7tH5MMN;OnOzxZvM{k|Ej$Vt!ddqJ)kv2Oic&zgfO9URV8ecr+6pJ0^mT zp)JJ!`lNh``TfXFzJ{1O#|2glkBnONwJQ{UuAYeXJzJqCTAD93P;WDBvIF6opweo= zu9)%kuId}M-*rkDJ_8=Y%se1w>I9DbNDR5jELn1Y2&S!Ms1!adslGV%G zB?A$xA%(HRQI{bpV`%>Eqh@^*l#cj*Lu zc+R} z8hJFZ_XggPSBTIVe0667wlU@6W`P+{Ht}l=&)huS^(D;zzpTvNkZlu6TX;Q6^otHN zQNbe(+ZYmFAP||7_z?sVTj7H6pWxvj!iZkfqdGKzf*pAJE>7(N34}3!^ojND?33=> z!zn&YIE6{mLOkA8k2w*xwf(P6ARG}mG(`Lef|la|Q8`RpyLL@YM+X?ibKBsj?VFNE opYdXKSQsyG`TD=-_{iPVdJ1K!(6aXm7)0RgDU*}=dYA9~Klo8l{r~^~ literal 16236 zcmch;c|4Tu`#z4IO30SALPBH>MaZ6nvSiCTWl33*J;u^v$&$5E$TB6%2-(dr=BY?o z$C7xe&gD3d^SEP8jrET*9cQAW zqdRJF_3{lmI(icLo6L9!ymPZE^$Q)HrlP@RZHtE^OWCndyYCtGxV=}sQL^SJ?Q9;V zwftkYOlj$GFV*B;c$#sWZS*^?uR6C4nXelaY)6&$zbJjdaOQ~d>+hmnPN~;)7!N(o z&Nn{!x~1~iXr*6y&CN>8q`LlqFRx9dv)!hlx8OwU&B|HAG$q@5(xz#Rt>7MQVSi;s zvZpz0IqZy0oy{}WcGGxM`Lx)J4h$Lir2_BFC+88Lo<*Tll zLNI*mC`WB1t90RzvSkF%G9zv1a5;fA-{#Q`rxv$rtcX|1U2odSg>#35uCY~Pj|K34 zTbaVdvK)OTDx%=)St_a=s>^vHH@Ue)2zHkM`F(Yh16Li+k-suB`{9y#>diGvgMsea zty;0Ias#JX&3irKj@!?gPlR*5QktKiPwDD1CbAtyT(=m-@&x+UZPLkKW)UX+&7?lE zlh1+evOh*?NtCnm+)8ID-HYDk zgRZ%8s*vL3Mk6JyCroGCPDr+J&<g+Mcu`&X(K-#d(Wiw2_Wi}46%z^l&d=+x z{5|}_{*lHk^g5d**u$U2*FQ9<`dmLU7VaA_w*FX2eMh5uH_QKu{}o07_%Yo_*j5p5 zxdID3LVA$^x$5Owb;6Wz*zd_~2rPW6vdMyd}+>bl~huS}Yv8jOekLFg%UKhht)h&)KMw+H;hx?JpVqsK2uM4tYePHL zlQ!oW&d0)V?`o?{kwd9_XLki@d92$I9~5Dfu!lD5a3GusYoXkllt<(2jT{6DZ7E^~ zjl1tG(3?HfdGho~hQbg9E;o{2S7(hvq3ShH{>cmTczm9~Ghru2I7JuV(|Z5dbV2#E zRf0$1Df7Dg3FtMqg^(?W&f-BUgV6R1oTD3Qi@W~W#Sn?74xQ5Y+uCgLKW}pB1FPje z!^UWOrr@XR;(eT9eBIP;uAP&*9rT1VtU(J2VIP+3cOx!{3aMD6jTGjjqr`T;RtWsJ zIi>v3SC$^6@vUJbzQ|#emj>DIO&0NXxmE@c_P#s+_myW!DMNWYs9?8ywaI#k%l*T_ zwU-z5#X?qGRYg=qz&XbXYvHHcu$n`#Hf*^jq=UP?<&rv;*|vWNnW^8MZEamMRY*}~ z`TT?KAe%*2Q7v3%k%ff?%etSz3N_Onyot=Y^uhWxX0QAOwr?!N$Dp~)3RR+`h`f=s zM<7h0ikhG+|?Wso|_C7w5ulqVz#Zd97L&Sn`XXCJ>xNtfZ=RXutg*lw#W#vN zc*bmywM}1O!$fHTTE=kxDg4ZTY!k_<$2>0XlvG`sdbtz1!EoN()G|@Ef5J3od(iWY<#)L$D;zcvz0%7?yhTyUd5 zUemec$&CH*Z*WIqr4)7L-VIq=BBO;yHWA*r<;;EmiB1rrf}#~A4ClQkoL-_}NRw*B z$^Xj}`hRQ?r)YR0Fw6MDEU6Hl15jgUg?{qZHZ0;2Vgs+Eeb!SOQa3ClXc`>)o?vnot=tYKxnEG#bZmHt6NqdDfeg}LUCmz3A1gM5-v_p-YSd7l;~b>V!y ze>LUw^;Yf5DLm-r&T(3%jNX`69kC(>MCnEhN=^8^3!m$R{&E~#oo=1Q_q-(%D#JB# z9+Wref_W6p%U4GOU?{Y2lwc%g+TY&V7C-x7(qs-m46_%0n(`@xO#G)Ck%E803%wiD z#DB$6JGy1lCteEnILR%bEa;ExTdMdGj&d$sdW_{z1(pWeSDX8^if zb2xbS)t>Wil8XNvh{5PF4%B2z8KP~;<_AJjIfdaoEn{(n#I}?f-vORtqh_A_;|p*i z+M9nj3wi~Av3zu;G_J6X+p_G-0oczXkXOKY6Ax7Ak5~a*Wp78H`C#;fC-Qdt(b0jS z+AFkA1m&kM#3gTj-&-w>3G?{AfzBLh+o|^CwT_{{r{Q-MRMZh)(X%PnET~?slxA zPz`ThF%q#7y!d^Woio9`oX;%9XhLccP8#}sWFG#WBI6HSF)jhvQcmOwp`8c4?zOU77m)*`|>%%(Lt{ZR187@I9D-F zwckLw>m)77opXE-iZA6pNFsLh$8%X)4gjbd$$IeF62u9FabN+<#!10F58*(t+3EZm z+c_&s1DJSxz95_;h_OV*HSLL8(;x3f_?k-xBJ6{8-(B2mO@^F`W)LG;73bKA`IXpx z!`fgKk`za1epXsc2W;gQ~`#+oHX|YVk5?oW*svbl-3TCHp)|Vh&i;=5B#(kds zKC(#1+x%qv;`@5w+|8@4)%#}`~s!ha7d78a$>wPK`5GmN$b$6v- z*qT1`Zl6fn8rk_78dr`<%W>96>H~Y}<6#OjIo{g#_jp>)-sU}BhrLE*5URl_Z6e30 zwF`vTdhR-K)4A=$MApHz0|$zKN-VrxQ(0f1U$!it4$UD{Y6k7BsBun`1hV>}CwA>v z_ozQ+?5zbW*mf2@d<*9G<`c3pZz(m%W%EI9)|mBh5AqQTC_?Z(Ux+o)qr&Ss3t9m{ znOb|V$&#k9@Kwv%A>tuA7vM`f=8w)yJbYV1*IAEfyIwGK!a9C+OTcLY8udyOW05eq zTgX>q5;q){is7dF%LKn3KlXVke7r3hJ5JizFtoP*pwRt~-#^D7?CCACB7@$4lM0;P zTx;2CadonfXrdzRWVxo3HlGDoUcycK#jpw2%byqcjgXU#iw-(6*Favg+`|&wT`J+C6SeX{R*=4-bSmF;JvdY6fyK9Pn%rynXC#Y#vPPb=XCQ0&97rfB;+x+{X|Xg zxQ(Uk#=wEJy{i7enkZ{qn-JLMn~SZVRHma24}WH*Th!IZyYNoyaOV|c+GN+2z4?Ra zGj)ceD#&N$iA`f@ucf*2IG^Rk7hkiXLxbv4a~bpsfU{r$H9?oyvBxd5zO)7~Tdkuq zsKt7nd(SOz&Y7n18O_CAZA$@P*nWzqMEG*&owcyMZ>ai~c&wQFm*X0tCW1?GEavyM z$50cfXzWw1cm43h>e5OBfVNd?MV1Fm?g=c2$T&%(KSIO9+_K9@Q*4Q_6S~t;wnS}q zRIAbjxuRLtrZkk;1CXkWNMa2IVRp?o{omOw8qL}G)Fxs0>9$Vky0HrEp#GB}3j7E2 z=LVGr%Hd&b>Zwjx-8r+)vfG0nxic0-S<2+QtEie_@wbYHn8!`K zfbRM?9@0s5i7%B|+rWhldvE*07R>`djpw_i`Z-KV^ez@4K$t(e13Z2ZW%|CX>7c z7pRQYhd>RL$nIxrS_Da0rKt}xINWHI0CAY+X;^=i3B@flKEjG$roOqw1=GfxN~ef; zj32z4*ffRlsHybr-d=AdVuvq97|}&Kp5N#jT=XMY#4&d{_5glrB;s9hn!3diw=G{F|%$uXCx3X0P*Q)70D zq}KF%NnR*{e$SC|7)o*scsaC8DSJh!6aEbre2CiCg{82IAXNRYqz&*FGn`+y9)8!f zR(s7-1=4%7q@*3$+_X4Xl?kDdggSDAdOqY=7BT&3LFcd8+nBuWy|cyMx!%A<)>WzX z-VtmUFii_S4r~)=0xacrTVl?f$AD@q-+VGe#oKP}qVe}~*fhK9Yi_G6f}>$|C5W+G zlRgoUAh#F|WM{xidnqp=ZQ_haxB%TyOd?j;w0RGgY;b2B*Sa^e!>m*t|lL65znm~n?S0!WB=HjnTC zmP$5#%cij_jUWLGb`cRnMpV zpc0`A>!(#-4!!K?i(ddt0j=dRm;tylI{RN&{ogPI~fuPGiirpM^s?O;Fb~l|?<^P>7P$ zCXVT0WY?3Nc}MKAuw}D`euC@~WBTLeVfYxNn=J_x*yC((-A*g{)*@*FV*Dj6HOpxD>H%J%@ez=r-TUydRHgjX@d}<-QvV!WB(- zeJ$5pg{Ur6$3Z~wkREbU$sIa@&2ga!rGXz z^rng)vDkdUp=fQ}|4|8zKbWjR4Ektv2f#K10)j$MRb1HfOe1GxGq)F_LrmgV7lfMY zd*m@?8s#p_v&XA$wLF{;3JK2!M9r;{C=II%7bRPI01^DWfo63g!+DJoN#%R3g@J+3 z%+T2E(n%$A>o3;>q-!9Kn%K`+g(_3kiyMJn2>{x<3tKAwsm9=Jv$zXjk2vdtc;iWp6RHf`Pd zI@A)G%IE-(9ws2rIIx^t!Gb?5J%?(x*)3I(3uE`Cw5oO$)TJRCIwpN#I&#TtuOa1B zXZPCzG?vsq(lNH^^Iv2NFB^&v7ER7G5s}!S#}} z?ydXDnlbu<*I!d=fIH|<7B$wI0)VlLX#aw0c0_;32WF?|2)U0AbmPD2jy!y?yhxRj zRr}%Xh$tiA|1FVlOUMWAyo@%4@6yuJu2)i@1Njl};0{@aHx9Y5dNvCp9{of>mr-40 zm2}U;1rPcIZIf>sh@cDYrb{fK8=Mnw_cODXKdD}n1-k` zo0Tat1udZGm5P@ey>2jL|1!kSZ?AXC#rxo->xObs`W?<*`;_=_?U$3PYY`fkEn$8< zdBaNHy`VmZu%fX$WGSPT`g_q13kz$Ru1*^{XMPV_=tHAnDS+yl2fW`<(CgX?1>fn| zhW*5lb@W0tn?)to+z)4&C>*3o)YiwF-utIUcibzL#Cc=}#jcikjQ!37h0pVTMwT&s zwNSCHO6n>$zj=la43DyWi+$f&H8e0`U+Y<_Oj^WcNxy$@UkcIyNXaDkV2~RTd13f@ zFAqWH*`eAmZ!YqXww%jDMkTZU1gMmgnY~&3ky7Qe|AMGvpi-C6!EsW8y0`g3Wu1{K zwqppggsC+ihvR(TIFyHOzJ^Sa2K(cJ0r&87Ex7^5b3p-2jofB1=!I+wDH(d-5_y-f z#yI}%Fw|e>;XhD`HvMVjZ5hrl!(rE0T104Hv9K);s;g(hUY>M_#-i`O=bk3Ek2^sE ztGqxs>fl`!W{8ec#7+|E=7OOgpWhZg$y#N7zI(#yyVC4d!#u+GMEM8ayy1rWB8KyD z)1pe^!T*7qq7JX`gh)n9`vl~2RZTCG?^-{XwDm)z@x9T!Ed;o?zq#HpVLg!}P1a$Q z9fMLyYb*c2XRA2M*kR0ANX4rYj#p!MzVC5QD-Xs6Q?AJ`Z}G|t?pC~#5}K}7rZGFW z;worme)G$wIN%^{e2WYc#MZ_>uV_*-`II9pI~tmVpJH}$zU*wjOJ9-Wn4WLp&Mww9 zjnoKXkB3ALbC9>@mzH2Oye@u&FD%4qfVXxA1G+lYU8GHe6iBN-_?T*esw z+#6^EJ*|kDIj5S{%_n;We0yQBuL4Yt>Hcto<7A)tv}SAQg_198wT=anvjGOZdaXxvWIrX(iEyo&#dqC5@LlJ+cf-RnulwnK~Ovq|0gXx z#kRgF*H0M5$bRYMZM*7e`g-q57YV;OUR?=JKQssKCU$=^&3NBYxGta8A3!l<{XO;Y z!Br$j`S4X4^$3kIUf-e6%h9?^XW^4U-12|z=IDBt=46!)m(BI*SNX0#}ni%ZNiVgT#I$a7-2DgLE zPS`?nIqEdt>&ZQHW>uHqi~dN@BuL{c%!@G?z?&hmZ_8=*ff+W-htke38T9hU9jhuy z?f&V4kcw-mOP#5n9){oSz5L8SFR2lp3EvtWxonmZc`TCk9BVO{vLqsWE5g`%m5dLO zTTENrctX`IAI%+VbzBz@H$FiR$RZcep1;cxA-;G>wG0G#Q`HB`Qyn&Cy;hokcu%~JGpYWSQ8{O)-y|=KC{9J@Bgyw^viS<5_$k^YBDCr zhKrX*qyuu|6b*DyO6N%pAHyHC&nk0YyzO&Qhgg&52=;&=gs?O zZdvxz`Rexw%)Kj0Mw8|yEZXNbz!kJ;c=J|*Ny^}A~N4gwzy8|Z++(v1H{eXfqm zHkFu32Spmdchb0Pv+4cx-G239X@F*$*|)Uu#ij0+@na*M^8ZiXfMXwSzX99=eMOGFNSpFml6 zpuC6jCmPmSRdzO?x=uRI>-ji4zJ$SWR!M3oH|X#r1aUp_q>9cZ)v|WNb)O?R?0EbS z(rsxK&cR$cgt)F7Oijwa{Q->f>fdrdYz>oAinh3ulM=G716Ur=3v6WDO1xB8{f>e5 zgIN(cMjdz(Po9q@3_4A5b+9bM^SO_|?M*`-u@)7PzH?N zrm)AAk!#6O&zd%i8W8sUJbC=qEEg`^i3;yyL7QarJxML|?RKEUrEc>vX%=$Qs(m#A<;3)w}9cnsveup1k3 zfq{X&JggR3A2cVg6E|{v45&k_{QHRixTQwz4d<)bP3dRPS2?KwwSoY=N+R)~`mQ|} zF6!TybaDNHt%Yq!qc_ikYr@{s3OKEWes>I&9ni@ zeTLTm%0~{=nX$8A-n@8YGT81GGxplW=;Q8{HstK8_kr^w(R5%o2YTkU2#|4t;D|OA zTCDMYWjKEu0^A`3C6U!3YgMAoUQ(da;vf-KN${@$blXJ zKX{l?w4zQZ5U>t?2;YNMr0C}0UM5+$ z$^^U)9Uo@>=uWQ_{&EA~ASAi7H|sb}P~d{6bA2MIRW<8QaRjX50XxgGTj`GIb6M}? zFffcJ&6K8BZvFzSA#if1=!WK#_?rAluO1lrBKDZ)u#o!Q?q0SFtDA=b8z296zwf=6 zG#pE7_COLyOiWCrdhTOh^UZF~(clJ_Qp`2Za>yxN!meI~po*>BT*mI-Uy@gc1YA?z z(Z*iCok}{c-3Z`buNx-)q0%*-UbDdit=(6D)|eQDe&z;xn}urE1<*BHS>0@L>YZd<669Y3ne%xDrDoJoFK1Wd~oZzV5f1mO);oAekbjU@G+T{pF7=z-4EnYj{|8JDLpn0WuKhIPn|E5^%dGl@clRv+zb zLP>wee+>T|w)gsdKg6bOQfR+!L>mX1Y4DS86&X0GdG~!Le(Rh}OK~RXoz4cryW8JF zbHvvC7Dq>w9Od5?)?UBIvNVy*O{OtW!(e*F`n%deY&w9#Nyy{H&;a4R^iI`V{j%_8 z^25Z5PX+FOvz6{NhC#C=uk<~9%aWtYRnK?BE#Hl)^Zm^~!%JYa(8!#w{+JQ;I(w%+>GymSob@&% zT4H|ldZ*Aln&X_*ebh_II0g;|v#qh@^O%(fmH=}`=9ryjslINQeC1mog(G- zuHnnJ*S{;eY<+D%t?PEV1Bj$Z9s~~!INRg1XUubb2h_7VdwDwsbXRxOL9`!!6zBG# zLGS0$!2%aY=wp^TOy)>;&Mq$@WhzLpUkV3^{I-lf;$UkH>xMwx#RAQRX}dylq%!Q? zWQY&I2(_sg5s-Y6#ePdxy#Aw0l;Dz-D-J@Jmvo5@cKcdqpWshV$qBllyPOqqPI65n zHniTxse#zv;>8Juo~t!`r=LErn{WNTW**WhmF`wWqTs+l7wt-)T$H$2qM`D*K&o57w4=~SLm}NC1y6rNk+vI6%Tz5XtPz^e#t^czy zh`==b$CvMyx6>}$<-kBMsxUMz`LvLAK`$glgg|m$ARKDYWbc!BbFp(?7Mk7!a1?D$ zOb@tpVSmfJXPV~QdW1KRH(;I7#lilq45|Su@v1%7VHw|Bfi|FO02d3(g>4NGTvon@ z05xM_(gVxFrgt-Yda>S(I@a~^7bZkkkG?FNV1K{!8a59>DDG<4j5Iw{9Q&>qdkfk&Q7gx%+07bFkEeLsgrl7rs z7ba(S_q}ncn>4wMw=+R%>Q5c)XgCnD4_H znh24HyaF}4l<4)G(7L{f*_gNZvbPbY&>L}Ty#@$iL$|!s`!VlgmSZMuZ=rq8k>($W zWZQ`$w;O;iR-^)r;quFK3$-Dq2Jg?jlU23g7))+0O+cjUK@$+caLWTh_HV97ByXLV zta90LHF%G_`g0U#G=r5ISSvnoJ3LFPqidGWeX z$dyB^DhebG{`=^@5crPFazo$Ei1rmq1eQOza(6VO z-a)3;=sw2yN@Am=h`zO<@3KSj;siP8s@>9-WS5Oi@D=Hb5!lKd73zPaWzvaC0hOV_qZPVDPF+xEokWsa@eJO2OCW;HdG-Kcv>EB)Hr zYQ}>BvrOt=H>oa;*5vCxWE4tGBY9=CU956dl0SsnuDWIg?Q-^zgMGZNnT#K47j0T# z;74rdi#nb?Q-3-=Koskd<`db%5i_&jVo>tYp=+7WT7EI6Bb#>z(;J$M+_8fFceI zRYs2kZI|!o+8Vkq$-k@9`o;bFkj;VR7baFQg%Q!4X4igL$&;aGW9Of06 z5od$JkFuZ@Xy{O|ihZQHFCt;pP7s!o>Sm;pnMq^&VhJF!1YMe{KVxJg>WKzAwa)J61GHJi4!_g047UUDE5mN^ObR*Z9!WF?=H)-H@i9 z_r}voIvwMjRRaa*Qum?srg0WP){d-^pKk-1KAhI$?LWlnEuJs>{s3Q06W0jQJ_w1q(E>KHjoBtuF*=g*Ap;VG5y~;Qf<%97!2HE)Q!B|MZr#JY6$rl z+8CAVt28>~HUNDgG?2idN~d*7CB!s=AP2%QCiCtw@O^L|8SgX$Ou$8fRj<^?c?z3Z zHm|vhsnEKiHM5sa16CHKphDq$yd!_@e+F6{I=XXb z_aATq@PH0T?7+*v4?X=qfA7@|Opzf1;d|3UlWdQct_;&?>6c4IAipOPY-iFoY?{#w z;2EqP!xD?`0H)CpPqUip^q(O>cp|v|#ktI)&iG5fwK)%Zu%TD1wk57*xie8%da$2s z7LJycMd6Dq;r9bdHnE|W?h+^s>#r)^vk?yE<#mp&!Xz*h5!78t%q47P)WUtJ(PB6x z5-GwqKUu1A!SJ@hrAZH!a*?&0U4sRqUw5zE^KW}2xlV30paV8CzhZk$NqCln{L1@K zHH+VjqXcCTl1tes+WEz(PTecbz#*zM+oehZ(YyE)PdYrR?_9o*xw7dI+*I)at#Z`? zT9B>aEn2xW1)j~I-A;i-gjwQvQ9KrRi0bPBT?zJKo%ItE2`O;GUU$RA?z2;&gV~e- zN5bJa_K+rJ$+v?OnSD968`z}T#wmL6q}sW6enXb`&7g6EX9~yQv!POR@q7zH9$(eh z3O-8_f($QCTIdi@Bst|T`v{axj_VFN_HK@b90WgZYQQ|(n2!@{*YR7`lf5K$F6KV8 zZ|eB3K&Ek{oHSAYa=2e4EAWz?I#XE;c1d@srK5(|(etIrj~t!NlRtB}k}363jwD(^ z2RrG@c&;|%2??D`c5IrOyA+SLE#aNzzUeNXMKKsY2 zZ3~dNg|h8```u8jmE{8uYXaZejE#Xxoi8`yoEs%%-H%x&?u4_E{l>OlaBWL47xtSX z;zy`Vc?tuLM`M%W;{F2v+l}|1-TPmz{r~H|i_9#alluxyjP&(mPHKeS92ps*%ll*& z_ITl##bh3tk!%>rN}n8TZT<43THrNQL&G~j&ydOF|8Y@%BK9MWL?T^vb1U4vc z-7S&2q`dCS_#e*wY7ZvD4rsy6E`7Jg3Ni|W%KEqb_;HHvk@$T06OCh_s2>K0LQ^UO z`KT7q9MjY7i<)WDcTg69iD{o*Q!g~m?#14t*V5>v{uOpcCx(eDg!)sm7}4fhb}8%@ zr(;8SML4=;{n(!4iFkq$^!mZ#mn~@l*iVB^o$9y-s&{W*(_QXm%LHvz->- zGjW3Qv|*#$Ma?L#`|w|GcQaGUPdtR}VRd^x*ia;Eny|m$B#(RqCD&=}qHe-A zq{9ic3fA4cGNc|BFZar!s^rPkse1hiE_wh%)|^ z4-Ld*zi7230V7{AyR?~I^V|}^5-bqnI8ekFUT^OtmlV8OM{9or*LkvSAUwtBrBWZ}evCWA@pEU&$o zTf4htd+`2z(N-m|)VDlI$gHhm#1gR#1XF5Gzw06Oa z-+*OU&6hws0mm>0pZK3lN1f;HC^)lps;X*lUj6%F^i4F-ry?C0e4EmjJQ~Z@q1nXu zWv_zANf+i$+^Mb6otPBad~W3PxvgTZRjlHtqQcGuXUbIll7{IYdo)5m4+{#4u)2x) zQ#Rj?7OtpUZG}q|lrBNij>q0y4ApSU7uyS)qEcuP^n)4B3z%(|5~z!mGh1w;pW}qH zrJr^%&aMN;nNQ6`G@1%P;Vq93MiK?*B@~!99`>1CQ+ha{v-xoC2j}cCqk6Y<%b0T% zG3Mu}s`vquz32M5CZz%rKbu}9T%g~7q#zeer>*tzBMY_N0_*&OjHNPdK zX8B5%2`vW{T*jwiB)TRzc&hGp4I^>w=fyQ<93!#f7iV%dGk>8c^4XBQ+tB*IE@S}0 z$RvBH5>kHl(doGc<|F7I>Q+C(B@m^Gxr25X590CUo)%cfvRr8SX<(Cc7j$wLy*v_w z9o%rPkZUdhjmRhbOJi$tT7K^h1-swWpLWUeHnYECaj@Kg+7$1#wp%;sCmUDwt#nAw zXPeg!DOlCi^a?&(HBqZHK>GA5zPnfCFe{w*I4X!SY!UkYmCO^p>Hs1Z-CFico>xwl z$V}T#mn~jhUETb~%{v|2?yZ4&E$c)*pB(0exS`YO%FLO0Saf)}?vV9O5!`g$RODf!z_0Yh?p* z+8movxlv_gMu1I)@j5Y|OcVZUJ4yQoU6Ga^KFzc)#}{bQb*!;&G()%=+_^Yd_&I++ zse;T)Yv7(-2JKq3F*Dt-6rwm@kwX<=i~U3lKR`#zWqL;(sG0)2boETvRLEeu@xvSl zO9Jv(!DE5J^pBXoD;H_#taH=iA~b<%=)py4TB?3*RKQ9%-(kOg1mSHOc9|>>zUK_r ze`_YsqI>wh1Q!Dxwn#&J{U{=dwgd*z^~ chart; + private static final double[][][] GLOBE = Projection.hemisphere(0.01); public static final void main(String[] args) { @@ -79,12 +87,12 @@ public class MapOptimizer extends Application { new NumberAxis("Shape distortion", 0, .6, 0.1)); chart.setCreateSymbols(true); chart.setAxisSortingPolicy(SortingPolicy.NONE); - double[][][] globe = Projection.globe(0.01); + PrintStream log = new PrintStream(new File("output/parameters.txt")); - chart.getData().add(analyzeAll(globe, EXISTING_PROJECTIONS)); + chart.getData().add(analyzeAll(EXISTING_PROJECTIONS)); for (Projection p: PROJECTIONS_TO_OPTIMIZE) - chart.getData().add(optimizeFamily(p, globe, log)); + chart.getData().add(optimiseFamily(p, log)); System.out.println("Total time elapsed: " + (System.currentTimeMillis() - startTime) / 60000. + "m"); @@ -99,49 +107,107 @@ public class MapOptimizer extends Application { } - private static Series analyzeAll(double[][][] points, Projection... projs) { //analyze and plot the specified preexisting map projections. + private static Series analyzeAll(Projection... projs) { //analyze and plot the specified preexisting map projections. System.out.println("Analyzing " + Arrays.toString(projs)); Series output = new Series(); //These projections must not be parametrized output.setName("Basic Projections"); for (Projection proj : projs) if (!proj.isParametrized()) - output.getData().add(plotDistortion(points, proj, new double[0])); + output.getData().add(plotDistortion(proj, new double[0])); return output; } - private static Series optimizeFamily( - Projection proj, double[][][] points, PrintStream log) { //optimize and plot some maps of a given family + private static Data plotDistortion(Projection proj, double[] params) { + double[] distortion = proj.avgDistortion(GLOBE, proj.getDefaultParameters()); + return new Data(distortion[0], distortion[1]); + } + + + private static final double weighDistortion(Projection proj, double[] params, double weight) { + double areaDist = 0, anglDist = 0; + for (int i = 0; i < NUM_SAMPLE; i ++) { + double[] mcParams = new double[params.length]; + for (int j = 0; j < mcParams.length; j ++) + mcParams[j] = params[j];// + (Math.random()-.5)*1e-3; //is it weird that I'm introducing stochasticity into this? + double[] distortions = proj.avgDistortion(GLOBE, mcParams); + areaDist += distortions[0]; + anglDist += distortions[1]; + } + return areaDist/NUM_SAMPLE*weight + anglDist/NUM_SAMPLE*(1-weight); + } + + + private static Series optimiseFamily( + Projection proj, PrintStream log) { //optimize and plot some maps of a given family System.out.println("Optimizing " + proj.getName()); - final double[][] currentBest = new double[WEIGHTS.length][3 + proj.getNumParameters()]; //the 0-3 cols are the min distortions for each weight, the - for (int k = 0; k < WEIGHTS.length; k++) //other cols are the values of k and n that caused that - currentBest[k][0] = Integer.MAX_VALUE; - final double[][] bounds = proj.getParameterValues(); - final double[] params = new double[proj.getNumParameters()]; + final double[][] best = new double[WEIGHTS.length][proj.getNumParameters()]; + + for (int k = 0; k < WEIGHTS.length; k ++) { + final double weighFactor = WEIGHTS[k]; + double[] currentBest = bruteForceMinimise( + (params) -> weighDistortion(proj, params, weighFactor), + proj.getParameterValues()); + best[k] = bfgsMinimise( + (params) -> weighDistortion(proj, params, weighFactor), + currentBest); + } + + final Series output = new Series(); + output.setName(proj.getName()); + + log.println("We got the best " + proj.getName() + " projections using:"); //now log it + for (double[] bestForWeight : best) { //for each weight + log.print("\t"); + + for (int i = 0; i < proj.getNumParameters(); i++) + log.print("t" + i + "=" + bestForWeight[i] + "; "); //print the parameters used + + double[] distortion = proj.avgDistortion(GLOBE, bestForWeight); + log.println("\t(" + distortion[0] + ", " + distortion[1] + ")"); //print the resulting distortion + + output.getData().add(new Data(distortion[0], distortion[1])); //plot it + } + log.println(); + return output; + } + + + /** + * Returns the parameters that minimise the function, based on a simple brute-force + * parameter sweep. + * @param func The function to minimise + * @param bounds Parameter limits for each argument + * @param numTries The approximate number of samples to take + * @return An array containing the best input to func that it found + */ + private static double[] bruteForceMinimise(Function func, double[][] bounds) { + System.out.println("BF = ["); + final double[] params = new double[bounds.length]; for (int i = 0; i < params.length; i++) params[i] = bounds[i][0]; // initialize params + double bestValue = Double.POSITIVE_INFINITY; + double[] bestParams = new double[params.length]; - bruteForceLoop: - while (true) { // start with brute force - double[] distortions = proj.avgDistortion(points, params); - System.out.println(Arrays.toString(params) + ": " + Arrays.toString(distortions)); - for (int k = 0; k < WEIGHTS.length; k++) { - final double avgDist = weighDistortion(WEIGHTS[k], distortions); - if (avgDist < currentBest[k][0]) { - currentBest[k][0] = avgDist; - currentBest[k][1] = distortions[0]; - currentBest[k][2] = distortions[1]; - System.arraycopy(params, 0, currentBest[k], 3, params.length); - } + while (true) { // run until you've exhausted the parameter space + double avgDist = func.apply(params); + if (avgDist < bestValue) { + bestValue = avgDist; + bestParams = params.clone(); } + for (int i = 0; i < params.length; i ++) + System.out.print(params[i]+", "); + System.out.println(avgDist+";"); int i; for (i = 0; i <= params.length; i++) { // iterate the parameters - if (i == params.length) - break bruteForceLoop; // if you made it through all the parameters without breaking, you're done! + if (i == params.length) { + System.out.println("];"); + return bestParams; // if you made it through all the parameters without breaking, you're done! + } final double step = (bounds[i][1] - bounds[i][0]) / Math.floor(Math.pow(NUM_BRUTE_FORCE, 1./params.length)); @@ -153,60 +219,119 @@ public class MapOptimizer extends Application { } } } + } + + + private static double[] bfgsMinimise(Function arrFunction, double[] x0) { //The Broyden-Fletcher-Goldfarb-Shanno algorithm + System.out.println("BFGS = ["); + final int n = x0.length; + final Matrix I = Matrix.identity(n); + final Function func = (vec) -> arrFunction.apply(vec.asArray()); - final double h = 1e-7; - for (int k = 0; k < WEIGHTS.length; k++) { // now do gradient descent - System.arraycopy(currentBest[k], 3, params, 0, params.length); - System.out.println("Starting gradient descent with weight " + WEIGHTS[k] + " and initial parameters " - + Arrays.toString(params)); - double fr0 = currentBest[k][0]; - double[] frd = new double[params.length]; - for (int i = 0; i < NUM_DESCENT; i++) { - System.out.println(Arrays.toString(params) + " -> " + fr0); - for (int j = 0; j < params.length; j++) { - params[j] += h; - frd[j] = weighDistortion(WEIGHTS[k], proj.avgDistortion(points, params)); // calculate the distortion nearby - params[j] -= h; - } - for (int j = 0; j < params.length; j++) - params[j] -= (frd[j] - fr0)/h * Math.pow(bounds[j][1]-bounds[j][0],2) * DEL_X; // use that to approximate the gradient and go in the other direction - - final double[] distsHere = proj.avgDistortion(points, params); - fr0 = weighDistortion(WEIGHTS[k], distsHere); // calculate the distortion here - - if (fr0 <= currentBest[k][0]) { // make sure we are still descending - currentBest[k][0] = fr0; // and save the current datum - System.arraycopy(distsHere, 0, currentBest[k], 1, 2); - System.arraycopy(params, 0, currentBest[k], 3, params.length); - } - else - break; + Vector xk = new Vector(x0); //initial variable values + double fxk = func.apply(xk); +// System.out.println(hessian(func, xk, fxk)); +// System.out.println(grad(func, xk, fxk)); +// Matrix H = hessian(func, xk, fxk); + Matrix Binv = hessian(func, xk, fxk).inverse(); +// Matrix Binv = Matrix.identity(n); //initial approximate Hessian inverse + Vector gradFxk = grad(func, xk, fxk); //function at current location + +// System.out.println("anetauson!"); +// System.out.print("A"+NUM_SAMPLE+" = ["); +// for (double d = 0; d <= 1e-1; d += 5e-4) { +// xk = new Vector(17.8, -17.8); +// Vector dx = Vector.unit(0, n).times(d); +// double f = func.apply(xk.minus(new Vector(-1,1)).plus(dx)); +// System.out.println(xk.minus(new Vector(-1,1)).plus(dx).getElement(0)+", "+f+";"); +// } +// System.out.println("];"); + + for (int k = 0; k < NUM_BFGS_ITERATE; k ++) { //(I'm not sure how to test for convergence here, so I'm just running a set number of iterations) + Vector pk = Vector.fromMatrix(Binv.times(gradFxk)); //apply Newton's method for initial step direction + pk = pk.times(-Math.signum(pk.dot(gradFxk))); //but make sure it points downhill + + double alfk = BACKTRACK_ALF0; //perform a backtracking line search to find the best alpha + double fxkp1 = func.apply(xk.plus(pk.times(alfk))); + while ((!Double.isFinite(fxkp1) || fxkp1 > fxk + alfk*pk.dot(gradFxk)*GOLDSTEIN_C)) { + if (alfk <= 1e-5) + return xk.asArray(); //a simple way to check for convergence: if xk gets ridiculously small, we're done here. +// for (double d: xk.plus(pk.times(alfk)).asArray()) +// System.out.print(d+", "); +// System.out.println(fxkp1+";"); + alfk *= BACKTRACK_TAU; + fxkp1 = func.apply(xk.plus(pk.times(alfk))); + } + + Vector sk = pk.times(alfk); //iterate + Vector xkp1 = xk.plus(sk); + + Vector gradFxkp1 = grad(func, xkp1, fxkp1); //compute new gradient + Vector yk = gradFxkp1.minus(gradFxk); //and gradient change + + Matrix a = I.minus(sk.times(yk.T()).times(1/yk.dot(sk))); + Matrix b = sk.times(sk.T()).times(1/yk.dot(sk)); +// Binv = hessian(func, xkp1, fxkp1).inverse(); + Binv = a.times(Binv).times(a.T()).plus(b); //update Binv + + xk = xkp1; + fxk = fxkp1; + gradFxk = gradFxkp1; //and save the gradient + } + System.out.println("];"); + return xk.asArray(); + } + + + private static Vector grad(Function f, Vector x, double fx) { + final int n = x.getLength(); + Vector gradF = new Vector(n); //compute the gradient + + for (int i = 0; i < n; i ++) { + Vector xph = x.plus(Vector.unit(i,n).times(DEL_X)); + double fxph = f.apply(xph); + gradF.setElement(i, (fxph-fx)/DEL_X); + } + for (double d: x.asArray()) + System.out.print(d+", "); +// System.out.println(gradF.getElement(0)+", "+gradF.getElement(1)+";"); + System.out.println(fx+";"); + return gradF; + } + + + private static Matrix hessian(Function f, Vector x, double fx) { + final int n = x.getLength(); + + double[] values = new double[(int)Math.pow(3, n-1)*2+1]; //points in array placed with ternary coordinates + values[0] = fx; + for (int i = 0; i < n; i ++) { //for each primary dimension + for (int j = i; j < n; j ++) { //for each secondary dimension (skip a few to prevent redundant calculations) + int k = (int)Math.pow(3, i) + (int)Math.pow(3, j); //calculate the ternary index + Vector dx = Vector.unit(i, n).plus(Vector.unit(j, n)).times(DEL_X); //go a bit in both directions + values[k] = f.apply(x.plus(dx)); //calculate and save + } + int k = (int)Math.pow(3, i); //do the same with just i, no j + Vector dx = Vector.unit(i, n).times(DEL_X); + values[k] = f.apply(x.plus(dx)); + } + +// System.out.println(Arrays.toString(values)); + Matrix h = new Matrix(n, n); + for (int i = 0; i < n; i ++) { //compute the derivatives and fill the matrix + for (int j = i; j < n; j ++) { + int dxi = (int)Math.pow(3, i); + int dxj = (int)Math.pow(3, j); + double dfdx0 = (values[dxi] - values[0])/DEL_X; + double dfdx1 = (values[dxi+dxj] - values[dxj])/DEL_X; +// System.out.println(dfdx0+", "+dfdx1); + double d2fdx2 = (dfdx1 - dfdx0)/DEL_X; + h.setElement(i, j, d2fdx2); + h.setElement(j, i, d2fdx2); } } - final Series output = new Series(); - output.setName(proj.getName()); - - log.println("We got the best " + proj.getName() + " projections using:"); - for (double[] best : currentBest) { - log.print("\t"); - for (int i = 0; i < params.length; i++) - log.print("t" + i + "=" + best[3 + i] + "; "); - log.println("\t(" + best[1] + ", " + best[2] + ")"); - output.getData().add(new Data(best[1], best[2])); - } - log.println(); - return output; + return h; } - - private static final double weighDistortion(double weight, double... distortions) { - return distortions[0]*weight + distortions[1]*(1-weight); - } - - - private static Data plotDistortion(double[][][] pts, Projection proj, double[] params) { - double[] distortion = proj.avgDistortion(pts, proj.getDefaultParameters()); - return new Data(distortion[0], distortion[1]); - } } diff --git a/src/maps/Projection.java b/src/maps/Projection.java index b3b17ee..3cf0d1f 100644 --- a/src/maps/Projection.java +++ b/src/maps/Projection.java @@ -305,6 +305,17 @@ public abstract class Projection { } + public static double[][][] hemisphere(double dt) { //like globe(), but for the eastern hemisphere. Good for doing projections that are symmetrical in longitude (i.e. pretty much all of them) + List points = new ArrayList(); + for (double phi = -Math.PI/2+dt/2; phi < Math.PI/2; phi += dt) { // make sure phi is never exactly +-tau/4 + for (double lam = dt/Math.cos(phi)/2; lam < Math.PI; lam += dt/Math.cos(phi)) { + points.add(new double[] {phi, lam}); + } + } + return new double[][][] {points.toArray(new double[0][])}; + } + + public double[] avgDistortion(double[][][] points, double[] params) { this.setParameters(params); return avgDistortion(points); @@ -368,7 +379,7 @@ public abstract class Projection { final double s1ps2 = Math.hypot((pE[0]-pC[0])+(pN[1]-pC[1]), (pE[1]-pC[1])-(pN[0]-pC[0])); final double s1ms2 = Math.hypot((pE[0]-pC[0])-(pN[1]-pC[1]), (pE[1]-pC[1])+(pN[0]-pC[0])); output[1] = Math.abs(Math.log(Math.abs((s1ps2-s1ms2)/(s1ps2+s1ms2)))); //the first output is the shape (angle) distortion - if (Math.abs(output[1]) > 25) + if (output[1] > 25) output[1] = Double.NaN; //discard outliers return output; diff --git a/src/utils/linalg/Matrix.java b/src/utils/linalg/Matrix.java new file mode 100644 index 0000000..d7c4743 --- /dev/null +++ b/src/utils/linalg/Matrix.java @@ -0,0 +1,199 @@ +/** + * MIT License + * + * Copyright (c) 2017 Justin Kunimune + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package utils.linalg; + + +/** + * A two-dimensional array of numbers to which you can do linear algebra + * FYI these Matrices index from zero because iatshuip MatLab + * This class isn't very efficient, but I'm using it on matrices of max length 3, so + * + * @author jkunimune + */ +public class Matrix { + + private final double[][] values; + + + public Matrix(int n, int m) { + this.values = new double[n][m]; + } + + public Matrix(double[]... values) { + this.values = values; + } + + public static Matrix identity(int n) { + Matrix identity = new Matrix(n, n); + for (int i = 0; i < n; i ++) + identity.setElement(i, i, 1); + return identity; + } + + + public int getHeight() { + return this.values.length; + } + + public int getWidth() { + if (this.getHeight() > 0) + return this.values[0].length; + else + return 0; + } + + public double getElement(int i, int j) { + return this.values[i][j]; + } + + public void setElement(int i, int j, double val) { + this.values[i][j] = val; + } + + protected double[][] getArray() { + return values; + } + + public Matrix T() { + return this.transpose(); + } + + public Matrix transpose() { + Matrix thisT = new Matrix(this.getWidth(), this.getHeight()); + for (int i = 0; i < this.getWidth(); i ++) + for (int j = 0; j < this.getHeight(); j ++) + thisT.setElement(i, j, this.getElement(j, i)); + return thisT; + } + + public Matrix inverse() { + if (this.getWidth() != this.getHeight()) + throw new IllegalArgumentException("Only square matrices, have inverses. "+this+" is not square."); + + double det = this.determinant(); + if (det == 0) + System.err.println(this+" is singular and thus has an infinite determinant."); + + Matrix inv = new Matrix(this.getHeight(), this.getWidth()); + for (int i = 0; i < this.getHeight(); i ++) { + for (int j = 0; j < this.getWidth(); j ++) { + inv.setElement(i, j, this.cofactor(j, i)/det); + } + } + return inv; + } + + public double determinant() { + if (this.getWidth() != this.getHeight()) + throw new IllegalArgumentException("Only square matrices have determinants. "+this+" is not square."); + if (this.getHeight() == 0) + return 1; + if (this.getHeight() == 1) + return this.getElement(0, 0); + if (this.getHeight() == 2) + return this.getElement(0, 0) * this.getElement(1, 1) - + this.getElement(0, 1) * this.getElement(1, 0); + + double det = 0; + for (int j = 0; j < this.getWidth(); j ++) { + if (j%2 == 0) + det += this.getElement(0, j)*this.submatrix(0, j).determinant(); + else + det -= this.getElement(0, j)*this.submatrix(0, j).determinant(); + } + return det; + } + + public Matrix submatrix(double I, double J) { + Matrix sub = new Matrix(this.getHeight()-1, this.getWidth()-1); + for (int i = 0; i < this.getHeight()-1; i ++) + for (int j = 0; j < this.getWidth()-1; j ++) + sub.setElement(i, j, this.getElement(i