From 79f980ece64036b5a834fcde7f5d675f13550c5f Mon Sep 17 00:00:00 2001 From: James Bristow Date: Sun, 4 Feb 2024 03:50:48 +1300 Subject: [PATCH] Adding flow matching posterior estimation --- README.md | 1 + data/fmpe/config.yaml | 8 ++++ data/fmpe/corner.png | Bin 0 -> 33579 bytes data/fmpe/coverage.png | Bin 0 -> 20916 bytes pipelines/fmpe/run.py | 94 +++++++++++++++++++++++++++++++++++++++++ 5 files changed, 103 insertions(+) create mode 100644 data/fmpe/config.yaml create mode 100644 data/fmpe/corner.png create mode 100644 data/fmpe/coverage.png create mode 100644 pipelines/fmpe/run.py diff --git a/README.md b/README.md index 7aff6a8..dabd7fc 100644 --- a/README.md +++ b/README.md @@ -11,6 +11,7 @@ The following model calibration methods have been evaluated. * [Bayesian Optimisation for Likelihood-Free Inference](https://github.com/JBris/model-calibration-evaluation/tree/main/pipelines/bolfi/run.py) * [Differential Evolution Adaptive Metropolis](https://github.com/JBris/model-calibration-evaluation/tree/main/pipelines/dream/run.py) * [Experimental Design via Gaussian Process Emulation](https://github.com/JBris/model-calibration-evaluation/tree/main/pipelines/experimental_design/run.py) +* [Flow Matching Posterior Estimation](https://github.com/JBris/model-calibration-evaluation/tree/main/pipelines/fmpe/run.py) * [Tree-structured Parzen Estimator](https://github.com/JBris/model-calibration-evaluation/tree/main/pipelines/optimisation/run.py) * [Polynomial Chaos Expansion](https://github.com/JBris/model-calibration-evaluation/tree/main/pipelines/poly_chaos/run.py) * [Polynomial Chaos Kriging](https://github.com/JBris/model-calibration-evaluation/tree/main/pipelines/poly_chaos_kriging/run.py) diff --git a/data/fmpe/config.yaml b/data/fmpe/config.yaml new file mode 100644 index 0000000..adbbbf8 --- /dev/null +++ b/data/fmpe/config.yaml @@ -0,0 +1,8 @@ +mu: 52.5 +mu_lb: 45 +mu_ub: 60 +n: 100 +seed: 100 +sigma: 13.5 +sigma_lb: 5 +sigma_ub: 20 diff --git a/data/fmpe/corner.png b/data/fmpe/corner.png new file mode 100644 index 0000000000000000000000000000000000000000..fb2fa5ce4e46f371545dbb0d1d338ef2c6471e79 GIT binary patch literal 33579 zcmbTeWmuJ6)HS;44ngS%EDk0aCIkY(k(ZNFhd>Zk;U9EVa3pZ5Jq7$n zz(ZQsL&L?|!`sZ=3i8&>!`0En!_m(CxtEnY)Xv43i;a(sljXUshleXvke%J>|D3?) z;%>wK(r7aXTm-{aP7ex!;F-Zc2;ant?I4iQGFdsaM?;njx9y4H zEMZm1p=AU4x!)pl(GRt{mZIxDL$u5F9j5gSyjjYpsdZtYMBcqcp}D`PF0$z8=xAzI>Z>bg7JURaX}e>-}+KfC8Y~pX2&-9ZfBlmU)nyER>PM-rwLQBn|2Nc{~9Nk zhY*1Sn$Hn*!C}=9WK)_jd@X2D|nee z5i3(#(AYObm~&u$H~DCByXIGy!5D^bOoO2NOllr8iaHoH?yrrAAA-gip$__37@phv4N=iGfR7%E3az=~jjdn8_jVXYciP&ge*% zki5LSz`?0l!@K*N^>>Hh9B3B`vFm1vPPVI~me#X|Wk2Q^Vs-{*<`}!luWrk}+Z=&6 z9+0h(tieeYv8D1?4e~DsoX^Bc8P%VL$2?r{nB03|y4F4|%4br#8uX^9;Ho)pexgJ< zFT0s~!RDTDW$fu$r&38wPQEwu&$=f4_}EqSW{39c+l(#0fQ?_Dl(Ga#YoxW-Oh@hkK56CCBB~^egt3~e$le%R^>W~ z<-!9_i}h-y-7BYsj-52v#anA^%KTg_DGZo!v`tJ1~?nJPtFe z+S*?8+bY{Vcw&VW74hUQt^W!~Wo2bg`^a5h`(HTr#gg{y&()-JT4A1U4X=Y3_?nv= zyw>x9%kmfga?oQS-s9Oa4YX5&=ULpa8tfJh3Ukxqsuh(M#%jsYoPoW{j*+Pn4S8i{=gYKoE*nnhmJ;le zS*OAhjB)GT$)Kf;%_+EhVsi4&G!6^??M$egoZQY0-huQ4647@(6IHA?%H~({IWSvp z!$cS_cAIoWUbw13$GdyoCU<(ULZ1>7*ElZ;vJheE>FH7WwTnngO9z1GS@hd0r*sB6 z1Pi)7vYGRE$9Q{r)NfzM>Gnqo6`SOZvNEm%tm0r|^e-03Y4$VE_!E_7Y@{0ehv<7K zaw2pN%L|t+|H-$6;TIE;ktyVe27W6qx4J&r97^NpDc=rq9?zEs$(7=BtEZrf@IPc? zBabJvGO?Z3HD+~GSzmwXvWSkpI>IIm1#e-Xiq!e-EX^xnH-}^-_^0ayg`MT~h=YrG6T7gSU56Wcm{7q;qMpX#_X3aeh+&!$VSy)+F z4Lzr|f(xE*{wCmvpp$hk%nPL~!igQ=53JNfr@6~I3l1N364Zuo_OE+Hz7yZa*u>9h0-5UNf7W|}W?rh9;S8_(@q0y6fR5SSj71EC42(rIPQ z$n4nJkGLJm1>?VjL0ypj1g>qb(o}FW9sKh+>G0b3Rnu=yJeIH>0ZY=v0ci?{F6{hdJ7$)IzId% z<9K1LBd&!$&BM#<4(9M+AKtIE(w6^NMY1|9uW_vSGJJR9p(5`P1gpW-!?q zwYA*&1qDq`-V|558{nrNUP1iq|yi*pWn&j;xtSUDwp227 zf_mlJsq7}WoyjjX=G)9X^>NgpnN)CyK% zq?99?<-HO}+-5vlrdPArUaxP7u74E^ec+>il@^&>AUO;I4t~<^Wbx$0g!=##+b#ob zndY=p;wpLCzSp|u$NfwywX%J4F{SJ4tL`2p8+=Z@0X=>N4w z6ea+$Cb~}XZbTk??@S=`v{T4R1C=J=m9TJTwMG9)vDkgo#DuydFS$$GIg!ImSzcKg zb}-kda1t2T)6LH&KU~)&FflRtZ})5L936!>n}lwfcXn)-eb#ZTN3(}PW=<3F=hgVC zWVYPuyBK)6^z`Y|*8x|q5!j^Nw?`dEb(h^5|Ghgi#M}GxrvHU%EqT$p^0^J8p zXA_uffC0pdCppZ2qF1zC=<_-+qIZJIVAQ|0x#_l*Zj+duJpwXQ<7TSCyhx=`WrhHh zfy-vZBe27LQ=#5@X<(4OL&?jF7epQ#xKu1DR|IrK;6qp#vdaAx6{Cph8kfKF59HrQ z>2C4MP`f~`VWOAvFT>c+_h1Pr%sLZOQ&Mhj9H8p-m{)&HN%5RmFUZgT_4~IO7+q(h z;(|&CJ3F((#iqN{oJUJ=Ot)N{K{@jm$jR`N4(98ejbJcPigQ$-Zl=-fjG(RV;EOyF zj}OV{w%@S&^RzXqSf`{P$DiV;P3_kcPs_W|hN)UVaHI&XZzUu+6!WB?}A%cG=F;aUk z$G0ACdmS$R4thAHrV6@qgMeinsa?8x4l~)V85Om-JY0&U5}C(33f=xiQFV-AbHqo> z8}_`dLUDT_7RyReUIaypL52G*`ndS{`ML1b>igv6WPC!x9uPhwNh+eduW~ii)w}Cp z#)_VvJpcavQ%ZXslJo8+>gvh|Ch-0HcYfE^PPgrYUwiRI|EF$bI*33)&H>nOtD@yF zzo;l21&gS0hT)buE}yBbqEJ;*66w>YPg5|McLvJB(o#pe*sr9VBlGz&(N_1R$o1?R zU6$6Hzon4Hcee+{Wf3UH-F*JOYZ-xtBCp~)3_oqx@m^a0KdVbOqai($3g=@f$cSK# z2?nrj<9mc)i^ot&A(|i2LOa0;Ezdn^kmg^~*0Vws2~zR%$&(-~ts~YtvpVDNAoX*O zWxCa&&zQA6m>f6?%FztxW+gOY&beg)D+o%Yv0d~0V6Vh&u#dc3cNT{suy(z*Wj!Kt zY1mnPx9TD3_NmLyToOl<+~8#yCnp{K1oHu4SuOzoA2;7WdkZ;>!ILbW_SwiJ{V&aL zfq<)oR=_+NRD46BRJt>WSTxPd7@ggFnZMPuu{-ttB}USR?JWOynserJ5kS)TK{Y=X zv>D#QRzY2D{8$fFXPZ|w z@Qc(=Ez2n%|6ShKuYJX^n5ZbXTONRzW2wcmp`#+7+1sx|LC66lFYDpK4GFp)G|_Ez zr2zmK71S2T*{YBX=jNzizvSNC@0m3CFHcQrX1}|!ySY3{$Y*K=6%ydP_rn={Nf{Y? zg8!e(Po+D`s_rPPX&5XOC-Jy46e5vs|D(3;$-Ut$A=m-zZj(KThLtr=I)`fbq>+L}&zZvzwgNxmtO}Zgmt#-Yal96#g zpVmQg1#IPq1B~sfO3p~jQG1MWrOiy4p7W!j{U;E9S|co*8uX@fj=X?PjNY8@EzQ3% zh;~C?GoTK-<0z1c`!m2?Nk>By0>J)hWzhZL`CX&i`qjmP%Up{$JH(*g(F`yZz>w;C zU{!ODt(VLGByloe=KL3-_qUh*mYmhTr`GoN_PQPKM0|aHr)FjVS4jYP4{`J^SM&u) z#9EZ3q0Y|E)?g&x^ed()DxzDL6oJxu$*OH|yP zqxHr9oKiYxxPfz1I2e|U&`?CSwlj4=3=%{xNOTO{!W)(`698Ik{XT{z&7Wz>T7}1oq;wm+!P*5J!x z;Y)aB7#%HcZ_mD53CaO4Q02ZkKp1)q(r>x|^y}U&NXb7u_q2`QU0eOAsqqDGEw7+p z4hAI?6fQs=d>-PC@xjV-Ap63q$y0}zZG^O+Oh&orMag__jrFOEaxj-j!b2G_4`uYo z_oz63-s6cvOdsy9g~3o6LWg-E*p!0F`S~v%?zY8pcp>E2sNg;bBm(Yq@R98$HrW8Q z$EKyy?tdy99B}QgbWoR04*=|g@W09$;}H8TVLB3#@%nyx-v=}TqMhqezq}BHpAmwo zyw6~A;fAfhpkg%+-=h?0O^V>|o(Dm;qjw(uw@{bLxh& zZgF!CFnkGqP!qpc4W|BEeUEo?wVsfap8gw)y>(ULhu`()XWxxavS~uzY|oxO14tIy zPw9QUQ;^U+qNb(Q8-j!uFXSREjZ$MX`n9Yq7G#{0lM~(c00Ck{*Dm>my#KohGTU+S z^QQqQ289`Uk55c=gLjjYm;VQUF0xGzdIi7ShzQgdw6rCal}Ef-^(_S}tqti3CyLig z(nw8E^$xEDF0#0=@Nk4fD{-W|&A}xWL$}~Li%H*rzK?7Tdi=jc{`fED;eR-0*+syr z+}#oF^lQr&yIaGns9gF)bTFGNz7;I(&4^Ha2?$54gjzB#o;ZjmD*ydzWg|w`zE2oN zhRl*<5|M|}KXSu$rbrCzY`rY$zV=7>90-n3xyqOC)}pM0<9Fu(UC?-qrSlWPGo_El z-UAi~UP3}%`=AS?65vOrr=1Y^ z`}|Kklx4R={^LHP@OtQa5qx{|0wTnf_uny@|Bj)UvRl6qYMMcWDF6l2gBnbGKM$|L zHx5BB@snN(*DnsZx)HhHejmV<-A1lgW>w$oc!}Q#0GMExMVevsUj;@3t9CuJXVi`| z^Kp9b9FD(yPq&(kq6waJ2&2DOj+w4d)n;!qr9rOf79M%!|DgL#RzVFjL8uwI!~$W~ z^z|`Fdv%krNqe3`DKrgLtHQ=UCHvaT?W}P82dA~}aAag;Q*JWf)ey8aLC;r&jLMoV zdi|dONZ>bN6L(jL^O^}um13c}(*jTn0|VnJ5fRm8V?LAc^UuM&w&G|L!W`XY3#R4&jt!i66%UWMN@=LPGCKkShOsdH!AOMeQGGBE@*^R z`(8rf;o;z_%|_7^zXcYJUMPMFJ2n6sCLNdc5I;@sNRC)eI2Q27|yL%0OTF5+m`vpF#6;!ThA9BHLr>Ccbr-k2R zcSoHLnYvD8(&EJm7Z)C{M}y$X5InOrEP0XV1U>kEuJ8uS%@ydF zgk)p`KqVXw!yr`Vw+7JzhaW%&D^-O3WdRdapkMB#uty<2MhudEA;|~)WgTPVc$tDV zkRUCuj@OV;P(Ed34E~)iRk@uX6?@F=?0hpKxL3AL^|-L|(6{kBiIl@UY_Iy~Pbq+6 z=S5(kxWJ*0m6cURzjO0016bCEZF+GsNEC0o-QPVzKO4GZ9JO68+bta7fcTQ7tl+ZH zc`jE;Pa`~L14}Wdknk~i_f8b@cp3EA50c3n@JHb7#b)zfh0wv-_;;;R@S+=T33au# z$8$D09MG+ybt)HUX9iBrA4^apkoY$?Hjbw>R1l6oh~0&_v|mR6(DNs`QYf9nqOS^w zVfy+39T{IN#jz3u>yfSebNJi31keHs^TMrT4r&daW3IU=Tm7(d=fQgQEnrt=Y~cz~u2fmZ!U?pS|K%Y)H}}+Zjn@}mC$jba1TiZEI64BL$@lJG|1}`( zzlMiBZ?_T9alMpG(v{k&MUntkQafZY^^KM(5co=0ge{4+hI4&=o$1y`+AmUXJxuoJ z&!70cKHW+_3%bMA7SLId0gR zmV={uQcUCRTO0`BL0O`*Hu8A7KANBKyv1*trgreXKJzIs)b6BuiHb@eFk~hI?6T4x z2oMjI&P*bb!)6~z&c)K54y&xu^vV=1yRdrN#5Vvt7c*X zP;14@QXpC+#;~`i0fg1Ckoq=**KJx$=ikjbFUfHLk6ZZlD*-Vv ztOVu_7qrwqd|z5YdHMW#a1(U9x&4sEhM^Xq*CjQZEz3YVfV*tD$(H^0SVx48g!A+B zea^;ZOG-=E*FI490qG>hqYO5B`cc$zS4nY15p`I^y)bPn|G2mGot~KwVwBZ>d(AGB|VTnfMAhb3tM_ps$Kr4 z8;gB^F=IlS6PZ<5MlHF(ifL&F0Sy)P1Yp=>@WCgUy&kgr9Y)mc4+4dh-}N~hSAI5m zCjf%D+q$M2Xb0e_ufPnE@;OsL=H}+kj#hF&mjt%=0he+h@(Hc(56h_k7B>;Q@erJ? z3=eM>=e#oqAozd*Iz2VzIxbBN))h9D&ng0BMEF$ct8&(WB%a8|#po1xtq;`iemONP zflw9c>JW2}axVH&f>P5Pq5b~Vd72d*;n)nOvci!XxpQL(_^gF(Hvl`7|y!QY;&5}-$QBh`qp)(3UfQ;bzhcOfo2@U|T z2?1A*2GhtfvbHVysfz?rP!!aGXovAYLqo&j&+JXnZgP0K=m~?PquYs!Y<8E|pcew? zMgW)q47~-~AFJ@GJxl|Qmui=l_DMFqlV?EC0?W8;{u~Hv32xSJS~7yIIMabJ56K0J z#@N^htC>5lI4qroo7Vsx&ppUNlICaBQ)x30E_V-?9WIw?b2T#58RHgu^ z%4^B}KYM%sfj-8^q@^@f@Min9s^U{#_40?QzSY!x3Z!!F7Sg2%_?Ugr{2uCAHdRmM8ufG? zU*bvtJP+*`{W8GL5$ZnV{s$YV0BlHLvd?&o!yPBgGry>IV67TLz@oX!gr`@2{m_Z~ z>y_YO$ktH$`b{&lp}a27*81N;Oq#F`fTe`9zV)%g%oypp2D4vW|)S z{#E-=kC4f{Iyz0OJOCWB-|BGIMp00A1zz=x0dT|sz|qgVFZDc~Svv^K(V^qE zKW(}5nG`-R{COq+GD=qhUF^S8#;K}qp`8p#QB_s10X=%1-b}@vl|J_4sXY4HC=fi^ z9O~v@mP3(+z&nWpCP>o^@Js$6E-by88NU#gxQZsPq8_h~K28a;v(c6z(jp(jPxjRN zU+tX$@qq^C&JRhR5hwKNQiMPj8wguDmqm@YNQokNbaC({XOLd#M)K1qVaCN>_Bu?Z zpN*^VyTyPYX|&~npj$lsaXu9zP!0~)-YL8?i(=awATxzzDfd6qF<#EZ{TNr3X9| z{c%Gqh(o$z1H5SpGE2dDq)dT$!MG{6-`D;xIGRwK%oQhCBx!ZX4C&tLcKbni9?Av@;KGVBu zt$aXVz%*oSeR(`Ehu~sj{sB7QsZ>&Y5*Z~>8M_>-;fkZV)W`lx(9q7JmmJHV#Qgpx zxZ*OntrtJnGo?TI3p@^K5p+GiKympe@wjrcGNW|i%G5C*c+uk7eavttY5rzGWG8Mi z3>+MI5w<8|%8GX$c%BY^`uGCcnXjR+H&m;CK;ogFC{l%$_~4UV!OWc%6#?PfMa?+tAjR3<$2;y8L9#{PKh=7jKh8$0mw`7q3 zu8;v1v!Gwbt6)#r_B*z*!lfHj+jM@n{5=Piy&U|$OMm7$|HAS0Z6b-t5%4sN+J+Y0 z6>oG+A^*%mj6_WU`UCZN$FFDpQ_{{3!jTYAJ!rxnKLSWm|GwX`dyJBoe}R5-Qh}vu0>7XeYE|s4YUr6@iov9tKq3Zf*fAGK zm-g%EX>`TpJ`S~i>qmwwSN;i&(0cTdHUbg$2`%qO7j2 zngV_5sQR?7ILMsjalxp|7?g2|D}-Hy7+fJko4zUgp4`dACky9qNdj0I#o_tmK@zOT z*EWdS9bvB#J(iAi`sW;0bZFpXf785TqLxQyGwPO$_8|vBiGh%ml+@DF(v&k@*zeHd z9LTB#A!HWJT5nm6TjKF0Ji?Go$IWCvKlkN@Fh@(MsXYO0W2Wdo<4|hm+1vt`Sc@G) zIPbIv+OOdL4ng5z9qN_teb#B`CKx9Czwi((oVMS7Q)F>i+ft4oI}*;bCD_7X7#hISnE` z6Zu!>{sIP|zR_sXiI}?SpEF?f@>qA)@Z?Va*++1BnMz|z6XpjN-Qx=TGiNv@S?&S6gQ=EnC(E%ogh;|6>Q z3J`Dn^xJnX#`bhK&W&iU{-_ri{kfCs-bxYR>nruq>zG%q<~7i(IYkUzd*E;3_*O2P z!wxT~a9>mLhEgZlS`jh}4~}9Cq9v6sdLAyAPNP1e(o)n~KL1W&^f){70E8y6dA#SB zNB?N;)|L1L4Qh>8b;*1xh*%sbs(+JDVoojd_ zyGJ`HmGg|ilC@?4hC%*iv_yL;8pY>bA{8)jHAOeTn(ByVvCzb>esTN4Z7;NR~Fw4CEEvzElDg6q&6Y^&$tX)nJ4ch1G{a)n2UD^YAKaSrNy8LkVvJ8THhdssarVGzhDYhBW*7U zeg^y;Rz=Mis!3?{rek(<(85o^1b{97?{I<{EKcz5BMmJSg9!|{ayx3bkJ+FvA99K& z6tO5atC{ca5Qg6-esqn+EcmnOp`l^uk!mmcQwby>(>$>mH~oafL~?c2(>!==Vff4NJkGMDirtwuM3OIo zE!tACA2DIM?H^v&L@Bv^L7mwmF)Rs z(Z}BdkO15#`P7&Qjb@Ih2yINa*OqM5Yda|HbE+&3Aj2C^wmeQ?LS~*#Rn?S*Ja0>N z=r1f;qL?1$r%eZcfxh~c-}aOLCM^%-6P+XPtz*p`29xX3*{V0r=^rsVP1+aZv{LCX zXjKKeKB5N~&ye|Hu)tmu1pAjn_gPcj+X_Lv zb3U_X;AKjg2FG7(phBu#dl4>lxJ!TmP)g+mLL)J7pM7{}i4#l;U)(@H6y!6QzToEL z<8qqk23;H(D75Cg({+>LNg`>|llRrreZw&@4ccCTFnlWbkM++aX;${Xzp0!&-q-5U zTsi+`UGoq4<4(cHqGnFhX>io8j*Q4yEVMSw6n8qU@qUHew)#eg9H@~=U>u~tINaOg zB?L7S2H&!gbaob^rxt(ck0j`4rK|Vos7OhTl80(N4bQhHEWOnExAo-(Fg_ta7k(W~Tq*MY6!ivXjU`j? z3D6eEmoH!LtXW3{SIKYV}dD5|MwcR%j<=eQH%OqSC%OvrivW6IJS zt8yxN%gM*)YKFl)ni(XRP{*zUsNO%<;I=p;gBHd+PV^Q~h*t&nY zD?!mOD*{-A`YvM3wn*)ykkQ68HN-XIwu4eDIMbCXD)FdA_jTxBs%+*fAtg zk!fsmD&S`YA$sUFuRaE%8X{bO1XK;R~hJ@%4T0`^E!`IXpur0MY&;A|j^oI(?ZR#?S8u{ltKkOc}`l z?%8Mg!)*~q4Of-$|_e|lp8ljk6cfwRDIi`I_zk6)cuKI%864y8Q@iQCXi1{0zy_ZR2p=Ro6KD2*$Ldc{h}|c)X(i7LlQWPz zCvXlflf%71A&`iq+FEBVf*2-#UOeRZj0+XMTu0E0FC>B+oX6kK9$;of+asiLziwYz zk0_ExaIG4U(}PgyOEbYsr^0W65#@e>3Lv z)dA>10Wu(n=18tmUpe=mb_tYU2rxRmeurs?1M}HmV7CQA8WzUUthoR@Ap!Ot-U>f@luR{nP2qaK(OVbgb&&__YsIS zN$ni#@u>g${HUj1Csi@Reh!&vD;Bcf<|InkCW;T?n5rCG;{~6(K(qL$k5LIvyS{c=UGB8CUyy{GC{JJkHSQCwNk}e`FKtx!NaOe zZsj>lQR48&Tcu5|X1Y&yAs!!pRTnc40r6ZAAUote&oER}e5$8#_R5im3Mijv5QOF~ z5{ljwWZt;94D>V>{n}sWp8gSqQm&07Ll-^gu?NI^_kEav*Bv9Y7?}lG5gJ7A&8p3J zj!d9pevi}St*3x$ia|v0Ty;K-k`vlL=XephANE@2EY93##TNB4YRC1W7+>V~I3k5C zv+3@+pEEF4Ot$l~^F)K#1TAH1K&ylY7VRw@lz^bcIvJf7Dp62v z1!!5L$4NP22je%(hIDK7<3``TqnY}H-6^@4_UP%}ZVeiHIbhnG%Sw6&U_yh}{2j?! zx8JMgQwZZRFvmFg1wzFGUy!6dRn`t50nr+XIRlFG)e4^pg6r}%adtgHYCh9c1}jjW zTnqEkmp0=Zg$zgQ6aD@U5zP=yRn_N4psXA)vNEJO$Xsp?d@7-nLiXLXVT zq9x4gqK13x9AAC6*j!eKb{M+zyr#;2^ke6XA6(T*x+NrK*}RY%1uUxZt4I9XPJ3Vy z0RFYklEGF68Ds+UF*Ex*HA9DdlYq2Ka#n!qE_I!}+85PEZpjsO1C?C(L>R(iBIr2i z`J#pzRkh~7j@trvI6W6zvR;n{hk;{2+bKG=<0~`3hAPY1R`1`^1Dq<@p#19vrxarK z4k>17EQBsh=qq$_v(+19UsYARWXL`(;5#3{9J4IxcEinHd1Og@|El$>7yY~#dM5M7 zZ0{T@cyue*C(M2u3bLeAU+DbXhxxz!#(8q`cyT&1FDdNt@40oH+`kqtthI{JddOl5 zUPlFF(Ti;@5O5Z|ask!h)*=IbLmSlRuXIPB(Q2d%(nXaypUC?Lcv09`_Z}}DDiBo& zzY$OQbo<$PO>8gb`gBurldxlF2Dct%xasT%eU5c-(U5n zwNH=}EhW~#zC2k7r*bto7S(+zVhUG`lwWl;#xFdyb*0R5IQU>*nxsia=G{IH?0D(} zO!He>GSO#$J*$klZ0&Ne{&Y8$GrQ?ULTws;T?SK_0p`44xe<9Nu7l8LWN;<6xpYGJQ4WKhvdsKEA{B_yq#9AIv`AjT^zQH!_i%6Yw_v=@uLwR9d zFG?pUfv3kNJFx4#j(xY5tU`z@3oxG3u^>`tZ^kLBxlpSRq5^mav^V(_!n<~GpN z^l=ztmr1ZuiCA=3Z*DYSNS^(SAp6E!udZvguPKhV1`;KrEAmQl$Lz3JhZuTioyA1s zV<^Hw#nL$)=vS0`FYcc2Wh&?^p;r!*_O9pbZ=4=%3ZJZbU8M4z^eaM(5>t;IP@JIk zCF3;6s1!LFdH>r2SOMY}7*Z%e_z{8%qxv+wSaf9*v;ETKV~%1IrW-Ru-+qe1v|7IM&SR*4aOUfKEQI z;8KYWJ3m(*mgUPt6*9xXy8)&-`-$qyJVz4Gz36sO(;JG&)&CD9}pt$gtF} zJ4i1WWW0R$zyO&z45Tj4yZrfD4}56f=uvCO_!U3GXX$I9k-^; zIZn*9+^YR1Y@$@0aII`-6rbRz!VWgKG2Q!s2YgOXF3U)9HfI>s*4>~{&}D`fHWQswHta{ z!#Lk5+P{V8pC7{?>>Nd#!)T^*3xdq~1)|n0e@1)+&K5)6Rkx_O)=A^ff-5PLWL$s! ze*9Rg>jcW8s`iQ!$Q($|-tJdA=8rh8_1vk84zbg-4RuZzOr*KVU;?ILbdwHYz>!<^ zCdSzqHQC>)9K0N37fp_dINLq(4L?TCFv|NjUyb~m5P@HW7dZ}0+ zd42z)4=TatKhmZw%#XuF?^z>x=#Lyyk6L`*IvLvZja5qM``RhAhs$xWIeq8i zExkQjK7I15R(yw>3{bT`9iah{S$7)2;@ls5^h`mAJ72LABR2u@@c*Q5!5nnwg9hrUa{eL zW4q_Ic7w@d?*%8-wW|;>3P56iM&#$R=I?H|5fRT8&s)SA+b;@%2A!S#tP7B)fFE_S zmwc)4T@l+A6f-4a`rgX>{d_I)s_G6~)XaEZb`OS;6|Za+I|uByz4a<#r$bz$jqQh# zLq2oKT|ni#B=G#IO?1DidvrS(b-$R+eq*YEq4KUC3#EU#xg#@ev%0Rq?|sOg^KFXk zh$kn4AZrIdS5{t$+TDH2{ZmcSrib%7VBtoOD3WRNWRna|T~Rt5y}N8U`8$AQD8vva z`pHP;-IIK-_`^_-6Di5}U7esYz9d-8l}4$BMykLZ`0F(5|SVtf|rUaP|N*TLn55lWf_JD^)OD&X1QeyhPK@f!V*22uZSu zq-&-}1|JuW!c#Hb`~aS&uLBx5wDY%G7m!Q!oxBFfzleAlyX7BV>S#8W|ELa!X$$-n zUf~p-EHB%Zzv(hb7#_*M@mF^i^}p@obZg~P#mw{rXBr(tXO7M%aKZeFG; zm(mnD$xr6h(>SkDf&dHYU|j*~k_<@A+jz?R$9qHNyoV+~fahQw+U{l;q4R3Jww@ch zSBBN(ITK!6@AF9L{n5Dpo*4Hzcm|@t%t?XwTu~a?Rm*q+Xy;@K=Rh~#TloH=?d(8wi1B69 z%P_uYFbbT`>@Cbr9)jRD`NS%}Mx73^NW@IO*CAXJ`GR-NIk?)VX6*$tCzYdv3=C&e z8NS}B)x85Dw6gv=F+nU+Ht5~f-cxz--aJrK0l%EXFbg9QzGQYlALD%N4!L{iaY~-% zJN?XldHmB#Qt<93y+`Bf^2-a*{MP%CmnM+)i-kD9CbO%$k`K|f7J#MW7RcTG<;p2A z=%vJw^QmiV_x-RMY~24bXgIp!vmOU`i1)Wy4%$s}N_s zuJ7Gs7Q~RB|5MlR?1i)iXPsHKna-I+K{&$IwfZ^Z=V9yEn?2y=EjeLpYO2m}Qz*^F z%REDCb`{ih;1&;XJzmuGDm1^ z@c8FX^{z&8;ca4_Ya@6U5u2PhHp8VoU9ZME>ntsXtvM3zH*d$oPWCHethwr!0E+us zV8P<6=4#hX1i-7K@US-KIRZr{Qjt>LmM!b_s-TsYa8Zj*k&~7cHy-Rk-T$Ct;))lK zn`%5I8zTBUO!O)P-V*##TK37vDz3HR5P$v0dLvyW-fp$gAY@n`+orGTspyc{5VhD` z*G{Z$zJ%b+dWVND=i}YObKvm`sIQA1j41XfVU_UEbc>@!YIZ5x_V1KMs+t4cR$1^G z*}(9qQLHSo$q^C~LPt-3{Y`plGS8kcCdi}gwr#B}{>VdfZnisUMyPv@m5>u-qQ&I< zi7py*#?@$qIVch>R9G{6BvXu;WsUIG0>;m7Cl!J0jATL;)8pTdxp21V4fK91oTr-U zS7t%a#P%!yl$DV>UuTwvAKWe7pYZyGce(+r#d-oNmElWEvM!EOM|$H10}P4YTe=%_X6?_c}6 zj13(itok0g026_YrXM#URIvTONiB-@K%DwgIe(EY5=R z@Z9ENvw?W?HPiqC#;+Y zt8C}Tu$hC61i+H%zS@aE<+novcRGFHaTwcGet209Hs_?K62UHkzT0kJ@oL%F;Ae9v zS?IYSRZHP|3*|J4f!=&7hX~Y1q4+a}6U+C_qQ`@RX?rwNU3wxfonC8$Bezq(KU5 znK5-!Ut)up0N-RqB`|#{fF0H#yIS-{KU2yUzAmi120Olns2{InO?cF_0A69Q)3Uia z-*8l4lpFKeD}JHxgsJD2m-Y94{s2>vkErEGu||GqntkMr%s!2S=@gwW*qET7G!DC*lw{_GO1@<=!3Py9|!wB|3tZ;HwxDDm-o}H*BVcB z5)KDfHzJjdOd{iVY!OqtOjoP>srIILZui`Z#uY`lcR7a+-pe_Mv0HFN&HQ4kQeP5d;8VKOQo-LOSHOwqm&K=BG*(uL@ z%57|JEI$`^LMXXD6uWki1@q|@dDBUd&^oTcSJfBIgO&{S^ZJF}u;}{BXrr1f|dNdl> zE;$!(tQsre9`-snw0Is^WKD`?*WC6a9M6A8-M=kmX&^7~R^ZnWCo&$fy)EB-vD=4d9zcK)K;P25&q_^5H_;JpR1yO z~6wJ zE&6>K>o+q=Gn@j>MTdpPch!LA zZ44f4f&w_c1`9{qZDQh!!aG&;J->bLY@!o6Sv~hmr*cX_We*215nFCkAXtG`CtI1f zw|i^X5wI3PGbMlVC5(DF9KtG?%I0-2I%n(iVDUFzIw-Pa`d~;_e{dzX{m@XiIgi4& zmRK-6yv3si4&jaTT;X|*uMU){CtSy-%cZ9O~ zye{lvc1Y~ai;`O)+8l~1xRJw>3Pp+VjhLfGvYVsDAYxT%A+l6)y zt^!imVH=m7pN17svek|ZyurKlVhl`5`Dennk)z63SM_J7Heel%hvPMDix**YjwH zVTuZ#+l_>x9m@C$1kwBZ3fNvcpG)?hgSUYHoI2jT5a1?|XE^bKFt;_z1&L8R9bY7n zde+3_O~=)Q6v{?-DSwoO z4g5%COUZXsaZg^~F6A!g2d6M(VRar?-APj3_YDz1CRFa`;o0>(vbpUX&fUH9=l(Ki_#( zOFJy@`w>7}`Xiz*VQ!V;PTc43pHg=%xS{>507c3$aBwjGZL{Hc^>eePmUgJu2+bR; z?br8DH^)jD2GzoY6ssqXpZOw`-ZxS%)5$lr3(kW zTMiiJZbkOy8k80p*jW6u7WAlmB@ue@X>k#+ap~ZBC9V8YK4V^53jzV!ZR}TJC8fow zxQ(_}rNs1{3SC1ZP}s>;Jo?yodAVe_Ba4Jcd|16>$kL+G%mqimD~ra<>P{O=#L_cj zq6@C{Fhm+Hi{tD?LyD0vq|K&W6U@oU(d4+numscqe}(jvp94Wgg_p9?MMqv_uw%GC z+|3@`=vd#R?1jqc%Gy4p=Z2wse2U-Mk?Xe zfwHnP*D(6KmOT!bwDj~+uVt`_gbBH>Pkdit$O@BH>-O8rEclY7m1$wA{5p(!a21V4 zJIFp(4OaQxlPvE!E5v-u=20aE4>yn%YKM*|SA-UhRcq|dPSJHccH8BuiXO`igwY9` z2gYZd`_`|Psv@50`n4D*o)|vrt*D*UukF-j_n#C3ijHbMPot!Zg#3c6dtHmS0Q zYobrMbn`1z^VXW<&|eGp-eO)SrKF@-SUYnl0${-qtRg5-(>JGEa69cj(9T&U&Z`jG z9*JETiL zFpVM50yBUE`9}Hoc6Q!3F}a?ankp$7r>XO$lK6rCaW`AgqD>W_7sa*TDMG?HKEt*A z`L=!nP)TZj8gqkN$)Vutd{}AZ=J?}4n$MTM|C)SA{uSn8E>6KZM%$XA(YbgAy|Tzt zk872}I495&x?M8hzTBmk~-;Gm3d%gRb%!zK7}PdvSICB#I^XpkGlE z-{PsrS)6OwPXE<;@BrU2s6#C-Eh&12op&Mfu@(qvQ?IGmqjS8Zm?#gMB?3aH?EW?R z%e?payqdk1*-@^?4i~=Ia=9)4`!bF-1H_W#XUmejI3y-AgI_xHzWPlH-Y2xD1SU8Hjv?^(C#{@MlpnrH-vVnQ*Cw0wXtlmv9al0?}O!3B1Mvrh_CL7z%+!urCV-4{i0;^+>r7urky$jSi-K8xMGWeM_IZ6Mw&9j3i*)Ge`qc-Owu~;N-}2y7w9A^lwm#N z9T(2T#bsfqY#+){tBcMl^MR!4-7z>EQpy?MzMcDri!1kcqfj^4)ivNkj{_YCvHyX! zE8^3q`9e0%srjE9+H#Aq8Fr@UXaYexwJcZeUFhDrx@!_nKqCM-;k{%dZ@inPU~2l4 zIh8azQapI@hfRbN)_#Af&``X_>bo7>VomD{GGq^Or?Wx}hAB1HR}(uo?rVV|C&}jg zP`nL8z&*-ml-a!e(#9aJhT_Kkxrltae(P5Mk2VA%EyjNd1>Wo5@d4=NwHTTF>UEmr zB{GWpumuLmE78gX;|chazDnV%>R6e%J-v3VTS#=|3L`_4fWQVcftt^D)p(v+dp?U@ zyPKmfltpq+^hBRecyNUzKBtIgl~?^UPQYmZF>t*m5-=2KOpQ+PZA+i8+QMU5Z`^eh9aaBnGiT!BPFQv4j`4xP zZ^KQJ-rGkUjX8Qou+@d#^ThBdB}Xf_>FQSwjW*L_dJFT}*ilN94xeyi4(dW^d0~_n zJJ!g>vC0bOENpL6+$nFr#M(2gSa#513k%xCQ1)nf?~(( z*!r;G5p>90`ndPFLvD6M(M?6g%B=0h;zj(BjX^Yy`Q1oX>T>TTiN zvTt`Z=?&Gv9`UT+P<3pd7x_}+i?l2|a8E(Wg-ni2HAx^J#NcQ~PYlHW(!SLQR?S%;ccSj4aUnJif1X2d!i zu^d$YP+L}@yHotVUMu&F-6b{giB0+$6e71Wv{gHW6UA5bpzNF2t)Jp?K~DSSiYED% z!`f+Jv&48TpP7J78aNi*Z{NO65KOiI4LVTV2{^;6@w$3jbm6|DKkKL&QVuQeJt2{M z6>%*lVQ!weruNH^VjN6Ah1BHl98+t~%!k;&D(z9qS65eDh?LZ9OeB;!X(txDENrWH zhjkx5FYz!Q#Mg9h>cZjJ``F;o-6d~cXlXNVA-qe4RMorz{{SO%-amW)7NVQV_teKk zprISq#n0mZ+6?KZuJjb=(TqkbFLVu1xNoI1@fDk%bMseEH?D7tzB4)1+Ld79)fAe~ zegR(l%;(FW3a5@X7XET4x#yxoQ|=_t0-aZ9W`DMJWnz;oGrppP^oEl&LAaNyBo-DH zGHrpW)VYEpHNt82K&g29xk#h2Q=Nc-0D|%X{M6?!5`hWt{dbo1&xTZ^FQwa`*5QC( z0qtnVTqsD5vn|N!mKfV+CpORmUA23DO3A|_598Gv@^0+SB|>F7@Z<@M=2RF;(+>#J zAD{ny8rjTJAXG&|v{9 za-FE-16!sJ-|A>a5;axO42nx)F7dT9M%KNbK3@{ew+?->)yjV`E9I)y3az>khaX58 znxMOT;HyuG&ldX$6d>E0Ydw{39qnQQ$&vu>Z@b zSsV-p_;qpUS_@xppMD_b^2Dp9e(4ITN7OEL814i>M-<~_oxg)P_ovCv{M@A1ZhiDq z30qEr+k(*YY50xjG`eIp2ufWj|Git9^Zn;HlnL*$g-mA{=$)LR0+&8}*tlFRZzhvx z_%tzun$iqP?i^aCHoAj@>&Cjy|A{kuK)d;MvKB?7a?&Zh?>rUbtZQt1&47A)djTE@ zzJAr(KiAGXadpEzXmr8BH7|Dxr$`}7=pDRj4xYLrG!V0#GsV9+ov$d$rO9gJ61d3l z&zb?rcI(@xsH$4D|IgXoy_}PBESbrXKECCO`M`6F@+)hPy#&b%{{2|ALBWL(?~2^lAY?pzKx;b_GYG=-FvN|7dAx?-psNitg}7U`z+4_C)M+W zgig3q?Lj1jKS+}r%``^F>ZIc0FrIbrJdMjR0u~bS??c%RfMB7+B}(p?uW>m2wo-Yy zf=wq>&Q<(wAoZ*n)%{)*#G-fX%2-?Lu9PmeARRzmwzJHE$OgJYm_a0!3IacwoEXKD zk00o0B&Bs&2{7-z(QBJxEK60biV$OJ(b%uhZoE)lx=|4es26C2DHtnJxzSl#qBoEE zy+Dih^BIQ6s?5tOhM)eMsgnb}AoROBHfzJow*G2ON7u9G5 zdYTW{rra3_>%f(21~zE0fa4L_AAcVQ5)Tyl03ouuS9=x@gtf%pER8$v>jsLt>T&++ zS(D#`@g06A?^IQmd;+<;Ym#r*ka}n#NVqWvC)noO5 zYQWI%F8kgS6&+(_cJ4kjTh4hX_E>CV(AN1yD8tv}x{ly6i@WyHS`<`73MQJCYSH{B z`Gb^thpvpKJ8jpc@HQ}Vj?jVyOwxLrC@PmBLfQy_ZX{qn!gT>ts>w)Jq-H_m>h)TM?lVEKWtOFff&buj}#3hL$AexWcs;_k0y>E zy!9?zI&XKZ*#up9*ZtzB8-~J8jkkPH#>{GlRb?F2ZvEa0uN+B(z1{E4!&{trBzJ$R zyGGs3&&qyLziAnlDnU`_E8jy*)EQ1eCaGdV3eC0KY}bC~)KjN#foKPmdSMUJX9BX? z|LJRE_=%8Pq1a?l9qQG-wuW{knEX1+LL@f#C@tIBYf06?k>N;g7wH4$)t(MfdSd~A z#bX(bDO;~6JZQ7?Z9v@G0aSYPnQ}qDKY7`T)eO4$E0%tAn8{Ep?Q2wM3N;?X#ok}i z@N)NBxmy2fTYgzm2DAIapO!k29s|3RUulx3dr5yb=JEtzC6QtnX>;&E|DrcIq@qkA zI^aC;WMdHszdUtq)m}^TImatmBB+tuw}bud7fvp8IN?<~g9@shQzZD?b@Cx{;g9}L zKyoi7)kqMgp|0NDIl-#HNJ%MWWi?0ca)|HdeD~^oPp)sklSm zTMQy0tiF<*u${DdB@vK_8i^Bwny~7!_gUeWPb7R%DEqa869GGem&0=&W(N#xwAu*8 zGt}wOqasEe_3)MM>nl=*Rm{~nfcV)Zf=Lhx>GlBqE3_Wuk(HGNiMiomf$4~2E2oa@ zpeYIg{KUnDat4znY^2kagYI5M(`IVsY{hM<7+Ssn!o>6(}?#=l}Rnzw%OS$^w!x0a{dyQymvwI~<7aZ8>-)tg)q?_1j@pD!QU(_i}JsyH(c3O&{%;c6=_{ zdDon0C{SOE`G7&g>gU(;7QYLV7D35MJ?e+IMfy&zn*ZDMT#==c{1S@4AY_G$d3_;2 zm^Ua;>?iNtyAJC~DNMlby=-V0uOGX>`5!NZ-$a~Q@S#{W1OH4C5uT<;YAkg8B z`9e#6R(`_K{wwH7rNFmEOOB2Gvo`Q<0RE&2&R?w)!`dSETu7l+212Oy?dZx;-)gX! z>`{?Sm4M|yaBN(&-bsHgHM7z4v&Of)jG?9mMqA^K57N*=)GK$DOn=E;H*7_+(UyDR zToLW5@5vhm8t;HJR_mpcG)ep;sr7&qkqxH7$SKr@@lLv9L>nH#IIf-NRGeDV zhlgu^XF@;}>Fg}w-KiHx2pxFDv_GK1RnaFqdD|^M_OWvD+te#^?i01cnqVR-bv<#Q z9l7zc)J`M zK(m-BXyq~n=21e$yLl$qBSqW^~UtJl~P@F!7`XU0R=!A0NFx>u&>H6NR zK$f*TRr_H#HjHdCKr=2fa;OTX;n~a%w~c!KqYP`ytGAin_nU=HrESE-U#xMJnnWwQ zBIdOW1Cg`OB0;3D=GYtv9M~Ru)%j99{}a0V-1MeOKb)`ahKk;TP5#Pr1nLIJpP-Wm-#TdKtH@T>U$XmiX@1DUCap2 z^!}XxA;(9ds%Eo26!MMr0xZn#t4kuk^eX~ZQi(xbWGf~VZ0Z*J%$;h^;_P%3>*LzJ zDcOq6K*#6FZh#pzCks`$_`pj73`!sdAdj$2p-^AE5atufA0zJ{DD;}A_HdEitnG6E zU_9YsCTC5R*a)<*T&GB z`uQjD6$~{nAh*guE>h7ca1qbQ6)5- zhpE;4c^XJYIn}`*0AjO2{jMJ`T3EC~X~r&gT9$zvrS@K-UGy%FXdig4gf|Ac>TU;m z&#>+hx}NW%gp^;+wa0rPJ8gt(EA04HfpM$fLYT&NnVy7szcJ{G*(&wB& za8oP}U}b_91nOWQ$rZE}s`rp9O}-mW)c_r zqg7t9*UnR^$CjWmV z5o*IW*1+9|rPqlFiT+`eRCIpx;*`J;oH8aLJUO{kFBq@s!Y0O$F;&B?4ffs`-hUFF z<>&EgtA}M2^&I`7e2_v2nD-B}? zIy3m8-D;}edKOAWm6+e9O~lm))i!keXE|MP#)BOrG~G`g)yBG!SGRWEtvihXgSif& zj?zpu<;7~?b~wbEXDi?#4pH8ep_y5m)zIm6iv(fiF2gXgbig zFLqwuNPrNBjcLOgC}SgcIh9F$u34YfK?OU!-s7wtef=SxGHkVG8?RIejvMpD+=}sz z=rlwRddyUwj1;+xq`I^Lves~LWMkwOpM zVypg3K$L?;(%jZobYUZTX{!F*!9=^mh>tBF9O<~-BApt=e5YSMX4bT>8MXgT59P4l z8rfCPXXI{G#fXy}P$#)G@jo^85*muCL@oK1-$V?JU+&=%EcZ#{Hn|~qZoMW*I!e3S z=U}vi3$$*gi#?e)e2`=h_eop^-+e;>J7Lw;o}JiTA}2psZy-VTuHJ(^dWe;7ED}rS z*1Zp2XAvFW*{&&imjv_sa!l_!6_gp+GrhoCf8!~wcf>-7%eyk_KkumBgvGMN{M}a_6v1f(sbQ#V#6Q^$zGHq6gn?HO9 zuNrKLgT6h^y=u|rFv^CCj*iBD|LK$a$RDWO<<2Rj{SOV=n1um1+1GvtINs+n#Q&UP zB=q=~V&?#MEKng!0S^VTVtpc41f}fcgY*RBUd!TliT0sOOS|x6F7|&GP%}2tmJvB% zc3Ec45#yRTcU>CC5n=+!4>7kz^1_zK*(;)8Eduqz%k3eqwZ+u`Kw`Qr@0#53N3Yn6 zKrQ(?3YsYUt=n^AXE9=cK(ypf!2lm+BG;*~+{s9ckE*}G#D@SD%(U#97!c(lH6|3K z0DP4SBqgJ}pQI2ddLQtv>Xr1dx}d&ls@kA|dS{hsrBY-{Dr*QQAHe!(Khp~_#mhUI zkOH%SNAS*TkkSJ=$wc78e911o!>`O;Ae{6!3Qm6qmb@mq9y>>{{1kjo=dZ`_x;s?a zQa+D&NeiPB5Kk0+NNyfVCs0rGEUs&HiNSmQ0mN^J{+0E8O>N$FOi1pvq>f=|XxgS` zsN{|AJ<8fQn6H(jo?0;1}=(;^HCs3~F!Y zO6V!-83nxRPRLLAv05|zw8IB#z(gk0uy)QK#08M~;pU^Sq=YaQZ{NO+kYoW`7F2wI*6vtX)%EC1$HbKwn0Dg%4peSJb#&r0;-2D+1LAbmV17aabuGSJ2g^CufnL z0HWG;*7@z&GQu<9_@ug$`zyt_<(BoKyuRTtw|TL13Q~g%9`xA_#b$?xlkJaBf59aA zRVG5AygCYSl^YnwIrwMKFw38`{z%sHagpWLla2wx;B2M#P@WF!if(N{*(=2-vriMSDEMc;ft-eoF5_<-}6y6HyvU}U&2_Y@z;A~97$^}FF$fRr6q8oVR1em}v_Emckk7rA1E2_|m*>pQOA7WLrC2@E zDpEu){G&(bfw`-MT=;qclUR683&QgP43sNb)%D~I}!A7O!eRsm^vvP+wqR9;Kdu8vKlBIdXnaxtJTv#M`3;r0xNI>uyxTky8jM>IRpBnF)kU?V
F$nRT^fixRFddFX93Gwvlj#k2J*&Gi^+`-QV)yydEKec}! z!mu3njPFLOk3uf8i110S%)0A9xDEED2CV&sr0ti(U@p&Sx8}R!xor;EAE2MBR|<%` zkM7xQ)jrz`+gcpS{=C;WP~u@g^V&8Mc!j-#E3SOj0~3c@ZxU>aA)#a#c*#lNKQ7Vd z03$`DWVR}BZ{QhJ^6jTx-+_ujL zLVuUuQVBr=(Inp}wei`ajlmbai45gQ3Am%=o(8*H9U*Pu9uum6hm@#MA$#`N+Iea6 z!@XmaaBkwCXFAojvL&=41b@Hjnizd{+PBPrH%Tbzw7s~)jun~+8VRt1%j7(DeJBVo zs^|w5>M+_S z)Q8B|7k~YS_ybg^?MJ3NX1?ZG=kGeX2X*8waqsn|Hqb%T({PJ4H#%@&q{}( zsg(-s1kcYv_g1E$jv$lJ2YJf`ynMhH{UZQZrMIVk@gxcyGRnmz`{aQ0zUof}00aCh z^n)FR()+o4ONyYU{k?UF3Mcjf{4v}ODHwf)dC3H@x@1ZmLEK)iPWMORliv1g@Gt_A z$q#r;kS>?ZdZB}39J{VK<@Xz+zN`_Z9nfk^xr$RmwI$e>RSljEz&(@uDw2 zO)^9&Uj|!%Hq>Jqx(H*v7T`bPkKU{mNCd2xK>wwnvg!knTUPJgkSliG|6QH>+j1Z3 zbSORjC-TdU#KKw`BS9@CR z?5di>q@gzDC?vsfJ%6rtOb#rX14`KZkDHKRN-C?@zT+UYeB$ir%X=UJMQMNi3uTym z(CQYswMU7M+HXiL_s5m3PD!k^x4$Hz!UU8y$33A~^nm9xEabH=ZEbDxv3w!GPzAM` zml$#I>Y9Nm`PN5ktckP@n%wy3JyUsmAE9p(kQE5x(mF6sB%B(R4n0;0@a84|Ryb(w z4QCx>Weo*;>+xnz2<(Gs;>p@1s7O>QWqf_?OmFyal(GkIw3hNAwm)bgXP&-j?Ekx9 z={jc&{h%cDGl5zHQsA7Sy4JZg-LbTw!+ZU$(2X7J3AZ)78O+HOw$i53X?{r{fEvJw zh3jelXW$pW8Fuq9PWifI)93Nf&u)__Xmth-8vQx!@84{3J3FQKGrjdSP)IXJ>Q`?LB>l0|xM(adO%&!`6u$Q1H;j5)Z+Q*6o!J6PpwjDg8G-secz$CY zPNvkXZ8{u!=Xt?tx;e+~ini3J)L#>>?-vR~R)F=Iu6N5<9-sZV1txfeS#PcRPlYYr zCofZ!It3=Z52*W5+Dxc-(kV+z@ZYtqeaj(ndp`lFMETm_EF-BWA+h7KpsiiwPVu<` zN?a+d$|?62GmZS4L4Rsp1E7AUo)rQcy1%8UzLKptR-Sv$H%_L?RlME ziY9DN-1h4Mp-%fqPA;Xa>=ShBieppQJh(#(2CX42B5`@y+F3BBKXN=IWxD#FtTF`z zEI^!}be4-)db!5M2LP8<@X=amW@2LUJZYQ0z825`5V~I5m*XG)evAE^ckh)-~TfAq4TR=e;2X)MWha_M_`!Vu>HFJq%Aw%yU)s zBm&!WXDGc8wTLILzV89E*y;W?2z;0YFg{Y@yuWt%2Xw=i+eMG8b8AmMbIYbK04D9Y z2{JJeccpfzZ7Lf50f}Px1$_Q~Jul*UpS|Y&jX&i(ecA{Iy2n-_>}piLCk?nno@nSw zLHGvYWe$`9bc+!|u%wjKZMRh&h$P8$WcT?S;zTD9{`ZKJ=zR&2KA@~}A0nwh!W!qZ zI|#LijIY2XLNZc-z=`Cc*4?(;h7yCSn%WP@ZfZr){ZZNrS3vz~Jm2{Ri8g@4(m6mM z0aaVZ!$Sa0;$WMKfkIPp>>#HMNm(>#k3m-Y8ldQPX#Pz=Y9onXXssM1n8|vuKx#(a zGdh{IsH9}GQ-TEG%q^+WdKc;-F>N2|q5{-hjs;4U(9VA>B^A(;Lum;}8ae^TAI18B zvPSYDp*(ZTe&Rk7ONpckz{m4Z&x!((B`GLD;2eup2WbaCL0tI@QBFh6vwR^fcA(0m zu!49c4Y4Y_R?ra zWc>%P@C?zKnnmVhKA|6KBP$vMmo#{~2hnCuL{<+8Kbr8VGVa zt%q2JCE2R(1!t6hU?3=UpEGzQ#w}OR_vbMxLUaLgsUZfmuImUHqlJ%-A zL=|88r|KEpMkFQ*Ne2So*U(NG9IL@>C!tJ%#9jW4#;9E8L4fh6M{BpKiOh>dmI?=( z#|T5DC7hh@K#|6%GwDIDar}2z0QQ3nv+VVnR7a8DnJ?rc$-!gR+uJ(=pCOm-EtEOP z|NHz8a9QY4GQoiY?x_)4(^Cf}+UwN?Sv8P=rLz~fRSOwUT6BCB?<6*GphsK&r{ zoCe<4n~=DtOwoAzt24x4utKHq2_6YUHWHiK0zW0aYP(x5=(~Xu&z(C_RaL^wqM~|( zMHUe)IeXc=V@Q(jO&I2CNVD4%7^}H|Na*e*A;Ge;&~pi+H=?)*lV1y@>12<9u^DLk zRmpJpGKGL4S$GhasA0$p;d1W zm)dFgo*j*|LqlT!Bk=iOY%p&l(K7f=kc!%-SVmICnJr`N1<{)r?2&kHb8+Z+B%h3Ews&`5ltz+7AdtO~bc6>{^cJhlcvbb|!b!mMY>SAC>%ME-m|0a6o9Ay-P)-ag+EB8yp}a+5oU zZf|J92Z-By*K7F%G|Z7NU%u4aNN@+vBf~1#{aPUEDZ?AyGbou(%s(U{WMK?|L~~iayp#CZ>RC%}Q!sD35PxS7_0~jG4pC65y$QFqn(5l9^I9X! zi)oPB5v4tN(oE`15)G-$*d)UDm6S;Bi@fz9paL!u5u>Qi%;m#&cxS_g*O?>YUs&Q5 z0+U$>R$Oz-XK7(?&p%?1jf98cU!*et=d1|V!Z(q0FK)Qea}{Bdb239(+@xe=FTNYt z`{Z?*yMr$ggv&;d8O2ueYfoWE50U^MG=3!T15B|F<1gesA?qQW%8=2TbxF<6;oDv^ zzmNSM%L3LO4_qvUai2YUxoEB#$ck9`K%ha-ycTeeCPW(WipDTJ5iRBP+}z&K^=ls7 zRgf85{R=D{0^a+#ATwfdC(&r8nlwEeGF^bm3DN}B%A4OJ_UeU9_`m996c_K584ire zBfAYL<*(mC0ffD}RvBB`bZ-uO*WkV@rE zab{%O3|L{Pv#eLhM6hgR>7N84t1~2*A+>!_7X@^gPp8G&GC@=)Y_%ZAg%Ab825<4a z30QCEp|tzs=TApCD1cQ*k@}^rwe>V)#$1C#uBgx9`terB9gFVk zNR)N|ZJ=KZL)hf61d&*fXaKAD02Z3e>}(UrrG(k{Lv%C-0RZP^=H|{Kxf{&P0q@_F zzJ^E)c(#@r-0l}4*|l%^$J}4Rx(;_5ljpm`zFM`VkXyAW`!nPU2CgllgesqWeFYvY zOa3S4!)~&(|A9A~fskM`h(|?o5@9$x6uiK|(+~!Y-xpwn_ylUI&rh<{(J()ehkC(p zLnG!Zj0oU{iY!L#8;Ft53}%r82r=o((-(sWc%=f^s&}qE44p=TDq)2hfwV{Dr6470 zUlOamM=W)oI~MTmJcji+P-)6~f*kiDX#=m{7*cvJB*9)o`e~8*hf6dx9wP!DIknT8 z^0MY_k(mr7YiKx0(9sri=OLu$yoT(I?#|96ShP$bbRc(bX>lWAxl*mm71NMJ?aXd_i_2z=~3l+1c4$HQ`X-am)`!p;0*T+{Nd&h zgMosaFj2c$(lMmX@iq@<;JV1XkK71`Bu z9OF-~h_>dt;&=Evks8JV)Ux3=q5&iaUgqf-MiplyT8Z4V`m}8AE^V7bn82JHWN&HZ z?D}5?ZspX8^o2=7kQ5>pA4y1roR}>LMM1nSkozj=(;hL(zcMO^@`F%B4`|mOK)?RLMF8Y*7yTeRZ?MYa4#2-4 zGmYL&9+XUwGzJl@l?+*VdC%6yslVM|&EYR$6K!zCI%Yid}p83gj?1);vP76-AD zfCFBKjiLiE2p~Cz9v+5l*Ft~$;h|zHkHyS%=xt+zi+>8*W86_uQAD)7V1cV|Xn=6> zi-@5c#;-9rz!A9_2t&A@o$W#33-W8j4E6z>8?b41vxf|Lu6R{gbLDb9LKwSrf?^LK5y_zfQwz zPWs@%16|7oWFHEf;7^c(1HNdcn>S_H><nmetZLXR@qAB-GQ18g$r^&Sgi&MpOL`Dh4V)N#BRZA*JJx} u{nEdut$!Vwki7W+^>~DwsQ&-IxbMgG^nz{*Mf3y;{*sYWkjxb~c==yEk=G;u literal 0 HcmV?d00001 diff --git a/data/fmpe/coverage.png b/data/fmpe/coverage.png new file mode 100644 index 0000000000000000000000000000000000000000..67007b732c0aea8b02357e31879624144d5c8d54 GIT binary patch literal 20916 zcmb5WbyQdT)-U`6Q7Hv!6agiryCejpQ$V^?KtZ}g3`9b@Q>3IpKpIrKK{};VC8Xod zwa>ZdeV+HZ_m6wV*yEhxKEhhxnDbL}1*s~_T*JPJjY6TW$;m!cL!rFj3WXn|5Paen^P-ubDuDYdJGqm#A09XA_48y5@p6KCh=PD1SL z&;Ij&V6%6$V!vmw$$9s6AFcAg8YM)E0$x8LOr@8_fTBJJ$Y-|O+UusqG`L& zBr`xj{07F=u36zk*&!jNx<6@!OTE?U4Jf? z?zzCq`C0R~4#Kw~w|Z-4isOof4yu2h92+Koi}TW#ydxL+iW)Nji)+XQ6+n$ybC)2U z8k1gu_JKIQVyvux0ChY|1nROlht6Fz>;U7khv-BY&E`QEjFV^UR{Xm@%*oAsev~BZ-!xt8=HTEk;6fHa zP1IkgSC#uHQ@Zp$4FR>WmDN@A3)QWytzFi94JIrcoNuF}qzVdMgP6tN$<30Q1ST!k zpUKP1e{E@z+P$x*e-|xaxomW1rr`lriDClR*yQA`VbcDw35?cpwt?)Yj?Fd{ky*+0QrkEnh$k%FlKje}M3JT2EI;ka#P=z{M5W)bxPOurl)T+Gv}-2tK%Dfbi?A{<&B|NO5`y8j(Pld)=-irAQ)`@8OHLZ!x9ADC8<^hX3sDM$u&_FKSiTJbwE0sjQMxmVULPwA20v&3rjI zxqHKWOU|N=Zoesuii)In+a7YTc*db>UAN>%iwtn|eWz<&zSP#1o3>!1Zt^^7O7lDK zuuAixS3uXu(@0#@%2T8JVZ$q*MR4oZtxS3UgoK2wnQW%VpPnE5GAJu63&~apROKNT zs!4)pE-Cr&;lrQ({gH{~m6f)|jwquV6=~@!)%eEo_p7x`HfM{vnI_TpUeT0&{Xny=%PT|eR0h{@&8c4_aT zB*(qGx~A=3s`$e~3J3|=Jack(roDGBXr(XB%+@x^HZ47!nuJ9=+}2B5o6Mp!dSJrR z#NiF)@!#|Go+gj2kO;Y3ZMgj=g-t ztw+FJ<%L=yim2dt34bf)c9Po=V!Q>spe` z-mYyq{x&Kn-*=ul`h7f=zmt~BTD1!O+UuyOW`+k03=9N>gcX%^^z^3o_Gv5QRQ?qe zydxtc)s=d>y35BX4JKmj*XZ2b+}+*XJaKOYCqrHoO_*C+zLCgi-(6wLi4<*MBhyXd zD*y5HK6&4YpZ=1#<31}MgB~wcHbx#keE8%U4jx_( z{BNOmvVOR?_;6oKL~P6s=WFr)4y6a9xA%nUGaqYL;osyDv30T(HW&rmDYH&HCartqx=!%?Go`i~0!BNQcb5 z7dxM8qw+KFOBJcH2{7C~IItSZR??mPq>+d24;OB-^4XBEc;IH&Agf;G-1_?QpqYia zd5PO^!`9GSfmvB}+*Z9mqJzT1dR=a1*bIMi@hn% zog(73&tcI?q*-w-TLL8BtlZsHulC-33vsx!=*T* z=PoWS)6*ieW})QVZ&L2N?Q9e`+R!z9dYnH$H%IL!iSg2b=hvMsiA+Jl3nG-)EgWgL zPf9%@8+#jnNpuntaLKt|!47Sliq_1~H7TbES=U%BTI??Oc>Ho&?oJ>Gkx!zP`)Fag zwn&*?m{nb1lz4+i`cDu4>=S}HSQBzlZ-I!2h~lu8kyx4clpQ3&#$)FVXVDfJ0~BFC znC5<3fmJRGVz+Iaes*Jd*B;#A=8kIcIVrOqqOboFSmcJ)kN=T*m`CqJk{xCo$?>Ds zPjgh0zZ)-$^tPS@m9v#$gre@b_ZPW?;#LMv4sbJq#)fTxEu!VnH9_ue6(f7ZXDa z?Ip)?sndg!nIVBGMji?rlmW~0(<5Dr>r|ozyn4lVMDmERS9*t5Faq037WSU~AYdfw zugc?A_ZEJhcc|j0>*nV62ABLRocpm-GnrHo4reLarn@DjJPTjFsZtEkNu|fjEb^c) zG+4aBqv}1}Uf|^BHiI6=&Keh`SXfYyYgq4@>IZv1-6HaU=75pz^KR>@{c{)o=%Og=? zCRr}x5G}Qu$YMMrt8YJTKU!Mi%{(*6Cm1^27?`RE&>YvRvMcpIa*SowdjpFA?e}*i z)msck)bjE&#eG*?FE6j|9vrtfNA?H5XBk993~cqX%FTB(~2i!TzL3{I5vZ~p`lWRn!86UX&ZMp(Qja>qOdEeHN=19Yq5`3 z=H@caxG#j7Qt2@o$j2n6o!9hSNPUkCk3N$rbd4!f{9734N0TG$MI-EcN@iqaWNK!X z2Pa@CM}>W>1V@VY0Zpt{Tt-0M)e_s$B8eEk#|L!H&vSt9?1e3@+liHD z+|m70DlTUV{V;ua->p-pQAi27gpz^RP3-Mep!r<8#mBN1ZgB_DP<%~2TgkDEi<8!Y>DxJ;$bd_`G?Rq;TT=*N#bYEdv3vi{FU@ucJQ=37vrdaEZY;JB$6Nj?l${49ld-j6_!N z`O(2swUh@`cvH+PVz0 zOds>Xd#&Xanb=YD(&4G7=gW|S&4F7f8;?ys8HpI?8#bu6MbZYrxh`{Dl5aRawk3tu6xnKQqM(Uu6MlJRZk{*hD4aTcwP=Tf?nJRxoYyP+WH?j3o=NIF zp;)_AqC9;JDk&MK>8q!l8v(NZm(x`lb2STZKF-%1E-)R*ki42nK9i-vBSsOZsSu3vS@oJ`NTLB(|t1lf`ZA!+zG5q;^vJdr~(X3-kxq^IR z;#!h~a@s!E?fb4DeNOgy?55(Sf=S4Gp5egpZhnvw0rJu6>oHv&{F6J) zFRu!E?0j5qFp^_2EVf19{);bv1wx54w8ngTDT1!X-%Iz~y{WiUq!5a2VdA>c5N zA>z45O-FYbkPFnuyBphoUjt`IU-`w5v)e+YHST$rGD{wv=#aOdljq|2M_!^A_VN7j z$?ryJ;tX_juX|JO4>kJvjc@dmBEK-zb6>yXSkmF?>F$-jm~r*~szQSx1SB1&WYg_@ z7Znn^Nkv7KD(vy~mVI4wapSqo51Mpj=|pyg)mxS``L*3R!mRqrbJFeC=bM}_JnMex zrq+8M#LCw5+f6l1PA0+q9{OypzW-}FzvAO3PK3hjJ+zEX-13{Azqo4(paso8&`dD< zc;s771`U9%EH$0GE;-rR^IKazTa#6JMvbB-CMFuWsy~PC-xNS9Oq+&HkWt_C&++{LQy}^HWbn30h6a zId++`4>J5kl~NMsw_ez04K_{z$EZ7)@jBg~w1kfughey9d!DqQ(tZXYz9CeUu8G#eDl31J}jqs zc##J3+X_UV$1qrrr^JY)baD3-Q*Lc&js$x>m3)Ky6Q0rH=hxTLA`Mj>dyb+zfVviz zm*0ejUMtkCnExI`cod6w{dzWZ_pynI`T6;;(f8D4RaDrso?KRs%~C7hdMKO{C$rO$ z+9fQfu9h#aGg3F1Y}8hj^-^YFbT&LGFz_mLQ`6P{jGF6G_N&w6bMS^(xVYv()IQq` zznhygK0Q5s8y1!gz_q&5e!T28Ik#nr!-Dk0XKP7CMci+1t|91LvxoqrI&YZ3i+##W zN*i~SDjQw3zV2F7y~yh56nU_zvN(HOh3)A4uU}z}jbg^0XQzfzzX|3#W0=W#ZHT@F z;$1BxH61Uv!n?(9XWkk@UfoG0;u&4}Y?90Uk1=X$YKn~K34YsqD)!+omxrk32qHM~ zv?*`)wA9L>G+1m~zVy1dSH~GY_Stma0{f_?{h}pZpeta$;y`2%RnzuI$c?*YEnzsUab_TwGi{ z#Zi%EjqUPTJ25YNoMzbb#Mf$!L!O}@O}eZv6UGL%?Tq&1prR}IR?aHzaP8Hj97jm} zZ)T~nDj2Y<_w$|Y4XK1s2om^nnzgFtDBrF>^Tyc~?3lh%FYk{>CHgQ+?dR%hZpHM# zKm@Z^(PDMqT=Pqh!#{VnCMqN&X@U2WiFi6UM8S2v6Uhpv$?Mc5;wRPFZJzI33b&J6Ne-s(q!U-PiCXw|N1O39mr=z{`vFV+uJ+D5=+kP zDIrm#0Y6pN4GM=C5mJJ1S`L2oKbpJiAsZS=HcKw$qOJK8d%&(r04f5#yma-tIc#8j z(ngp(+0MB^+cl|-?vh(cFP}Gkf8ls1Dmn?e{V^+RS0;kp1TQL-F9ULrc zYGizU&z3QyR7v&Fm5x5Eoat^XmCDRr1y;2FOb1?cAK~1`QeyO&J^4y-RX+G#oSU1Q z_t(1FH=OP3y*w7i*@3O3)S-Ob1Iq#HdV`!iufj39s)~QI%AWm?11LC820mQAapOkf zGn{3=lVbQB2^?Qwk$h(&t7bfLpaa3ZuvG0((J~jjG*w|s3cqJ1lz*mweRprdw)h4qsaCGp z-oj{@3gzSJ4RV>_bgge7d%O93!6 z?;z(qz*n}X-E<*sRL5`iW(TT`k4i?#E=+m6M#3#-_@pizb1;*X4oFo{;sj>XU+ zm;*YqvzZx}n5es|oPOyFE*#&h-6{9uUcJIdX!78!niiMvWh($<`deR_ZA&QEz9Q}`z$fNMd#uXUJ{`SrckrROL6&uX04 zzQUW9l$7-A6JuRwP(auJ{N#o~t?SZYmV&v3MXAj&6PrQJ_km0qLH9rD6pt|_O64LS zuC`#r2DUu^>x+#F+aFQfREu-dtufAH$mGYK2e9hy<>j={dTsqj-Z!{R&bK4QBk9m8 zZAZgEw<~ksHWd~T0i7|^V|N*<)6ccFdosiUsw^@^qUOK@dZ&G__)C&2v^p0kJ-EeR zf3yA({=IwmoHs^7S683t8yG+%b+=~6y5m>oGLmml%V*^CH%32@5JP85@r}QdpTpjS z8LCq|Vb}M*o3u;r>-x~2B6Ld1%9hY$zJLEdx4PO6kLcjyGVG##MIjdBu&g%-K!~la zEu7-TTxoo47JAMEi#uAFjZe{pTRW#hnhXEIaE~_~572lqB566F^i!Lfn&xVk-bJ{K z-E>XO^{(`eNYrf#iX9(Q%}nWQZ?4^Ht@GI3*^FnwbgTx5*SFtRveBKgd2cmhb+S9* z!K~+xT;nnB3`Qbs*ucIZLWVKu44`ccy2xHpiN*Nhv1%DM~PGt&^~r=XEI|-{iPOs$$tLUqASQf`qSq1KMmLN2fyyH(H!$~B`PW9J`5w{1?FE%+XDuh~9R_3j`nVI_ z23||)cHb>wElU;g*&=0&?4#k#40WCTU-kJ5n87lIn%p2MfoPgnJ=J4NmZQp3M%$Dr z^S%sLM1L}g*ETz$^%I<>?`>^ek5Z`r{CE5)YcoIJ6zKc+iFyu~@!^n{T;zACe%j@h z*Yql%-l3<@{y|U9NnCcyW1z1OTJI$wE+Z2SCf}|=8AiG9t==GJT6KLy*8Jz*W1Gs9 zx|Wq2rBf=ssY9i9!;Xr7^8@WAWnPlO*@6865&%5yr9_7D+F%xdh@0MOqx|r!OI>lh z0F=t@XGIPD&RkIE#~bWL`qf|i#4d2n`yTNbH=zVS5h{l3vTkow*B0QuzjN-epZuJc zX^g97@}+1`GADodz^%)FVadvz*VIrW#YTPU;xC;kA{^kr@2?LB!-;J;-BI|ESyixqD`MD6nxas|muL84t0>Xric zn+gXj=Qdq|XwF~Ce*rqw!Fo9Vno0o7vLQr~$uQ}x4JgKMt>zL)JYhB%G^7XAuVa34 zZXTUbqDU;UuXt8-Z1>b?RPcu$Z82#U$eNqeR#sL@?>>-8gWJKNlG+XWgP`xJK#Guy z=}K=3Ej|5C2SZv3Nacc`@I%I?bNX~D#A4G_W~H-YRT`N%4)6`e}5eq zNOSHmc2;`Xta_Of4;?dlGq}?>EoPq=IB!Wc*9Ipz|M%`6FS5#6O z?#498khn}K?CuW`GBq_dS=2{2j0uA#GP!7RcTcPzF?qBe%}Wt^RP?Ck7#tOs9#R!IQc366yfcn*uX z5*+Aj^2D=LwH@%dlsT4SjJAtU`*<3UopbfZKmPJC`0^Z0OiU~?Dr)UV9^rl0P5g5D zHGdYY=Yt>RrfTDS7CdL_JPuDzcE90?8CzJifF!yEKh#)wg6`X*f|97$M@{7R%uvZ7(8qsF8PO}3y^4^?=CCBr3lE+ zXO3akip|7sM*SNETL1{A1Q1+^fZ&b6MvrSOxKFS>tHF+mhNcl7A1@n4uLOEH6=;)C zBF%~Kqru5pA1&@fl+2e`uK91~Xd@Dh>+=8;P(fY-;)PS)Qac{CqzRB0A?G{NN-vWF z43)26za9cYu6@JNo01X^PEo_QcPS`VR#rI`l{gxbX3XT8qjdMiJL-`; zE)Jdl{ReStiN)9*f6%`NMIh1&tLR2HKgE z&zHv2onVS#OwZNjObyL9ZGWI*$P&aS-E*&%x<>Dsr{Hqc!AuN1s9)2yyDMpaj?f@| zPqz3C0%_|WRPu@+AN*!-%fPr45g*+PN?zI7B5TSqF%g)%WF?s<~@_*ZoUl7 z%8pQ&vL9oN{1sfX2k!0yfG8|x8+;IKG-`r%ad+}5hmMZULuqLO5YhqOl$Mt-z_WtC zr(TpSK*QmjN3Qc!$T5W-&wxfyIy^iC)3Gg*P66z|cBgv$W)4Jd@)S;GzplSz zBEP>HdaUdU&)k;uH$U;mYG}ampKo|T@o`8WPhcUfU%T`m^6m|-nJL8=4-TANBN6Bta+5`gwqr_!HGvMXRmWR}y z%CQ*khks1Mz>=_K$xKSR<@ReTA`?3!S{6Y%4>ey&5qgPT_W$Ji57I#^t~Is}0Ho6e zkUKt}Y+$$r1@~px#w9f>b*8~b4|kW7T4sT=ldZ(aXNnDB|HL&t=Aqlr zk@7VQKj@TOIl5;-F zL6ubSASXIxT{U>8GJEnxbviV)Z&2Y&D=L-%F@ov1Aw-+uZ!YK0tXa_1*Y_S!2x6qc zX2LZ`2*lIrn#Zs@1o%nu-8&^uB-rWEEwvMw-K!oewYn5Z65qEq2(-dSKtEU*SDS!n zKiuKwD4Xifdz5(v3id+gdr{NQ=zDNzLV#$25klVcFznbaMa72;C?IxOgETuH|(Y5>c-$Ys8 z%b(5Ku7kXPy$Rv{i%}}jOI5YYlF_nra?H%kGI?1wA3d9GNC8TQ6i!%4Xp*h)s;jH3 z!co?imh_^cq8f#|G2ptATUlAz&C~_Gd?{IK#>s~l$*t`($0SD47^D^>;`!}2uMB;? zPSRN8^+%ZjWHp+q;5S+Hr;CG~E*(PFj99ic806GC1jfUk9-hb-dF#DC z!h|O=wZA zHl8fxf^BdA+YgCSfH| za6s?leq3;$Oy-2{>fwCT<3X;+JPrsAyK7WKJ`Fj|C0XKK(d#UcYX8w>Tm`dAe8VaVJ+x?RC*>Wis}4WQ}(|=P}BEGQfTv zVBdJA$$h40eLfYkUCyMKWPZUyCl_fhGJm08;~eBR>l4>_abDYGW#;yKrUe!cfT>Mr z3iNOx0MLtIA64&f|M^2LB9gK(T5NkxP^@1a2{qArKQhBLtHxPnL74NRw` ziv4^ruT1F;Y8+1z>8zDc$0cNCF~h>b6qEQ0n@+Tf^vRu^oc{5n8R+Tf7Z;;fMB%V~ z8yr;I;QA>3@(L|2Er|XJEAo2auhMdGsO}2gyN5+eO6u+7V{|0>m?;KK_JxIo%<7d~ z+elgm>g79WuiG5yIHB%mm1%}-X#}IO+8RHLxL_|3t(3d&Sg`@xN60PyFrcP$Yik|d z3EUPTWzYY$f{sNaivWFauU&(6U>k3JWn#r!y)lP^D2Trvl2wUX$T?-6PLy!lQOUj(7k_>_PQcF24h3BffSWKerI0gYZ*MQg<;!Rsmfd6x zhx0g){Xhr?xM~&yAEW>cJPONRE5vb2EoI?J1muC(NM_o4{J=XyCCN)03Mdfd2NDO8 z69K;pRyz3c_0SL&EDlc2N0pq7|3d}%RGoa~?fdtmHCGqr=OvyzVN}mn%y3v}6}W$& z+j@{z^z@IU&inxXDCI}Qw-E8Onw|a3q0TYqLJIgGII})9ER&lO>Shac0}#nnfk*ED zo>c*h)8%$7HMb5yI%Is9V%`C&bU50vK;7hfDhtAL@DC~D)5C4{`o++YknV~?FjByL zEa?AB7U|WOT5jXVB9>AUD$kms`k4ujEE2$J*>cn429bXt^}W^FAH|&- zjh;NYnx{F&5_VvoDy9l2jE!l-x%uei3(e%%`r|O2bc&`B;zv2 zBxF!VVh8Yp)zrAk4gpj*h*ySq=i)bhq1j`{NU+Z+#cgisYsadiEoFYQ7< z1>IX!2-hvNTx^d(;t9wuY|l$=*Yn;0go9Mh=Brdd*RRyo)u9A{lJeS@;6;%XjBgW3MuM>ZFae5q?J-vOPIG+?+USuYACXF0Hpea=u%N{oAFY^PtV%642hm;MD#n8kt?RRMV zStfuDyti9OL0{mrpHXb|Jp&m=^&M<40BM9XGv0MS#Ni6&vZg$-Z%O` z5C67-%*$nEexe)l^mkZFsKei1qgK#M2Ay+_YF!Iic6G4dNR5;>xvc-VA;(UoMG%R| z)aP9XXeDJF3-$aqT$Sq2m|_UaIpVSb+;iKGu%Keuje7mgcd2g6(5?o6oCOM_yNAcp zO`lTW&C0359gt?AmEpa@eE#Sr^S8`#pI>Z|7?#-XwXbo&f>xUr@BUhV1Do&e={Zd$BG=Aua>&6L|f$DcBFym2O)TF~rQ8?T=oI*Q7=FLNv!>Bgh>J;?}z2hF6ZNAkkz(-ZD{114(ih$?B zK&wyy&%MrPAQph19`Emuda$=kXfE#$tnle-Nql}Dj#!Rys~k< zH6l4x#51!x|I*Oe#;$asJp7Q)Ti{SZ`AUNR2hPl!;^N{WNrF9u$n_=*?rn^||MEp( z#(P%@1rGaIjZ41V)J0dfTKV6H!f6){AQ}QE7$hQ#3Dtk7`$HHR1bxKXfG2VcML_MITJKk0;$sM=Z5(@jos6zIx5a;o{t<#$`hq{vNCD=edJfKT;L%*ZPpu z_&D1)yCD4zT~C2q`m@u@%j*^Q?3tzQSChc<;{c}_%%@JEmTXy`3g;q;7dgnR`=)J*Aj$Hz=0w=cm<{ zNuiM;K* z_7ayQwEKJB0f^*#cI>d)e`j-R3+vjorHWC*As6j`{evADAA+DG%zE!$BWBV- z%>g{br4mgBcZ*a^USy|H@Y?s=4@AOaVp>3y!n}4Z3*wra)3u|k=X2O}3bBFIKoh7L z8D*nRQA`@SO$=$?rq1QfvWxoY!9=;Xg!GEWuw~jhI;CC*cBQB|)hAHrfCxid47WW& z1aW1@mL zR>Ame8xmX)*kZ`8Ab6)3%iRP$YoNkb3y_l{s5nSc1%7@jq)ot=feSs;UX4M5C+=(p zED-30#I7!_ScU*t11##b$yn3>iDUJuB5jB@Z8!^ju}u*&@CVw zRDV`7%bPQR6^ufFiy+ zn1w@0Nr?#1aPHCmfRA=-I)FR!%T?|5hKqAA>Gzb^nn0U4T+MtBbj6qjiV}4G=E1?p z!$bGNdw~ey27`d&@oK2wo3BAc+BYHA+5xR1oHhsy=<+86cOCXw+oeHN0N?x7n>T9E z`3tvTO_rc1IBI0@e^>o7NJvcF()8jI^e{*eh5*kBT4ZRKkHv^jOiWxck9tof##ck` zF2x@|h`8@{_4W7kDp`Ps^IBS(?%WkOEc4EssIu>a(v$B|sLH6CKg?oK6AkHrE?D1w z#W!@l59N~~TnWf<5z+%~Q~KSNtnG*t4J8_wsS2c9JG#5gkN4JiY<|3jSWc;VN95hx zPJmCoR9BmrnMp!#`h4uWWCgVl`VRab)jxU~N$I1<6+mcc7sNVS-(H9}dK3w=VmgM? zq4B}b0R|1AKn4mJAh#@|M&D#Xr#J8?q{r^y`7m&BgaZ|JINq~P>vR7W6q?>pYkR(W)?Iz>?1djwLIO)3fzsQ-q?ep6ISg`fZjSy^8 zgFx|&-^H24u$K14)FqYy8^umIh{yp0F#7?nJ8Zs94+||Oom4FM1G&f@Z^rP*NVQd7 za!OI}xAIwE2L>M5)@EkXz}4??N`uNC1)-?t`fxtAZ-lkHo*YExYFcTPELkQy%e@uQkYO)_{O; zEWK??i@d+!(7ZzAcc5PtsUx694SbX%%vMg(%{9~fpIj;5TVJ07*mgQZT)}p9+1E*Y z{rd92z(6(Scd%{MZIz<_i)njY8t5yCcTnDT-)%!bNA=l|)_+XglFuyDykTh&c>yAi z7*oDJ9&zNu7Pa#4@~ZBw*eJ36>&d{!PKJBeQFizl(tgmIsPEpD7>a?A-oKfUQ%o^r zS0Fe&I>G0){}iNJoeCRbHa52E{nERtBMnAmM`lIGdKO34`=P5~L~e zS$ClIHaW z;CGplgwo0(;a^C1n_64Tsj7yP&CstVvsn3qH0WW74w3FuF+VmLVo%8-Jk^{>WF&Ht z@mv<~EG1q-^cabp8U)f|?$AB}x(8J$&GrPKLWo_8W<~ssntnryiAo4=Afy$bCg!+# zHO4upND!f4hBF68|FS=r(*QW-w6)`*AYa?Y|M(n{(V+REo2xv(^X3cI@k3}xU_TQO z6aN!7hPrjMRW+Ndo^u`Uu3M0@Vmgkg{OQ-Vfy}E=qb&figSMfGOkR*cHMrxo(+N-7 z?jlf67u&)@p{1pzoF>XzQ-kR#S9#ABsA0T_7az2@5k z)Qez-0`l=V-rq{H^@1V_`HT!$quUPXxvWe~?|>|VBQ6e=dc4k~e07y&RRkz0sB5W` zLjYZ0mY0_!EiPHXcpHC_3dx;AwDcpyx)v6`0VWBxOK&hvHot{)`ZJI!GCVw8=f0O?9eeBg=d!D>h=^OJ4k)9OJxA{Tg0{cq%$1G)q#FA6JegxC8LB zPMJBDxVSjLq6}5WhdYG0kxy@D@IM`UkSQGs0_cFE+VS5}qb?WsdXx1oQ&Eo{3`pNX z_3HTIu;6?VxkDLn(j5I`bm z(1B3!0WFTgMN?$~3ionb8_#D?yvgjm$T7Y8$NCye#6RGskpKXc7Nw_O9O2&PUB zME~Io+*%SbZ2sKth*hA$WDQuUXbYC7-RXmdWHvc=&VqovIw+^`axi5>aR2^&2?>cm zdwV|sY^pwm-v@cJ_wYNHCt#kybH?Yu0Qkc#!MLXEJ$&;YxsM=ZgQusb2TkZ3MC^+; z!ZE~2%0w>DfM*MECUH9ro`XSDG!=H)rm#6_bmGpN??hIwd@s%p#%tZ!XZT?*%&Lt5enG z1HH>MnAGa%dwY9`{!m@XZqz^lMFn}eR$nX%3W@^qMiC&dOj<<+3vV$p5DO2g?;6Mu zAoTu(CIwe95R9%807#j74Bz8r0w4qzaZf!PRbt$XUKpqqU>pLE346x?_XoVwA=s&# z{B|KAj3K2Yl`!ijm)X`8Cjj|yVJ|ho;0Hu{!op|>UPUrYYR>^R@6UZiV>k6h9J)I6 zh^CGXxuS0y*D7ul1$)djp@HBPK^p|fqa1mLhPHc{g76Dv!YB~k>3eJI4v?pE*DYN* z9#GmZx$i6n!5zix5qYp6hYwZiBT&}5lV7z3hV=lco6X@Bflep^Vr_!>?34av)@j?`ev4>4i|segchfEq z=-%IrS-2qZZUf#8T8sYf1o-bbI666jL>To9UHlDZZDAgE{BZMXlGY+`=%X8f&k0us zs8|SeFwpjp>+AeM+sYr5_66{w6GuVtSq5mIIQpJaE__)nKt_8u0Af!X0 z2g#>kVS$ls2+mkQFA+1E{KZYuli#DQ?W%Dtin6ll-~+&lEd2Z_4o(`hLgcQbt9zcb z+y6H-4DsFX3SGWXb3WV8reYH_se+*l(0rHGh(Zug9Cr%|i z$h<#yaA*eg5HLK^{a8Ql?{*vz0|Rby1-uZ@6L^>pumbX!N7RkwA`!O^{w3;TP7b%{ z-jkmBG4>fMmJr}dTsA)#|M6g;vm&XK|Kt@fz=zX?Nbs?(>*kpF!-p6^0$9~>+0Amb zj$j*?02_Ym?NuJlzzmpS93eM;@BjsPDa-q4Cz{EVkw#Sqsk?M3lFo1oRA zxcRSG;Zlvv5=467vElFHiJgXn)nH3EKnmkYPf}9o>FChLgYVyq>GXbIcs$EnA|La& zv#bW-5Yi(W&-QP@UHJ*%>1Zag-8!(LGAB2r)F69?zIyVRv%v)n%ya`101N|niN22u ziVTZ(3H9hcxV~V6Mg87ZAL>TCoO3SJVJ5QqGgB>38vq|DQeT1mK?S3cgP`z+R96cEozR^mBqTHvtnJ=e>e8um=Y)M7Ra{DqgHmAA;}RDg z{n6w3&L>}ou40OSJ?A>ZgZaO@$GkeZ0Hki;k45GqOiY@g!vY{_V?I#-k14>LJkU>W zR;yUAB&cQU2vb&1EG!_`&Wc2eK`BCk(*dl?xS9z~B%?}BI7 z^D>E4M_0FD=Ti9j=%xO_L5bo9|47>H*^BcyU>(LVL<7~uXtkOQXV6#DJY^;<%?_FKNU~H3%1C!2 zyf~@kRN`9n8@vV6eD=x1r`a8$Y6*qO;_syyy1t_1Fh+I7v0H&PfG{wySO^FRK6r2t zP*L>(jpA7XhN=LK%)b6F3BZ@3Q!j@~!s^ux?_3D=XF;0kboSzm3Y3a|7l8m{{JkAq z{Qpll|KaHaJTMY()$crkz^GBp7&0mXe&=vu6>_LFB|95iZJ8R#{;=xEW`T{R@iPV*$5r=0!TF)<*g?UUL=V0o1CI>8i*$b9nogJ6y zdapPLLxYACV~~@O0H~Cj1`))vcAQg;?hq}ugw?Kp;?u)zvh2&4si&H0BL1td1Xv~bj+Fi= zsRk@Sh0-wXf@z+1q?97`IIJH-oDm&Pg*fvCvB?E8_6Wt8+G}e8;6@ig91Ka*)HD-) z#W^aL#S@2oIs&J9_FdbQ`$g084LlsNnDUf5BjT zMaupH=NNwx=05=kN$=JGJ~q9LC*Dsuk6VIoD13>f-Jl8aqJd}u;FoTBM+}wHxZ^Sd zrbI>%chgz6Ec|5F#YQL9144Oj%bob9sN!`k)9Z1wC zfy+Yvj)XWqwj?}U|F0`iZddAjPWTLJoPW;Ldm(aCcsM?1rzHL!j77XopT0(aa027S zNRAZJN9!Yn3eaX_m>x4#9K!lgAnSAXE`l0Ju*LVcH+~f8L;-kRSX|V2H-HtOstt=T zAzst0m?oM6PX=z`9dz^$MjlN3ahUKQM9BoN7j$=otkH^zb8vVjpsMU=mBFxtfbPQH+7P6`FyI9K0z<&+ z2!uv10Q4Xa`2d-cgrWDdBEyn`JVz4y`ppF)Q`X2tg zirjy=MUb6~sGtJm2>7es&h>v}x`;F|-Fm}1yNS=Dz_ zzr#>wvz~o@jJ*fwRaeZObOrzKZd?e%@H}AWmi}YQfBX7%=XUvtFpiKGw`5G#*YWY{ zhm@aJ)!f~`6iG7PXdPoyf2OXhn^@qoZh#I#Fc%DD($dr0JX29r^j>SnhVZ;ZOd1Fx zFJ8R(XOjNQmoIZ08z-F(we`{e9w`U?wv1Y7&`X}?DIu8^(nTHFsq$ek_dhPUIBXIM zbX1t|=I}i|6h5$?uHl6gEx4aL<*dTO$-wPabY=kW+S+dDPeND?=QN8+QAg4Q<|-DR z#^qJY{o>m&K!3G)`6W6!`jhc8IvCEtxOz1zF|iW_NHntl;$nk=EabY`4#C~knfm*n zur2@)Glj!f^Yjfhrl&GH;T!_+VIC03CEcAnFF>>Vr#<$=5cT_naWTA-NhlXkhYprG$^sj9Rro$M^>guzSGVxzjVRldUQb_(ws2NgO{UXy)%FsBDa^ z|MM#qgOh>OS@*0$?8Sw2Am9hn!?N=7T#3Hgtz+~Yy{-wKNshiUbcMV%!OFA`s09Vd zQ9$Bv%iP735T|8eP;c-)GNqEa(-9+oEnc^983~<@kJme@4JDg6fiB9GxWpD(s$`tJ z#r2mD!slQMfs{D9j-#(Jj~l+ZR%ru=Ki3CblF((G_-|yB$d@S8?V>2GdDnK(r9IV(p+`REy{H`5)NdDTWG=H_Qreq%}7B z$?;1nm`r?zXVTxE7#n*M!}+-n8FH?k_^8SVJ6)R{QY4lG7;m85O68MQaTtV{VbhAI zkPl*jJAlka0#X^P^KdHSdVO}Am80W3c;QMgQ4`5$H^pdB=l%%f!01Jod53QvEbE+u z<`P+5UHusdBc!1f2BH7b!64BZ(BwoeU{ z&s_5$;Oi_%qsYdw8^uDop^srVGOY1AsmxY<@W6lZ_iuJE2vv4}gRzk{NGZqU4WmaB zwBwVLwk=5@`aEv!7tMizctt>B3!& z!2Q%9hg(qzX%zY6KEHDsxk%c)VVe%%(9oC?md|oj(&T{x4Z=NDNWWZ7T&Z*s8xun} zK0Yo7EI8Cs58mhp0FY`)){?K4)_~g7j&_!eJcZ2;lSNN9xjTUD$jiwMJRSPGviBi? zwh2h7LUH|m;z%Xhr7N<&zP^d0EoXvDZ6DcR0KpY>StpvAnMr&G2?=OzVCclGh(3Ds zNOuVh2GDqxKs;8KkOx_S@TKk+ep==U*~z zOI$#Nk-YnH_iU0t63YeG476q_#7VG?HZ$CB`m7;AgD8W5vu&P@4zFks6(9RVF04dT zMJfM#EX$fiuM}LcWItTp2iph$RN||Ia2q{{&9r|I<$ydgUB!HjIa<+19d& z<$Q1-|E&&sMqip(+|o%h6(@N>y2bW<%hSoqPz$-8!_2aY=<`Up%W<%-(QmJrM zNex~l7Q9q~R3pvg+k)vV!K*4eu0u zdpDmbM@f){mDPHrAc5Xc@;n*13|POn0M?u{Z`|wY%VRX$KQXtre_vQw$dcyct0;t~ z3(my{DIA1wPr_0_OVG8#W4D>EsnGQBhrux=5BS(Ya&q$LPFgxT2gvUHYo`gD;3bb> z^bEv@48S(U1uL+E!k&9?VKA|Sz|HsFL^f%j22*dESaHy(su)B2uQxFf5zXOuBn@>` ztbLWGwm;1>2#TI~6Q6LtdjUNn$^6+f)@L(yIxy!Gxd-_aSxA2*vJoOPV>~bwp>zcD zjx5A#2jR$N0k={BXwp@F3zV0apJ6Z8;4Nc4;lniR)6CdI*|16qprH=Jd_=*p&s{-5 zok&_ad1d8T@MxH%q@?5^DisHjC3ar>8U0V7tHSqMisGK>mKeVvx_OhSrKLqKM>!?( z%5@6H@v$*kVAnrvs_ytG3nFml8W{)msH@LmGT0{n6d0@ChXf{?r7++2&GQf%-OM^ z&P>(1y$2&T#@!&m`o=CAcCj-g+9xWX1%qk_@xgCHLxfpsa*%?9OSJ&m2#;Sm8cZ-) z%;I_mXe^6MD#qCJJ9+vac~Fj|GE|?>f5(Ch{CgO+fC_pOfEk3?6xEe8R7*g=1mY(J zy7dt-TQp#>{ICygt*vi>(E%Vthc)9lX=!a8gwGBj_#C`}pFe-<{=GOM5NHDe_+;t} zCwx5siLUB@d46ZswT>xF>rwtN;K`pbtqZ0%C~809%#74&f%T4j<)oqy4jl(K_X6}- z6f)3&@JN{D1C-%`xG8SG^k10aFT$6(!izy`BqSjC3RKRqPy}tUrjHM3B0{%eImtOq zF#_?ZEMfXsF_9+-bi-} ze3Xx2+S?q5g-oy-*J^o$588=0AhH6zsxYYW$e0iE%@nA{$hS~HUjQ~IWt)b4&I)8J zHDLXy8Q+s(h%K!X* ziH!bkpB%_TgH~w`00e&FAn*t}FRv4Qr`_Qq1HEHF29WXn2xv(TcFetI}5mIk^lt+tBuwg?t@qH-Ysd2MxggJbVWX z5LkkMgJ}`Q0o+?}KMyF4$}m?1d5DhPvn$tbKHnH&1rVL1Q#IdmmH7zKqu@(rC}EGC zSCF=VwyFA2z%UDRbGX!mR8%T3>2IWK2iOjoE(Gkkz1t^N4lJ^e3liw(|SSV%L>)8{U5_eR8Uk6+fXR#94;shR7HucGy zH#)$kqqh}sWz{kseVZK>fxwyXE5K&P$!pibHf`TN9q8jFz@_=w_xk(#rhzIf;Q9nl z;JTn7;OJ$eTodTTg-hA%_b!W6e+2aLtmJ>c6+Iju{n#!p-{cYWO-czkl^+Cbwy6Ny z+AD$CdeZ#)^GiNFa9sNrIQ^{5cnLTLT`X&?sp&b(Ah9U~I8J^NxC*?euy7)1!cEi8 z4tTB`aIShwrz!Ysy+!iWIzkY$H-~)Y;Oqa4p94a-yl{Ic3p_!Q!PC{xWt~$(69Dve B6^Q@< literal 0 HcmV?d00001 diff --git a/pipelines/fmpe/run.py b/pipelines/fmpe/run.py new file mode 100644 index 0000000..cafd10f --- /dev/null +++ b/pipelines/fmpe/run.py @@ -0,0 +1,94 @@ +from lampe.data import JointLoader, JointDataset +from lampe.diagnostics import expected_coverage_mc, expected_coverage_ni +from lampe.inference import FMPE, FMPELoss +from lampe.plots import nice_rc, corner, mark_point, coverage_plot +from lampe.utils import GDStep +import torch +import torch.nn as nn +import torch.optim as optim +import zuko +from itertools import islice +from matplotlib import pyplot as plt +import numpy as np +import torch +from os.path import join as path_join + +from model_calibration_lib.gaussian_sim import ( + read_sim_config, simulator_func, load_data, DATA_DIR, INPUT_FILE, + CONFIG_FILE, make_dir, write_sim_config +) + +def main(): + infile = path_join(DATA_DIR, INPUT_FILE) + obs_data = load_data(infile) + config = read_sim_config() + + n, seed = config["n"], config["seed"] + mu_lb, mu_ub, sigma_lb, sigma_ub = config["mu_lb"], config["mu_ub"], config["sigma_lb"], config["sigma_ub"] + + def simulator(theta: torch.Tensor) -> torch.Tensor: + mu, sigma = theta.cpu().detach().numpy().T + y = simulator_func(mu, sigma, n, seed) + return torch.Tensor(y).float() + + lower = torch.Tensor([mu_lb, sigma_lb]) + upper = torch.Tensor([mu_ub, sigma_ub]) + prior = zuko.distributions.BoxUniform(lower, upper) + + loader = JointLoader( + prior, simulator, batch_size=1, vectorized=True + ) + + estimator = FMPE( + 2, n, hidden_features=[64] * 5, activation=nn.ELU + ) + + loss = FMPELoss(estimator) + optimizer = optim.AdamW(estimator.parameters(), lr=1e-3) + scheduler = optim.lr_scheduler.CosineAnnealingLR(optimizer, 128) + step = GDStep(optimizer, clip=1.0) # gradient descent step with gradient clipping + + estimator.train() + + for epoch in range(120): + losses = torch.stack([ + step(loss(theta, x)) + for theta, x in islice(loader, 128) + ]) + scheduler.step() + print(f"Epoch {epoch}; Loss {losses.mean().item()}") + + theta_star = torch.Tensor( + np.array([config["mu"], config["sigma"]]) + ) + x_star = torch.Tensor(obs_data).float() + estimator.eval() + with torch.no_grad(): + samples = estimator.flow(x_star).sample((2**14,)) + + output_dir = path_join(DATA_DIR, "fmpe") + make_dir(output_dir) + labels = ["mu", "sigma"] + + fig = corner( + samples, + smooth=2, + domain=(lower, upper), + labels=labels, + legend=r'$p_\phi(\theta | x^*)$', + figsize=(4.8, 4.8), + ) + + mark_point(fig, theta_star) + plt.savefig(path_join(output_dir, "corner.png")) + + joint_dataset = JointDataset(theta_star.reshape((1, 1, 2)), x_star.reshape(1, n)) + fmpe_levels, fmpe_coverages = expected_coverage_mc(estimator.flow, joint_dataset) + fig = coverage_plot(fmpe_levels, fmpe_coverages, legend='FMPE') + plt.savefig(path_join(output_dir, "coverage.png")) + + outfile = path_join(output_dir, CONFIG_FILE) + write_sim_config(outfile, config) + +if __name__ == "__main__": + main() \ No newline at end of file