From b1b6d4f58a7c626bc3ec8d20a157bd7fead63ddc Mon Sep 17 00:00:00 2001 From: Jennifer A Clark Date: Sat, 14 Sep 2024 15:16:06 -0400 Subject: [PATCH] Moving average (#381) * Add moving_average function for visualization and convergence testing * Update versionadded * Run Black * Bug fix bar_.py states * Update Changelog * Update the docs * Add tests * Formatting to align with Black * Update tests * Refactor moving_average to align with forward_backward_convergence function * Update tests * Update test_convergence and lambda tests in convergence.moving_average * Adjust convergence.py and tests for codecoverage * black * Update moving_average to block_average for more accurate descriptive name * Address reviewer comments * Update test to align with changed handling of dfs of different length in block_average * Remove incorrect popagation of error in BAR * Add tests and error catch for ill constructed BAR input, u_nk * black * Updated version comments --------- Co-authored-by: Oliver Beckstein --- CHANGES | 10 ++ docs/convergence.rst | 25 ++++ .../alchemlyb.convergence.convergence.rst | 2 + docs/images/dF_t_block_average.png | Bin 0 -> 44332 bytes docs/visualisation.rst | 1 + ...emlyb.visualisation.plot_block_average.rst | 6 + src/alchemlyb/convergence/__init__.py | 2 +- src/alchemlyb/convergence/convergence.py | 139 ++++++++++++++++-- src/alchemlyb/estimators/bar_.py | 34 ++++- src/alchemlyb/estimators/mbar_.py | 21 ++- src/alchemlyb/tests/test_convergence.py | 132 ++++++++++++++++- src/alchemlyb/tests/test_visualisation.py | 45 +++++- src/alchemlyb/visualisation/__init__.py | 2 +- src/alchemlyb/visualisation/convergence.py | 122 ++++++++++++++- src/alchemlyb/visualisation/dF_state.py | 3 +- 15 files changed, 515 insertions(+), 29 deletions(-) create mode 100644 docs/images/dF_t_block_average.png create mode 100644 docs/visualisation/alchemlyb.visualisation.plot_block_average.rst diff --git a/CHANGES b/CHANGES index 993ac8d8..62fd9cea 100644 --- a/CHANGES +++ b/CHANGES @@ -13,6 +13,16 @@ The rules for this file: * release numbers follow "Semantic Versioning" https://semver.org ------------------------------------------------------------------------------ + +??/??/2024 jaclark5 + + * 2.4.0 + +Enhancements + - Addition of `block_average` function in both `convergence` and + `visualization` (Issue #380, PR #381) + + 08/24/2024 xiki-tempula * 2.3.2 diff --git a/docs/convergence.rst b/docs/convergence.rst index c71e65ee..a4a96056 100644 --- a/docs/convergence.rst +++ b/docs/convergence.rst @@ -75,6 +75,31 @@ is, where 0 fully-unequilibrated and 1.0 is fully-equilibrated. :: >>> value = A_c(dhdl_list, tol=2) 0.7085 +Block Average +-------------- +If one obtains suspicious results from the forward / backward convergence plot, +it may be useful to view the block averages of the change in free energy using +:func:`~alchemlyb.convergence.block_average` and +:func:`~alchemlyb.visualisation.plot_block_average` over the course of each +step in lambda individually, the following example is for :math:`\lambda` = 0 + + >>> from alchemtest.gmx import load_benzene + >>> from alchemlyb.parsing.gmx import extract_u_nk + >>> from alchemlyb.visualisation import plot_block_average + >>> from alchemlyb.convergence import block_average + + >>> bz = load_benzene().data + >>> data_list = [extract_u_nk(xvg, T=300) for xvg in bz['Coulomb']] + >>> df = block_average(data_list, 'mbar') + >>> ax = plot_block_average(df) + >>> ax.figure.savefig('dF_t_block_average.png') + +Will give a plot looks like this + +.. figure:: images/dF_t_block_average.png + + A plot of the free energy divided into block averages showing consistency/inconsistency + across the trajectory. Convergence functions --------------------- diff --git a/docs/convergence/alchemlyb.convergence.convergence.rst b/docs/convergence/alchemlyb.convergence.convergence.rst index e5512074..c40e3ad3 100644 --- a/docs/convergence/alchemlyb.convergence.convergence.rst +++ b/docs/convergence/alchemlyb.convergence.convergence.rst @@ -15,3 +15,5 @@ All convergence functions are located in this submodule but for convenience they .. autofunction:: alchemlyb.convergence.fwdrev_cumavg_Rc .. autofunction:: alchemlyb.convergence.A_c + +.. autofunction:: alchemlyb.convergence.block_average diff --git a/docs/images/dF_t_block_average.png b/docs/images/dF_t_block_average.png new file mode 100644 index 0000000000000000000000000000000000000000..039f4a7af3621e8792d2bb032fc5354e87304069 GIT binary patch literal 44332 zcmdSBbySsKw>G?KR8qP@6i~X{beE(^N-N#n-Q6G|f`F7XNVl}KlyrA@!@IV>=Q-z$ z^PcgJ?~iYM<9qiQ$gTUn?-g^cYhLr3a~b+pUK0H=@nZ-Cg8oKI>>UJxKna1stD_== ze@S)nvjM-j9mLfgl&p;$oIltZLgYR;*jQLQSeSmKbTYKFH?_86XX0XFXP`83aImrG zVP>}cj~ke*?Tneb*VUE4OCH%sy|;%zFh9V4;PM6YOd*h>ls96+$}VZUbIvYW<0*oN z2e_t2W9|nP^|NJ>jL`Zw>LFCPI8PB8d>MS-{#;Lrq?(E`j2C*Wm-?YikJ5*Vkoh^> z=p4!GPnLLv&V?5|zx!twC~!<2^702~O&5kvHp@5q>(q?9jyIbx(oUM$gs~8@uwkD! zLi@XLe=lP%^dKT)uut5F|F|G_fqkUd(1iZ)yO1`<^nvwySN|)8%PPTgwIn$LJp31S z@@M8;2>E5A>a5qSfhmAcTu8n?pZr+O1 zQsc32eSi=xfYx|dU!K9cy}uu%Q5+;o<1b4iDlQ&a8O%ZJezE_>1%1?WZ>~YjVYib& zmv2R%x%u*AUSo60GK~Xl`k^_K2n`ExQRJ!h8S*G;d`7n#RmingIA7mNu}iqSoD4uS z$IWQTTaJIxU)~)Qi)mUt=itC?Zf<^^QFE3S9gX?(=g(Ze3x^?I5(KQjFZ?!%P}+o; zkXGd(PlLeW)q03UiMqTIMGk*`?K89QkU?p1bG`@haX)Ax7~(rY({bP9yd8SD#l+n| zhA(n?xVzZkP%3vjGYg->c5a6RfEnyT-}CBbyX^Ueq*{t#cMMlY%0lS#``O|l+}EOx zvC>tBvK#Hipo2HO>vfmwpRM1jtB1tIVBIcZdbbMRA7Lu1sjaL;3wDjG>F(_=jgZ=4 zqJniF`?;VHo4Mim3~9fXr1pfI(mIy)5)$FP&d65%aA>99Vk?OWe^wQT;V$20Ecn`T z8nH?va^N#Ge45sMIh2+%z88y(4E>9|@kyE=B_%gE5nBq{!-?R7{?YAJrLirTIBFUi zWV0$_Q&ZF0i+T5vi3ytNZeja>4QYT>fYQrjhfK(6eyn?S74Gu^3-6HEO;q2!goMQ5 zgr+qd#H#6RB2>_jjp+XVq9wbkN@eWfX5m2q_c7A+^z>mbZN?wIgZsN1MO|IY6xkq( z+eMVD4X1Zy-W{IwNwWk zA|)f!K3!#=Il6V&MPeOtEh#Cf-yVbx(L0(oSQ#qTtigEs@+HpR`1tr`8>W{}c{!_= zO&b!I?K0x^%?&FDN5GK44H{TNE{AO>d$+HeFfV*y0Q5%B1Pg4GbHRgc$I2buG7LI; z&U@`rg!{|Q+wQ|2eIcR6+aY}`kCgGt-R@dV>n1!D4AR!~y_uuRjJy8jx@f_BLJE&V zI5&59Y*JDNH!lMN#Dk`b^zoDXedm_;;bE+-tgOQ+;|w{2__$t!+S}v4BQT44N6Vc~ zu9gdT`_7g|^=_vllC~RPoHac650~0QLPJxIXEkfAI&uBm@o?F#sx)O}Xtb^kSq}HpANh;mjxY=W`gE!Iy#m%9zRADE>ovLf(Q!>*Pc(C(_YKvj+oHE zKQZqs>zV<*pi3nyC@GOpP{4#NW`#3Va5K5%UEORKmQtenjhCX3yR@#`5@ekf&l+2SZOxn$31kthw6@k(pUt`> z{quu`==%El(DUd2D9Gw$wy{W(xlK6~=+bHM-Fw)oZa%_lAH<|Nd6z_i7l(U{5 zS6>3(nF-ZIPb?)pqhUS2|5o9zBLsbz5~m;xEHJ6ngt^F{v0d(SkMExD9eoh5oY zQ?W*o?36Iy@-zk*uxjI`%Hze?G>0csd@l|hqGY#AmH#gL0VkqA4n^aa$zfM)MZH^& zw9M_q4>uRh7yNXI%QuILii$Yig5LLy0|%Vr_8jd0c&FJAZ%CO-%9ic>M{hOwD)nW; z<%r54dcABbS6Xk}-FXFW4@DuZ0QQ3K@4Y_jc1Kg7KuCEUP&Wp@;CwOm2vvFio~Zf$ zdM8wnRuKALLxWgEL<9l>5MIUIT#%lf@{dnmHR@M7Xg)eV6_zqBFVEQI+u7q~a(7vY zgtEMT9kjTpr)59DiHeF^c2)JU^Ra=kF+y1yQg45M$FD?sI7n&B{dLE#q?eb#CnjyP zBToP>%{Oxg@DLrBB_s!DXU)@s{QTC2gJ$!|$(_i5**|4fS#O+pbY(p~8jiY1;0XQEzR1(MS@C%u z41VD~N@}?kA2S*d#pH)=l8Fgbi$7&>dGcT2^=2=Q){=CWf6Z3dKWfS0o8-dc)QF=c7e`%&hEe1Tr$Rn7*}30IAj|T3cINE_Y{>M@r**ppRc-CIO>bOXxb>-P&&51Wr$gcux_hP#vreuNBY#N({ z4?af-q=jo5Ja*0bUAyFQ<~FJ2v?fpPb@|~I28lUII$?kh{G4jHm1*Rmfw3|BH!r(N zWuOpm59Qlg&IQigX~!xHwf+32CNkz0i`m|XVr$R7-fUZ4EkALXeW-4WaS^qgsQja+ zT&}nBj{fLJXI50$eN<{tE$oz}@C-QHnanOOMvsk+&H5S4!J7t01)4d=g7a<7NVh+y zo7Hr5FJ5nCS_y8GVghz9S|JM>-ud!{4eZbj|k2=%~-yr z!~N6ZEtr~`Iy^P?2&|;BF=Z{+wZz(o`?1diti5kEVSSQ}CKuPyNf&=eNWit2tH*;Jf<6C-uN070S=k6r6q&!xLuOr&4cMA= zqQCGuf3&f|LPA11;hAi0g~Y|j$Hv6?6qYs<1KdchXlc{)(a8(e$;+6$cmboI$H-Q^ z2=EZ&p>%O!VK~60u&}Uv<{F%-p-|_+jlRLbt|&6zGO!s>&(1#VOcbkE8llT)2t-Ci zpcQuiP^-xqLBar=zOTp|as(V$z7PElN?8Jo>MW;aidUTwY&)Z1Y2c_^kCM zt?ZPxbWId%9xaD6Pp{qD)sfPx772+5?zM_AB_Y@ev7kC#rT?9c^nU7dp9ma8IE~_o zk_~zW1~ChZqPf90urF5!W}T2%GlR*0*0kJNKBu9{YH6Y1b>8ooZ*ng;9wLXdeqz)d z9vcgo)b;jUP6k7PbVaW3?d{PqFa(sA9(2SNFWjy})vL_#*x1;HCntmRGQ5Zb0|Q<5 z>o%Tyd;Zgb9OVgerG=H_;oO>u|H$jAtL@OO?(x%Y#YmdjFLToV(_WX}}VtKoV8kIVD$ z0!tC!(vRg#Ad)u2{p&VDLZ1H<3H*O^k?n8}fvxGr5fD!HFn3r+^C2T5@)+}pz59pK z0{R{F561qF7oUeiL;7v>>A;v){_%VNt=~N|=olK+U-&T4Fr)b6*wRD?po zMzpwLeLFfkBHFcb>ed|{3mX>#fbhoFeP^Gz%H`@JtrfPyG!G#r6g1?oKTx?mAawMs zXY|}pcD)E{m)w6xrNElT0Zn77sZpn29iAasH6ip2;#|b{`8JAomHz&?7@6cdqPSuY zfap-24Om2B&mYPlw#z&yQcc8g|>y01N)bLZ|)p>DnW~a134yrf%?3&NP%U`2 z$DI{Ib*%R)aqFO6v`e|6GINq7`ZHI?w7?2z_Sr!vB%sGXFd*fg3AKxdqEb=AfX`~# zBlEw$+G6RW5O=|Jw5#;|7HP_uZUIL+%t9stV3WDBJb>tqZFHsLseWu*3LPWzFKM~x zBnI~XatqhVSc7#VpC=ZZ&LeGy^s{razul%?JJk;9`SKzUzHkB}iAh8e9R)v=0mUa` zP&byXqE8}$<|n0bGBu>~oxni~C*Z+UV}1ec@?pVNhRqz!*U}+Ey$E`GPOn}AFtT>h zW%P}B!wqS6`TtxRxQ`dBd42my`!*2a2#J4?c|clq7vQhm$e;qSo5*28N`MG!rI+-J zn6zhBz(h{XtT#N*V8F6oix}P`^P>iSEl261ZoQ%PIMF9V@Nim2XvB2Xel4ffx6~(` z$$pnB6ijm^*u6b6q}+5~MMsspmAfMP-?dZFYLDC|s;B`M(L(4$DGVdfG5=F-t~ zz`eWoAFVAyGFEFfZ5>l~?39j7wYf%iSkW zBqR_f=jp6Sh!LHd2ge%0yhN$0AvE3*A9frt{@uITm9yj~IZ zokU;^hAa?CjF1_8X7=34FZ6!H^6pZpgl_ z7S%G+6R0;~^z=eYjQ*`8t-K1fePQeiVj-x9!Zhm;F06OTrmVwrFRV1#!v{91zZpZ4 ztgvdTYTziEp%3;r019oHyH`^%b#)h`8hm*V;EwTAKE3N(>*{u%N4Wi$ghY%pQu_?H zbGOT;)=vVSQ?gD8Y|I}rf7;-2Qg0N5MG))jWmAuXG!GyZC%FPNw3vT9f~Q>1H4k|w92*aEDvP7IK%Tha%{!cDro3TP;i7Pd`fB( zKnJVA^w&mJ#5HbnyHU)^W4d|l=|khSmXJ8W(z?mr(O`j1C=RX7R!-8#m99uSCZ_6X z-NNx_Y-~?41#emKnRP6xhGMBD%S=aUVDutSA){-3;Ir+h9CYbVxqg(#OED)-p2r@Z z%Ck-%d|kBjn;3LWiT)_7Ky*Eg9u5HPuA;lq>@6hxTaa<-^rtV8wV8perka#g7#S@Ce#8h54{vk2saSZY^zL1mrGc?jgIWpzZcd_;{%1}p4*L!`Pr6(Tci1$$N-7?BYUCkk6bo_KA7?Os>W zdqW=_Y}8WpzLa1gM7dy(oypy92eJx;7ZCpYNay_2cAg9T59VubYAWpD@FSm*hldz2 zvDeY2rB;%?NjY?a>Z}_FG6u|NRYFSswFBe@&Ux|mklJ0%#fSUt%L6$u$!|n{ z%6bzZ>VJuzqGMzm--c&Udsz!#gCSH=Q2`TNK@h&&EqKpr)v)(EAt6E5dlA>b$kY$0 zQYN1n*chMIC9PB1TXr=zGsp5LTMJ;Ke|BAuxT;^C5wiPUyVs$pkBk8WsFc8PxV3NZ zmeyNl^4RN8%io#uWfa=(KsqKT37+6Hv>vbZ9RU%PJe+_T#+^4jnv;_W!>6LEt6zt| zkL_Kxq+tGH9=V<;hX?^i&Q35P_zKd8 z?^HiBF!AB!bS`t;%6uXk+Ay2tfw0n|!vBbO-tT`0!U%yI*|g@3Q3b)y?eW45UXEwM zyk25y)e4y$MlEsgTo7&UXl!SCSd>nl0vrx<-JXN$qzj24T9;Fqi%|BplTeop-?QI3<-9(`J=o5JiYP*i@N*>Qyl~qHNlI4 zKnpdE9Zvg%~eef1tx6JZ{(3wi-F$HCO}SxR4r)Dz@}a{m~5uT1iim! z<~G-IGSu;L`K#v~XAXtTAX63>7yoH=Wl#P58K|le07^a#rSsP=htWd_pI}Z*PKHOl zuXot$SWrmk3$2-Vd(zw6n_pav1S!&KB8CMb%-&Zrx7%thB5rOx^R63YMOt<5DkjA_ z2sWgBtSSxV2OD>6ftSzK;LBAIgvU34ducjI0;of514u&*Pa31KbR{Mdw-~VLbZU7*%BdMajoLc=^=w78m=l@e zF~NNLt|AZ)8WklQ)_&}4pDovu4#vb8)vbOn9K95l1`cIZi7uk8!=Z{VL0a!hda?Nd zW%iJN5<>t2O$RzhHZ@AeXelvznM2)47yFc<8lV_LmQ@9NtG4>K9ZyubthK;9_Ot5k zw$0WjJthe~YCl%fYl7_1YFC$pZ)tScj4Cdqxp|N%y6$MKab2cle&LaW6PKYBc+4ti z^K@*8xna1;57P$_M@@P@(>mlIdSA|(2cP0&Xn9xGS{iiJQizp0|B3+vo%L&D@q{c# z0XA1xE_%}$&2I10L$?uO_(~>;Brlk@uny?OAm zjbnJR1M(wlu?ec#4&)Y(KYEss=kOJJ+% zOmKM5JMrWtj7Jqt;D2pWnf){XQbB2Jm)Jsi?EI>9zE9M6t@$c#dyy#_MQ5G15f2A! zVc}Is4;~&~7zIi&#P(>AF!Zh*BLw9)hS`T6T^$GVs zU%5&kLxj>OoWFS!K7M5MtocTpD1Gkx&z_J1jpdRr@E@~&pgwWFR}NucAhA=pP06fV z972~{JB2PUoJ~#fL|{a<-5Ph1u|kemriB@QVXAPs6RrwH7Q=1s%VSq{qD*8SZypv6 zvS+D4da$d;(~1H$UdUO^CO+{mJJ|Ua+MS)v)>uC0V_jP8bp4LaxCOfEb;&Mj-%(Cqmc`S|8_i8ADcljv;W_t?ekeVtEjT#5IcN@Sj)rZ+}Ve`NbWmFC2PdrK!?laIN334m0$ zQzo(t1Uz!U(L_g*3fZ*g$4ltzvJsV?`$^<`IAoo#80(Kf7k9oXg0a3aw~&oG*@Me^ zt)S_J1emtn6S6rohtm+Kk#U|Ak10pgQt3#h-kbP*xjZ$!hhD@@bgTNlbijS^Q>Wqe zSdy8B4B6lwvEG5+>cbvKFsxe2>b_4^d@klbkMhn=OOP$Tv+oy=`|Vqor~=o}ke-wA z6vmb?ps$xkZzQUNto!P+galr`DH?v2A&s-UGEZ~wrwu4Iolsh8P5v|ZfxY|JQItTO z^>%AQlQSwf=eNDWYXt*bQyC8y5;WyeWn!nKo??HUp}CTcTjup1VYQ;y*#jV+qy(P) zUo55ObiUM=gAEI&7^Jg)k6W{f?Uk~5EvPAp;Y zCnjEkx)=Z)g;y11T5(izd{VqT1cO4%=Bx2HN}x)c9UQlZsCv9lU(M_*g$x3mDe|SF z;7GNU=#$^`>knu9eO=(FZW)`Z7IjC9lR4==LLw|CMn2$irD8IS>UYb3Apyf00z(Eq zllq=>#qHSC;gugpTr%%01S{Tka4OWhadLISf}T}4#5Oj_(dwES1U!)RZL;sJBJ=!C z8D#;9fb_wvZxde&QkySzR`v2PsWfYSGZ53;S`g2n;=WkerjJ{!YyDeZiO4Zth8GKjb?kLo59_ zejsKxe){{s-6{HIR-P`Snaug&zeZwTLZfy1M}a@6gK$W$6G{U1f{`ip*elx)xuxMB zM@I81Y;ggao3rV9<)}uwir0Grll5o z5&o`^nvbz}uW}bwfTcmWq9twUNq34o>Mc{(@6*3(MH|(_`9chna>nAlhI0n;LV^UZ zDS#%iN_LJcv4*KQbA=%<&Oxfg;qER6F30-E-BB5$`bxY!yf-XtY++?(LIJ@@A4Yz1 z6xFkkl$S4!W<>PJpu)dZrJwgm$DKrp)G%H_J*eg@MPeko&#J0?p@JO}PAcSa>VqRG z0V+uNPo^*_a>!X@{pE=xf9w+bkn$-9v^3r!oqHT}O8n@G;Gd`%jZk1TGAWGVP-=pn?K|K;hTBPkE*{+c{*cUJ?r^Bg*KGgh1;?u< z97##g66XozIaMyj<2vVWjz44|$Gyh1hcn&tBCg`ph;VOTt@RAaw8eEk!6L&88;*hG zk2i1b6A39O2>zlcw%vBX`}K?Ds3+Tc(K)&$a-c>l3WcJU9DAv4)LZl!rn1o|Hz6n?WMWYdC$=x8amfgrl(Ecd7rqYbx!_6UkKy}tsg8d6)-^$yR`^I%p4b?0i=H!% zYGe^C*_PEcJA~!zM&ql)MrdXRuUa(9K?_`7ix)Ptzy+k)b1jSri2@Rn`qV#T-kJa0 z`v}o#@0OWCCpjqYcl{McXwP24WRS4WmVb|AfuzRAQXI)u0RuL@v)?rND*9RHIGwDU zMrJy3Ax_F4MpVbK{IH7+d9h=!jtgt^;V$XWfr(Go`4Nr>$}Le*hx^?r%_mp&_fkF3 z(9ni%w8ZwsN|Q3P7!$P-5=A7SW3^;&qFCDZ=`d@NAnRv z<*Pv(+0~>SGa5Ebd8EB=&$$H5Wb-?u2Zso)RZW*RvyrKPnL7pT@FaNrPIJxi!{zqQ zb7U$pok4(z4wsZ;I-lxBs2qi8%y|}kv%957W`769#wyGHoLN<#F`xhSnumvw;hG# z1t^n9{1g7OI;Z>vYWHAQrrwcjsHB`&Nvj@&xIb>tMd`fI7UPK!Ps8LCAm;%vr6r>%r-Vglv%+k(dZ$q>`xh@Gux}k;YRAY?Ii}=BIKt#6jGVQd4 zL?V|1*a5}g^^I6llg^nUP1f5d32M(UD0vlT6_`+t@rhcb!+oJT$P$rj=Gf^7@O<9kSjK^b zh{!kb7C#m4;e1cqLz{ByAmP>D--$&SUuPOVcjmPfOBbatf(M=#4XRIQ#yE zfz+xnD31t2jH2Z6 z=y{WUWFccejfKx}&-bbjUDDJ3Yb13@5y@zZm*IELQ;8R)(1HSUCIre^@?$p5yQEc}V8+?N9|8S-UOup-emak1s1dzjtLFS{vh27Te@cr(oEDfz6R@fyowO!cl zly=;JZQ`9R4{Vr*@*f^8eR(U{C$%l{&Z0EGxQ-w4LIp5(_ITVGl)Bkd@3TYSDKc+^ zciwfgi{s9My;8hpDUXW3Qw0$H=8IW7E)*!$xmt`lLwOtsZ{1vt1RM79g8*N|hR&tU za>m<1hjh8Y#^nz@ zN_=STb(ZucT^)IN?x)1xoY+2)uE!gn50ZIMiAn>Qcobe)^SMxisMLuDW6e>CAd7== zHQx+XTt$?QExw`5Zp0P8z)>IY#aFsyS?wraIEW_~56P12xO;mB-Z-7sP}#X^H<|r# z3{_IHKNocav!b3n;#-NAyPWTd@G6ZiOuQ}bM4ptW&cCPZg8_}RoR(13e@nkGbI1%5 zCuz&2UqO8O6PUz;j!UX*SV3sgEH9!TZg@ukN^YnuT~P6Fe*fC}$?{EQIS%g(O=LE~ z6%vOvcX0ObNfu7@d91xrB?a+#@}xo=PhZ@wYMS+hq;d31S~FK|!sV>%SkjFqQ*2Ws zRVJQjs^WJJRPR;IZd12Ko@Ig@G)T5pkYJ`B|n>!74i5sun{slx0y`J z8A+8@4IsI+1!0-NjvGWo)c4g;GM*;W%8l5cE@Ad6TkxlEjbQU;nt28H=rSEsRx0Hm^wLIg&f$FZIXM zEDTcH>$@qRO-&v>HEdRXbqWZm&2osAe>;54_wAF*m5hbXB=;vHlY)-EoS;r~z22FW z6ZctOVG(AaymkFH@b3mR_E(AQNJtSM*RqUAc@kMha>e!S_h+>t)|(p1!~Rb}|+~#rC5Yj{3PIE3+F31wSOv z?*^UhsZSbXNm~{4me}8)8VgaV_*mG*1#8ugc@Qe0BZI{N&7kc(Z&S#&{~cnY<@OOX zRRcjh#RO2x*Ejb1ba1W7@^~`qbm)o1J zHzgm$3Wen2x7cauov#=5wdB)ny+3N3WQ#r%lEtrhG+f&I%9KWuKJ~fX zFRx)`VexY`KRMbgwU(53ilc5ukJU7>x~R$bN*zGsViKl7VAA>%5cElCvv_u3z3Xj7 zZsTjTX~m<`p+KTNd^jnIcbKboi(-nqVw$2dSoJr`{azd5+n@c#Dh>Z}Opdvb4!xM{ zW^UA!mvc^}RNT@^nm(HgIv`x8O1(jY1CCHPgTcS0g}Bu6f&2DuEb0bW&sl>o&$&%R?Tbk=X4Qq)yDKcA7mJb1K}=_fgfX!ke=lD?RE3U@f4jG-uMkVLnySAZ!ujyzhlS(j?F%DA~$lphSn3}m%0mo?1vtoTIz zDQywvrsXdY@FBpru_2X`nCzUkhrG|_L#F(QkphPuFy%>v#?T{@o>Qn6Sj_~qHYXUB z5}g1w8UqP!kTlMj8*Cz;faH)V!^*y)O)P$LG9`6}(dR=6iOGw-_5Z{m_1i8^0l-Z{2~~UQf%JBK{jeTNKaAINE(W3 zr=wT<)SN{01b^Vx#R}~V0~wA>iTs7x6HFozt+Cd$!5WrqZCoETbI-TqSU-P)+t;#A zHg^@XDhlWam-yT>{8ux{U7hN=`B!aEgDa+UD`LPA0bAs^%|({d?=@c3m-oHx=~4CW%w zJ5SqmVQ|}n4;IH3D!fguEZ+?_F-Js4@)cuV|Int8jAU!a%Dph&_7Doyi*(^4+Coz_ z;f1fYLY1;KLRh&`M*awrL>{%+=6$8=9)IFQ5n?9qkBQP49jC9^!aJOtoIsxvBPB2ilS>6HZ2MtJThv1Jw;c5 zAt|pxd-15ww;4`1dq!5x&;d%MayZ?(bb6J zC^zUZ`y!7Nv^oUo88bb@*1*jNrhezBVy@P+zK0Y~`DG1Ccda zB=AxS<9~?hQZpgzH5ct#@3biVZmTvwCukP}A@|6Bc2>SbrZj0uw$6}v<8l+uL=q4G zuH9b}Nc15S@=C8RfU5fK+PE6=Xuxp$wEJJaKF#(S}`+-4N7khhHQqk7qEel zg0)!M#_#NjGNXHo31ZWv8i90RP0pQR9A`(*O*RK!?@7tqI z{u;#+ptKxT?mcqC>$K}5OH+ACOhlA&cG+_S)R)kS<5zfY-)PUdtc5tJ?YM1v!`uvV z)tRzW&GhLY&%4iUEvr%3_V2KCS&*;HW1lXV*BnqK!MPea3IlR5JfBkzEKqn{|2Dix zN~S_HuaRW{`k88HPc6Hga~QjwnwN-FH7k2YVXSLSt(tpP&Di!0V_bZeZLfHE-mWwF zGZv^rm7P_U_vDWrJ#yJC@3vX)c$CIv3!kLp_7ZlQV)^4vqNkaena#M0+K0J%M=r;0 zTynQfs#)g+&p%b>N%^{<h!eRyD_z?}?W}pe zpcxMs2?8FDO~QTld-gIPE&wsi+6!Mz;*Y0?f#5nlrALArQ(L;9C_GlN6Ixr_{irh{ z7GZu{zSH^DNon0Wv6kbwN>6{kn2`~U*UikrGfGMbD2T0^eh4Jc5qkITX}a6?OHd^a zJLJ&XayQemq~LuWpP%thaqMQe3D4PpFX6;>$aF;=>(riD5XN5tYp^pBe~59&ybEMw zvZ^z&DR9+}vRdQSKi{Dql89@FM9;M(iBM5G`?}zN%Y1>#@A-;EGNyrBjl}4v;e+Ym zSS5vz?Uix*)R2UEo}h7$R`S(-bVst+#6v0=NzCOXY8sI821`JR1P_n9qUwat<>l2K z{N&ErpE`S5WE2#|w{P|SDdPChj;uJ7sHOudKIHW*`f(xXXlSsrQNWh9*V~1*+hYZ; zk&%s}?;~rEU6n5>PrMpf6HC{+hEF@lfup0EXbch+M+}lwA01cv;c%HHs4Du7;|9rax6gSvY-=c!u@19^+)z&pn(cP zYj4fWFs`JUj>Fc5qqac}s$gi8{EoTN46oY}>@Op%WX)9{I@6kmz=17nEG$T(bYR1N z{gVdg0|_gulJRCxClL}AgO5g1Fc~q9Q7>#AR z*LirL^~mqxR;fUP^z?eVyRcW8qzl_K$?z7^)Z^w}udLwN!EtXbzmKV3rv#$@G?!mm zPWg+sKTwpZpp;BzUT1!MUqr2D=A4JBE?Y1$mxg0MGHAL&OX}(&DMXJ#B2XpVC7FT5 zd{omy(H7!_c4WI$d$AvBe=YF-{d-UWEC$X4Iqgjc&R8{|gHw5xL;*eh>wJOU`wV2( zx}YRjwOE7m=r)zp9Gu@GSWGcaZv#aTmS<^u<9jBB)kXFq^X-j?Uo|*`zH!Z89V&ll zqgIa@^Z(^Gnw}9>|JESS(Qcw|pvv(&^f4lQSC(&-jznzyL@2x`?f=O2gi_T6n&Go&C)l(*2dXo#N=NlD%tc zO14rE;8#Ml4{nK_?+obWGYYzNo%aZeok@r2DjjTjs`=doMMO6~O||_~Is$<(=r)t; zdf)v$Z9t&wF_L91D~qO~(Lu3pX=(ZPRU#@TX3ZuhNplk*ECH8RS_FJ!E zC&6RcY<)Xr&VZ0<5Sq_@oXgY5>9!;Z2YKb1AI}@`U@6K-T^5m`_)gZ^#;=3D)Y(^D zmalIz6(gyr3ucxU)QW#oQx}kI&@iqsls;+C>US!1eAaW8s4On~m^i|S{f2dLvfHO9 zfeBY7)$29=j%eBQLOfzl@7G$s^0bxtOYLiFqxOG(uKgHO_qdMkE_+yf4-N|xtBh~3 za&f?W3v&J*UAF1%3@^7D5NYsI5XSmCuF-<<9#sCf-0fLsRVxVI;r1ml2S-Pr6<7&l z!HRYQpIAmK{ZnjSdyhmrr3)fr2e(mXq7>P()WJl)@D_>SMkF#WZ6gi0@iF|6+-gX4 z=?Y7h%HkJ(?yPFLg;$O@Pu9+@&C!`P+2DV{L;1w^+~#>^aDs*B`62BQxaOc5TbQuMCNecPT3&IS;FB|CixukP_hQ~ zmFBJ-eu>|iWa_xCI!EV#EdD>MaVLY=OlmtXfm42+{+I$Gya$aql2N2>uv3EPGu8@e zC(xEFbg--gLBT&b0@`c}wqdIKTW+TU0=;aZy)ya>-Y<=m3$)kQf8P?JmC<2=-N>$=L=@vZKLb6mx2RvrdVdD0a&8ia4*rz4XYGOGhCiVT z9{MNA|E}pT>JpwWEK&ue7x$~-2o-aD`dH-eZ5xM+6!;OP97U1R$Hj0*qHCssIai$C!FQ`HkGBF-HGTkE-NyzYjG7IRUAVeH^4;ohsNp(tCF5 zVJ4|MGQolp-5u_pgg*Td+VC+qUe~if#2s_lXF%#a1vxkoUp%-%iD%>P!DlvH3z?zM zylq4B9nEx5l>l8y#NxCkbYw2CLR}>G^ZeEih}{Q?20924`H`@rJ}O%rSqoA6y0Oil z?vGt&!tbDdw~~8y_qP!kiV!xg-SFv?HYu4hoqZqDp<5oxX9gjCWyU>Ba6QwUSC!3& zzazJ7Z^ejC-4K|kd-jya0z>V*=<$%h(VhX;waRg%h!7IJEjjuPyLyGLq&nfV!1 zPv8bn=$iaKsvIEY0jC#c_5@EwjFP|LvwOJKY)XX`AzJPhYTc}lc*?oVCNk4XoDn7h8@`Dyf|%5A`x$zsA$Abvwq?OQ9{oCX+NS~`J)%$Tfyo1O;+_Az8e7Q zTc;0wPuQo$;iMM6K1Lr9mU0@s7D-R#wybCt)>Tg4d-x^0FX6M097%1YK!b2h;GlOf z`3A=l&~}cR8)~MoKJT1~`C{xBRB*srK?d|F=*S38_m0euCl;a~E{6z!7Vo=T^xVSx zNbvBg+dIaLQ0zo?fz;qAb0y+zd6fsGMe$_}UhRkuj~7-04hp}r8@x499{nyl?FYsT ziqvx^GsuUwWj{wy=0a;`61&5#gm0?&$UTz`1UxUg&Z&lV8#9~!6fIJ{frWBhdWdbl z=kKjH&`HUCQ2l&b_NSvc2?Tvm?pIvL#%rk9so|I{Z(g2zozrT0j!R2o`u3Fcv!{s# zq7QPn@VZ95m;xS!W?>|8NlX!qnOIwS;VZLQs)(57RDt%}D;{qeB`^u>9R#Ic9vVmG zO7dzJ2*!<&Oe;rumh!3}p32tsvuRFWoeRzH#w3rgsaf%eZZTk>3&B#gJ4|;k_wTkU z^Gi&H-84w6`RysrJ@1V@Hn!mN;u& zNFLEfh7@gWt7J}C&@j?Ajb;kUNC*1ta@b}C03{JrS^VMh+a2&YZ`Wo2m55boFbP>b zpdf4SRJ%nJtiO4akO+PDlJ&7>MlDJ3xj-OEsV*jWyh=O?2}ykRFN;)u>XfUiUai!)ZD8tx;0+n-5Eg>>L} zUpkXcwctS)Tv!n>k5N@3sfM@Ykb(`cVtQw94QLq4yFVh>F4dVtL><`=(2KQY43|7 zv@|VNDM5B%m}`qb<&6p=s6JRK{V=@E&`u3JQ<(L0gXfM74sBAp6bKbcytdWP%(9B> zsR(5%W-cpN#frapWi4m>*bYTS2XGp|-1O5)`d{m+8e~C@rx;A=PLwBWYN4!m3{QDb zoE`~1U8b%*kK{`Q`5a1b#*ftRFGDRf(`Wp0Ct5=r;@$j77ZSu^xgFZ*es??-2~cx& zE=q&Lz?M{(m39iKOiGs()+$&P-C?)*fA2{Tf6A{~*p+6d^vYT|82F8m>#yIb5K~@i zk?*3mQ$PLzIPo+HP^@-=JkNLC$S*It8W$9li0{k?21@8Kn&5{Z8JjX2SfQbGs4++A zH+zZ^5X_<(HYV}oarjg;Qs{ZH@*W*=BKecJ2F2G~5>l>~bofTqqI)|$DGiO{Ai=7^ z)3H@WU|20$kl-~z#bIZW$ zq29fDz`>~fbZlO29(%=)pK6I}xp%tJ>DvKfxU|CZTW``;!R>F@13F(H2;Z%^i(77d zZ)eS2v9x;0S%vhqNwX0VD3G6(m_=5ep*ZJt;s+Not=X5K6NgOk?EuV{V(V_EGZ@*5yi8OLR`~S{2?CFARI8i&)uWQtygMO3 ze}y4LxoE%(@iLOfQ$mZTl0jo?elNJ1NsvTxSc{_yS;JsBJbK@GIb?RG5g{f+@bu|m zM$J%HA*fyWSu%&RYsd7=okgLEgovlZERTc!aiGr=yUxx>H^?eM?xDwr~_(;?5 zmcf5ku_*wi6L>JG5c~yO#*zr1`PTPXFM%qjG3KovH;$}E=lBUoEwWK(OvMnJgqAq6zHe)bh<%CDnh}Vvr5M z&&(ezC9Q>&8`$$M_XR14&PdDN#u-<9Q9BC?U$9N#3~OFVSRE*M$4i|XpZ0w2gV@)j zeZ@8T&_;JpD|F2tA}ogO0=daE6;G*`Gg2Q_ZD7sm`~1n_{ea} zPlxmCMzSD1y%gUfkL*JVIA}gG$bOHDp6IqF<2UP{y?rNiL(Z#Y)z;aC-MmcFxs z^bPe(UO&E3=V>SU^%O{W?ru9!3#zrk;sVgyuF| z%@7$Liy7i;c)sCTQ+=Mknwj;1@#@{8vlI-67e8AGcCUoS^Kx|9I(Hn1u;$}}gHDv3 z4cCE0M-lt*nW?ryDNaoZtHSs4Rh?`VMbB9~?gmCD2H3gjlW?lHUr?#{u1>sUY1UltYhz6M9hMHf_89{oS%m=JL;8xo zp+&QzA;GYZ&pRHFnw5PWZp#yT>Tua_@`>$Ylu>pl}elL<4@h!(P&_ZS6}fi!YzKW{k+YSM>4&o{bo^%f^^Z5-sh#O zd=8*CAkqW#Wep|1lpAEd-r%n?a(*lGp@U?)c30mUHM!6Gebm^#i|6>9^~Lyr0gZw9 zPxZN8$^h@pvd-*Hm*|s@2OGte9$v*lfqT^JQg)wH&Vob8rQTq1I_N}wM=QJM$r@Ke zzyC;`Dv+T>rc}{M*JIx(#jXDv)&SM$nfX=@pu- zG6(eT?+Fh~3AlpamHdMve+N~8Ig+Srrq=iv(A=A9odyAjd$7I_t*S&ZIeFU!ghS4wk_jS z06UNaN0YPLrEnC#!8L}-?xkypo}vOmHBMAAxw_SNbA}%M{4H2zuqPU`zHL_wq6;OmpPY{(KbFb8D7!> z)X?};_&vpym#yro($HzCW!9v^UuVSA#_ z*T%gb3Ed?2U0m)wqTj_Y<{VdpV1H$W&<7K67G9$qh^qZog{YF1>8;Z+pAy$!h1L7> zx>{(Lm1OR7nk%s?JiZOcK1kK!UXit^8FR14dQO26@6>aoo&P(1Li6O$qX!?J8Z&?vIeTtiAr~ z`c_ad&-=kylsKfz|F)}qRx=>~< z8dX{HrdllV4Ry6QYPVQD~a{Bx- zeyzWq&oaa5alL*;8$s}rECGAGKC|orQ$wU8&Ry-qkCWKH$)+vDwXR;NO?&=1XQ1zD z@Gpso4>r;9>u20#q2WPz z6{mtm3GA2=7$s}bM&gP6LA}0Tnbq$M z^`p7`RQ=4t$U*G|YMUqCq-0#95FN9sv}3UIDRznqty>jc^H&Nlp;OH)pLs!et;KCq zdn3S<=wNQd(z9;jFZkM;cUzum^{PD-hR z%7JsWO~I|#oA=`gGvZ^em7Jl>HB>YkyQDXri<(iikx)S-tv_EhtvT5DneCy|3n1LM zzVGSY1R0m4TRWn9FVD4`xBUWBuJOF5Q#F;_K{`3Io!0V0{NU`CXn3MMeQGMYpyzp? zghs+x73IVlJrhi#&`-AZ!fI zpcH}XPURspNz}8lI}ytdFYVxVPodnR)J8?3(Vo9&zlz!TO!!(n)#vLaUZG4AkIbV; zw^dCQciET>X3F;bc~ZMFwLN)D^o9K)_IoQOS^9i;#M65%-6!0tzN9r%oP&c@TH5y1L=p6JD9TQ^o-BYEiQ>?mKTq4@9vif%^CyL-+Ys__X}5- z`Let!`9VLTenD;ldJp2#%%VElI^~r|8W_QRZkBKyQ1KPJ3J9n3iJ^*>d;hl zD}#(p-xQr;Stn_bjqSm8Y~B?15C=ecMIPm;yqLZY*^(4Kax$hEtzQ_cEBhrQnu@LI zSsiOi${jtJw-N+t0hD-^8(*{_*8cutHYIJ(Yar-wKH0V#A8xFO7 zdMrUs>9R>ZzTn+cfn{hes~dbnv}{n`j?n_+FoHqw$?48W-&OxOJUuB?DaD}0mS+-X zlCsDsVurP6=2^LK86C^ouSXf+zG&z|JsvuyIUs3`6=L@oGTR6lZUfM=XezP&5~*qF z?Twpt1jk7&mrxL9-&#uOynphr&#VUp11K&#Y;z&q4p-JJa!>7Os$D3|H0f*$T5Xma z!J~O?c*5>G$cUL^zPcrtdkOPeD>uS{7@pn1Vb*nKYEwq^)7J<=$+84GO*6#ie3Ku7CC_&%cyGOGd+UyW=b8JG1es8q<~~K2}zO=kt>KkViS~tk4* z{rKo-R8yDmf-jVU2&n1P5>MbH%yl|QHme%|v`nLqICXsC@lQf z#aPdUfA94(2u|*R)#bSJGw>U`5J$o8loo+Bf-BBpJ>yD$I*LI7AYQzIhem|N3z3Z(gefAAUq1QsoWRTL_8rBDMts|H zH-|anDmi^=1TuhHdd1OaUgI$qx!ws)5vd^n5*yh~MNtzBl=(=21i;!FGRMadOhy9E zw)9y-98GF|u#~(b8NX$(J|}PD(G-4JAs(bvB!o$jnV^SXRMafA-5-w<__k=-Z0$7R z-bV;Ib87EuZsgu}Njf;Aj=;*9w=lzy4um0qh|=MA`#UpO4%h6_rBSj3++)Or796c0 zW?Da6S9U$f*T2Aoa0z0ZpIf#sI-+A8%{P&Z9&YMS2XaZCF%1kD((s`Feh9t!KZftB z`{Do9bs^*Tch?1Z%`zr=)g_ zlB%|upY*HOKkcDw$_*epdA}HbtGZ7I)kCGx_uBUF7O1I)zj7N+6_S3@;!WS*jLrM zPmVMIi)-R!Q@#em7j@0fu$ZU?&J|fV8*+aiQ#7i<8!GBlOZ{akUCsv|xw)@Rk+mS2 zw^JyW1O%VYB`!-g$i<7VbUlp~!IhW7UbINxHoQDeBjWV?v3VJVJ@~cR#9_spku`))GnRm74n0A#+fC$#OSyw+pnQ zpo!C*FZX%T_1wwa0|JbrNX8ovJd2)ali<}(d0c&xz*YZyIGm9^!QzPK{%?e_fjznC zpD$#?)2Gvw20tbbH;n8u^D5T6kAvzhaG!5zaWuKCW-!Tsds4Muh(Mvzl1gQXz$${8hdWi(CXg1%Ogx z$KiM9g2+7xr$MI``EO|E*8fNJkzk|V{L+-ESoCe~|hmzMTR}if>gTLdG5tw>88_-aEncRJC&9c>@)h#}r}fEuLE4 z?H{hj?qz^cjB6c|dbWEQGDd1vyO98Vkq@pU6*2(L{9kxI)f)!_4XPTH3sa> z5su{_Fzb)I@z}M`H!By*r;QGjme4tY zJMVk=!CO@aDyD6$(*)=5@m)Z-4t* z{&@8|y`S?_HnOG%1&DV*Wyv!HaVP2Xs(9lvo{i&cP^`ngL`~>bK>S#qp69(!Xx2nJ zAmhmgv4rFG?>0f9tzm}y{5&%2aIN4>$g>pPmtKnC1+)O_o^j>cpOsmuer6k{oepjP z#FM%z7^3Z`^~W%}6D5LDDp9;I_fXX_~XO-T9EeEiUxOc39cgc7CY2m8|} z6krnLpgEGcD!b4s0FV{nlmuhd_dBBjm>nJ0w-&uNb;Z}@$IeUx;k#?Dwhgz)04ztk z(D%_lq`drrPw8_D0t$ik=4j!btzcA}CvzrlOH{u*93UsXz*R-WG^UD;Ixz8+ssXEC?T{G92^iB_UbmtlU)tmjf3mO8H|TdOD3 z-?8k@safy`$}8POgt9T&cLZE{Kpi_zzJ?PLqth&Z8PXwot(+BZK?>nN82R|-wx`Ac zJe5{Ur=eRO-iUR#Ro`eV;pX2YSLpqL)g8-Bj&hwl+Fb*mzHwMmSC(Xc!zM=`EWBPd zVNmoSb%_O>?Z%gL{8>;K6t^8L0_-d&Q7>S;@tG2n!}A}uV=N?@jOE(l$!xwTgm5Go z{As{|voJ?Rfj_(6g%0hgdj0-BRr^nMpB?!atd&$FeTT?FqWd>*lq;*P9{iuJF%9hG z3%`*KoKfskvN~FHz6E-*G4jUHDs)fndq@c%G=G0p-k>NE?O}=0GHtfp+%M65HTR2+ zN~T}(`uU&F*Q&dj;*&-dmSbAfBI&s zy8I>@2G#|xlgxz$K*x-&59@eBqQ$^EfDt{pK6SANd_SZ;w35Fd4uDYLWN^(y=n;#z zIAO-^JHh)kCZ_ufYevuqUCHZd;)NZkAZ8_}2phTOvX$re24V$>-E4AfylrvLPwH)N z&kl1LEbG&@mQ8D+m5q_H)_-i^qH`lrBG1yu1jiR=2&I&^9R!Z+xuigJ*&s| zQ-()~T~DNxzd9BUxu}?-OI2N^nGMh@P|7k9!UVsvr~ze>JdSHT)j7w}or{$VnJuof zqbI7mhrHzy4uV5@Jwf!?BR(_i-p7B+8GKe#E6-j@&iU2T1WNb%*cV4$^xGu+2K^gp zE)Q%+18xcLRjv7*tQcEV)A!Iy`kjx@JId^v;^qr^n{1(sWbe{AHF>E*!86O(9q8;` zyLsz+1&ol3VAJ0&F3-d{(6`MR?67aO7GM?eBq)lmIT9xFrlc9@8sLrp;qRv8awk{iJuAJ>hq|e?gb$&cUPkB`U5~WHV!29OvntmyuVLlgX9C!=~L% z#mRpoS&hP0f5~fo+*V}adp|})>-&S5?GlfdX4O`We-?9+WcI@b4B`CD-Nf26OSPdO zM!uHXxVl&7BtpQWcLtD3=@)cnlGM^4iD6Jfn(lVsS|~UHLld*aR&ib*l885++Ir&%_^wTOCcK2?^X29tZO4 zC#SuKp1t#j{WKa3yu7lX?3DFnII%djR3p2HJ5AJ(1kz7A4gph+gk)(>l;bAE-e6MXKp3LA%~)- zW4%^XKaX%|x9>0i1PSlfURAXb0>;}ddx9wLwO@OyKr$f!T?cZB zGYF`w5~qq6q-vMtv0TnFnyY>tR_r=^oq67o>d}_S`$uoP*co2jr0Qo=u`KuJ5-uPE zl%l|^w~2ScH*QDOAIeLqW_t{U5u$e@{7dM1kffaLe(XKAe7VNfN}X$;I$2tu@e(IL z)_xGadfS}9O!6HuF?Z;K5Zvq!!oJl>CTBGV6aSX0bA06nIrMDb zoSUCf@AWrCvl}F93Ec=H?%72Pf68V;T;X^LBVGXAV_~;O!`x<9h?gVoSuBEV=H|sw@yK1dAy>aA!TeHMg$Hs;FBi zl<|P;cJ$Hv*%b-@Ap6m9KL4O_dd`QsQTdg;O^e>PU8(t>+4VZ;hrHk5pp1T;JEntl zAvw8YBs$sr&t#4g{UU+z9iXCM{LOtYpa_yI&_4KP_{&hzIWp!*APLuL?$US=o6ewT z%g4Pi_E7^;V=AWmPPOzOlUI{B*hR~5?*-hrN`juKV{&h8<3Minh3hDJyuL6Q@MtjN zze(x~+ahql>24~%9^YY&g zJ0=LblA$Xw@5!6Es2BB1;uf^DLD)YI+&3z!tW#O9qX)W6KjFPjp1N^wCx5VWlI4+h zs@r}pIi&6dpynY<%f@1e9m2De7}1nk=Q3`9=t=7KZiey#bQQM}Udsk>G&)xB!EC}G zk1F3m6IBBA!N>1pY#TB!Lriw>#WeZDId?gm>|3*a)voVLp~dUjK6u?-e09;g7YpID zi5o{|bsXjet^gXkr2uz{8RR8r12-*V%-dN@Z}a3FvYaIj{=4ocUaFG>jFz21*pn`A zc7`a%-D&b?(q{hDbxxxjakS*5Pxp6kC}rEl;>ExU*EatA(E>->qoVc3d=VtgDFlKO zu=NNkUjt4$OvzG?q#(w<>@zzKs?zhRSe(w?L1Uj;2BzQ$YJxIAi#5z-diWzc&*%Uk z5w-Jkcy(ew-C(@jMzhqM5Xhtnh>87ueNjD6x7!}kme@^8uOP!R_p1l!Hc?PnmWKKL? zDjh$7ZxuWEKl|eT_=61&VcFt4XPNu~9myB#Y(6iA`qrq_GBQSQtzC|QKG2OVvNE2N zSS#TAyT-tP5-GkuZ~z}e-8 zEwKfOH}W%>RrnvhS0N`L>*%G@c#GtN`?8c|B^&7jeE^0qxI*r`*MI;1&3;_J)r5QJ z-o2(r%HkM^_KuE50Ml>fNS<8+24``^%)b#8A`s=)qrP4z`?JXZ9Jj0A7$E}^5MXhA zay|D_TSJt^W()8GK=adWhWMp#89+;MjEQWr6EoVs{a3BHhc@uZ(RGZ4cg~@lhgL~& zR9BO{Jn!k?|1p{+e>+UbzDK7ICAKh3Lvfd$H*s5B(w^5LGrx8mAM)s&@pm$(*3Foa zj7!poE1nMpAm}4-~fwsR7T~1n}_s+*ibV`H3<>8;jxwbuU=m@Gl`5M^5M**80>JU zyPqPF)LEz-q=Q9|BmkF^La4Ek{Q9!Sxjt1(1uIhhq`b70UPSw}k>od zOM`tZlz|EK?eJ>Gu*1d;!m9FxN(tE{l0?xbq&M#J6deQlemCqUA=lg#emt~!Gx25X zd{>8R$t>Y2ZXe*=Gd|w08HQ2|TAf77SyYlPUZ=PWhxE#-(>ZO11HhUy8Dn`98>nbG zoAc44^77b$)o*U4KkINi(kC2=y=P`-7C7nb@*ceONB6T()ucAvPI<>z?zrl2lN;ag z)w>p%2?+1Xv!DShRM;4l5J^vX*YAHWK47|L@Bf@>Qr6ZTYk_xAu5OB6!J=&Ud9`vX zQSaN2kB!>zEv^TZ7+-_di8kvhQ$>xSof(r||NZBpGJkV|v+8^Y}#Mzg+oAl5+;Xd;&@Ms5tX|zDGoX{XDCj3`yf>{R0eo4@$VR z&?IzwH~A&*LC&sP*5=2RMo4>~?;F-i2eztjTT|G%R+0j{Bxoui^3|>Os*na}f7I;X zVeGwcyg-h|X2;{(qT}^9IM?3{C@T;VcDEEykhU+r7k!D@=+JU0!b1{@`qC)c*V+!C zeL^>yXYbRhp^s9Ilw1C@lOQ-#0jvZx0L2dIcbQqs2IorNm5M+Jz-wu7(3W8P+c~iT z(cA)rroa1cM9M-1Y|0S8+Jn>n+f~jFfZ$cJ)voV5J7tkFWEQlCOFt|4`$LZN2!;6&c+A6!okv z%Eo26_?IX57`4jkF>nxxxHm69dDD+4Y$iEz5BMKmcXQ6%`_z;%EP75?{E-8(@nlXL zFRXa%tqQocqx_cQq}TR?!OzZ4r739@YjHJR;XeH?LRH&AF(H7ku(4?YcGYMgweRGt zz6@;ZU(>`r>|01QK;bMo8Ki_H`z3GPBBArzCEWit8OE#_H`A{`*8*K1;?cfj2pI1O zpyZzeDWw)91H4fiUd>LAsAiR8_Cg^%Kr|FfrA^AcPu2`iNLfqark0C_EJFDY&s=^g9kPuZ@CHq;rO^*Jy4wy=GQ+LRI zcK81qH9Sbv5Y~a3$BoG9*+44L(Gl!2GnbX+1h=yxGtb5I39)GZQ2$2I*pBp_h>BjW zHNd@>CQeuNSO-nm*m=2C0P9?{Jl@-?$x<>TOoig^B;JXjjrk?noYC_PNEyHWQIN3; zCn)2Ma(o~(BeDTXsrCnJj^FAZdBoE}S9~UT2n9j)n80(V20ehNXDug7bvHta#)~O; zx-@|G?Oh_rF#SPspEaqmK|7iAvY*Xv8lCj#yc5EOLODvaKHUd(TFl2tG=zZ2(z2n* zQfHwIssaHqw`G=~FfP(v;>wP~)r=cwv!o*bE%%Wol&IY0Sj(8>TkpzDs+Axa5l=Aq zQXGjTB!n(cz6chS-x$es>3Y2`(DSzHEfO3aa3ypsl#OP`iHAh`^JS^m^Fk-Hp@-86 z!R=wt(9FQ5>p5%xo*jp2?8s!FZyGYWk;o1e00BdKi05p>c+!h9EOc-Oh?N6O6OlI; z)wFKGd2@Pw8M;IxSVaHmD&6={lBWepBs0%EQceAN;tr!e3;ztvY2famb9AhdyF6>p8Wt4Dlu$|6=DuS2a+@~?< z8Wb;Z3IC_LrWy_c1^WhUcPyn6bqEjce4g)Zc?SCy^;8~uW*J8hauVl3l_@$t44bCB zUFk6ZPo`d0k+AZopQ+@wMkRJ7T}hK+p(?~R900zba?S{Vpexu#%|1&6NDg$f0#eCE zSW-9;kU@zX7#6G8QVJDb-khdQVu-t;8#|S6e>_4u#*DgAQ2M*hf-L>LO5Af3@%V0T zO`?wu{o5~IljMru(>CM+Q<7Pp8G47R)-0!3qS=t5dw{w(`D{zp;!4O%NZyG^)#J|qA_;ISXdXMp0 z0DlF@?olB_(_MFpo`i01*}yA|68FBCi19xML`=%F?d3aXfOojkV2bSbt60sjh5r~EN&zlr1kgc;NkhBXkjFThT zsUCl#M-`eu=x`+{d(t=aL&8{&?xYk04SkE!|!J+BuyWrd$j{*dsrcI{rdy3M3YC?ufU$Vx-tC00;h2x zMF23C_Z{9j8;QV!am>TvBxSt-Ts`e1e?b!>Ub8*o)*dAYMWX7F-~{34x0&vhJ`XXo zsyh8vb|I+10_-LX*QIn*wjg`>Y>=3B9fA;>7?@(#{EJh0%Xrix6?l3)fe11@u4T|D5j zAS^>h-=&s6!?2*m;1HFbHKK$r?ZHN_LZDV)2}Dbmdy8Lntt3ns9;th(NoF=x1y0i( zS+%mW4W;K^@_VqHa7z`o{9M*a;{JFGU||^kb=n(5*+l?x~+y#7y98aC{owdBuez^zN%h{o3sDNnATS&~whclSXBLn0vZ_a~^Xw+Kav0BB06IZ#%y2VcstBv{#qn z=I4PRF)hd(Fs`ZGH9_AUyv9X-5!UP8WUeF2kk*E_R?(0Ua#`14S2|@?xrOjY9?A+i z@!*~5{`opytY5X}OgFgYP99s2lW4_j(8qU7Cn$t4h^*vr*1|bO%X4APysKU2PCX)c zL^koA>)-(@R%rSis1Pi8YauI#?=WldCSIL`e# zySaJzg{eTQAnNv_Wc^%vPZOdUt!fQFD@(H)B?s@j6p{Jp16dY#@xCeBGhlSmv)((As>lF=21y{97;Jas? zCf*BmaU=+RPt0czlE}q9H3Mq0AJGZI{VXodEf12m6AfD-451 z|3;%v!9|KPpY9eGzHw0uoapzH^2&erm4&sl9dnMT#pa~`^ZN6doAbSG#$*MUP0Q+E zqjV=iQ5$;oC(1W>j#5j$ppx_E-fBlQdtI1r+#vj5x@mW;s8qo{p1&#CN~Y`P8HZV7 zD9of;tV&+@{R~s4zt{Fj1vTEMv&e{Pd3z6INcWnps59!yYTAi#ZyKrOgFC1)UzfZk zUXxzzEQ;#zw!AFD1UYD^C#KX=0jHdPg^f``ULeed){Lm*W62YN5^2gD4HLc{)rCw zAA)+6@lxT1!vf^tk^1gKho>taRudjuy?jS}H|YL@CU|lkV5L8>QxceD4?emdHQ_0h z@p|8$L4p(qH=e;fu{l z9AIYuU1e2=d-}7Yp}WGdH%QQE4mByPj4aH~dgWO|hSwlnn-=c@AX0;aO}nP2i)}^& zYK)M>?nFe5%C5s;D_WzUqSZa2`%pZ_YW70qwip*yYgOfnF#w|>;|#P?r3y>d zbJQpWQwsHZO*GK+>Se{81Y(c1XA&hT;!mGcguKiZy;vc#F7a;+Prp1v+iFv2>3y>E z#At~B7Y$>BX0YR*_@VS<%hf3mF=F;f5qFYAx)=}vCRDz7Ujiy6@9s>*#aeac#l_f7 zJhSB{$5p45Fi_eNE4Wa*ozJ7HLZY+bwLt8A3lwqc>KgK&BSzikyd8z3D%uZkrPaBK z`Hl7a!R2+oElXV31R%+t+d|nkvT6*fN)a>Ip)p+^wVQT`Ho9;_7#TbMyOIcHufngF ze~Clg{RP@Zi#z|JUDPgqvESj(_nsx5OT_CuzuOzo ztMLYrW22T(670*DW!2P(3u^bZ-^9k=_n2_%6%T^uc5{A8(5yiUWsr&2L;j)Nm^Z92 zI+{=cIO>tF!VHi)Te5@(^vdNlkIx{2y7Ej!>~!nv$mgeM_wL=B+1#|-ORyO)e+qOs zt{_W;{0gK^Azkcc`k6q0+=K0z5ADUzV8&t%*?Ce=c?YTmBTn;7pFXoncOy!HMnmyri zDXt1cfDuxfnhD(!pyy}EayzSk=5TwSoQf*hoQj-W%F&Vc>eZ|Edp|Vhzkhrr&};E` zVsdr)2nRWl!D8tO==8qyGavxc3^L5o>FIVQ(@4^lcztgQP~Y>X$3{vv>bWhFghxbV zfaWP8Dm68NM$~2IE52kiM$kScq@oJ@14AljU;w=LjCr(7@ts95Psx1rK+h3{<}00WV8zgejg5Z3n^K|x zZ11wNFQ~X0ucD(M?|Ta}ZS$R6^|F?hKf@#Uwzs8aWKcmj8!5sA(vf52cSa5+K;tkt zBBE>cgCH_C=tdQOt3PccCnQbt}L^Z599Yj;pN3)EbC#oKgNkT)s2vj|#; zpviPg>)62z5<%P#A3kI`1{KN~ z1<7MV5TGRhEtYJ>c*8&EXO&PepgDVue2QSPG*E2uYU0zsFN_f-m<-Y$z~*IOWyPY5zyxk1q^dK=e4v5%OaSI8V@BDKc$re&%HS>R= zSN5=m9l3fKBUy;hrS13W*E2vk8T7(OD(vWn3ia*EdO<;|$h4dC6iPoS<{mRIv5*93 zU}OmpNC`9~B@L7X(pGlz?YwX^GBR9GHY#kqeI`HPWqb67j5%`l1P;W9v|c<#Yt_c( z;Qu&ING;2>zk~3*{pq0t*nSl5o)+c6bn!mKfRw`Fk z%u7p4HO0+5oJH@|vCADb;+)hX4K|eg7YWz=ynn}}&vYg75p{mBKKfDU9S|2snx|cc z<4aEuGGxG?C3?OUBnh(5O1m)~#+ma-wJ<6&C@}@U!CsMfzASU`nD1xz`REU(Jh+ew z>Yi{e-nnxpp}0^!tMO)Gjc-N<1IWx~*45G8xN!sNG|0-LBj{OJS`EuVce>k9-8e;; zGf3v+3z?6fkHs#mUm}!q7_m_=?haJ7+Th#_`TqSo1>M^GydQ`j!;ywmkJr%my0zdq zuLJaEFM2@&?E9a;6~Hz`yf2>7ZlPT9Gsxu9*6a=G9Krw+kdq7XSbHgu)DuM2U5ZE*4FFI!@z6A*Pl~WMX|THM-;i|E`ENxENQeo z*Mb!00_hCjuR-`g=x&mEvEH1PAJhQ-8nS-xvOe+>JRPi?8U=5DFJ26X+qKA#Xmi#V z;UrQsqayV}qnnZO;QyUMEp?HLAxCnNt=h7GG*Xe0@b2x|EFc(6{=suuYF$?HeT8|h zknZZgJ~ZNV-u&A){cm=Qi70O)?Y$Ns`DJn%N|Lf}lSg>Muk>;hV2v3Ubi-`g@7wEQ zPLquvPK~r4w2IVmm2cM;MQldcc-rMTWgi42_od$9jK59u7TR9Sad@k?C4&YufTkv| zz<)CgzGOY`VMO!K-wV_U4pp6;G%;&l>&`7;P!@?ZE ze8Kc#@FreX7To)dUe0YESL-~pb*@d}edJ#lsNMsEbjz=NX~4rOa})CP&z-X|8r7r5jhGr8N{(M!5Jw3pjzJqc7qGX4cFP(S^Z@A`*-Yv$>-@?V^)Qi6yw+l zf5D~>lYmJe^WHB>zhqZ`R!%@f#DpCAS6+m8yyFK;;G&Yz`kXqIkRY*Eljh ze7N)B<43|N*ZC`ry}in=+bv|Q$4E=rbh25`1#*nov9Ym8KnXLLHl%#t!&>(0gbfeA zf}9+MJBo@#0@fq+Fr*`v>0@Ks3=$H0u$)pqJIu?U;7>K0o=hRHg^I|?%Tod}vb4j2 z)n`OrrP}$QJ&zn`Kjju&#iJFw4ptZk3y+!r2NxG^(1j;HJ$ch}S+oh<gjpVf1<9Yb`9=m7M$SN8IsX9>*;t{OZ`^1=H|fV+vKnR=XjXN>-=|o zJvLkHJ=dLX!GNSB%GnrY2|`07qW~;gSE)=56S@*FJ3MM3de}fnMRyrz=Q8`51{4M^ zF6GMoF)P;NQ?QL>b8~b4KBx0^GLIj}Zoy8(g`;v;wCjf#Fse1d3S|e>k|Kg0kRl%6 z&VT&);hUbGZvW@(02eviid?nV4uJ?0!j>u7bg-?sH~}(;lbhbGZ{p%`>Q5#&(X$`k zpcZ%k-dJdAO8z2B=L;++B%7onbm66PoX(NeC+v@eaeisHMDS`kIu?WfvEydfy12B2 zQ-9QLXDv}Uxa!}~AT?Vyte3nalx-S$u=r)IfXI8xk4|aSY;_>J34|0og(Ld4mX4G5!VIpqroai4Z4-D-lZ&pzQYZevUFqHDjF0$&P_8=EvIfEUSVYbSAw zh*;?z|J31B?Zfa^O-#IPltiOY^>pCdqbXny+OxiT-2GkZtq7n;C>Vwir zScpy8@c^4XfI|+)Y8be~UNG~+B98Pq&5DYPW3lu&Ud;*c_xBHJD=I2#0&nW(&k|TZ zv!K;kmnzSZpNLK``RAxN0L+i-ICA^pLXec?>IpKhm!^-0DrCJ1YJT6Fg|%NkuDR@| z__=)}{MD)AJos^aq1-M?4n>>mp1)JRv|Al_~#C9=oqsz!!A*7%Pfk-yU zg{tt)(vDNfv(<;qt*vv#ot$661z>jn1*=0AR6vm@Jw(1tl9F`D8Gk#xwsB5LZNj

NV*TYV8`>#iZjOelR?g11ke0JTJ4qn^z;y5D4|8klW?=b4mnMJ~)lhJOly|?B zX>&kYVnqmSxgQq(*dMvdU?zoEnl+ISJ~JcnZ!iHjPAb1SzaAJxqfg3z^$Krg&mHbN zY%D7)s~^-@-lPAe@=Ih?#B9Gqi|nN4#RspOn(MIGAeBOe<+2Y ztq>p|zNtyZ17C~R$;}Qn?KX16HOs9jO0_wV_Da`f)C33(H)Rd9pFD}+sz1Gobnv$U zbdr&oi3(0V{xcC8u<+-dvxS|sOTVYoG&JUC8vRsqe~i?62#Sb^*n{QG%FR_Q#W$YV zo^MlCQ}f>#uMl0hCDi08>6e_GybX$}^Z6B1SIee7ce>5(?Ajn}!hnqm7xo5&Fi%#0 z81STR^A6k~65j$;g+!3aX=-oBhHupEq;4XKPZL+17GpK$LRixpj`!B0)6!l+ssfeH zbL*uJ#%{?E=mXsI`o8*XWo1Rm+B&aWqQEh=8R0(>r^HflegoTAF1aFU&0`0WCiBq> zw2go4HtVu7MaCpw_>hqj5E9(@!hVN~R83ZnUD8VLmV5-^=aNCm`?G@xSp-idwh%OA~KFKI8BSDnkg&KV_^?j)x0tE=|8sSD@C@~^H5 zKC_;N_Vzo5X&@TQ1VqCO3=D?Ksgz&FN;|ug1yzCmGdLLA)x+cK*chcz+bfdyR#sL_ zY-~7?7pZjLV_WJ;#Z^>P+}#-0AbhNnkW`MGU%*rbnRLcA4h^Z#6@mXTGdufI++)81 z3>U-VV-~IBh}2Y?_4Rds$?D2Vdhm55BqW*;?~04B%FD|m=S%JJfwZ^`5sS%P-Gt0|ZW@cvBu3aRnmOFwLrpkfEeOLqq%U=~JlO_UP!SVgGsk{U>8f?a@rg)7pN0gy-Y5HJlFyuxBeB z7jZxX{`;#wf*ANg4zAP)dg?RfY=K<1v0oRQ3K7!F| zoS2w+2CObPcg=ZPB?z@iP2SdqMMtFQE)aUp9PKQV3s@)E*MIwV4}pm9l+)Fv&dbXq z*1UNXHN2Swe~H`{P3dl8v=Ga@1rs+h^!`2-Qu!K0pj9+93O!~?@&s(Y&Tls??*uer zsXVr8B_9|VFfcM=n$UwF7bX!Xa&>t~RHDA@QnZpIojhe_(2Id*mWxj*x0W4$a>j~nZnt#GFrk41B;JX>B~TasIag4xgALTUKOLo4-5Pg?RVFP@nKfD zZ{2#VzV}mB?~Db+$;D`u51qhMg6t&y`gjg5^f1KGGRpAtzqjT5Xl zoWF@OI5KK}>z$rXZ)tA+Oy_<0(!#>xlE#(o{e9E4v^2&OS$X+l*!D{c3(c@w+ZdUc zv_SWpySB7cEg`6F9uUq;MezCl?&_oLq5J}4n`Cfj_&DP`J-xi-;(WV0JBgOGAeeZh zp<#w}U^^bEsJz`h{t=hAQS5c%vcH_}Qw&D{DIHxeERfD;&z_Am1>gp_f$8k__4Ct$ zfe4F;FrgB(HP+>f-xw`XSl!DF8`D+(aaRW;_e<*+oVYxj62Muz{QQ0cc@eYXVDzgEzlMoVe$;sK-1SXfA zpa0G1$jIj#kKPBlaRbSWU;60N58JxZ9(nqoU?!S(mjTjv-i1ATTcY6msI&$UBQmNNBg6 zsBAnw*z9FzXP5YL^~r4*8v2?)hz0?hv1SvIO)xt<2M1(2fr9Mp z;U7PJDlW4cE*=@tM4tEBf$;=xAhBm+YHF(Vi+QI>t-=J+8my!Y#Zb!MN_a!h?jWy_+k=rjz#t&Ub1f#gOdGF!4uJ;Qj16L zv$DQY)6gWQr_<$r_&~6vRcX&qQeJL)KZ*aTySsZ=OG^L-CMH|SJx8XsJ+(9V z@^a*a!-!nn57Dc2FWOn|4Wr=`6H5^m73I5f<%+3=MF&tOGyg7#YLyxCCZeM1#Kfa+ z2fgw{n22Q2=bSe=IY~Q+V7d^WAzf~wCe2`M6AzKx`;Q+xTsEgTwp}4dXN`1X>r5;y zznOL=^nb7|#(>2KSmR_?uQrX9T72c_=TE4sli1naHC|m^ z{r2gVE_U2z-2QXTht=R~aw{k(ye}x|>FVwtg}|cul?NDNQ!}$R*e?mSwYBYjQOB)M z^P)YUD_}GD6Ubm&?3&~K^^QT-;yU8u{}}wQU%{Y9ehUB0&OrSAE|Z`5#b5s?^Cl4W z?{~oiDu|2k$zZEe;z<2{se89kuo)Kg+diILCQqznV7^sd7nVtIg0Ck(-MiHI-I6+7 zX=3KPSD@~zliFby+tcV2XzTCgF}BwbDlPf2Id^^X zO{z3$p-At9CPF}@gY-^{|iL7H-H@&fSvb!#C255ba%lFyn&UhQneU=U_njJ_-M zL&e^w;$m6sG@K05RPRVbpswX0{5Zr5^FY zYIbzgJt1q4dS)iw7P?>}I2Y^;0#(+D8+b9^X0N_J{C!L-}qn4-3i=BuL(mm z)YXYVWSCu}VQ$?FAos%Mpb~Cw`%9~~K`4b$OJ|!0S|SfVMZqRHs;+*&bCF?A@bxnm zWo;K^sK9VYd&y@hAeHOKJjJ-&Y@Gw#iwq(iE!y|^;8JWs%5tO`3T;~3m*lnP!RKGE zcaZ(+KDYA%0SE}K?MALhcSUi6J)P1ULRYU9kyAz%T(3**?vkA$Z*{e~6vr*!0k)j* z+7CT=QbpX3{o2$u`g?IKHdIIV?^Kx6FOH2AY)>otA{8J(ne%%(u>stJC5PnP2 zL)f7V&mdVX&7;Fih)81BU0b6imT_a@)mcyi_C8Yc*+_K(EAxzT$6PICDA{%-;6d+Q zfj!zW}$h{8&AEk+iQ-F_FQg=>S{NPq6=f+^WbGWEARt>2bNCrB%Gm*i_OIr zIa@^1xWI$-bcSW~V7A!IhJd-tjXP}FU~+QF!+&f3DpyDTjmDa3vq$6T8x9AMn4Cp9 zb)Df5)w>=wFeo_A^`@BpTakpKGq$1WVuh^_=t4%e3p%vsBsXzfH~}bFtvehe@O?5& zT*AT<75dn$l{Mo0oFnyXZFm0dAS^1s^X0I+;y(X}f3jiW?N)pHj{0OEq@$Hb`B2pfglk%@uM;uzx^e34;Y@Z^ZQ4+|1}6;lTVIdg#Osm^GY^9_2f3Ai_wN_VK+t(5Cq?}I3$JPv z2;5s_(|h@nBu)e7DuEFLbxed(3a0=*$sZn}-8+47M~1b50e~Uyf9-B4c84Pfx!zr6 zBa{93_|$+O8=Nj@PRvTeXOG)MRvRCFme=$N(9KI>hpcdqjW<@-OQ>9awF121zB~t* zU7GdH?0OpVUi)*F@{W~Ssye+Rv4c=BO3?=)g)84X^&qfXdD+uDs+r!Am7W|GkbM;` zeWFrXD;Yo%7Sw=$D2VT2bDpQ^gyHzAe6L-Hd*BKLCG$aj?VLvdeHPqZad5)Gz}+_~ z74&WRzm2MnoT=@5+aBy$yf)(x8WMRF+o_a8FX*=z8_8nzykmtHLv6aK$7ic2VhYYwn5DyM(6`=%2rHWScx{yMn@oUjVfNu!tPQ z<#F~WmP{Zp5^6;Z3=Iqjj$H+)T+eP}9JJrqTZP+bh*ZLcz2ze`d#bvkRwrd@-?MK@ zC6~*CAP`3&8lEayfrZVRg3ferH@La7JkPw{wfyePa&R=v7{*d#mTo0x{p#ipcu|CtCyN;gyiYb~(iE9%hpUJrY z+PEVDCOI@)QY$xGxy2_SAYc*(RD1k8PEJmTe%gGtjlpF#y%O^A=q!MntB`c`-`^N$GA`@-1_<-!VS7z%ns3># zIM0OyRZT?IiS2qlTE{?Fx09iC+af4T=uIfxP|68uZ*2|9Vo~%~Zm=-%p+NtE7Iyp4 z(|xH#?2cYG|IBXn*5lIpNS%W8IX<=mM&M@KqG}P+>{o<`Ece-0^Q?!D?jrImE*4S{ z%Xb@zaGu=&mM&3IQH0$Da3A!zxy~Bz;Gc-KA#neZjQ44*+2UFyiiyR`^KW3eEd1vAYYXbK0<*6}gp_`34>< zB_&8_v;zavCP_UiLyD|d-uLahR@Kv}6kuK!RJwk7S;4n&fhwfEc@d#Z|J1k_kI;c8 zne~PC?kwXfug|Q!y}9cgvu+6tv>SfgUVa+|*Z~3|%;S~4 zt-cYE|A66c7xXVEnrYfJHQU+Y{&}lAgzMR=5kbJ@p1dq+Hs!xi{egBmV@K4`1dUuB ziXg13a|p6QVyQBSX6pAp>afS!`{$GsC(KXmO(>|TCEh?$6?;`{maE}e z)P4Uvw{2s>;y%Y`>>rb|>W#W#vAWbir74D2T;3r81u3EfS(=A90PI(ss3gmKeU?-1 z3v)P-i&-2veqr`9GAYR+KmY9t0CIvB9#rT_tsWOf-fzdd8+=j-3w?cBM~o@wN{?nV z^N5}UYdiu$XHqYvbUyk%Ki^te=NA>kJK^Y2NmLWA2e&gAGWpmemZhBFael=Zy~d>D z*@t!HRa|2HKdb^o*7&eV!Ea9l~GxRPS#q zeO$NAvd3@dTyhqcy<-95tbBUf)-R2MAL)wPa{0~Sx0R=U+`Zq zf5Qq=>`pYwHn%9o95OqJPyMr=xY=IuwWi4swX2b`<&KD%F;nx9insTzXjETLON-Lm>_DUy$Ul~P{;y!eosv&0 zd+f?x$DL1e@!EQN{G}&j*0o>LlexFNc}TdER7Ps@;X{E<27XLe>)PWsMlW7)WZV}) z%^qtK4?NS3i5iNP+QirfIrC6_RBw|;(R_& z%^VSVP2XV_NjSMDA)NDu+fn($EvL%ilne)CIRqPLdI)M84ZcFE>_NN5o!1n@u+Q+MHa8Aj%#qwDb6ixH{ z?}G2_3Ky%#zT^ncHO*+i(=s@wX496#rt=nCvC>afx%t*IvXKmlc$A5&E9&OY8FMv{ zJ)P~SUpWQigeNB|;SG(n$#jt{QqR<(s)@gEvy-Hxn;(55fx(~=D=P&t#`Je zF~j<^1w!TX3xe{#={!fpm_+N`Ay?bV|6x zvunwooGf2cCDDU!QDTO$K1hMCuhD*TOC8y1*6(T++&)matO1}IcGaYr#wUIi(Q0Nr1u1+j zo!>ma7!+B&8!k$><*s)t0^r)zI-b}sKnMt8m(S9|Y9m8Tz7wD+kxJUIY~te3$&kN>u_qRSj{z;B#x^(TFxm|0ci1)7O(ju8PQQKy2cBmbv-#`UAYG2lK=FJ=lN0N0IL z9NMwu<3cg%dEfCufBv}Xn!O~+v2R@DPYfF?#&}3VXGEDHcOn6uW@*e9>`->%bx=sm zUo4vnInc~_bT@RYCWDEkRCsbE+!Eq)vh90b^M(?4!OF#qthJiv0(6~tpP|#us2F1j zA#LXQb`X@ethZmnQpQ(`dBULgx6=fT_esdL6GTQnD=5RLsVl1ja@14)wzP)~O^-~l zvq!wk*@d}4`p$f0WVv3mw7D>>yB(tkJU($*$BQ$|z-a-}-wiK!gt6<&6L+(nBWy4W4}iv!&^7k}njV9(e(Wj^eKf%yaU#iwRc z7YCLg*tszgE<|`9aWb8}S!)gJEL?{ZprFgkqdAD75P)}kMc*TytneAQh&YJdajyiJ zUhuqEI8j_yZcEFJLY`1JS)6ntI8W`k4)9z8=`J=-VPV>B;BeNr@;u{k43D}QZ-U<= zMykNR>%pWp220yV>gOUHV4y5PuIa=?r|0RZ!J`pVkO7FLy@n1Z< z2UAcA^sV%~;x;2|iopjz 1: + ind = [j for j in range(len(lambda_values[0])) if len(list(set([x[j] for x in lambda_values]))) > 1][0] + raise ValueError( + "Provided DataFrame, df_list[{}] has more than one lambda value in df.index[{}]".format(i, ind) + ) + logger.info("Begin forward analysis") forward_list = [] forward_error_list = [] @@ -170,14 +185,9 @@ def _forward_backward_convergence_estimate( my_estimator.initial_f_k = result.delta_f_.iloc[0, :] mean = result.delta_f_.iloc[0, -1] if estimator.lower() == "bar": - error = np.sqrt( - sum( - [ - result.d_delta_f_.iloc[i, i + 1] ** 2 - for i in range(len(result.d_delta_f_) - 1) - ] - ) - ) + # See https://github.com/alchemistry/alchemlyb/pull/60#issuecomment-430720742 + # Error estimate generated by BAR ARE correlated + error = np.nan else: error = result.d_delta_f_.iloc[0, -1] if estimator.lower() == "mbar" and error > error_tol: @@ -262,7 +272,7 @@ def fwdrev_cumavg_Rc(series, precision=0.01, tol=2): float Convergence time fraction :math:`R_c` [Fan2021]_ :class:`pandas.DataFrame` - The DataFrame with moving average. :: + The DataFrame with block average. :: Forward Backward data_fraction 0 3.016442 3.065176 0.1 @@ -389,3 +399,106 @@ def A_c(series_list, precision=0.01, tol=2): d_R_c = sorted_array[-i] - sorted_array[-i - 1] result += d_R_c * sum(R_c_list <= element) / n_R_c return result + + +def block_average(df_list, estimator="MBAR", num=10, **kwargs): + """Free energy estimate for portions of the trajectory. + + Generate the free energy estimate for a series of blocks in time, + with the specified number of equally spaced points. + For example, setting `num` to 10 would give the block averages + which is the free energy estimate from the first 10% alone, then the + next 10% ... of the data. + + Parameters + ---------- + df_list : list + List of DataFrame of either dHdl or u_nk, where each represents a + different value of lambda. + estimator : {'MBAR', 'BAR', 'TI'} + Name of the estimators. + num : int + The number of blocks used to divide *each* DataFrame. Note that + if the DataFrames are different lengths, the number of samples + contributed to each block will be different. + kwargs : dict + Keyword arguments to be passed to the estimator. + + Returns + ------- + :class:`pandas.DataFrame` + The DataFrame with estimate data. :: + + FE FE_Error + 0 3.016442 0.052748 + 1 3.078106 0.037170 + 2 3.072561 0.030186 + 3 3.048325 0.026070 + 4 3.049769 0.023359 + 5 3.034078 0.021260 + 6 3.043274 0.019642 + 7 3.035460 0.018340 + 8 3.042032 0.017319 + 9 3.044149 0.016405 + + + .. versionadded:: 2.4.0 + + """ + logger.info("Start block averaging analysis.") + logger.info("Check data availability.") + if estimator not in (FEP_ESTIMATORS + TI_ESTIMATORS): + msg = f"Estimator {estimator} is not available in {FEP_ESTIMATORS + TI_ESTIMATORS}." + logger.error(msg) + raise ValueError(msg) + else: + # select estimator class by name + estimator_fit = estimators_dispatch[estimator](**kwargs).fit + logger.info(f"Use {estimator} estimator for convergence analysis.") + + # Check that each df in the list has only one value of lambda + for i, df in enumerate(df_list): + lambda_values = list(set([x[1:] for x in df.index.to_numpy()])) + if len(lambda_values) > 1: + ind = [j for j in range(len(lambda_values[0])) if len(list(set([x[j] for x in lambda_values]))) > 1][0] + raise ValueError( + "Provided DataFrame, df_list[{}] has more than one lambda value in df.index[{}]".format(i, ind) + ) + + if estimator in ["BAR"] and len(df_list) > 2: + raise ValueError( + "Restrict to two DataFrames, one with a fep-lambda value and one its forward adjacent state for a " + "meaningful result." + ) + + logger.info("Begin Moving Average Analysis") + average_list = [] + average_error_list = [] + for i in range(1, num): + logger.info("Moving Average Analysis: {:.2f}%".format(100 * i / num)) + sample = [] + for data in df_list: + ind1, ind2 = len(data) // num * (i - 1), len(data) // num * i + sample.append(data[ind1:ind2]) + sample = concat(sample) + result = estimator_fit(sample) + + average_list.append(result.delta_f_.iloc[0, -1]) + if estimator.lower() == "bar": + # See https://github.com/alchemistry/alchemlyb/pull/60#issuecomment-430720742 + # Error estimate generated by BAR ARE correlated + average_error_list.append(np.nan) + else: + average_error_list.append(result.d_delta_f_.iloc[0, -1]) + logger.info( + "{:.2f} +/- {:.2f} kT".format(average_list[-1], average_error_list[-1]) + ) + + convergence = pd.DataFrame( + { + "FE": average_list, + "FE_Error": average_error_list, + } + ) + convergence.attrs = df_list[0].attrs + return convergence diff --git a/src/alchemlyb/estimators/bar_.py b/src/alchemlyb/estimators/bar_.py index bb0afd4d..52353064 100644 --- a/src/alchemlyb/estimators/bar_.py +++ b/src/alchemlyb/estimators/bar_.py @@ -54,6 +54,9 @@ class BAR(BaseEstimator, _EstimatorMixOut): .. versionchanged:: 1.0.0 `delta_f_`, `d_delta_f_`, `states_` are view of the original object. + .. versionchanged:: 2.4.0 + Added assessment of lambda states represented in the indices of u_nk + to provide meaningful errors to ensure proper use. """ @@ -88,7 +91,7 @@ def fit(self, u_nk): # sort by state so that rows from same state are in contiguous blocks u_nk = u_nk.sort_index(level=u_nk.index.names[1:]) - # get a list of the lambda states + # get a list of the lambda states that are sampled self._states_ = u_nk.columns.values.tolist() # group u_nk by lambda states @@ -97,11 +100,25 @@ def fit(self, u_nk): (len(groups.get_group(i)) if i in groups.groups else 0) for i in u_nk.columns ] - + + # Pull lambda states from indices + states = list(set( x[1:] for x in u_nk.index)) + for state in states: + if len(state) == 1: + state = state[0] + if state not in self._states_: + raise ValueError( + f"Indexed lambda state, {state}, is not represented in u_nk columns:" + f" {self._states_}" + ) + # Now get free energy differences and their uncertainties for each step deltas = np.array([]) d_deltas = np.array([]) for k in range(len(N_k) - 1): + if N_k[k] == 0 or N_k[k + 1] == 0: + continue + # get us from lambda step k uk = groups.get_group(self._states_[k]) # get w_F @@ -126,6 +143,13 @@ def fit(self, u_nk): deltas = np.append(deltas, df) d_deltas = np.append(d_deltas, ddf**2) + if len(deltas) == 0 and len(states) > 1: + raise ValueError( + "u_nk does not contain energies computed between any adjacent states.\n" + "To compute the free energy with BAR, ensure that values in u_nk exist" + f" for the columns:\n{states}." + ) + # build matrix of deltas between each state adelta = np.zeros((len(deltas) + 1, len(deltas) + 1)) ad_delta = np.zeros_like(adelta) @@ -150,13 +174,11 @@ def fit(self, u_nk): ad_delta += np.diagflat(np.array(dout), k=j + 1) # yield standard delta_f_ free energies between each state - self._delta_f_ = pd.DataFrame( - adelta - adelta.T, columns=self._states_, index=self._states_ - ) + self._delta_f_ = pd.DataFrame(adelta - adelta.T, columns=states, index=states) # yield standard deviation d_delta_f_ between each state self._d_delta_f_ = pd.DataFrame( - np.sqrt(ad_delta + ad_delta.T), columns=self._states_, index=self._states_ + np.sqrt(ad_delta + ad_delta.T), columns=states, index=states ) self._delta_f_.attrs = u_nk.attrs self._d_delta_f_.attrs = u_nk.attrs diff --git a/src/alchemlyb/estimators/mbar_.py b/src/alchemlyb/estimators/mbar_.py index e0ab594a..fa399cb1 100644 --- a/src/alchemlyb/estimators/mbar_.py +++ b/src/alchemlyb/estimators/mbar_.py @@ -33,7 +33,7 @@ class MBAR(BaseEstimator, _EstimatorMixOut): .. versionchanged:: 2.3.0 The new default is now "BAR" as it provides a substantial speedup over the previous default `None`. - + method : str, optional, default="robust" The optimization routine to use. This can be any of the methods @@ -79,6 +79,11 @@ class MBAR(BaseEstimator, _EstimatorMixOut): default value for `method` was changed from "hybr" to "robust" .. versionchanged:: 2.1.0 `n_bootstraps` option added. + .. versionchanged:: 2.4.0 + Handle initial estimate, initial_f_k, from bar in the instance + that not all lambda states represented as column headers are + represented in the indices of u_nk. + """ def __init__( @@ -135,6 +140,20 @@ def fit(self, u_nk): ) bar.fit(u_nk) initial_f_k = bar.delta_f_.iloc[0, :] + if len(bar.delta_f_.iloc[0, :]) != len(self._states_): + states = [ + x + for i, x in enumerate(self._states_[:-1]) + if N_k[i] > 0 and N_k[i + 1] > 0 + ] + initial_f_k = pd.Series( + [ + initial_f_k.loc(x) if x in states else np.nan + for x in self._states_ + ], + index=self._states_, + dtype=float, + ) else: initial_f_k = self.initial_f_k diff --git a/src/alchemlyb/tests/test_convergence.py b/src/alchemlyb/tests/test_convergence.py index 2bde94b3..7852217b 100644 --- a/src/alchemlyb/tests/test_convergence.py +++ b/src/alchemlyb/tests/test_convergence.py @@ -2,7 +2,13 @@ import pandas as pd import pytest -from alchemlyb.convergence import forward_backward_convergence, fwdrev_cumavg_Rc, A_c +from alchemlyb import concat +from alchemlyb.convergence import ( + forward_backward_convergence, + fwdrev_cumavg_Rc, + A_c, + block_average, +) from alchemlyb.convergence.convergence import _cummean @@ -26,6 +32,113 @@ def test_convergence_fep(gmx_benzene_Coulomb_u_nk, estimator): assert convergence.loc[9, "Backward"] == pytest.approx(3.04, 0.01) +def test_block_average_ti(gmx_benzene_Coulomb_dHdl): + df_avg = block_average(gmx_benzene_Coulomb_dHdl, "TI") + assert df_avg.shape == (9, 2) + assert df_avg.loc[1, "FE"] == pytest.approx(3.18, 0.01) + assert df_avg.loc[1, "FE_Error"] == pytest.approx(0.07, 0.1) + assert df_avg.loc[8, "FE"] == pytest.approx(3.15, 0.01) + assert df_avg.loc[8, "FE_Error"] == pytest.approx(0.07, 0.1) + + +@pytest.mark.parametrize("estimator", ["DUMMY"]) +def test_block_average_error_1(gmx_ABFE_complex_u_nk, estimator): + with pytest.raises(ValueError, match=r"Estimator DUMMY is not available .*"): + _ = block_average(gmx_ABFE_complex_u_nk, estimator) + + +@pytest.mark.parametrize("estimator", ["MBAR"]) +def test_block_average_error_2_mbar(gmx_ABFE_complex_u_nk, estimator): + df_list = gmx_ABFE_complex_u_nk[10:15] + with pytest.raises( + ValueError, + match=r"Provided DataFrame, df_list\[0\] has more than one lambda value in df.index\[0\]", + ): + _ = block_average([concat(df_list)], estimator) + + df_list = gmx_ABFE_complex_u_nk[14:17] + with pytest.raises( + ValueError, + match=r"Provided DataFrame, df_list\[0\] has more than one lambda value in df.index\[1\]", + ): + _ = block_average([concat(df_list)], estimator) + + +@pytest.mark.parametrize("estimator", ["BAR"]) +def test_block_average_error_2_bar(gmx_ABFE_complex_u_nk, estimator): + df_list = gmx_ABFE_complex_u_nk[10:13] + with pytest.raises( + ValueError, + match=r"Restrict to two DataFrames, one with a fep-lambda value .*", + ): + _ = block_average(df_list, estimator) + + df_list = gmx_ABFE_complex_u_nk[14:17] + with pytest.raises( + ValueError, + match=r"Restrict to two DataFrames, one with a fep-lambda value .*", + ): + _ = block_average(df_list, estimator) + + +@pytest.mark.parametrize("estimator", ["BAR"]) +def test_block_average_error_3_bar(gmx_ABFE_complex_u_nk, estimator): + # Test if lambda state column representing one of the two lambda + # states in the df indices is missing from *both* dfs. + df_list = gmx_ABFE_complex_u_nk[10:12] + state1 = list(set(x[1:] for x in df_list[0].index))[0] + df_list[0] = df_list[0].drop(state1, axis=1) + df_list[1] = df_list[1].drop(state1, axis=1) + with pytest.raises( + ValueError, + match=r"Indexed lambda state, .*", + ): + _ = block_average(df_list, estimator) + + +@pytest.mark.parametrize("estimator", ["BAR"]) +def test_block_average_error_4_bar(gmx_ABFE_complex_u_nk, estimator): + # Test if lambda state column representing one of the two lambda + # states in the df indices is missing from *one* dfs. + df_list = gmx_ABFE_complex_u_nk[10:12] + state1 = list(set(x[1:] for x in df_list[0].index))[0] + df_list[0] = df_list[0].drop(state1, axis=1) + with pytest.raises( + ValueError, + match=r"u_nk does not contain energies computed between any adjacent .*", + ): + _ = block_average(df_list, estimator) + + +@pytest.mark.parametrize("estimator", ["BAR"]) +def test_block_average_bar(gmx_ABFE_complex_u_nk, estimator): + df_avg = block_average(gmx_ABFE_complex_u_nk[10:12], estimator) + assert df_avg.shape == (9, 2) + assert df_avg.loc[0, "FE"] == pytest.approx(3.701, 0.01) + assert np.isnan(df_avg.loc[0, "FE_Error"]) + assert df_avg.loc[8, "FE"] == pytest.approx(3.603, 0.01) + assert np.isnan(df_avg.loc[8, "FE_Error"]) + + df_list = gmx_ABFE_complex_u_nk[14:16] + df_list[-1] = df_list[-1].iloc[:-2] + df_avg = block_average(df_list, estimator) + assert df_avg.shape == (9, 2) + assert df_avg.loc[0, "FE"] == pytest.approx(0.651, 0.01) + assert np.isnan(df_avg.loc[0, "FE_Error"]) + assert df_avg.loc[8, "FE"] == pytest.approx(0.901, 0.01) + assert np.isnan(df_avg.loc[8, "FE_Error"]) + + +@pytest.mark.parametrize("estimator", ["MBAR"]) +def test_block_average_mbar(gmx_benzene_Coulomb_u_nk, estimator): + df_avg = block_average([gmx_benzene_Coulomb_u_nk[0]], estimator) + assert df_avg.shape == (9, 2) + assert df_avg.loc[0, "FE"] == pytest.approx(3.41, 0.01) + assert df_avg.loc[0, "FE_Error"] == pytest.approx(0.22, 0.01) + assert df_avg.loc[8, "FE"] == pytest.approx(2.83, 0.01) + assert df_avg.loc[8, "FE_Error"] == pytest.approx(0.33, 0.01) + + def test_convergence_wrong_estimator(gmx_benzene_Coulomb_dHdl): with pytest.raises(ValueError, match="is not available in"): forward_backward_convergence(gmx_benzene_Coulomb_dHdl, "WWW") @@ -52,6 +165,23 @@ def test_convergence_method(gmx_benzene_Coulomb_u_nk): assert len(convergence) == 2 +@pytest.mark.parametrize("estimator", ["MBAR"]) +def test_forward_backward_convergence_mbar(gmx_ABFE_complex_u_nk, estimator): + df_list = gmx_ABFE_complex_u_nk[10:15] + with pytest.raises( + ValueError, + match=r"Provided DataFrame, df_list\[0\] has more than one lambda value in df.index\[0\]", + ): + _ = forward_backward_convergence([concat(df_list)], estimator) + + df_list = gmx_ABFE_complex_u_nk[14:17] + with pytest.raises( + ValueError, + match=r"Provided DataFrame, df_list\[0\] has more than one lambda value in df.index\[1\]", + ): + _ = forward_backward_convergence([concat(df_list)], estimator) + + def test_cummean_short(): """Test the case where the input is shorter than the expected output""" value = _cummean(np.empty(10), 100) diff --git a/src/alchemlyb/tests/test_visualisation.py b/src/alchemlyb/tests/test_visualisation.py index 8be24867..f2535f8e 100644 --- a/src/alchemlyb/tests/test_visualisation.py +++ b/src/alchemlyb/tests/test_visualisation.py @@ -8,7 +8,7 @@ import alchemlyb from alchemlyb.convergence import forward_backward_convergence from alchemlyb.estimators import MBAR, TI, BAR -from alchemlyb.visualisation import plot_convergence +from alchemlyb.visualisation import plot_convergence, plot_block_average from alchemlyb.visualisation.dF_state import plot_dF_state from alchemlyb.visualisation.mbar_matrix import plot_mbar_overlap_matrix from alchemlyb.visualisation.ti_dhdl import plot_ti_dhdl @@ -147,7 +147,7 @@ def test_plot_dF_state( def test_plot_convergence_dataframe(gmx_benzene_Coulomb_u_nk): - df = forward_backward_convergence(gmx_benzene_Coulomb_u_nk, "MBAR") + df = forward_backward_convergence([gmx_benzene_Coulomb_u_nk[0]], "MBAR") ax = plot_convergence(df) assert isinstance(ax, matplotlib.axes.Axes) plt.close(ax.figure) @@ -218,6 +218,47 @@ def test_plot_convergence_final_nan(): plt.close(ax.figure) +def test_plot_block_average(gmx_benzene_Coulomb_u_nk): + data_list = gmx_benzene_Coulomb_u_nk + fe = [] + fe_error = [] + num_points = 10 + for i in range(1, num_points + 1): + slice = int(len(data_list[0]) / num_points * i) + u_nk_coul = alchemlyb.concat([data[:slice] for data in data_list]) + estimate = MBAR().fit(u_nk_coul) + fe.append(estimate.delta_f_.loc[0, 1]) + fe_error.append(estimate.d_delta_f_.loc[0, 1]) + + df = pd.DataFrame( + data={ + "FE": fe, + "FE_Error": fe_error, + } + ) + df.attrs = estimate.delta_f_.attrs + ax = plot_block_average(df) + assert isinstance(ax, matplotlib.axes.Axes) + plt.close(ax.figure) + + ax = plot_block_average(df, units="kJ/mol") + assert isinstance(ax, matplotlib.axes.Axes) + plt.close(ax.figure) + + df = df.drop("FE_Error", axis=1) + ax = plot_block_average(df) + assert isinstance(ax, matplotlib.axes.Axes) + plt.close(ax.figure) + + ax = plot_block_average(df, final_error=1) + assert isinstance(ax, matplotlib.axes.Axes) + plt.close(ax.figure) + + ax = plot_block_average(df, final_error=np.inf) + assert isinstance(ax, matplotlib.axes.Axes) + plt.close(ax.figure) + + class Test_Units: @staticmethod @pytest.fixture() diff --git a/src/alchemlyb/visualisation/__init__.py b/src/alchemlyb/visualisation/__init__.py index 6955dcaf..6aa4c358 100644 --- a/src/alchemlyb/visualisation/__init__.py +++ b/src/alchemlyb/visualisation/__init__.py @@ -1,4 +1,4 @@ -from .convergence import plot_convergence +from .convergence import plot_convergence, plot_block_average from .dF_state import plot_dF_state from .mbar_matrix import plot_mbar_overlap_matrix from .ti_dhdl import plot_ti_dhdl diff --git a/src/alchemlyb/visualisation/convergence.py b/src/alchemlyb/visualisation/convergence.py index b74ad336..bf232b5d 100644 --- a/src/alchemlyb/visualisation/convergence.py +++ b/src/alchemlyb/visualisation/convergence.py @@ -56,11 +56,10 @@ def plot_convergence(dataframe, units=None, final_error=None, ax=None): Keyword arg final_error for plotting a horizontal error bar. The array input has been deprecated. The units default to `None` which uses the units in the input. - .. versionchanged:: 0.6.0 data now takes in dataframe - .. versionadded:: 0.4.0 + """ if units is not None: dataframe = get_unit_converter(units)(dataframe) @@ -136,7 +135,7 @@ def plot_convergence(dataframe, units=None, final_error=None, ax=None): ax.legend( (line1[0], line2[0]), ("Forward", "Reverse"), - loc=9, + loc="best", prop=FP(size=18), frameon=False, ) @@ -144,4 +143,121 @@ def plot_convergence(dataframe, units=None, final_error=None, ax=None): ax.set_ylabel(r"$\Delta G$ ({})".format(units), fontsize=16, color="#151B54") plt.tick_params(axis="x", color="#D2B9D3") plt.tick_params(axis="y", color="#D2B9D3") + plt.tight_layout() + return ax + + +def plot_block_average(dataframe, units=None, final_error=None, ax=None): + """Plot the forward and backward convergence. + + The input could be the result from + :func:`~alchemlyb.convergence.forward_backward_convergence` or + :func:`~alchemlyb.convergence.fwdrev_cumavg_Rc`. The input should be a + :class:`pandas.DataFrame` which has column `FE` and + :attr:`pandas.DataFrame.attrs` should compile with :ref:`note-on-units`. + The errorbar will be plotted if column `FE_Error` and `Backward_Error` + is present. + + `FE`: A column of free energy estimate from some X% block of the data, + where optional `FE_Error` column is the corresponding error. + + `final_error` is the error of the final value and is shown as the error band around the + final value. It can be provided in case an estimate is available that is more appropriate + than the default, which is the error of the last value in `Backward`. + + Parameters + ---------- + dataframe : Dataframe + Output Dataframe has column `Forward`, `Backward` or optionally + `Forward_Error`, `Backward_Error` see :ref:`plot_convergence `. + units : str + The unit of the estimate. The default is `None`, which is to use the + unit in the input. Setting this will change the output unit. + final_error : float + The error (standard deviation) of the final value in ``units``. If not given, takes the + overall error of the time blocks, unless these were not provided, it which case it + equals 1 kT. + ax : matplotlib.axes.Axes + Matplotlib axes object where the plot will be drawn on. If ``ax=None``, + a new axes will be generated. + + Returns + ------- + matplotlib.axes.Axes + An axes with the forward and backward convergence drawn. + + Note + ---- + The code is taken and modified from + `Alchemical Analysis `_. + + .. versionadded:: 2.4.0 + + """ + if units is not None: + dataframe = get_unit_converter(units)(dataframe) + df_avg = dataframe["FE"].to_numpy() + if "FE_Error" in dataframe: + df_avg_error = dataframe["FE_Error"].to_numpy() + else: + df_avg_error = np.zeros(len(df_avg)) + + if ax is None: # pragma: no cover + fig, ax = plt.subplots(figsize=(8, 6)) + + plt.setp(ax.spines["bottom"], color="#D2B9D3", lw=3, zorder=-2) + plt.setp(ax.spines["left"], color="#D2B9D3", lw=3, zorder=-2) + + for dire in ["top", "right"]: + ax.spines[dire].set_color("none") + + ax.xaxis.set_ticks_position("bottom") + ax.yaxis.set_ticks_position("left") + + f_ts = np.linspace(0, 1, len(df_avg) + 1)[1:] + + if final_error is None: + if np.sum(df_avg_error) != 0: + final_error = np.std(df_avg) + else: + final_error = 1.0 + + if np.isfinite(final_error): + line0 = ax.fill_between( + [0, 1], + np.mean(df_avg) - final_error, + np.mean(df_avg) + final_error, + color="#D2B9D3", + zorder=1, + ) + line1 = ax.errorbar( + f_ts, + df_avg, + yerr=df_avg_error, + color="#736AFF", + lw=3, + zorder=2, + marker="o", + mfc="w", + mew=2.5, + mec="#736AFF", + ms=12, + label="Avg FE", + ) + + xticks_spacing = len(f_ts) // 10 or 1 + xticks = f_ts[::xticks_spacing] + plt.xticks(xticks, [f"{i:.2f}" for i in xticks], fontsize=10) + plt.yticks(fontsize=10) + + ax.legend( + loc="best", + prop=FP(size=18), + frameon=False, + ) + ax.set_xlabel(r"Fraction of the Simulation Time", fontsize=16, color="#151B54") + ax.set_ylabel(r"$\Delta G$ ({})".format(units), fontsize=16, color="#151B54") + plt.tick_params(axis="x", color="#D2B9D3") + plt.tick_params(axis="y", color="#D2B9D3") + plt.tight_layout() return ax diff --git a/src/alchemlyb/visualisation/dF_state.py b/src/alchemlyb/visualisation/dF_state.py index 09d685a9..c050cfe5 100644 --- a/src/alchemlyb/visualisation/dF_state.py +++ b/src/alchemlyb/visualisation/dF_state.py @@ -1,7 +1,7 @@ """Functions for Plotting the dF states. To assess the quality of the free energy estimation, The dF between adjacent -lambda states can be ploted to assess the quality of the estimation. +lambda states can be plotted to assess the quality of the estimation. The code for producing the dF states plot is modified based on `Alchemical Analysis `_. @@ -266,3 +266,4 @@ def plot_dF_state( leg.get_frame().set_alpha(0.5) return fig +