From e51580561a6cfe52fcf2b4439bd863d65384ca5d Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Sat, 1 Apr 2017 17:30:34 +0100
Subject: [PATCH] Added implementation of weighted histogram method
---
regtest/analysis/rt-wham/Makefile | 1 +
regtest/analysis/rt-wham/alltraj.xtc | Bin 0 -> 60388 bytes
regtest/analysis/rt-wham/config | 2 +
regtest/analysis/rt-wham/fes.dat.reference | 57 ++++++++++
regtest/analysis/rt-wham/plumed.dat | 81 ++++++++++++++
src/analysis/DataCollectionObject.cpp | 4 +-
src/analysis/DataCollectionObject.h | 6 +-
src/analysis/Histogram.cpp | 10 +-
src/analysis/Makefile | 2 +-
src/analysis/ReadAnalysisFrames.cpp | 27 +++--
src/analysis/ReadAnalysisFrames.h | 3 +
src/bias/ReweightBase.h | 2 +-
src/bias/ReweightBias.cpp | 4 +-
src/bias/ReweightMetad.cpp | 4 +-
src/bias/ReweightTemperature.cpp | 4 +-
src/bias/ReweightWham.cpp | 122 +++++++++++++++++++++
src/bias/ReweightWham.h | 54 +++++++++
src/core/ActionWithArguments.cpp | 5 +-
src/gridtools/HistogramOnGrid.cpp | 4 +-
19 files changed, 367 insertions(+), 25 deletions(-)
create mode 100644 regtest/analysis/rt-wham/Makefile
create mode 100644 regtest/analysis/rt-wham/alltraj.xtc
create mode 100644 regtest/analysis/rt-wham/config
create mode 100644 regtest/analysis/rt-wham/fes.dat.reference
create mode 100644 regtest/analysis/rt-wham/plumed.dat
create mode 100644 src/bias/ReweightWham.cpp
create mode 100644 src/bias/ReweightWham.h
diff --git a/regtest/analysis/rt-wham/Makefile b/regtest/analysis/rt-wham/Makefile
new file mode 100644
index 000000000..3703b27ce
--- /dev/null
+++ b/regtest/analysis/rt-wham/Makefile
@@ -0,0 +1 @@
+include ../../scripts/test.make
diff --git a/regtest/analysis/rt-wham/alltraj.xtc b/regtest/analysis/rt-wham/alltraj.xtc
new file mode 100644
index 0000000000000000000000000000000000000000..9b30b5aa4c53cb53febd8168250f9794b161ce8b
GIT binary patch
literal 60388
zcmdp<c{tSV+yBX)RAeVH_I-<Jkv03i6J^c5Ws4Rm8e^AaPek@8TNEKXSwgn#WZx>S
z^gCyY+<$zJ=eXzj=lTBTIOaCj?da?^@9VtI_jz5P&rCQtL~S@YIMndOJfy`N&+fcN
z|J)g$Eo-DgKc#Wc6v6+dmiI&dmM|L!2OoihgIkV+Lu>?N7Z_8)FVIg7UlI+ZjJGsB
z>)fUos^jcS#Fkv5-?~Aw?z=#Dx`=#ZoO*`qCgI~scTA{D^5o`A35b+Dl~U)%*tE<0
z`6Y8n5cWpH5@X~SKjj!cCq@5te_5@i)Y?0C*f76>|Llza18XX<2>@#hu(|>3Eg0{I
zKZbtqntBm!NWX-vA;g0to1rx|%80iu)TsWXR6XIxoox=fw~yOd4vTf}!&m0LAcXAS
zC@jVq^$yABSoxYIWLMH(yqfkJ8AvZ+z{@O=ob)sv9dcLKW~ugn&1wQ{7QosFteU_I
zym(_ftQO6qzS0)1^nxXfhNx_o2}>4*E{8GFSq+h<&OipY)pv(KFE?qEt}$krFt^cD
z9uf$C5%|!<VkwC7s%k93#iw74OG>$Aa~QvjAZkwIqNj6*6_@WUT|(lHj_?c2V`uyy
zbPWJjI$-q%RvBOog)v$HOb$F>dy?%qQ#Z+yaiTwTQ_FZVF+w+Y*1f4`tk!wyVAy9#
zMcP*Hm6EGgO@7C1TApRvA7}3}m->1^T}j5R;fTmkUi+E#8zTqF9uF&1h{E5rhc)W#
zzs~^!VB-f?Utkpk)=bd#BLqeBYR=W0SGZq7&I*ch%%hfHeVdF*c$9C3WNdGL+_IlL
zJ*@c9*`b#|`G~1TS$qf0?Y$W8srqD}I2A8AUa)q+fYVeYbG^)GAA+~}I88CvWXWFE
zPiF)E<|q6D^Y|Zhl?PT+V0{3rdVg}D6aFBYRb%qFd%S+xhsj91X<-5dt?C53vOL~V
z-YRd=f&<Ye^Or1V?xo6z-wg5%-nyF;c68O(V7usU*?evuAI-K1cV~0yYF)C^@GwL8
z!|1Q3TWFk}9N1@cHsN3AfFZE50PAC5H3n8k7<2q#)l>{E$ac1*L+}Yjl!gAh(J;Au
z-XM2V^8VO)&LE~i9PW&7R;=%AuEj@ZI9WR3)P1R5i|;x#FTNZr#Fu>GoWNLbX?B~*
zd7q&n0sf@)$j^IOUz~j=D!QX1`dy$OJLCV5uP9&xZE@>?RRdUcVT{%klLM!6`h2fi
z1`Y8@+zt+9l-xqR$#nB;ogbuYBsw%quT|akYKrBZ4Wm+SPJCErexB4#MOlH-lN>AG
zPV&V^&I;SlmJGBc`ZbvzjlLH(n-B+KcXY*LKKt+M0L(wW9k3Px>v>=egfR^S#jLIB
zNYQNdnL0}Z1J%2bh8*QdH{&h8_b-lUX5Pvyi%Y9iQe})59|&SVWs$s(Axn__ctysF
zN8sZl939<RXWgR<3qM!1$n;82Yv5+ASp-1TzpS^<&f41U=!pIe`mr<q54ys<<I4hT
z0<c~HR(DuiU*M0SbyaILV++Actm-kAjlb29Y$8zVV^lDfga7@a)D3NgbsDwe`51-y
z7WSkmw%fNqKj8Kz+)$)niAkwTRBBiF$dtWgsE%V4EulxxR*+w@hjUJP4;Q?|JpKnx
zTi}G;#4P|$9pF@lF<MW|{8v_I#A{6-9pz4Vd=7c{9d0_)+YhJ(ZaSQq$3#vf+U<qJ
zlgL3wh6$P~)j4OLP@)UpmLxBVmn|Ha$+srhWba8UXPg}|@b+hLj7hS-p$}nq=Kr5;
zzeoS0rtP31=!yFfSe<Zi2-tvi0RAjy{#mzr>|71;<M@V;)rT@Uu_FlZ=}~j6vSIh`
za9K4MEq{1Nev{Zji$x~kNc8<Sc_pW7&pLaumhJ?4nzP<^p|>Zo^A>PnM&(>J4|jYO
zfF{~y#l~Orfi)Od-vBEx6OjF3#mxW8B#Axqq<zlBQQw4{q&+-SEXXIgY4+#mOmVH(
z4*Z;VU+zIsao<1d??*PT)0+2{?E6ae3Qs_X@6V40xgXs)jy*flRrmdC{bx>6u05RC
zcq`01UJG!p0%tLB3Ik^iym`$0D<z35sn4(V3)(iug{%=4q~EEWaad)hYq^@EX_s^7
zoZ`*N7y3B&y<`FkyioDiQLVw6Ok~B>b0y2=0d5SemVtdtZ6CBL+YUF3cW1VP6uX*Y
z<*jIcdWD0t540qPzJ&l?XJ7BKW}6^2g{+6T-X7~ip3Wk%Oo$;uwb~i(M>V~coO_tu
zTRi1_EBBMitCYENp2gMB_Oo3Za*vso<*T)*ueaviSa^N4l~cxslxkf*x`NDmFDo|Q
zdKXyX4euidR(MTt4OmNeSiR{Mjm&P{W8ypBdII^Wz8Ya(>yDZ~pH!C=e9JXt+*n<C
z%DJj%bVG@P(50a3#et0`<~Yw<fd###&Sn-3ex3|@azx(wiH%lD>!EZzjlHbccx%!h
zU5@}O^nnB}z}gH!G3)9WlPpr^o<XX>g&#&CY*yr!&Cm5wpA1p&C}l+lDfT-}bW$~R
zIi9x5i;J0g9B_dux<*%-r*PT7Z-L6;_4-|pmBMT7E`@duQ&MMjdIKQpZVq7Ot$Dx-
zI_?8D;xJ&n4y+vz6qB#VKaUxy#q3YzT^#5~y&cc8^c=%Y+?1diROW1b>8qdbzhF5=
zW*vd52{4L85zjYV$lPcgDau;0Tj>&Q$`wp`EZ00~r^}>Z)<kJyYwxj_6&r7T3an6{
z_Z<RO@GS!WKe}RA3n>zP<=d~&Ge40w3MD`)Sl%`|i)!UG(}{9Y8zDW+<yU7lNMhYx
zCXJksaOf3aB^BqT-<VBZRw>h6<o3(UE0@D(woli_kDgU4JqI@SS68gO71lmpHLwx_
zE95Fc3a~apP)rU;GIkp|)W@arG<-@xYS^(%mU#0dZX_x8>+chej|}hn{PyVZRot3p
z55){6GWA!Xq-4d!<pQ6cUPa;Lo~b99y;~~V^g#B^m)MFHo6P`-x|;)d*mx_tX2EkG
zJFs2`R#;zu>%g(&k9`%xHAe*!=8+-a2y!fI9N(jwloO5E>{?2%aJE}jT3OM#d?T(2
zK}8E{-SlK=*%)M=5LhBFQaiwI&RtwDIyJTEDdw5P^)>$&{`6i}Z2T28#pLS&V08di
zGhlrQK`}Y7%%;}$BsyR&SuN$NVtyV4%c;++=c|7V9B*-F{OCkAdNO`9T`_C93qg5r
zDCcO_tFO3Cfu~o!%k{EJR#T|Logbpy^k{u-pQY9w2vC4o`>QMZ^<O{zJ^DNT4?RmZ
zu>H5}Uwwg90(|Q}=;{o8%N5Xh4c;e4*Ovl^;(fiyQ<>NOeNa3k-z?W3JiqD^U=vKl
z<IIEV?HCVZzb<`L@+#F`JqgW|)Oydt=$`p$R8^V;*4DjEY{*q5pCGhXZav_Yzu!?S
zYXHLT_5#>=tP-$#18Xj@8Uw2ru<qo5nuZPE6)M@91l~`+<z5{_9NU471&`l9+KksY
zb~7ZZ??CZ4ZI$U$<@;8KZg%@Wd#y6_Jyc{RcFMB2QB8jK+ksm5JM!!Lyy5dWq^w4N
zIkECrUEov%PVfL?HsB10F~g3g&WrWg8r*&9Y(iOvp)#r+=|ktda(hn=w#X{sn@bWJ
zEflYmI+C9LPV(}|*_#zNN2S+?o;JCz7nJKd2BN0YdKGNU_POZM>|Y`7e1h}zANy=<
z{1vUK3b5J(t2(eI|6XUm=l^R~Xts`dhb^~A1wD#NQGasAFdMa|%1c`PA&4rhPsxvE
z$^pOah=SXlFO9s?gAZJbVphzP&rS1`=P?V3lhrO|TjniS+k0%>jh^~+x?|7$<6`5l
za=^+7tapGF<_M!J<&Li0zJbkd_1hc*mbsE44^E0B?jGLGSMGRl?;VouUF$=IU+Yc|
zX-HpUs=iNOOjHQVI=%KMO5A107k^ZgW-?HI0o7B?VK;G<YHmM?E`Gh+9#*XU74}v5
zLcj`pE6lnI09`kBSmRhCvZWKb()rz(6Oe^cgXvvQ`IFDf*dU2=T!ML-GshD04|?%z
zaEZ+aJgG&vQx(1Vx@~=R`02^nVjYjfm)@#}3WJ@7Y(>30=3fTxgL1RGuJ&Q$uV;W2
zPTtUmHpFJ2s|~QC4;nBvpxqq$mR?yV=7u~U)RPI4boXm4P8PBs+nkAQmrfVCoO>nq
zblSosMy&2fxXk@&imgk&3Pui9g%uUW+Sdn4by06$mOC)cb+}C`$g?kh+QW*Kzrw`e
z125)WNCa5X=Ry=atbTb9vdw2=5rX~cz9_YW%IVu3uahq-BoA-eam=gg>SW`8zWzL+
zmy+s{s=m$_k^9B$bJg=kKOPmt*3>j1qnUEMzc4wA7WaB`%RGO60G95qu2}i2Ca_5X
zE9i+?TTY-WdJ$lBP4f=ORygf3$ssE74aM+{Zt^glFlxGPm3EXkGCGCb$cN5y#h#o%
z-R-VHa65bSonBkFEHx*o*-*c{$hCB;+U>0FM(1f`zs0cX;n4@MbpNtq<F7F9nDdJW
zV3h}5k$ZJD<Cw`-_7q6xIr736+4-n%a&<^3S^nxJB`v%9Hk5A0*JD!q6Ff?Yslp{~
z`tV)iwjKt08?O9F<X<T$!Wpn2+&5mwTRIZt8Oq=lrLV%ZhZSDJV`p#ud-QkwAN;i)
zbo_7Gzrs4gy9qgfJ{LmoU%5foo%1%$SiV+YpU+5oE<+Qs(Aoxd#OzWVs=2PWg<ye?
zPBAmAz+<x1B)G-XB)~J-PHO+HuI@V?q~sjivfU!67TZw%rRL)klJD$Jzp*^;N#O;y
zxSInw*!XJ!u)<yo7Zxe;6JXT`R^a}fuO`F;cD9G===l^j*-$?{ohJ=>a{MN|?)QBC
z5VoxpxVBMg$$fHdOIVJqx|H94^KyotuCUx&>bXx&Cz6{PjtwCDisrtkzosAmemQf1
zV0jNKRvw!ItZmS{k^*Zcu%gefnjk2;&(=DnpJitlp3B5qbo&F6_MnPoNR}UJy`YVn
z<&~3b)J2sKZ$4bQFtLL3%0GP#zuBIV0GBu_VEu@9JOvk7ow|Xc!o!dty>E<9OW!cD
z>%BP>@K;w{Y&;ePtV6(x3#>uF>I$r{cUYN`4ZeOMaia)7b1G!WWF2DgkZ*`5n=;}U
z=bKXMe)?wX?uOMqsq6cwVq@MKiU!6pw43(&3RL9`&Yaa)<;oqK?g*-UEdNDBpe4Ge
zE(sRbPLG9^$EE=*=(vv_SmS^dYQb&}s8aQD_3WRzDUNTSP9zn%epn&v`*I;NF7Z{7
z{lOt|A9I!jA5z+4eBGdp<5eZE=5)=4a!sQX$6GHxe*WnUKIQ7b4+Rs&;HVj86(}>i
zoLKp5K5)iDt^zCOTmbgm_>jApT-9ivx@@G_x=Nw#AgIytag?WG&FBe1uZ4`X`1;bz
z1+Jrx*JC3Y&yd~bs8SN4mJw(k`AGBYh5r_rPD|os(&p7g&Vh^xoi?vvW-{w(!8wTf
zmlYd-?E=;;V2uY>n0o>VV4eKKs;NqCAfdf5PR(6W>5ba#3Q3<IvPU&H#2jGc%~qcM
zJgrADH5HLq;-Pha^x29-Fyj5ZnJ{&+z>ixrg@mupc_GeK5#2mV?!xxU{E}zhIJ9!R
ztXTOgTo*w1)VNc?S_Q1=bF5c?Sk<bgOzkN@eN>M!^nDm6(!<Rd@AE|5DBSGgF>2yb
zT|IZ<K}3yJ9O{wqW@H*ISMk)6Gtr<+VwJsyT7r#BELD7su3Yoi7g8z8T)caFY^=N$
z&Mxu#AWskf;Zz3B2H?c3|7-Q_zUS&Y68R7s4Iv%wGL}D2nq6_tL*f*d?JwIdd!lb6
z?8w!_Z{ST8t;62Q)x9tlpg`%UWBT-GM3>R|g@TWHWYkodccuB*%`)sl(kk|7ijB9z
z-U`EtKEr~%B`^Tiom!=N%=!-Ua))Fh-?_ZmkPl_^mZo0=u6Ws#zV(PW<gZOtWpI0=
zp6{eW(ZSBRLPa9B&NnU!KMWNwB$#;JPAF%6&E7{BA!r<KBKlnFjjEHwz@GVsm+<(P
z-kJ?;2$=tJV1+sQKYi`>BCw)+!9-w%bw!9i7b1hVjP9>A=zT;&sBnJ>^Y4#;%w(9z
zJmG^Zyb$}f)y|yenqPI3bzc7)-52Mcomdun^99F(>mt1Xy{Jpoi>3J5yT?r$-eer1
zt*}7x9`Hc+D<7EyCI03BR^F-wZ0H{QIk0;FX2q+8py+(nINhv>bS`Q(7jzZoLh?M!
zN#E{$U@|Oa_J)^`>FbDvpPuSb(WxJ=QEQ*c*W5xQsq{iMn%6kixQ4@2Z51Z*%v+W+
z+~yNJGO64+&!U6=WyQ){m4KB6Sm8_%)3=7f7|n^+Ra2yVyTVK1K^AwU)*x!6hHf%e
z`A5FtKBl3M`vp-~*onK(pJYAUPC-ny`Q5WqT&`=q?BfEvCPQ&eb<Hm;rA(Eceiuo*
z`5U2EMhaa%{OaAyijB9b0UH!OTuWe80A1r?t&RL))d;;)pY1}cxDx4HjB`Jb_}#W|
zROe7SZ%B?;PfR>sy5QAMBopKdyiAkku|dUD9b6Wnny7^nEXso>xQzMFd%a2?@H$VS
zXwpj}kG@~AgQ&Z;6)SIrT*c=F)=R)D4Xlq~O#Fvc^Xd1>W|c$tauGlEhfpgIMGzHZ
zMwK759?4Abx_&y%()u~Z>PvH;7~ec)5U!xjla7PiWSXr!2D46*y#7BXQsrXH&Q7Ha
z+3a6FNH=}FBYY1lR{o0KYeUZBg2xh@1FP5X8uq)sUQif}m-oq-<ToBl2#xP}Kk0SH
z({J!4t-5W{MG0f!KAeNI-X*-`9=ub^pN53h&purw4v4l%@+zmMW<HijVP#}0=&o^G
zzaU)Te#<Euzkt8GVzqxo_ft;53VIS>1y+4v-2wqHYm2jo4ymZ~(}agJKnZoeC(}}`
z>5IwxZF&`6k_%rG5qK{YDKAjDj()s&iAn4UcX{@x2zki|rp+rT9h{5TU!;a#npEy~
zbc`#dWu?>T1Jb{&Sov!RuwDb!NMKC@R%;l8hy7-Kx9_8`&sQcBK6$-6s9$^umejVH
zCL>bkiG;K*p0iZV9KLf*=0v1Ja;Tiaa&Nce$rw5w!(r-Q@5&4n98L_|w#eH`8eSgl
z!+%by!!!9$R&36P?*SXEeO!n{oC>T~Fg~`!TEqF#PBoG>k)t+FDO9$**fM~Y!(^;_
ztVt;P!gC87!)thy7IUZW@l`mEUMyJ6U~}d^(*xT99z8LKN}s#E>xeU%irN-NuY-)w
z4k$v<J-TA!u@1lnPh2>|CQboX6Bxr8{BK=5gDjC(giuTfBAI4nsFar_<zrQcF{$KM
zLDEUdqal*NnvPx~tYPJSIxRBd+--Oaaaff^Ep@8h{bN;a^+OkX4mO=OCuut7_OJe0
zkN?4n8T?z1ZHN5)Z`r@XHB7uf&=sz05F?>?VFz7b!MnzI!7l?vU#F|=ar`-#rbGIY
zy%0)a(WvhlI<5D(=P$pb)c9F9#pn|Dw&)$*y=#ZIO>)QbocXQevbJP9MNGPa1G9#j
zHw8~DT2nT&S$jHo?uW2Ddl#%cwh&ms{%~P)N&FI6!4L4e|FEhBgh&|aapL#u2OiBB
zFJt0)nQZuk(8*o@(ZNll9<p_9VbSzbJM&N$iBDnQPZB8g#IF|IbhS+?brVC_QI?bk
ztpr}{<0;(OudS>W|CbXhe@z3<4(MHpfD>vKX7AbsoalO`RaM09tGztS#3@xtf+X}2
zPJeL12erHr+WQEXSEy_H9kOS>QoSTG`>uZ;BRzuUXD`f)Ny1z2#<*<Kr(Op4Ivyf7
zkAK}LE&Xhed42>W+0_)Qy(_#|yis611gsIj3NDBb`vMHBCaIr*(RIlT6Lx3Zc%%%~
zVEU<{OuzMImA<T-!o1xf^bDPpl-t4g1LCV1yMD%6rQL7Ujy}om|JAJBE=BFh&ud<d
zRWGlp4fEpGbKF?NW!}q*&G|6YC_K2<whs?jUjeHru)-O^?_5PBJh3x4`XZ6rhOq*r
z!-tzr?e*Djantnu)MG-9=8tD@T3QXaNavQ+@88(u*O(QZGjjIlxXZNY%~n2Wx{%r#
z&{tIGO@<1qqcHVM-h|`wU0t#A*Qdam3%NQ7tWaYK=zw+L535EE?sMZaxi>#h+W9`a
zF#xv>G}9E$@~M+}N~I%nv?3IEC0rg@W`!O0sK(`wnGj$LlA8(S@(-r@**ZdEB2r(^
ztb;c*l(k{6Pu=Sugm#y|tk~>dTYxnWSf2r_C$NGC;P*jLbiS(hxY0#DpmkOXWi&Z%
zCil^WMHO}IeEXd7!a&OT=Vv{-2^w0t*9?P=j76$M+~*4vEGuptzU<~J&|gQ6NZaBL
zihIM4dQ0hlYbE&B;{XWwmlZ32eFCgez>4m(Gk_KP*xl>e8u!Of*^}!vH?fK$pL7rO
zb2k_?7p`a-Q$>kX%0JwyvvGWII^-l(6kgAEq{6FuLALsO+#hcK0Y^TyI^LmYe&Bel
z=hnDw^4n2L@^_}XZ2w@z#$W&JwK4e$^%YZFck|VrZ&AsBqO_mSa#Jy6=23@b#{)sX
zmW%OvdLQ^cX1^L~k+mFLYtL{}KG40%-{NJub)Ur{|1I^rPR@zjO1uoMxBH6ny^mLJ
z*}e`i#%UeetLvOJjQ-`mFb#AC-@>p$?h-hGt}w=|Ej6{`aeFDzBLe|-$h-sIaeQTL
zUkf+8!c~a(Pd*wGr_&a;n69_3APip?YnPEyHWk=snLe>1_(7Ezhs=jk=Qu~r!ay3;
zc<B?Gv5!y(cRc{Y|7*@{;QDXbzrvmcUlno{eeLxI^eoV{<HI@D?|Q|LGn}nBXOYgA
zoPk7{=ub>`ZwsPUES-Wz1e6<+yNhRQEjM+E<WJq?o5_zHqWw^N%)!%2x^R__L&!cZ
zuK03R1EFBNLqrk-%h?zk64-M7)fJnw;j_R7jF@Xg&<7AguZ8yuf}(4H#*=HB$O|&j
zlN=w8<sh4%ew%d0A2;d!kxjf5YEqGg47k}zM_bQx%;MR#Ns?>!l0@~D_WWT2v*+=s
zW#28y2h*K6w_!|1U-~mBL0GH+0{&&i%3t+?RRvfpKvzRxg*^dAS2U}pgQG>ft&q7n
zM>NZw&}v;=%QkaG)Ob}-WA1@ri3ut4gKwS|iiS3OkXJDr$Yvuv>L2PqnD(7vgE}y!
zPGQ7~rCXkaMsEKM`Cy|(YMGw&URG@Oujn2H-CyMct2?kleZ@m>iZT6_L5%}ascP9A
z!NYPTBuiR((jsESWbwMp#QFAxw6$v7QaY=Ix;%0rYkD916}bz1(!^a0oUsB+F5dY;
z1hzRXqop2)7A!lb7<&?8j=e11%ZiP^ngOd0u%dfG*uxUSJ_aB9g5NoSOd3F{+Ul7L
z$nF0U;!YNBsT)X(n(!wNSs=0Wxo_0^sqLs6<@Z?2c&{XuGiiec(@GkyfdZR1ZP#wl
z^sO@Ki;~oyJI{OK!dTRS%Xk40bvIwJ@>leEAo{!oeU1e^452OP`fi8ygt!z^xi7&S
z;moZRsx(HIPIElUWU!i<Xepdh)R9Jxl9b_XZ+3=`bbUX|t-A^$=WkuvJe~FH3M=Jp
zBiogOipMkOjPGoBZ(cK|6#0fm{mY7tzXBUR@Zx3yE6hLo9#{0e@4s`P@klJvN@;SG
zy_enyc~3UnQZFjQZ!-J%p>I(d%#2=M3v^SrQ&ULi-$%H1IJLU+yHiTk7<+T9RJj)2
z3pOYsu9Z&@j`xTVWz3X_j#(<(!-|!^8UULIu!7GLhXO0CHM~v;ie6h<&jgH+t_wk_
zd>+Y7sPiEq>FRwfCWB`t74y~UE7HDMji{6&@HFj(@4oO$4)DXLsPi0ZW6#a03a4}8
z)OOGNGG2@)Q+98>ddOZ)x9%#9^j}@EIwOV}g>M6_(7O;rttCWz*1<ol8jeK=kuJ9?
z5?L=_qz#oLr?T8W{?%dl=RKl)GP7??Kf>ooX*7Dz(o$%9wwh(=o_(R@EMu3MHL0{v
z{=t`HmEwY4_Zj`B!(Ya|rSBQF1who@9KgfIUoQX~)KyH61${R9+;#VU0kyL+)W~yr
zh;*$VS~sW$hgeAU*$Z5IP6xhd$U|}6>o`#SV^uD5>(Nu($9Qp`=BD04mwY!|Y6Pyi
zx?VFo#8$+4P=4yjTV~15_?P0)4gRnE6+HjHX77q-y#<;+fcbZZp5+8+4882{b;apy
z?yGXi(SZ5At2Gjb6n`?wt2<0pUvFTVr=#oL<mUs)sEIGqcf!Q_g#+7CZZBWnyhk%|
zGU-c%Q424Q$-W2vEWJr%pHGP36c_gq?&sPw|5)u^!K?7V190C0EA%h~qQJWA0~ZzI
zd_ArRQ?qIddm+7i_z^@l>wa^o2Krx=Zc!-4Z9Y-4-0Z2|T6{yN>*LB6*81%RjqBnx
z^;GM8IcI%bf_hY6<OfbOJBmS~xmjN#Xft+m04slmz8bF%dO><%y$7snz=|#`XwTAk
z(KxhLbt-A@aqFWc_0a}uwr;UQp_X!&IkVSBX%Go{EWI;{u7}MOJkrlSweh5&P&!4g
zGhZ;^-D0Qugi3F$8+F4(;})Wf-S`LPxS<n7{L6}szk;Tiy*3`O`T{HT1DJbq(X8s6
zSEYKGzF%b)%}9S>RvWEzNbmcKGqFy_&E@=v+ABFn-}5CO#98|kAa&F5;@P=o97G`1
z3ICG&-!v<8a??1M56;l1Bwbb${>GyqUp|0F{mY7tzh(g|^uL(Bb{MdN$Kv-xP&BI=
zlb|zFpH`>;yyFGewICB-se(I&HEfxL(gHX5sfWYw3E%9fpxIAy%In9u5b@AZ5B2gU
zcb8{}?knTcTs1uY%);jFz!_8hQ!l9%7XWpq24Ll_Xie)OS7B2@{1`Z)hsEy!MzjZL
z5dFC4YZuvS!s8^@i2O`=6VV#Oi&_v=)HJ{4z!w%4=yR85T!lh4U}&cPg=UVty}SGZ
ziFw>K!-()5ii1`s?D>P0&(w;vrE|TQnXC_hh<{nJ@zxSxy$h`1QN)SB$_lLA5EQMc
z=6F}Duj!?qaoj%k)~LXEZc7PYbyO$)-4p>5grwm{%{Xe(_EpL3VpsG%(~h=hsLB>j
zqm5*y@kWcXJO_iVM_O$X3eRiBM((W}URHUyxQ7)hZ*2qCSYSo_74(FddNsJiDpc^=
z_k5TnHTwznAe3{UJ>rtNYuHuF-tM%}*6nhRdu(I2gPZEhv|<aNLzM(=)Q$wl&ZLO@
zZyj-IF|}m)8pOb!eCmC$kBT#HTY}u!@x83r>|J60@xZ4r_hN!q5io+T6FaQ_A*sH$
zJbZC%57llXe;uDe=(S!%bu8-`<MXuQ`_L?CT29{Vk?(ssOWRx3=}X=dSH`-^a)4u7
zj^e4T<2l@iubO9XeluSiCY9^&e*BNmJ7DFla9)LBg|k7-Jvb+U72UI7av-3}9VvU1
zX^8j0wLa9?i}#jKs1m|F<fy5ia0>|C8OMKdt}dOXy0=hlN{PJDgCn%(x<aH4m;d%X
zzt*crdgiJ&=2(HTLtbR<KNW;3x3u=?3NPXDFZbC{`~Pe9ujrmd26Esb=z0!%7U;3?
z`Qa^NbbTBDQAu-T!hkE^SQqJ7NS+?k*^)T)Q2+KxFPi~`KjlI0N&8oM=dC{y4aYmi
zJKD?@=hGbW&02ZBCgikiQ0!CKz~&f#OMut9TuhVb1K7-N9l*+8J%P;*SmEq~81j`6
z`T=}l2#T&*T3Kvd5=vY>ah#iL2B<sVqAhEb-28@K#K|}&j7nY%)7FV7w>?sLLwlo^
zh*^4;kRr?(-$Ku=smkuyaaL95d>O}X3g-Sv^$Ut_6~?5VmH~fRvGP~+x#crptpQdr
zYfO(_@`qLZzD`L|OosPBxN_#F3-9{r1+I+{85hnRF3?`6Ztp`7L^m~2?K4h)>ToI2
zSz+{(V;mlZ%@en0!VbvL>u1kwdF-RA(Mh-4_Osl2rzj0U_poB)ub?SD7qDIiR(MSq
z3}YS$iphZug$yIbypT!$D?E}Rwa8J+h?uIx9vKIvbK)8b`Tcfw2&*T=gKWnwlDlFA
zYC7~K76^MBFBh)(UoO&0&MLBbRG-~S_esG_5l>D}Pz+k4-5kK``YU=r5Cp8~{i`dm
zYCv5W-eJA}g{kYNPMci{iP&k~GU94(om<Lg<1Mqn^sz&+u9o>Hj8et9dkZBAGv#{b
zI5W<K#n#M89P<9wYb0-OYOri?dQAJxaP~_0i(lb)7sMdyE-O|Zi(Xsky-Oyrngc6(
zZJ~1&{nW^hQ$l+4tD19&p6Cmyd+vZ3cRJ!X_VT`hbM%F8=GV_Hk<wA<^OX*HFy&6f
zmy-DHm$Aq`o%Q2Rj{B7x?6<p4n2cHsFT{3>)A)aW4hsSvyZrz*9*fQaI2*uRD@Nz*
z4H&~d<9806R4dDNv9_}078_LzC0+bH8O!4A&}$g%rLLbgesd+|bGM}d{YMF!jU&mt
z{KJ<zz6nXBaI0{v7zqi5G?o=RCck!glqB2b%>GE*`FR&i;Vvsy9(x|xynr<abcOB`
zbN?_C1jXpuTU+ewKo~N~wsGl1sO3ZRbRmI7hrw#5+j)mao-;q<2#uup$oh@!;qZ{<
zIE&W&+jpkdrw<o?)@abGDC>K@B;&6)(h;+I`N2jw^JPe`y{y=HtOc+d0BZrT`T*-~
z7*p=BcK%rKwU~;tWNtI63~?|Cnfx`m<S=yo;<=?aDUugVn{I0|F<Z5<=L~1tCclc=
zT19wKCrK8rA=0$p91Y+2QE^i>!-e_V{iemxt)=aNf3RY+e?`|*;Kk&t8?eF|0X{B7
zMb`n1(}ACn9&Vk8(4ve^d-<%YZ$^1B-u<i=WqtX?Nztmwo5Ty{ny!fiwy)ejf2jOv
zTD7Eq`V_~?bAF+aGyO%@vGw0?CS4{AXKJ*x9u>FW%Zh&euRXRM)_nr7qJ0Z;;Q#dJ
zUuS?739R>lRUKHtXaDwDpH&C!E~e!f@REv9hir4>Og6r!3A<i5+QIgQvuiLm#q*dQ
zp~zRiH2ooJ-<(W~5SOp1&wS+lzr}jA(7O!aKfZe`Osb1*RA;8Oc%Zg98<gCs16X-%
z6R^U50JjoYLxB|mth@EKnt#>K@U%)l*CZ`vXxME!1mRHqRozP618sR+(<h=RUnr!t
zyj9^*5lb#3f8+VMbI$D1`))VWptQtcbw_$-LQ#X4ktxf;C>3tJ$}@4XdswmZSm?3w
z;(>JnSfLhR_Av1Oe%Gu!v)V|H!Mq_Jb?<)UND-OkkhFW4%57SlJM)d%MjIyIFTKr-
zJDIXiVwev{PkT1f=AKQEvT=1_?CGXos}V`%F_}H-u`>8?WIKk;GY*XJWyQv0D}c2K
zSYa<f41EBB3h27qzaAG_@wHuuA7U3t%?gQ|eTbOtJ|1TFm{zI3^=4j0M&W*gfy_JY
zUWruc$C=LkKrN@R9se-mXPO4z_aQI1a5wY^h2X7B4hvO#Qqq?a%Ky+6n>}m>u$BYs
z3a|n*fgG^D-_ccFY#w>(VSqV%^mY*H7n6H>z8YVchbSF&J&&ciQ9qH)_0burQawVk
zs2|AJ;f~7}7Q&a^n@zKO)5loT^6y6eI8rGbZ@f%(&qP~A_8~0Y-F(H$V}TP7_JX+6
ze^|wU_1zy<HOJ&WV`ixtPERx6$GD6X2pk!j!dr~DClAQ_eI4{a{9L%QB~vwSe<88g
z1BXi07!OJIS0SzogKv2&9xP9L+#~&}|2-@%8KrT9g%kX7XI)`+e-qRx^yhJK`+*be
znm`&jdw~<<S%*6qjeI-(O?Y}`Xi)=;nUh{W;4{!}P3<-Wxhrp~B(EJ|j650KUa{u1
z+J)E0_HE$eo7A<qrAsG|97*Ts7gb~!TjsV%H>z9uCTTMn08w{YvGP~wt<dMexLd%Q
z53FE=_^l8W(+_@l9cy%6Xn2TqU~&N2ZKR#<kcYVH7djT<s#M=AsSw<9<72d{jobmM
zH~Xjt3Ec`Gjc>({r<y*}eL5<V#&^FEPfz?*%1us@8!0Yj47c{?Dpu#iS-?64tYd$4
zl>}Dwy29u>Cex&JNldJt^CyR4$UPA{%d1}Oev3uOPA>TZmqRlnw~rjI$bFcl^o-E4
zh1g6dNm8Iim@97L7`=YYQ*qX_GDJGA&o1sW*ElW{P{x3k^{=kj_$zu{p>_QVtdO?^
zdcX>`_P4)QHJ&lLvFbLe#kuygxIa~}M(1&1l2HS~n&&b*j!t$luO4&0R#!+Y^%oL6
z9vPlnwjLzj<E9ziC8IKF`Ti$9(lnkNn0=r-uvSsVe(?2PR`lzC?XO@Y|22D9SpWDY
zpexiWVm0Vnj9`opZyB8fnst@4LUSB1vrDQj-<CPy<HG0V&dp1!T75g`Q{gy1-lFTt
zj~yth1z(dDm8o}}C;Bqh6#2Gr8Qy1}J^#_t)YxcYT5{M&^|DvM1uu27CLsv9QwOl}
z*9Bmc2UbsDy$G!6dZh?K(Q8XhiIdY<tgY*8)JvsI1}1Y(^)&KA+)G{hs)PiWujY9$
z7n2+^n8%&!9XL8IrM_%CORCaL>+<mETgrG*l}{1#NSt$(kL0sH$ZU$9Z2rrMmAAG5
z7aMTq0;d9SM!{GL7%}T=d#Tgtyw<fu{vQ1kq5aASEfr+O>4y?u)?XAxBtMsk<58t|
z9XvU6x8#8PhnPkI|GJ82T-{8|REm8=0=z-(M^9|XJiIL`ByigF8bKli-NTBFw?_QY
zGz?h1ffag!-D|~Kh8MGu^0x*!xKyW4g;c3Uq?5@^o$R$bEp3xSdqG0%h{yit=f3lf
z&Qp?H&vhiJj7bpskRzYx{+ynBW5rieN$FZkUPv7a&uqmC@620C$ACSo*mx`SRhau+
zPXjBQ4-x`1KAgM#uC01FzmP6?Otu^hybhsRtJUc*pFTVJEuGir4MIOhl4ty!(&SH{
zVvnkooYgoj+!)`))5<gYO5$6(<SI`j8?IGbaE&gaIPN6zPYFjh+%?#vD^}hL^$PR7
zJ+O}@z5=YUFGPRN9FqfDUW4!Boi4m_;I!`x4e4+bL>Sj2l(onGqw@k68Kq|}1l}I8
zI*Qc4EQ+ffTmPv3QSK8yryji#7oBrXWf!Bh2VHawncfG@i3!-TZe8g%8vm;+R{jco
zALere6~O8Qtgzq0Kf0r<HmR-ejTUlC#Di(S(BV{X#B<eSCmTW+Iha4J_uDJ6wc1-U
zg?!PNKBeV!v7$~SZQo;d!9y)&i&Ni~QViIcQ$KQ@T2Zcl(OT4iC}miNziPJzVCAnr
zfsF&!Rw}S+0c#|TRd;k{KS*n2?8`JHC_}6mTE6XqxVWlEKS=#EAju(x=F4l@Z*^sX
zEnB>Jd<st;-wW71{SlILdT4(`Rc1n%C~1dN>ZM{m9jZ5n`NfGgMlLm`OYdRDV*g6S
z3~b!Mngpz}z#0JKb2z_qK;7>q<JwSKz(CXwvt$$PZ`6pEw-%-j@8VB+RXw{w$9biw
z;mT8$KzAC0o~X&sQ9P9IUzoj?^BbkfrhXpbdMKcBK*hk__MS;XXAw<0#DvFAZN=)0
zcn;X$8~_($VLtbg2IDgj6jQT^D_r9(%0DIYY~D~n-Jmy4FMJqI-*0pp_u|s%(AVDK
z1tMy}7AbYoPc)Sg48)W|Dzc;0YJaF^@cq>8UA<q>B&UdqfH-Sxp(kpzvp%ZeHN
z+cV;J=mq{;_OFnynD1|bx=MTk`j%6m>ngl!v}b8>DoPgR@sJLbUSsn#r2gT+9l>0K
z5-*7=uVk@#tMx<YT2PGYS`EFv)!S|7JF62hv3W1A>_78U^`-UgekVNbmPonyIgSf-
z4)$8hHTAdATJGckR{q)#tgt{a--`xkZv>RUIu1cGtjjN%ySnUNHpkb(NA+K{6(VxZ
zo{8+g%PPw~;=L+BoA1^_)6_bXuI{C5p(UF(9J$^kA0oqF8`bHoQKHTiBAk2gZJW)o
zJ}GZsM_1?2URJC;_6@Kifwc#8bpuwYt@zN#{`P{$qX+D)%cJN8oGqF|>^CwIV$X)J
z+B7ka3>T)G4jGIGPV2Rv^(L?un_0(eD0q~Q$o%|iaMJX4YTV!!p$XxeBbSmXMp-B~
zyzb=|E65**usgb9bq_K63?mI#V}KQUEzD;MAulj%OMS7NzblH@N-5^CqLcQTrwU^1
z<=OKZLgUHCkEmMg;sS$TkC3$c$X+KETjC<GjT)VsQGC35{h5hz8ij2`^6eLQ4QRX8
zLiJ>;FFZ&0Lh#sS#lmBW;LHl6E1DH@l|USH{RKhMbwFL0Td`>BVn*}*yPjbt;y7hI
z<nA|2HP+Lt-G7u-dM>i}zkc-XISu6_$yU<!0Dji>CD-`%Ro`I4?9C950GYTvf)^Ud
z4&sCfAyM6&00{V(6)TT@3#`yvVZ7iPutL4UUx%RReAPTZLujWfm~X)BoB1rnbpUtr
zSinq}SF0Ww^AqN~#Br8*eJS3m%$>1fLuLD=P4S*+G5x$xi}DX(h!|{_AHt~`W*eY>
zFIK#G8TX^Vh9+3!uC7>lYy+^OYgQHLng*=Oz}mLM>f}{s<gKxfp8GNTbcpID{mFF+
zan!biS@W4U>BRBtp_cVlr(a30<-gmH^XrM2JzhcSOXEUe>gnQVCx{(Q618LcT(dY>
zG8jKM-tj(I563yXtXO&M2(ShNYbmhW0qa>{eZRwcs#eraM{>1au;jTYDki~Va%$qu
zRZs3?1QToA3e642>Pt<tT#&~b^%qoEb06WX<ZyZ!e+n=SJKZ=exJb`R=x7!k{Jtu-
zo#*W&N4W67URG>8mH=1-ffcGB`Z@t-ZGC~Dn0%F^U2}^Nh-gOXoo<)0HuPpEZ51>v
zjS2`ODS@@s+ac9R@!mc9=SHC?+1=L6n5dK9w-yiiwBYAR&sQpaVv8UOvVC_>k(T=M
zTokNf^kZjjVdJqxWWeeJtf9cF4y<S|_zFSMtZK9=ghhk651Y%=^V{`*9292@$XX4J
zvnNhgPFl{<(o*~>dPphs0-`WgG@+%xkk_a`-1bw!%=L$lU%vawv7wRl!Tu5p>Bx$I
z`yoNKznp01|F8G2Q1j7yE8G&$6n>(w6{35t-}4`ii*U<4sGnVtd+fGGsAvRdsraC&
zmMO6c!%txkZZCrF*PLEuzsS6z54sm?u<@SeI#KsZ&bwl3VDy$s=A+cYz(xTHw*o>r
zYq`aKX|#~L^N*Fkf>&Wa{|LT?`A!@$7@vlq=z66_l6AoN7_IhTl8I6d&Oj<_XbpM6
zt0X18m2`VN{#o^(>*MM2b+4V5!s-`jWbfHnD?gK_b;)Rzd8n7)l6%Sd-2Jf&jXy3-
z-yP&{;QGsnmACc+7awTa1)Q+o#(Y1jFfgL?RD-1AjM3$=JPOX$_(`P3bVK@4zRi=v
zbsTt_Q&vsqpP9|yId&9@%5wUpa5ppe?ks74Q<m*b!YP4$DX9#xU6u=Zi~Zj{UKn=Y
z7`Ls2q79FoHHDS8E&{6(ur>lK{3H~Cu`C2dv#Ogu9dnzTd#@BYdzQ_3<!~j3bHy{$
zZyFVT?7UTqRADFWUw?}J=y%KbiA1ILfp_aFr;{$GSzV!0t2D5&y%=m|I^uNPXc<Z8
zFZAHq^8g6AhZPHNCHxF*FgF+vfVn5U2xIt(nSWJ{yo4Beg&faKzez^x58Z5s8hC~Y
z3hO#2r-`gew<`z*PPrcrrky*_);%2&6L-)jo*$24GQ(V-puLAyS0uSsJ+<v)mwaMC
zjLC~VoLKp4A8?%j&R3wRK5$yY7`*kjueOLJDA`;hFkoZ9ToTIIcXN`L^DCb!NpxYg
zZG2cGZ*XpLj`Ixdo$=C{<7E90CSnA<KSz$^Go(25jvZD>u+fT{xIE#;r&4OhsjT{$
zUT7~XHrHM!flUxtYk^hz4=ap+Yr4;8Du#r%acl1B<nz#d>ptsWd4}W`@eWNOh&it$
z^u1v{6s2q4I<~(dM=NvH^J30HTs(%2Mi1jtYr^{@?k1{))TsHkk&d2Y91PpTiIt~<
ze`5C7uy-L=1x{E`_|RKn@>G)}f-vgco#p12=NG3H3M;;|;koshy3KOk?2RQU6r^6M
zElmFwDxGqV`s9E+pVhoh_!1S<MV3n2Qw2ZwJxn-U&p&;M*HfO=Nc_IuT8JD(-Ca{y
zoeNI`8?1ZG{6jSK_oVNhU#a7?4@TubbZXA>pBXXxP<NfZ;-Hslz%@4=HfxEcf*-gK
zU-wn>*bkcCel%@FLF#&gllR)Vocm!W9-hx-{hu$I<*o~HI}iA2W=ETZ1)x#?vSM>C
zjIOEBpWwndCg%IRffpZq_4oX@wP(g0trp|l-kxMPeb{aIQ-fy=-}$x%8-mp3iLp`~
z=QXDnNe7TBZ{l&5Gv4-wernvukT{KdEpatQ_oSao!~Fc0@u94iMb&J!J)CIff9<I-
z@Bb}(Rw8m>y#t!w2Th@#5^#d1yL(m?yN6qjSbTE{;p-CPHu6?}4(a8PeXqyIJi^^L
zjodl>zLIHv5DU|HD}AMP*|pn+dO30NIsf-CDeK$%YIy3;BiR_{MDFL<4#}?w3&x<8
z+^PRqc`DR@ya&LV0<3DlN)N1XP4joHn&Bnvifv?VzQd4|EcXUSK;Yhq_eSdH{p>xE
z-2{GX+Y=Q9#zgGLPT{T5hHVpVYHNJeX*}TQMyWe|!p&dd-B-(1LXRQeZaIgN%hHaU
z0ee`n^3^3^4FOi@tufz&!vm}ve^}M1y6)TQ*vzJCeY^O2`{Y5+sWrm_pPb_&&g6?;
zj+gw-6n87zy8Yl1#U;IhXquILEMR`x-bFXXij!@)Gkp27^TErm5~LFuI)}J@=KgAm
zm8U{Y#Y+axNZ^F)R|HT`@jw6JR4r4Ned^om8|12M+1$wLPLI;@d%wB=wO!=HHFesC
zPH`5FgB_#2287J#aHGS0a^!sMHvOM`rhEUZWz9tAMRw_EK0|n}xm(o}PMy7aV&kXZ
z|KurJPoqEm^fm;gf~SV0!!hKA+^5w1Cg#&f=Z-;ysluGYdw=r4OBAD@?^6G&8#^pt
zX5Ew<cquvfY(=mZ^PPff$JA}hoYGFo5f6HMt<L#p7Z+4yGYD?UM+ZQ}-8{w0Q(^yr
zxxYOLSkY$yB=Cg2%kMmOINTQXWcs#Im5%2pBl2|(j_zr7)7q=et~@v8h|*NlLn=af
z$nwkvS@{X~u?kc#w(WEB;#GBrJmo)4DocCKDn5Xr)|cv%ebz^JV{q<0tk`%etRK8Y
zn19$?6WanOKg|Ef9oEawFWV_Ujx*uAWNC-giTH+i7ChqQeZpS#bPwKIOVmfvC%33&
zZ_K9azxg#ULOl6qz<f*mZoFv#bK+pq)DTbGB18MFvEh;!BJ$Qgb6Vm(tXTQ#3b4Z3
z#|;2h6<~$K8q660ItSDl$RBjg)vPNeypWjIK{D*;vOf`NTHSInwxw+)QKYr4tjS^`
zQ*<G1Vwf~so<E<RW9bF!GvALgr!2Y=s@!rxOV8hl-`flpiTpJL=S}FxPOZYqSHA#j
zJg~x(7&OJ`x&=Wo=M}o5Bwa!B3<G&i)Sry$Cn5w5o8B9p?$Gr=m!BgV#qi5xAfr7{
zG}>%7@5EbkHm16tah}CCPj1wr=opUj;zbTw#VDL=9S>KQx?<3z`t%>H*!b!!u<q0V
zV_@Y4U3WY{bv-HQsW<*d-qE+|uS6Z))wUx2#=nN&b|k;xdAQN{rwdQ8K%JEy#S!}p
z>GSGQcHt|VPfhPG4Xfro_egzDW#Lp^T^f`2fcnB2y8U}JMKk|vU(E(C1ZY|gnnDfu
zKmC18=)I~ku)^%1=b!Kvj1RzDM$f-$M#N>U^;WjTqj5Ye%lmZ2{aTfXg&&!Z_OHiy
zACkW;dDW0fdxqDOep>ioXKzN366GsaujB_6uY_o(KX*OrFi9<qDr9U;er#DU4dHh3
z6e~|11THS%)CW$pKmDH)1VyhYE$7h@38&_Tg92YOgOHm`l$LcR87lL#<1NI9xwfVK
z>OA{va10qEH+@ub%#FUdY~1wys#+N!Q{nqdY*^=7CjN3?d>*@<)wzQ2a!c=q?Dw!@
z<Ecakft3bWU4T^+SRo)jx&uaQs?mT`W~9=5pHP?e5ONg}f!L3qUg#@gm*%Gvl@Qh&
zdXe{4cF?A}<)TpY(=X+r{5fA57+(-F<TyVaTYqTc`2CrQ#sR~t&pM7zX$#T!z;5`j
zuGsi0=4=(#JLdCProal?;-CD(s+P>2Bzf2#HT-yx>kg{TgkSd)MWK;+yi;}X`8Pkt
zQzYC`{mA4C#^)M|y^LLs5{dTmWYfClJp18c`@UY`$Fn@XapIrvl;vfQ=J8hS;l#pM
zi6AHOz5yo^G=+02Lccwnn%nHnN_MMb6o~ODOO(Ztt@InL4UBHC2G31wgvmT;zXb8>
zh9@&e`IhG0rk9V$G5&EZ=WM7EL(!+W?(t!x3mb1O6l;_zYSb12JcLPx^!8|qjjs|h
z1M7ESMf<8Eu$lm?<R48{PrqODIKg!Ere*u2%Y~S3{HsBUQ+GT_l-M#+`8qdc&ar;o
zpo$@-uVn7LSK8knfl_=G#t`sOu$j;Ho$q(DH}@94j&imiHOh<=gQLaWHHFRjFA*Vd
zLjA{l&KdSA7|s(roITeEJY3R{scEQ5ycY`Hc8$~Y6Y>QO9XZyWA{VqJ2>Dc>n(CV{
zzFOZ7cPBb`-Tdb%yi-waN^L=^ankt&M_L^wu9T>;L`|Jz7~R8(m7k*b(+~smZw?s%
z=S>)+H$UjLuW{=Vj-95aV;tvacTXhGk*&!v*?aV3npbX<J06QZ?w(ylbZkTN8qTtq
zHy6>H)f1N-rRC*ash12%C{{_-6i|1JCyxk3g!=@sxjP(`41kEcnquXrkfZo`z?uVd
zuMMmrFg^o8F|4MZQbw+Gn)IAMHaC%3p&STK_gea)Rx!Lt>$I>Z?e?{JR4j?tZI6_&
z-nmWU6T)lFZd<Iyv7V})CrM!S&6qLGtgI)(`H*g+XV&KZf3jlZr;wxgB)}R5tdLBY
zdztBg732L;6Qp*2EF<(B+a;T!N{f~VC%x13KUnQJ^3PoteeTz{FTjc|i>ilMUGv@u
zGK_^tQK{>yh)34ac6GMf{ZgJ^v85J~@fz=jl2<b`4saFj(G_08<6nLTryX)N0azab
zs}p!X%>Vz@e;<cP09a#y^**qIM-jk*68;LjWpoZ``tZ_dB`gU@JRPFtQ+@q`hwsdK
zm8qbK7TfpXH$6zR2YppjwdWrzAE@_QwXBkOL~NR@5m<ib`en^Q#-GDD&{MVCWjNTR
z<LhweRL8^K00{ZN@l_&N?|3l(m}?bMz{(0wi1|AQmR~0ve*JkZJHkpb*)U>XA!ive
zsWF{ri{>Kx*K6kPU!GPVzxY4!ef>s-bic)`%U`#RQvJ*7-E7K|1ld30MGjiIH!ysQ
zJ6ZOFiG=k16?E7gR;+xL;}2^D=&B5y&|_itsu*3X+RhfCq7*BNW_a6G4yK?4veK$d
z74r@&=sc@l7|l}V$WKu>_N3X*gM$bl!U^$OcyL2;;BHO#s}5Ar%LmzqMToy{1uoxF
zeL1DVgSaidmlcb9n20EW^)~2=0#@*T%(dVV2#V<ehU1QC6`$~zK(YoL*Z(LYz<g&o
zB{JT0NVTq_HV417_Ti(KI)ha<?-U;n7rZL6z7sM|=g{z3;!=k&#Wqo-@Wty7Sc3d?
zZoRGfN@VMSM*XWRHr{#+SiOKX3Ruqps~G5tF4|~Tb<di7t!dvwIhtp$9h&SZJ{V{0
zFf67b^GKI*BU1izyx5hDHV2#v-4(=(M_l+Sg>7%Wdlgf*q7FE@7%q{yQWc~gQCYVq
zr2l+RP$4uQjk?Q<mA8rmE7X6?^(rl3r3cn|2#U!8aREeBYE)ZxRK5Lc!;lk=JW~Fr
zOyhEt_%`yYTf?NE6g*+HP|JSDr*%dsHWVksW{k>>F^#om2)7{g(aQ%p1|r05Tdw-@
z>47~}7s(l<_p)NKXC;ya)?2`e?pZW|b?=^4BWZ9!DAYlm3%PYf-e~>iQ-1LUQdG%p
zL?%bCq1^GIlkxXi`_8K}sauT_^`6)Y`;pxxINun)Wtie*PS}1RT)FTxQnhl?;gYen
zsEVm2pzi95jmHwfo5ipKFXlQR=#2jZf@1PjO4reC=BZS6!bC0p`9%s8!sumQq(yIN
z_4jYZhnOzEiYPblHC)EIFH0)d{g|>Lxq19p->u8T7A+^e6wlT?cQIN_`ShAF^k%_O
zu;XJi>R(oDuCJnV;2P+916ZNI#_R*XLQqT&bkhw9Rdq;l1%!3u%YAwwz`-WIWqM9%
z?Rc)E3ti>DTh|LT`m_(3B7QX*27C`7@z!d*eEu+AO=Vy~So8W#9Uie^H*Q}`k%#Uk
zr&K#8{=tfc#}YwJ#q$GJ=mXII1|7jk(DmmIYwk#y+jm!Gu8PaJ_*zr;2bq1@u9-?P
z&oNrwEp#BPvo<IgQC6$_D0|rUTwE|-h^2SOnOl+;HMc}N;IzF)L`7WiR&(&Bs*~<{
zLyjWfr6KH2eMNKs|Ln1ZIiPDC=!!mLj|5$DLDw#L%jkU7{3XL><d&1sFCev69xBVc
zg_!eD_50RqTJJ=5uF>I{d<nwx72&|(hH7eI!n`{l!<-Erp{^m*c>2eu=?}A_B(|GM
zoXhLEU&MTV?}<y{2qAZM#pZjO(t))GSRqG=ftP>}SnGFKZ;{X$IX5MmunQJ_LDhtm
zPIBg(J51l=*;mq@OKEn5s_x39LDK-rw4SKo*+Ir{1zp2BtFafT7vkp1i=ucrE^fXK
ztFzbnNlG52K6r}3We+Ph9*e#w6t4=t?|>0ltAW)3SZjA!)1Mf{TZ%ZCbGY#Zhx7!S
zO|JIy`VD3+lM5*E*{Dv2yOLNQ9SAq7vIt_DQ+iqcsNiJ$y2_cXrEX3wp3e)pr6ICi
z4)!63d2jG@p?YZ{^scVh>|F_Ae}y*(tgtyC2A?H>b%yyKM05_QNpUG5?cXO+8$Hn!
zta@g`F5AiZglyc#xq?W5ZBtjOsF>JCFgS&dfHIpn#?52Gu76+E<Vvmd_h7HMTK{)k
zxoK}`M{yU_KOg)&WVeSC3xCCY|8Xa9LSsez065PBXFo8a>y?IF?=sS$U4e=9<q^V=
z@xH$F6RqWbEA0o5M8EdX)#D%S)}$qQsC0sxc0*vWyD*eWGD+PpjAPlF(}1x2%yzc+
zWeI~8=JJ@lhL)ym^oiqNR&4wg!-~$;U%(3M3G;miySW-BL+2}V4wp%Abt5a(hvGV-
z^^SqV`!ErC;pKFVqrH-WGLsqjxOiMm`@XF4G2=PKdVDR4OkrBO0yk>DQe`vw`uZ9(
zw}R2oK#0PJuJyvbtXTLfVF}DXu;C5^EBwUtgK*A(Ljz9@{b$a`M>Qz*wYne4JJz{U
z=(*i8AoY<g<ncSYYC>)|9Z9Y}M5x?@FF9ACAmF2|EY0A&zoeW+hVFS%a!JFNGGP*;
zR<qi;4;xh8UygD?#NAxQ#$O4GfVC8I01sI6fEB%ec>zJuxvH*s#o8_CHJ-h$LTw@M
zWHjP-yTe*9iGEYtCu+T;c^Aw)^%I>mi)M!eCIg+MRmAg_YVbDs>-PKE%zE=5)_i&Q
ziIT^%V~7;bpz*QAW;E(wRxJD#!wPF3_dT$J#so+A)Bp{Y(qqPieb>km^^A(!cs014
zZDya~k9RX7-x5$sbbB~0It}$3CzCKg-w2jvX5slz8PsR?#uP~(JiMev&y-|}co#Y2
z%drtqV_BMl?vnSgV&SiZX`t(CVBG@NI$&i6R`~h5Z#jE8)y^b9F4aW*aOJuK1&`8#
zSK$kV6z2HKfzrohN(UDEowiRKgb6UL86u@wxQnN*U47?insFw6?9k2g;b-%X_=RU0
zHT8HE2Uk>_KMYZKbw$7a*Z!IfYzSa22Ua84Z$q*E`}y(zWB-OawCeZ*z#0#%p1|t*
zoBemqI#nouG`@O;nxoe5QHY(C7ot�JWOaq?=mn#NB@J9Zo&v2<rTPpE_j3i|zMQ
zzQh#vovEp-?QI*M8fWAa@B0|aH=UAL`tV%gS2bM#TFL*Fzh(iO7qFtwu+4xq6vpVx
zLa!~&N$){Bhi$w#mNz1Op>J%?PbLSj`E?NJ^?gC|*!eyZpF3=o#B{9nOyj`ZUB1|=
z#y(dm>Z6MX=Obo&O&&L<ot!`SbEc)}b8cDX?GmAHuvOk=#mZm7Tk%0l+-P960oDi@
zGw-k-BMq>-9-H37-#c;(bzJm(dc<)lzg8{1H&+&rL-CpEuV}60m#7l?%bRsSw9a=t
zQ7%y`JFIeWg66T|jaM`Wv_7hNTPocN^mpNAS7#FM+slf@9+nVlDZV$b-UZf6z?uf*
zV-OVMTPgEj?2P=jba*d4pbW8gTbeY#FzYZ^j3mnoUgT9%Ckv>jY5ZPD73h>^P;`WI
z_0hN@57F$wC98O?eK)t46{u8J4!>QZh}on&?cAvbLHFp2jmJWr!Vd&im>|q&Frbgw
z)m3dP7O6~l@0wEImQuzcQx$}424`WlnsnVo^GUVm-IND}BVT6Muk*ha6u?nPJz;NU
z6=*&8|GGQxc&h*R|C4bjgd!_t%bq1Gdv7w1y)rU0LS@T}D7(lkTahhW8IkNrL}kyU
zG=JAQCEkC&&+nhl`|Eb=)_MK&yxq^^T+i!qUDs=^CNrd~A@i0J%9Ex}^hAEdJ82xF
zdP@7`9w$cr3Z{pCe!z9$yat+v!7(d*d33BQbd~#<?=VoZoIRQlU61X{7<rvKYFtzA
z`Q+lNnaD_5CY5z655=z^cltvgX?{>UCftk9Fh)eV`23xo!p)~8l!J2R^C?E3X+8VS
zsx`}9-KQxg^I`NoH`btO39zDafDf>e?Xn&d>a<q~IIY8*>GL?EheB;e2w&cB<F)Rs
z0z%8Mn#=bLNv$l!DSC|%eP_DWQs0{+>h)D`o0BXSy&poZZRc`rC^Xu1HRj4}gp)28
zQ9*0`cVA)TuT{Y20Ibkn62SlXu)lDbc3EXAI_!<d#c#1Donej`Vt+T2^*Y{fa)Kf%
za|tn!-`aeV%DR{5P+sS?$1TT}unQW;S@ZRzR7mCpY)1zl8ZgwoRTfnw6jhh!Tz+RO
z`xIJ_Jys0;6&IQcbgW(hR)||%s$JIlC3&CoWL^5K$ST|jy(6SEn>SSa1}Iabk6(QI
zqvqv81eKL97e(qs6!I~mPoSi<b0d~lag}8D%=^*g8aABa6}1r?f{P!%`A|Mk{V659
zpA{2-g`A2DK85WDy7~hv>@W2FZRi}JKJ?Sx6-lDcn@4y)qK->xMqPHrZ_&Oz)<NEd
zBc|m%RvTr~*U~2GV^l@#n*6m&+y^-swqn0EpP`7AIdlIM3uz01k>1CJtj>1_JktXK
z^{=j|&;QzApF=J1ADI!OYb@vm(C^BE9DoNs3@+R|@Oy^A9Q?k&nV7?{q_)#Wt9|nl
z=d{cTvjdg_zGv0zRImg6)aswmKiL@<D7wAEon#mkI}<7sMfbw(<TkY|mxr9*hK(-P
z-RT>B%6Hd9aUE;}VcGvLe+5l(Mu7Dj#K0`DG5~8CEQ-omY6m%s?XO_VkHjqYDMyJD
z#J*N4Xw@F!bnwiX3=ItLj+RR?*wl?kh&)w6?SpeKI6CoD$yryvaNmR5dF!UU2EoX?
zE=SB(ChiTfWRLA}V&ttEz?lb}sNN+7IKi`UA>aPy%uz8-aUG*dGcA2#Z<%Vv5_zw>
z)aUAgd<6aUf>NSS=%WD!WUX{b><l}xns*u&&P-@f+S*a`x)B5OH+(duqE$jd5}y@B
zYH_UHVK40JiBaEL44m1(2?mXxU0eoE=nsE;Kz@#HzAbeGBWth5=g4(){_KTBy|kx9
z#*_Jq4(?}9megrme+)UnCKbfZvdo#}*5vCv<Psh%IqUoO1yc)>PODdF(QaJuv#za0
zopM^rcUbl>D<-}g39OLtuzvz;GH7ZJtPy`TjVkj!#Qx)Ut+jn)QN|3Ldc8quX6Qy_
z3jWZOc+sQJADO-}>q{)Wr(`}*k^JivK^S>X_S)2xUXruK@EhN<7mpq=*yZVSE)X%D
z=iB4N$X8Kwiw@|op!XvHZP9ZJ)P6$y>MiL{5uS1IY$gub>k(?tUS=Pzj`y2sV?kQ0
z@~XwKtv-HfB`LY!nY<Rl{L#_jTE~@p_b=0UWvn=E_0yOwmgLmR$|)|~BoMsD#}&<)
zf?9G{Q;d8S6{r2cN(8KLfK>uC1&{q5r#~ou+8e6Ar8N5dlV7daj7RQiSZRaATao;X
z<*qEv-0zPJOwKeKq=!9sBtOFEefC2CGfF~NFVpi{S6F_1*Zg??;EF?ULQO$ji#4Oj
z9w$bg3cokHzg_~)Uf?_loVz`%X1!>my_TEoIMd7e_UQXgK{FYH*p6#NmyAOB3%>IP
z>m7b&?a9Z>+Msbm^!8e2bmA3-*jDMNj0NtjGdG*r=;b#%4hK3Y3zH8$-x(|9;Mu1s
zM(17Af%O^0DFi$L?0Iy54I2Oc{mF4>eeC>b(m33oc|}jQKFPkH)M_G|d-|rNW?z}B
z&@;L?hY6E5aTkiXI+@Plke5gK+Zs#8jUW34PkBri=d}%zSP<p;p3!lQU#z~Y9Jikp
z6HkR4g#)p9;5e|>0ILS*x|^p|L$011IPc@$C`Vz1zd*Y|LEws<3rRi_hgYgb^Im1S
zk5TU8`KMKa-^A@Vmgp~8WK$ozbfuDjVsauqv5r0OCQBCnhu$iFRdcIXjQcp@BmDi#
zbFTlTeiil)>bVQp&<moUB?z@OE-idzRGp%Ag3#H=re&0p_mawN#3XAv<GZsJe(PBZ
zI_WJX)7|rAQmyrN2a7AT505f>(hvpf)SqL#+Uk;5Z$Bl(V4YqhC;XD-;At0D`i4OZ
zF_O|~sl9!L(X~wAow$(m(eKQL`Whe375{y{fksScw7u_x{9AC*a%j}k!zZ)3_(T1c
zCUIWa;qpwXi9$_9=Kr<v(S(!dMJ2(YvghnHEPS>kUvmt#6?Dv3#IQMS2$g>oP?G6M
z&|QBAUw)4jBTuaYHg8~s*(!ktu!g`f3oMGxtFHoNebi}iM_3hpT1PZiG|aHq)B4RC
zQU`ThbfB5f-RY;YVx*^i7@W}|t)Mn@t!%L&BbMzo$=mJu(aU4U@buPH4dpt+f}>rS
zOVv-vdF^M#q+dnpY6Pr(z^VtVp>RyU%X;KPoQTc2&JiX;zUt_eJ28xwLrQ+{2!_LO
z^ABXD$iJB;wc3$6^!V%Ag9}Hv&Lgi^f4zH6?lQ^3qbxlS#D`3g6=mhtuGh8WrI&&#
z#qFWg{ObW2J@*hbvycMThrsF!tf*RLS68*DBWDdQR`kh{awO|(b!J?KKjun5_*k;E
zpP(_EH_B4t4|saB{f7yA#8Rh%F`3AQ<%qYtE2ZS;a|iJJjh+b~8d8(|`ZG_Y;p2sN
z@an%_Wq2+%H4TFtAPbz3pU}U*J#do$(Nsf$Im(DItAWiF+y3F{yjZScd$uyXIA{KT
zTD<AORxT6H&cJ&EhVE7P>A_8ov+UEPBkP$X#*a!q$>hJL37mRDcPZwi!y@O6gA=!o
zxeCFmd#squgJG_U3-O7*?#c}`h1zAWui{RB%I@m0LC!92H5{S29F~1L^Os|vW8x9@
zh%Q;S=Dvrd*0NM5`q_lloT^m#4z3|rS8;!QC7G$oF~}qpBQ$ugqKqK$t`<AdpA@I>
zFtU#oBX0#gaiM0x_5odEfi(b*+5YIN8JF5*FWV@t!$7uC5g~DtpK%+X+;41zBB{*2
z_4end=5C}`;`FpH`f@ug6_l;K!sSLX3Ry*ye$l24lHPF)Ht$kr<&3z|<Ty+tai=@w
z^nO-M`c<f>(7iR(t*G-X=s5uOE~`TYmydm4sy-{F+WBa6^N5)niS&MR%43R=&uq$T
zE}K!2S__I$lm(<3p4WNKai}^bquZ^4GHLE{w%(z~21oj4RZ`dp{Bv4IC$rN+BmcpQ
ziMK+niVHO=b|bLb0xR^{xHNxQHK;!;7|}0ZyK>%1B>w7*K2xYBWf@t&vqb0#mZ4R$
zAp#!x)%3gZ8v>OckqJ~2Y}js*cTMH0U3j}gJ?bm9DGPbS^m{3>N81@(trhz%VAa3z
ziu(Mo{S}hPf2LpU2G%Cf73wThzCz!BjVgpt{;KImJsQ_bK3kLLLNz4y>zxLRz~>jp
zt5QMKGJ`UwefTnumUUEE(in2or7lWW<`a4(Jf3oB{M-cR=woA3?}gteKOgUZbsSHC
zkSd)#Wgrli+|5}S`RhEeJ^)r!Ex-w^kiT#tZhw0gZo&M+p2TgX6t|$es`&zx%nO|p
z$c(4c>8_pAc{NT(*HbE(UvOlSG2!Jgtxccc#6wC7+%59gS$&09z$rc~XtQ+kPXEj`
z6^W`5Lq691tQh$#<Sg|2@m7Hq>Hs`(U<IH3%^EGr<+d3T+F54X^IADDo05T4O4Hn+
zRb-_u)jh3)J&VC7nZcvr;U=!i)K{Z;_W3-g!ZOl}?$&M-&LtNWjlSA2GMMyS-&Qwg
zsUspE{s${2{t7h~&Lhwj#To{z#=xrnhgE}fi+DXIzkRH*-IiHl#=@Irqmj7LkI3SB
znqj){ddVq!qg$JzMn?=eeA$fL$hb-jVlKw=Rja)Zdct!@Zwe2;#E(qWFm~iIXJ63T
z8f5Ri7{JJ1q5sE$9DuzJtdMKb&yu{b%Q}FMV?<p2lid1=NU;Y}hmDQW4o+1YZUz1n
zs+zhOBeHZbZ-@3dMc%mz!{^NmoQV71jSu!2XDOAw;m<7A*v2tE%>Vig-l0psZUyx8
zLzlJ3ic#-c2dpq>!=3?F8(@_LT`m8xs{Op%eoL;~MocT!Dz=^5{ah=H<1^Nq(Tc;F
zT}Al?TBhC$X7*`}LQ1a6DnnRaTqY$Jh6~$5@(Tl_u}faz2A8D;f@fW$qj|6ntUw;$
z-B%d(t@nVl0yyV@vlciNfm7xWr#jJm+bz$f@6k#jjGDKuyZc;<{QB~}e4Ijme&mIs
zmOfKoy#S`^){~)%TE8A+sYc#8#G7%uq|?~c#e&U?-FGLRBY4xiP`<-?M{*X24khVd
zR*d{L3s@_G^)s+S9e}3`tmpo)YA~M{?+>%{>r6Av$u*VS;W<g3K3Usa?{_QP=!;pp
zw+#`MqC0c_a7MY<J9$zvZj;Uv&5J(O+C7#{L}q$xbWii$4vYG-Q?^R}{K~}sk66XH
zcSZgFi@*x|3!SeF{;;YKCpX-ZKKbjyuap_{(FS+Y%cuMXKgh*KXueZoE%Xdb_m!M6
z>AOwl7(w}98B0~m_X<&d9(T{1*DfjHm3-lNWr#N}-YM7l2zRm#L`+b!?&*q&zd~L`
z>$(N3kf-pZK-XPXjU6Jpho1{^3?mjKi!~LsG`P>6@2!2gRnJNOP?E^qsp0qq*-)ff
zIH{1hyek$q*BPN+sx9$qR%hp$?1`ZylaU5mwkLQ|Hhs6$Gv$*k_Oqfs|7(B!FZHjG
zqfogNeJwuJ*!baa4D*5C{p<Nyd>_S%_#|$Rfr<!np#rPnS1<h5=<&4EWw>AW@Qz@$
zR-Gt(LfdkoJj1~O!KYJON#(22Zd)&XbBHmlqX^&YrBvGKOSx1Hf{bF!rJPcGF@TZ3
zJ_a_(rRZyh-GMa@j!|_0>YqlDr;EK>BkL_T2{F$|MM|6O#0Ob^D^_wE??ZVlRXbUw
z+DM0yB{ekPrFav`Zbqsc#bs2RAz5YmvOzpUf~$(~I(4(3L;m6SlHMt8d{}Ug6(fIr
z2yBp3(a+<t0#;Q2dIA=ugnycymx+Cpw}wYpX(&S?Y`=zCy}V%Jw;7*!gM^^P9&wo=
zlgbL8fHrxUP2<)a_aemzx3^@|=?Ri~Z8p8bHi}GqD~jLV5rm1#x~%i5lMH{_&x)bH
zqVI9<guMk@72g|JkHVrTUA2t*TYXH2E{vSa`IQ@mxbTYch90Zm2DM@W68q7y7j_ka
z!L{@|QxrGL%S|6?&K3&HHohC~Eyb&MXgn$uSyHHc`DKyi0jpm7Oyt?(<ujXs`*g)*
zE)01J{VsA;PXKG-d%`gpaH4covk0T~@x98@cYgj{xFS53L&~15^wW2vvkz%$_@<fc
z3HYOOvo~|DBAPGKH~E`ZxD88&S*&3zyB;s4Tsq3}h|-35Rr-f#OHhqvJk)Z3Jqy$R
zRUdx;WZ(on(fbO;i2A2-dA`RdpgWx0l-xCC{OAa;PqAGYcCjtR80!qtBL^uW9xSB?
zZD&sO4dQw@xSOelZ$$9q23DVFqimu(AkxhirP6v;5En`L!u0j)&Mt_weXN-Duc-M2
z#3*{!0`nGp)C`O553A-c%W)s=EhkD=E(76cGdFcc<lXmvGyX}53-diI`N8g_R8|Q>
zWN8|VD&x^H4#Z{~dD#SU?_=Lpu^BwbS7hY*q*&ZW$m6NbV)x4BV}ipzR*d`=>Q+?V
zN530O4OkuF`0yWA^{(6Y<v4FP3wKmz46h!1c~MDjag%~@Xrz@dSyRCpE6kVi=(#|%
zTs)i|X)>&ykiJ8A@ID<M({0WF1n<YN9f?ag9@}YH*H_b8o!QC^i|*-)k-ws1RT)@w
zAqG&n^%5MTmPN;aSkJuDwXDJceUrOw>Js#)s_4SYh`-w!g&w14pGnRn;0Z<EeH;3T
zQCKi}km{we<$A=AdHZTBW&9kWYr@p`kYeKm${yJT{HH5?XXupv!HS8$7W~l_#p(~N
z*Ws87Rz>H4Aoflr&!kO##$VC6QMqMrXBu;!_^nOgYn!(6P7*csHn&zNik6xdNJ@TC
z<}w@5X%D9U&V6~dRz5OuXNIC(@T<<NS5_+i{x$1gIk&ILNbTv0;{MnE`d{i_C!lWW
zfEY*wU0tAVf!`n34Zbo;SB<$FK5kRz?AU9+MeFHVdtKo8WX)}~wB8&v7&Ys^_F{pv
zgthFQjumr)Y8lt#7LpL!%XAyOlk18*$1BvX&whD18df7bU64ATKh4EpBLyq|%~=@v
z>o%~02cVyI06u`eW*j^g&8lg9p4#oc`+bpkX^VH#+j^{=KP4j#W32g7+lLB91iOPm
zYpHTh$-L%Mx2$swbvdRlBxm69z099qEg7L6{5>|xrpNORd1Mt*j7p5lM}9vmM*a%D
zAdVNXqI^IQSam>Gs9k>ht7C$n*2>RwHM#E>a!uKOd|;6fRX4ehN6VmPb!Oej{7QDs
zZ0#9#PTd({#O<f~4_B92Oj*ue2?$f6U|Y#<I7(*n5lez@5I@=KWwmS9(Lh*s*9$N@
z@45!8sQKU%V1+!3W)1$ss=CAcgtW%9wm{*9<-?Y6ZKX#?9Xsf?Rl?_4j(Dkh^sn;R
zQ@kJ5u8PXFj7f~{;(z|Beq-XtMXv)a?+zx$4Q?dxzC3{aw7cQ-xRus^P7J&ifd`zI
zffM!@`ne3K-sSc#XBBIQj`NeS>H8<M)2JEZHR}Y%+p>ltx*p8wd%I?5t#TCVL>#g@
z`fIConJYO^GFkq@;zZ%4PqZhRIY$WZNWYCFSC7&w_~Mzc0oM9gPmKE3FTfe`hZ8a%
zdLH2Ohg01^vR%ja#B`LQq?P+OJZ%y>8)OG<d@tUGQ<`qCG$OeNIXqYW-EZkMWHH>B
zbDn0js;hRH_$Ge3VC$gm^M0DuyAp5xX@o+LE`;1KXop4jv0~z_>%e*kSW)|54p<dI
zQ(#BOsrn}YVYkNGH_BQ)VU~Rf{KV4-7TXz)97I^##cr?(MRJc{)+^4wgq8L<2h08&
z?-vCtHJkqWk9V%-H(4Be)6TeP6geMbUcOnFHRy<7{|75Z-uez$n}HQIYda6DFvCJ!
z-}rk^Woj|<HJVQ5KWGuv)19D=;vM)AZ7MG7n~z1$m428)f6M2!Sr5DEBZLT(_0@&~
zjt=72E4H_<wR9#U92mb3#h7%*|6q6_ka$n&SpF;aKv;HnPhsS(sQv#CSmS^dvI%-F
zxV!%~$S3*u>LL&|sr@fsBe{!R*lcprZy0V_x8LP#5l`iKDJkr8Z%t1P&rOTu;Y?46
zo6m+uv$^_?5~a!U+~_qE2G3%G&*|A`=}2GM2j+|b>K-d5-Z~Ae34db12w0^-SJ-F2
z>wu|93%=NK%9_TD0eVth^(4&k>MM~nI9WaRqFf6skwmy^)DJv7u;gtsw)^Ln2mHCR
z*d!3nUv&iZ6yu7Pfwh}tOILVdgN<&n@FHZvJyz7`f9<XLz{Uuy)u5{()B)gGdp*|w
z;I+akzzQ=V)U#ReAy@rn)f$tSu{ZHcP2ymykc|{P#BbG!JM6c5FN_S;yQ<nHG`ChJ
zA<Hb>FYEkFO!FMqDGXxTlC=t}=QuaSa{Farm*C6f@p3-~s=AkL-{cSXfRelMijlY0
z0~^Y(?gA_1Rs1M8riDdOwV>u^qp7`T<J|$a&E(w3M*>^fKT3NXw@wcak#}(l*|py4
zr?PJArjCdjFmPVsm}8#WUdlgc#2-8JY0)4$f{pC?nsSY0m*ZnRE0q8fySRh<STXWe
zR6k$`tR27#^FZ{q!l-i@6!1@r+B8DMv5`cd>8*QsWF>#4m3<7A-}-R{a%3w{dUTzi
zRPgGEU&(us4>Ph(G=nj9&S$97YPJyGcuo1rsy%uG%h~gZH$=@J#SsuF$%cX~e|5#k
zTTyi@FR-HO0Pt4yc^6bah+@@P&(`zt(`X)N<DagMe12?bCTbzhZ#|J<=t!G;B3Z?A
zxK{YCipgZ4T4tKtDV81FF@{r4KPlI09LV^aSen>XWo9%+o&~J&%5(G}9|gjSd#o7v
zD}DvA%KV9ei@*v#uy?J1tI#o@^Zhm>Y*%t;qkl21WdE$4^jkMkAf=LB%~GvTmug*n
zT33pAR4gng9Pe6SXf;U{H&Kf4viZ2k$<(H&=ah3yHdLU^y+&(<wHP}VXm{fkBaekW
zg`Tx!gRZW?8vTbA9j_Z4h4zM0r2`z)nyt|VD^D0bcux3jof}FGDX?Q!l*%EXvbKLs
zQAaweF<6)1{I)5rm(J$VLz3CV0)uA}iWE=6jZ{b!%bh8VyQXx-44^IAW5vW{%Yh9t
zAG*)70akZ7CWS@Otj;PDO6G5I^*Qt{KSz4Lh?sGev-0~yJEX^CZ}~nY_g+7#RlWG3
zPr|q5rxlhf(toZ#TIEYxk6l&?&~j`V$-BRyAel{I=RfFT%dVc!j2eyY>WY!aqFA90
zK%du!FOHtI?9PQWgmKc0*qz!rKg$@FZIxx5WVoDMMwsBma)S-~dzVQ95ogT3?C7#D
z%#|k-9MxpGzZ+eTcx+`{p7Tw{=PFJ1v*R>grx-4%uDrS)CLrE61gLwgn4Ett0#+oj
zLYxu=0xK$Kq2@wpU4L#X+3V-N9${XW%ZV&%s?BcYeBt+vXgFa(>3PE3N5>i}D>D+>
zhO4dP&9At|RZeFaA7C*182h1OZgeV6k)mfwnHncTL{?tv`X|+DK^WugvCd1u(Z9T2
z0BQmBI}+i4f*@dpT42u$h`-v37zY2);mM$U7l|hpmQ75~<F{rk|B|k?be!mR2v+MF
zYj^&Mmg0wVy+10{W8^MS-9*|WJgIdG5sjR&rr&Zav2C{KOA9-53LIq42JW!}DE$3P
zR;aK4Gv{BYfED@{>@?67hBtVSS8?ybS4PEvI<-2hTTSv}fxNw?IFA=UflAI|Uut;J
zkvGG$-b?o~j*rXIJ-C<oRwO;%&8aUUftm4Md%~5iaPOMAywt_eud_N2=VUwvL~N}F
zso&f{3Hdh$F!5N#0bqRwtgsYn=7z3cb6`=F7ig~i`T)HuMgHx&)dbT|H?3GMV$~Zy
z?<No}q5U|=X`4n@sz@^3CB+)OJb9T}H4^E>L0T7gsHq&^@%GcV(!Kr#oj&w|gP(Y)
zqtbKTTcu#h-8_Ji$9@D>sI}1d<Cp>~JLn3vELvBMHGz7qKGFy5u~Gec`V!tUygzEi
zVq$~}dPJTQHUDyGq;j$1S5~3uwX>bBB7Hl1!|fM!^IOxFQCq!guSG2P_gdCHo+Lxm
zUsP0<a$Il0s(Y-M%!g5ZK_IZg{wFX7R@8ZIIM4Mv54`w(NW|s^Z7iom=1}w(mQE|;
zGl6~+4;j=2Wj)zzxVZ1k%&3}D^?x43>y^_JdzZ7oK$m%Ir}PK&Qyv@v;q6Sly30IY
zD3Th!-K6w7FTb##6%&u00ap0!v7z2Z-H(H({6|-G4B)$mx<&JnbT+inIcX)fcXJCK
z=QYynvWPDhH~L@}x%s|yRG_7QSq=L$*`|D8q{B_68N(ij#{Q>0K~?JMc4j;=<z}AZ
z+Dv4k&TPv{&HGp}^4M8m?F3fnT?x#AmFdskLb0l!tLWv6FdbkoiXAkO7C(KRwV@*2
z$RtQf=csDA_mc^8+GkPKYLYjk%BQ~QUD7^xUuxO3rQ(|R<eZh}50TQ(1~*l#De>+r
zx!k$p{}R@NzukDn=p5`fU`+y6@B~y%h<*-cEG&x30~!t8JZ?MTWBF+Zs^xT=ZWS_a
z4~s>bJt-$7Sf%HfZOo;hl;_aE{}JYF$LwL|SBG@xKNVYL`jzyIo{N_1Hx21a&sEm_
zhiXTfnAe0m!O8bnG4fbcZ+IVAQC?sLtY<;jfL+#8N`6`e2aDL75Uz6STV>rGiNEf|
z&fRIAeSDoH`UF|P0n3$~iTG#49(LF|<E@2D<*5Ssr%k_URq81`;@mJ{@G3K0f3tSB
zlS60fc>BOVSTXU~ufPg*Av$Ln04q1J!an=0>-|$9TFY{3>=lNR?Q(--g)H>F(b3_D
z7H#Pa7NWkTOCQ$TzSr@{w8LlU`eutB2eb6gcKRz%zO>hwP*iemNX2`);jddMo^o0>
zT#^Zb9kEYWOg#1ju+{<V^*_3z;`PcOR`n<3+k8b7yfi1a8W9ssx>|=++P^ax;7%Ff
z*l5sA>*UIclMmivm$o7LUN#u$K@zIQCUMxGOnu}><oFI5yJ^Uf_z%xmnVA|SpYD4U
z=^pFvqkrkKkX!#VeJtuc0mLXaFr%)+K=-lW2furmMSp7_|4Z{pj7f>QQJNH~R+R#l
zQ7gyr`j5$R-R}LejMYYsPg5+PH5`Zh@$|*3kQ+=44snzp-)@ayV?7uv{_JdWba<LC
zN^&WtNdIskO3U3mfQiTAKLu9U|JYu@3ciK!3&#)#zrDa8|AqbeN0vIACK<DlRG!ag
z)FgTRmXe0_Ckqu7k}nNUQ&~UAIJ|gdRJJ>@;9?Q0dN?J6Zg}kbM?3~Sq{k?mf((c_
z!?m_5j}s`koEAE@j};@2ef)>@0kFCOYa$%e?y^!{!}f8!ADzS%fOjPFZtq#E&Ws$t
zkKVU?T^v|*<~F*&P+86OQ9mv^Gj%1Q)V#<n#FT@;iHT&6v}j07oQ<(+PHBRT;7Vbe
zCnKYphR&z`teAK#>?_o>y0II8RTNmE=i0l5cxZuqJ;y+BtbpSUoytO871Qf`F^0{w
zS$bYJ$F8w{JB~vYb~Dt$dU847yQ<CVs}I?d2$K9m-Wzx4FTB7p&@jz@!I&$NclI&z
z<pWwQRY3h40~lT7gwhol(a%IP16J@`+>@{<$_q4?pQZX-G8xh59{7SA{o^e)<F)1$
zzZND6!pD(IVgv>)>r_@2a%A6s95q(zS&g(kJx``VagTs<S&2bU-|l`^;V&h>?(2dN
z#N-*@s!CLP?qkKsV_}AY3v(>=ISudv{4_Y`fJIRrt91_dF5E}xHo#?IR~~JCI-Zd?
zD%EfGHD1hE`{IadmAzE!sscyp@n9$Rww2Grx>4<l>~|=85pH!=oJ}R)l$oSTalB`R
zdz|ZU$FNd-Js7yBD@Gn$32fj2XkAh98VtwBffKE(NE4}#lGRt8lc_Sy5jWrD&V<~|
z_WPE3OLt=-@KUG3fE$(704{Ch$C{zZzT%3l2j9I6-LKcg&bck>A?VH=7vE4Ow1}LE
z7uJ7fqSMB38<yS80~mR13$RH8E9i-SucHSX|IGtdK@>iE4`V2~uVfEK_+0C-GX5;?
zw~}%<l)PA$V5j<s1(g*$J#D=MqU&Y&^HJL-{;VpzI%|?;cE+Kg4<Cxjr{t%2@L}lU
zD|T2dCU<y0D<&R`%CFkM3cUgP`h1kGM|X8iU}5#KZR|?osXWOXwYp-&Xgp5qw{s}*
z79U|)vElVmK2oa%0gA$h(St>lYe$36I2z_3e|LpsS<!OHA3<4kcU0jp9zn`ySvtkr
zaTMM9`&lvZSX2x^oMPtys~^O`O*kflMN#ppVfyoxy?w4AIrpQPoX8RxJjM?QJ-?M#
z9fxT{R!?^SjF)O7pGvDgJ4jjo%*L|gp2t&C`*UouYqLoHusSmhX)<LO((cx*x(ffu
zUwLVPu;SksKz;td+{Ypmp>}~ig`EnzdO@89`3pB5zA|cWsmn{tk>*UQoR!-O&T>QY
zANCgG=&x}MeItl0e|Y%_FZc7KB+_%SdU|~pTG@9){Q|6ttsc`ZvmG6MtvUKU@zAN%
z<72TTLYd=K*AKh3qlDbk6@%wCA^3n5=EUfn1$8SP0$8E`{@uqGd+WGWVM%o+u%$ef
zdjIGa8^NZviEy=Bo<3bz>4RWaiqgDm4M+3@KUx<#e7f&oz_BKKc$&W{1HU(Ndg0z3
zveQ)Y_T%c_9VcGhE_N%E+Q*8C$0Bfm73x)N=miK+tS5jK>Tk5J8r3w_S}AO2Yf^&}
zob<cs6d2=Q+eHxff0rh-){MC5?(|-nNp~9W(~7<G9nF;(uHoWl12&;2UkHk<w~rOA
zsiY@O87SP68j(vY;KF8w-LR)CCbMD$>>-?FV9f<q;KhSl;P3O0G^L;3;Y<Fhk)Qr0
z_MO6I0(?$KcSU2~$2EEB6Ti5`bzYg5Qgs>83cK)|%SLphopO3%e5TpBqlvniIQsgN
z04446sY1gq*+(BU)197b{wTGd6@xw&ApooezzWx5pyQPuc(GwoR36a0Er#Vbz!I%g
z>Zc{7_~mgQulRU{@iTmTuRKppqRAbc?Xo(Ft3hJ@s*Jj?xeuMXD9C!rFfFL*=%<q4
z`GSn0=T%N$YQn{OPQ_T=z1I^LxQ`VRk43NoD{60{a#j#<qV`rTEQ;P+77^<@7Yh%>
zGCX1Ykh00oH+_|3^t%Jz0}uKyY?a+_uf1XCVn6zQ-I?unbIXVEGP>;R%GbnJ&YT`y
z!8K~bRdsG)iOusCdRzBSW8$h3EV>s1n0PFLA6R36H5ypKv(WGK+y29<c`QT1ZJV_+
zKW!@PwOqbMF!%M7r;vtg_7wI#ij#MFDY(8Azu%!DC?^iMSKog&-FYV`S;aZp-@PmK
zFvH@{VPPAQk6%d#?|Zi2ejFPU8n}-Y6OTol0oFUf8Un0P58?@du9JUQHNxUWw1&w8
zI~z~B78`dC4)EGs%!w4rF?Sh8{wRFH5mIcy(0B2j=8AK5>Y-=bPQqBB2bRi-zhp2v
z9Xs!|=0aSN)3LOq+-+!0BTq4&4XggmS(u!UMQ{OY60oB3E9_J>EA)%M<MonAgjS$;
zQB6c=-E*A~mn!blxrxSe*D@2$aOUWO2bB<|y_I(ndV+L(x?R!Cwfb>yd}D@m8^_~&
z1L3hrB-m9AcZ$a-59@@BKpzP%zK<0Hk43NnD|}=0><wxG^xS1+msMdRimzTKB|qzX
zVW@HSnOK%B>ty2-?d)NdHx4t^$Y0ASz5U*^op391M|=nG^(l*L)+-M!Dw~u~7vBgj
z31E*em-~9nxdWMtQ&-E8bw_F+D}02%|HFN37sTgD&=utcnozsIoB&r4zA`EgsEMoi
z+DmD(BXup0=u&&Duw26(Elr4f-7q$i@<i&q7Ae1ZMgH1w<<?>Uu1iCco7$H>im=vN
z3bF*6@qGjnt2Lh>U#B{C&&{9NB7=H<w|2#-hwTO~Vc-k}P9$)`|9^V|4HAY7L*(!k
z9rggRq-Yh}hHTCWUBCIBsO~|+#TWfl{n$7(V@PBe)7F#Q76KmW^jD6X=FQ?SeQFdQ
zlm4ZB$M$TVSv~DJ<)9RNW_FYj?y+LxujpqXL7$CX39P7|0K67`zb48Dv^*Qi47Kwc
z2Y5L>tD=P;bY^1-nfom>h4delFP&CRzlqgaJLOoLL?~f&TATO4x=stv5j&w(oH`-u
zVcdGW&Sk%W1M@QuNwNl8cZNEzsRiz1#iWOQ0c>Z1H3C>cU-bQYFwa8&{+e#Wq9T?D
znUc6U@<t+r%@?xgwHJ^x0q)w*C>#zKGwQ=R?u@Y$Ax4i>hw94m3tYH5%P+Nl!C!b*
zqQ}bXeV0GlKsE2F0Zp`<IIVj69LVuE1~BSjQL&2Zw;uv4{Exm5k_r|@>8h^e9uW8N
zOMqdti(iA@dhJ)<E&3f3CBdv9%O6IZ4aR{~?dAieM>8wWICL=MdtBHFPet0;`SAQG
zs&;G2Htj0FeL1>#csyg+X_G`W5Ek5H#l&Cn>w#4bSRVi@^j!GPaLf#gqGP~Cc$?ik
ztNa$X-oeumNOvQ~VM|%R<phPB8iXa<ai^mKTh|O|O1$#b&-LxZ2(~$G9HB~)r4IKx
z@@Y|!w@t~$KR5Kc?AZ{qimZfZR<Ahsvtn{C7P2$CFAN7(J79(S<?lNvwGzU-4PCY?
z2N*6WK8+asdNTV9YrWr!en;=^XIw+Q9z0mBGxA@H9!7c+-n+qeEhIoanX2+<Q=3O}
zRneye+8Y_aa(?KC2+J-wf0lo~_2$q%R!saAA96nGos;PIW#|Iy6*y*rMbY^x`c{Sg
z#dTwPQ-$W#mRSw9W+<=uB-HM{TMUitKC4E`uc(@RBeP?g@AIvaTG9JC2XY)ETw0?L
zZ~X6Wawa`Cit+dvn~-3x*jC33akZPTF!5G&4qyY$M9|a{I0OH1qJDo37WOnBwV~AU
zll(-_BABqe7|mGu{62J$_r`WRsSRK9#cHDx?k*@Ge{$fnZ{^|JNs+{&${i%1?9O}|
zIICb_DdA@7cQUAfi~4#RC)FcZc^@kV-iluZY}UXE{)&3;F#csYhWZe#sYYgAqM^QE
z$+!;DI&FHV2JgjN56TFAt`QSm#8$_Cj!VSD@+Lp?V9YC3zEIDr8HYOFMjWdqy==#O
zAbUf!lAAZ}R0Mf1yV`Ku&J#1~_ZItEQJ?>{w?2nD0O}R=yu}IX0N7Liv(NU|0@gra
zh1!+C9atrSb+^7!53=qcO_jgHeqR{Ltcu6_lvfZhoTKSZeaolC*NtR|g+*Crt_Q<n
zI0sIhr|3-cOOq8X<bU-VDJr^o?UU|kiiNK7<uAuwTwke0vY=uL{&wRP6MscW04vl0
z=zCwltMT}NbsiQ)?JLc^RRXOwWpeh=GOE|6@h5}1h|KK_EfD0(%;opOlcfAV*R0n3
zx_dAButwY=J%vM0Wg0tJuFeMNSx+2ypT*X5w|<c$(5Vw-_ta($5v91F6%&6&2m>q3
z0kC~RS0i8r|HAzNtD<z!9Hi)TdpxSf9)yrF(IS+a;b6&~;N&HnJnQ<Z?22%-Hi2`k
zM&X+vtf*nCwx4aAMFeAS65T|_K2}SOr4y%qpBA(D(sc8F%hB4-;;3i*a<JsCt{C;N
z=D>OtSfOS`_b!b4VnFk{XaH$c+c;HRat@u=nt}{taUKG>7EM4bJL@t~cDU2K{<(xY
z?pu4d9i*G!$VgK8T?1dKkbwK0cL*yzB#%cr*;$a!MP<nroZEaTKpD876_frI_B{Ii
z4lutWa0FeMfb|P3iq=)}!HQPm@&;Ayl{g~ZuxTr<D>gmIIkGRra-90D>D=Kwb?lQy
z%MvE7Z3b5+jNa1wyi;sJ*abZ9Fs85#!AkjMJ>Nl7*C~-}I-mPEv2i~uCjN?mS{41?
z3{+2e0a&Slbr@Jt@v1&R<tRX{g&04(nayI!O~%wbUN?Gn0BdR5)cYtQzlC$_yeg}v
zb5+)Df{Cbct@dQby7c3;pX%hs8dhI@6KS(>6OEav>d|1RH=0g@%>aM9F@RD3ikb_c
z_7)7rQO_B`gJaw!SQM@6*Jrq-Un3*(>$fwWYfyf4Vdj27iyY}np%A;G{xax+e_#DO
zZFTlTPR1#EN4UPs&W_v2H%Tk{UznR9t8A0JpX_w+0)J**)(sr*cCKSQf%{l7@mB;d
zu)^<;eeF-Yo&?s9yR5f7Zn*U#vGeOcwczSk_21`OG)_g1%4i*Tvi7cv9hu^(ZOAed
z!|~@q=;|3d^({SMHV&Ir3=qC=OdM$}VRq2z-scU9SNOM{n39Q5bpC@CgBdY`<BzVW
z7`O<!Lhrpdb4&J_*ZO8x&7RgBglnjpuW_<SUKKe-L&}{X=w0qf=zFtvxzdQR+UA?4
zk*()Mz=s=Lio!pb0-ko>Ig;Ln#d6*WKYZi+K=9>ks>iWmL;qmKz+=(#!35A1W-aJ<
z(l7w)@~*CvZ`9o~ulZ7?akt<q*s9-W<<S!3XnZ@7%iZ$w7sIw%Urg$DQOT{aAocUf
z#N&fq9j7Q0UaU{W7mT!huz9g#dj~Hv*SDndewatd`6hBHSavsF;UoP0AMRg4L)7^e
z>=e)y=9uW(1?od|9#H-Cg!AET>2pTv_=h-TYxw<Nry=KXMj{VV<oD<)G7(=t^Z7Oz
zYe9NSvnjWROpUi=uh#+EB`WM+w`ATRw>uU&bl&seb~OgOp4ZoaCI4!Q!F(7WycB(I
z+X^_LUPYgOlYvE1IY6VNS}(#}IzykW_PJcdJtqf7d+Gzo55dQkqx4vxP4i#Z!6C0X
zHRfe1IeDemD$F3gQW9Yvgunc;#YBqK|6t?;L-Ulgiz<<4jdLGfH4cOo|FUA@ulQ5I
zN&}i^0xS3x`Z`jW!J%2z9p`VYum6a%58h#IFtfW)&Z5BPYZ55#({;t9W=n`e>*O;T
zo}{-AC%rw(NuPfHdMPEaOgO~zyK|vW@~xLO4f*e8D6&6NzQQSFlYbfri|%8^q<`J`
z!wMct0D9s>O@O=0s*#a#E#itr2$G}nfqsO1hcqM8pu6Z)$f;+lBR1hy`c2PekEL`W
z-x6i0Dy}E*C<}e%xul1Uzy0=}pu*I5ypl?h4iX+628Z<m@3FIo|H+Dhzv2%88>%OW
z0#?XX=;!&KgJsb%z{X5wXus^I!?A@BhzfBv&EDw1^80|V6`_)3`z*6BS{jGG^{ap3
zOi~uD9)i2z0xLqf!t(k@V)0EpUqn+uGErucV6(c*m(*g4!mdCR>~6eb;;;B~zzT7S
zzHUVHkFL;j|7I1PPg8Q>b4}uXv@H^4J4?v;=&GU0Tj`w6s3dzOvmE<_A;;-!8>n-)
zKU{QPs}HzdYa*_MjXyUV=))xyn=)F@w0V<D-uf_?T}-|+6vz8nG3j4Xdy5KK?|`mc
z!0H9Z+^{IB7SQli@L$jIC|40q6VHjBcRa`;k2U7Z_)5Pdu~ORX#yZF3)uowQCI|Cv
z^I5*QL;1+lw@y52^TzUU6m`6<7q*jg=CSqM^*jR#?D5DLZxrfZT`}=j{Apkl2G%fO
zh58D8FElC!(5w{cy>4$mr5Zk}alNW6O}xpaeU>(=PouzrrnqI{1M`Un*HgUoFPT0t
zIffi*RFZu8Eyd<EMZ<^{p?TwlXOy$H28Vc>?Z(e@rrMehKDGY`D+d0GKLM<Cz?uxK
zs65~f$FRqKvpNLKu6KCoj1|3GTQ=z1h&`ENd&uYuFM~{Z-tA)%?NrUQ&F&s5gE;r5
z@j?dZUBeIzrB+9RY3-ye#dtT~SEU;MoQmNu&8w0)Bp`Y9^*>lK@K^j9U^@Y<Wq)D-
z@&PW?GQU|*%wVnG7aH=Tw3Gf(kJnS0Bo`AF$r!TH>NVlzhqdWA6dC=juC2=ud$g!m
z_ZvTNQwn2$j@R@Mhh%WQZp@XO*ih@`gm{COhmp=rm4W+NQJ?>{zdnb$B^6k+fYk}=
zmjC4*Rv%a)U!m{WG67Z*(DglhWmL}6JTe?bn(dm%{`hhnTz>0{&oV_n$FX*)qIB$Q
zG1<cjnrv;2$h@w1$I@O~CSK?+j&!WFdCdQv_R}%%PG6>1+lkgE38Tj($2^-}%VkS_
z*xzem+`~dGh~8U}tMFh?q5DGg-qLU#Aq^y5WDle>Oi&;%)H!8ubtv|i$ztIZBk$t^
zjuPexll`)qC-5xC4j<br!IK{~<lqe^?#Od;>YIuAsnrpuAhr7BRT7;;_{qvUkY)Dr
z045%bFalk#0BZ!qfHAN_{@T;kTE@;T9ch+wV<JdSdr3rwIScy~$2|c{v4_4TAG4}O
z94$jk3OjpV9DH_%?Y>U$p`4$3Q8K}_JM8RZYKyNn<5d}D?)De|+$_p<ZkcJy-Oq}F
z$0Fc|!SM#xNMJPtR`~7ro<H6BPJ}cjv4-k4eeICcd#r5cm*hPh8SWNcF^`IgUwOpT
z>B=q@wZ(?3%!S2_K9*<OF+mbO9H#x5!#H-u@zqx4Afpn)@LTdkLhHO!iuvpNSuyZf
z^fd#(iM}?|6j;FraKX!d=Yee&TG9yo6{_HGUvLcxmq<92e?`V-np|3Jd-(V^JN|rp
zKT=QXgqW#QD?3+is+DRz;gQ4x#6Ow1ob0r?za|+~Fm*WAF?cl~y6Ue+KxW*F0Zcp=
zaSC*W-yb_1bVc>BkiT$0!=mWD)vS#}x{&rQ<*qy}oz8U8hm&HQMI0lR;*Do7m9R`+
zdVx%(cPM=NlfHY@nu$X91$L20a*d&`*yj>~v1{B^I5!VjFtJ6b6y54j@}3-(Bihf3
zfybiXQ;F(dLV?v0SZRS3YOmj{G9kF6kE~azGX032>tx#YoeGbOMsB~#BxbsJf2ONG
zXS>c{Tl+jy!iyu{MJk12QWLmjh9l<!I`E9CueC|u!F4JXil?r1Hx51G-Npgk@t&@j
zcq~F3boByOly8{<s|2v3=5FXbKr2y3`f2+a`^}&9HU?8<7r5$LOm7~aupGIxY$+?m
z^zH6vmh-NiSi_zQOruF-g5R7Ezo*0_{yF#1sfN+(NRw~LD<xId@(#v?pb72hw12Q-
z;IU}dL}0xRta`x846LYr8_mjEhfDgN$tNZEP0|e+Iigrj*UT%imx*+io@t+AkZW0r
zk9XCxxOP9!ElcH4%!z4LF56$o>X@I@_niEt)Nw;3%9;>8X05{XWQpVa^Z#VUz+=&@
zs2D)a25f*8O5(lw(A@F?(qvi5l!R`fA(I^vMsDkyM#zyMLt<KMGwK@ksg=@p7Bzt~
zpQvLOGOO;Ky}VVWVjke$)GDJMc|Y-9)Exo3){=nRoV-SaC6}+xH-iMb^(#u_f9<jP
zzy>wTf8w6lSzv{J0DVsZ<OFmNBMV;{6$5IA#DkTjcZhWKOz_raq!id*oT)5z<GyRr
z5gYnMz#zwo|6xqx`>X2K2gCC0Nh_5HLLB)>U2GqpdJu_t$4y&lfYnoR&O_h<3l?ek
zmwlRI(!;{c3Uz%eHuwShnmnjoQ0I-%IY5ii<#>cg+Kmx*t#6?bqMH3PcJGHo7i2P(
z!_C>d*=>5Op4O&R7ZUS?b%ow;_8i*WJhGy=l_}j+=Uw^ehC>sPntsK}EW=-vqZLhw
z26qzz_gFFUSNsKFQv_CsQS?2g;9ux-TIhX+)0c6Ig|Jbs_AOPlqed`0ORrq1lZ0We
zS7Suhr3~+c($)($bwAf9x7#)l7k6laLX-K!Jt!K^s|#qZ7z&0o%=t{~B3}2t<J~@_
z`-~J6_#3O3%!pC%cR)XjSO{35FF?;Q7+_H}t5!fGfvPh*IkJz;`YeMA$9!Bz>D{w&
z&B%kE1-i#pZCyr8IQxgUwtkYZi#50&+%}Hxlf|y?9`JV$rteYcJa|C7pJ~cBp-#So
z=HU5%^aL26gVh0D`G6H>ZuqDe0Pf#0FdKtM;4FLmAh(j;W~6SjJR|J`Z{!F0TPja2
z><Ou{hR?)4VZ?Kgys6RVCiz*@Lm*6_;pup9<ddgTLRY>I`&LeAru)5ocR^@$;Zedu
z=+)*uT`}@k8er1|R`4qHz2;DF;i4)xl&+eoodX2s-#b+p*3^$ieaI$d^rma)*gkT(
zc|eL(#`SWlmb=v;?vmB}dgoE%(cml!W%livKCJ6I<y|8ujP)<d`)Vk$;|UlMjU;$F
ztQSFR^fv}DnGxfU0^50Dg=j+gD}LY~UD2!weE0;WG@tQ#-X7?Wx?JuvgS<nA+z{qV
zNO*hcbA2Lp-KNzgu8k|{=e&M#Z0FhCi67uSeF?Ahb84FUuZ>S1hWORz*6>8~ipSgZ
z?4t_f_pxH)ujsu6`5zk(STFryWr0P}tVy}{O2%C=-5lKc+R+k3q>StF>?&_Mgv_a%
zoHmK9^~4-zlBIr`<TNa1D-9l5U6LSquOM%pu_G4va$V-tstLE+h3;p3*K)@x9PK~T
zKq>cESB%cXegz%jyc;%n0s1*<;9a<=3IomBaPu3LQ!Gt4>)jhdkvDygGHTyBp|VK-
zankk8gI-On5k_gNsBrE8hbp@sf-^Ng7exgLb$;Pceoh+G{uY1V`njU!eK~~G@jP||
zUsF|<{V{-pfxqHU1Dgu4LJff4Tc}x<@E=wU>4@(HuI4XwqFDG`Raz3BCW+iUM5x2-
z60KWyV@Rx{{dQ4Ro{JyJu<f@kmGbzZZK=x)95t-lV=dZS2U3?!T_vj65JnG=Sik;t
Lo;8OY3+w*@nA-y$
literal 0
HcmV?d00001
diff --git a/regtest/analysis/rt-wham/config b/regtest/analysis/rt-wham/config
new file mode 100644
index 000000000..f5872ebf3
--- /dev/null
+++ b/regtest/analysis/rt-wham/config
@@ -0,0 +1,2 @@
+type=driver
+arg="--mf_xtc alltraj.xtc"
diff --git a/regtest/analysis/rt-wham/fes.dat.reference b/regtest/analysis/rt-wham/fes.dat.reference
new file mode 100644
index 000000000..aacce3031
--- /dev/null
+++ b/regtest/analysis/rt-wham/fes.dat.reference
@@ -0,0 +1,57 @@
+#! FIELDS ff.phi fes
+#! SET normalisation 1.0000
+#! SET min_ff.phi -pi
+#! SET max_ff.phi pi
+#! SET nbins_ff.phi 50
+#! SET periodic_ff.phi false
+ -3.1416 53.0597
+ -3.0159 47.2084
+ -2.8903 44.8862
+ -2.7646 39.8037
+ -2.6389 37.2465
+ -2.5133 34.1986
+ -2.3876 32.0535
+ -2.2619 29.4245
+ -2.1363 28.9321
+ -2.0106 25.3416
+ -1.8850 20.0828
+ -1.7593 15.1329
+ -1.6336 12.5434
+ -1.5080 9.1217
+ -1.3823 4.8916
+ -1.2566 1.0762
+ -1.1310 4.5409
+ -1.0053 11.2837
+ -0.8796 17.3947
+ -0.7540 23.3114
+ -0.6283 27.6347
+ -0.5027 35.2913
+ -0.3770 38.0352
+ -0.2513 42.3201
+ -0.1257 45.6324
+ 0.0000 44.9696
+ 0.1257 44.1181
+ 0.2513 43.5800
+ 0.3770 42.8887
+ 0.5027 37.0370
+ 0.6283 32.2551
+ 0.7540 26.1036
+ 0.8796 22.2865
+ 1.0053 20.8150
+ 1.1310 19.5912
+ 1.2566 25.6362
+ 1.3823 31.3547
+ 1.5080 40.4811
+ 1.6336 52.0285
+ 1.7593 59.3373
+ 1.8850 70.1980
+ 2.0106 77.5370
+ 2.1363 80.1242
+ 2.2619 85.1814
+ 2.3876 82.7074
+ 2.5133 83.5724
+ 2.6389 79.7823
+ 2.7646 72.3409
+ 2.8903 64.2447
+ 3.0159 59.1219
+ 3.1416 56.8935
diff --git a/regtest/analysis/rt-wham/plumed.dat b/regtest/analysis/rt-wham/plumed.dat
new file mode 100644
index 000000000..d53eab062
--- /dev/null
+++ b/regtest/analysis/rt-wham/plumed.dat
@@ -0,0 +1,81 @@
+phi: TORSION ATOMS=5,7,9,15
+psi: TORSION ATOMS=7,9,15,17
+rp0: RESTRAINT ARG=phi KAPPA=200.0 AT=-3.00000000000000000000
+rp1: RESTRAINT ARG=phi KAPPA=200.0 AT=-2.80645161290322580646
+rp2: RESTRAINT ARG=phi KAPPA=200.0 AT=-2.61290322580645161292
+rp3: RESTRAINT ARG=phi KAPPA=200.0 AT=-2.41935483870967741938
+rp4: RESTRAINT ARG=phi KAPPA=200.0 AT=-2.22580645161290322584
+rp5: RESTRAINT ARG=phi KAPPA=200.0 AT=-2.03225806451612903230
+rp6: RESTRAINT ARG=phi KAPPA=200.0 AT=-1.83870967741935483876
+rp7: RESTRAINT ARG=phi KAPPA=200.0 AT=-1.64516129032258064522
+rp8: RESTRAINT ARG=phi KAPPA=200.0 AT=-1.45161290322580645168
+rp9: RESTRAINT ARG=phi KAPPA=200.0 AT=-1.25806451612903225814
+rp10: RESTRAINT ARG=phi KAPPA=200.0 AT=-1.06451612903225806460
+rp11: RESTRAINT ARG=phi KAPPA=200.0 AT=-.87096774193548387106
+rp12: RESTRAINT ARG=phi KAPPA=200.0 AT=-.67741935483870967752
+rp13: RESTRAINT ARG=phi KAPPA=200.0 AT=-.48387096774193548398
+rp14: RESTRAINT ARG=phi KAPPA=200.0 AT=-.29032258064516129044
+rp15: RESTRAINT ARG=phi KAPPA=200.0 AT=-.09677419354838709690
+rp16: RESTRAINT ARG=phi KAPPA=200.0 AT=.09677419354838709664
+rp17: RESTRAINT ARG=phi KAPPA=200.0 AT=.29032258064516129018
+rp18: RESTRAINT ARG=phi KAPPA=200.0 AT=.48387096774193548372
+rp19: RESTRAINT ARG=phi KAPPA=200.0 AT=.67741935483870967726
+rp20: RESTRAINT ARG=phi KAPPA=200.0 AT=.87096774193548387080
+rp21: RESTRAINT ARG=phi KAPPA=200.0 AT=1.06451612903225806434
+rp22: RESTRAINT ARG=phi KAPPA=200.0 AT=1.25806451612903225788
+rp23: RESTRAINT ARG=phi KAPPA=200.0 AT=1.45161290322580645142
+rp24: RESTRAINT ARG=phi KAPPA=200.0 AT=1.64516129032258064496
+rp25: RESTRAINT ARG=phi KAPPA=200.0 AT=1.83870967741935483850
+rp26: RESTRAINT ARG=phi KAPPA=200.0 AT=2.03225806451612903204
+rp27: RESTRAINT ARG=phi KAPPA=200.0 AT=2.22580645161290322558
+rp28: RESTRAINT ARG=phi KAPPA=200.0 AT=2.41935483870967741912
+rp29: RESTRAINT ARG=phi KAPPA=200.0 AT=2.61290322580645161266
+rp30: RESTRAINT ARG=phi KAPPA=200.0 AT=2.80645161290322580620
+rp31: RESTRAINT ARG=phi KAPPA=200.0 AT=2.99999999999999999974
+
+REWEIGHT_WHAM ...
+ ARG1=rp0.bias
+ ARG2=rp1.bias
+ ARG3=rp2.bias
+ ARG4=rp3.bias
+ ARG5=rp4.bias
+ ARG6=rp5.bias
+ ARG7=rp6.bias
+ ARG8=rp7.bias
+ ARG9=rp8.bias
+ ARG10=rp9.bias
+ ARG11=rp10.bias
+ ARG12=rp11.bias
+ ARG13=rp12.bias
+ ARG14=rp13.bias
+ ARG15=rp14.bias
+ ARG16=rp15.bias
+ ARG17=rp16.bias
+ ARG18=rp17.bias
+ ARG19=rp18.bias
+ ARG20=rp19.bias
+ ARG21=rp20.bias
+ ARG22=rp21.bias
+ ARG23=rp22.bias
+ ARG24=rp23.bias
+ ARG25=rp24.bias
+ ARG26=rp25.bias
+ ARG27=rp26.bias
+ ARG28=rp27.bias
+ ARG29=rp28.bias
+ ARG30=rp29.bias
+ ARG31=rp30.bias
+ ARG32=rp31.bias
+ TEMP=300
+ LABEL=ww
+... REWEIGHT_WHAM
+
+PRINT ARG=phi,psi FILE=colvar
+PRINT ARG=rp0.bias,rp1.bias,rp2.bias,rp3.bias,rp4.bias,rp5.bias,rp6.bias,rp7.bias,rp8.bias,rp9.bias,rp10.bias,rp11.bias,rp12.bias,rp13.bias,rp14.bias,rp15.bias,rp16.bias,rp17.bias,rp18.bias,rp19.bias,rp20.bias,rp21.bias,rp22.bias,rp23.bias,rp24.bias,rp25.bias,rp26.bias,rp27.bias,rp28.bias,rp29.bias,rp30.bias,rp31.bias FILE=bias
+
+ff: COLLECT_FRAMES ARG=phi STRIDE=1 LOGWEIGHTS=ww
+hh: HISTOGRAM ARG=ff.phi GRID_MIN=-pi GRID_MAX=pi GRID_BIN=50 KERNEL=DISCRETE
+
+fes: CONVERT_TO_FES GRID=hh TEMP=300
+DUMPGRID GRID=fes FILE=fes.dat FMT=%8.4f
+
diff --git a/src/analysis/DataCollectionObject.cpp b/src/analysis/DataCollectionObject.cpp
index 2b9152f27..9950d7097 100644
--- a/src/analysis/DataCollectionObject.cpp
+++ b/src/analysis/DataCollectionObject.cpp
@@ -25,8 +25,8 @@
namespace PLMD {
namespace analysis {
-void DataCollectionObject::setAtomNumbersAndArgumentNames( const std::vector<AtomNumber>& ind, const std::vector<std::string>& arg_names ){
- indices.resize( ind.size() ); positions.resize( indices.size() );
+void DataCollectionObject::setAtomNumbersAndArgumentNames( const std::string& action_label, const std::vector<AtomNumber>& ind, const std::vector<std::string>& arg_names ){
+ myaction=action_label; indices.resize( ind.size() ); positions.resize( indices.size() );
for(unsigned i=0;i<ind.size();++i) indices[i]=ind[i];
for(unsigned i=0;i<arg_names.size();++i) args.insert( std::pair<std::string,double>( arg_names[i], 0.0 ) );
}
diff --git a/src/analysis/DataCollectionObject.h b/src/analysis/DataCollectionObject.h
index df44aa377..9c5c723cf 100644
--- a/src/analysis/DataCollectionObject.h
+++ b/src/analysis/DataCollectionObject.h
@@ -36,6 +36,8 @@ namespace analysis {
class DataCollectionObject {
friend class ReadAnalysisFrames;
private:
+/// The label of the action in which the data is stored
+ std::string myaction;
/// The list of atom numbers that are stored in the object
std::vector<AtomNumber> indices;
/// The list of atomic positions
@@ -44,7 +46,7 @@ private:
std::map<std::string,double> args;
public:
/// Set the names and atom numbers
- void setAtomNumbersAndArgumentNames( const std::vector<AtomNumber>& ind, const std::vector<std::string>& arg_names );
+ void setAtomNumbersAndArgumentNames( const std::string& action_label, const std::vector<AtomNumber>& ind, const std::vector<std::string>& arg_names );
/// Set the positions of all the atoms
void setAtomPositions( const std::vector<Vector>& pos );
/// Set the value of one of the arguments
@@ -64,6 +66,8 @@ Vector DataCollectionObject::getAtomPosition( const AtomNumber& ind ) const {
inline
double DataCollectionObject::getArgumentValue( const std::string& name ) const {
+ std::size_t dot=name.find_first_of('.'); std::string a=name.substr(0,dot);
+ if( a==myaction ) return args.find( name.substr(dot+1) )->second;
return args.find(name)->second;
}
diff --git a/src/analysis/Histogram.cpp b/src/analysis/Histogram.cpp
index 53f16ebb1..6b289cbc0 100644
--- a/src/analysis/Histogram.cpp
+++ b/src/analysis/Histogram.cpp
@@ -28,6 +28,7 @@
#include "core/PlumedMain.h"
#include "core/ActionSet.h"
#include "AnalysisBase.h"
+#include "ReadAnalysisFrames.h"
namespace PLMD{
namespace analysis{
@@ -162,6 +163,7 @@ DUMPGRID GRID=hh FILE=histo STRIDE=100000
class Histogram : public gridtools::ActionWithGrid {
private:
double ww;
+ bool activated;
bool in_apply, mvectors;
KernelFunctions* kernel;
std::vector<double> forcesToApply, finalForces;
@@ -202,6 +204,7 @@ Histogram::Histogram(const ActionOptions&ao):
Action(ao),
ActionWithGrid(ao),
ww(0.0),
+activated(false),
in_apply(false),
mvectors(false),
kernel(NULL),
@@ -357,7 +360,10 @@ unsigned Histogram::getNumberOfQuantities() const {
bool Histogram::onStep() const {
if( !my_analysis_object ) return ActionPilot::onStep();
- return true;
+ ReadAnalysisFrames* myfram = dynamic_cast<ReadAnalysisFrames*>( my_analysis_object );
+ if( myfram && activated ) return true;
+ else if( myfram ) return ActionPilot::onStep();
+ return my_analysis_object->onStep();
}
void Histogram::prepareForAveraging(){
@@ -522,7 +528,7 @@ void Histogram::apply(){
}
void Histogram::runFinalJobs(){
- if( my_analysis_object && getStride()==0 ) update();
+ if( my_analysis_object && getStride()==0 ){ activated=true; update(); }
}
}
diff --git a/src/analysis/Makefile b/src/analysis/Makefile
index 02080ce4b..937de0085 100644
--- a/src/analysis/Makefile
+++ b/src/analysis/Makefile
@@ -1,4 +1,4 @@
-USE=core tools reference vesselbase gridtools multicolvar
+USE=core tools reference vesselbase gridtools multicolvar bias
#generic makefile
include ../maketools/make.module
diff --git a/src/analysis/ReadAnalysisFrames.cpp b/src/analysis/ReadAnalysisFrames.cpp
index 43ea1d855..478412346 100644
--- a/src/analysis/ReadAnalysisFrames.cpp
+++ b/src/analysis/ReadAnalysisFrames.cpp
@@ -79,7 +79,12 @@ weights_calculated(false)
}
if( wwstr.size()>0 ) log.printf("\n");
else log.printf(" weights are all equal to one\n");
+ wham_pointer = dynamic_cast<bias::ReweightWham*>( weight_vals[0]->getPntrToAction() );
+ if( wham_pointer && weight_vals.size()!=1 ) error("can only extract weights from one wham object");
requestArguments( arg );
+
+ // Now add fake components to the underlying ActionWithValue for the arguments
+ for(unsigned i=0;i<argument_names.size();++i){ addComponent( argument_names[i] ); componentIsNotPeriodic( argument_names[i] ); }
}
std::vector<Value*> ReadAnalysisFrames::getArgumentList(){
@@ -102,13 +107,18 @@ void ReadAnalysisFrames::calculateWeights(){
if( weight_vals.empty() ){
for(unsigned i=0;i<logweights.size();++i) weights[i]=1.0;
} else {
- // Find the maximum weight
- double maxweight=logweights[0];
- for(unsigned i=1;i<getNumberOfDataPoints();++i){
- if(logweights[i]>maxweight) maxweight=logweights[i];
+ if( wham_pointer ){
+ wham_pointer->calculateWeights( logweights.size() );
+ for(unsigned i=0;i<logweights.size();++i) weights[i]=wham_pointer->getWeight(i);
+ } else {
+ // Find the maximum weight
+ double maxweight=logweights[0];
+ for(unsigned i=1;i<getNumberOfDataPoints();++i){
+ if(logweights[i]>maxweight) maxweight=logweights[i];
+ }
+ // Calculate weights (no memory) -- business here with maxweight is to prevent overflows
+ for(unsigned i=0;i<logweights.size();++i) weights[i]=exp( logweights[i]-maxweight );
}
- // Calculate weights (no memory) -- business here with maxweight is to prevent overflows
- for(unsigned i=0;i<logweights.size();++i) weights[i]=exp( logweights[i]-maxweight );
}
}
@@ -118,8 +128,9 @@ void ReadAnalysisFrames::update(){
if( clearonnextstep ){
my_data_stash.clear(); my_data_stash.resize(0);
logweights.clear(); logweights.resize(0);
+ if( wham_pointer ) wham_pointer->clearData();
clearonnextstep=false;
- }
+ }
// Get the weight and store it in the weights array
double ww=0; for(unsigned i=0;i<weight_vals.size();++i) ww+=weight_vals[i]->get();
@@ -127,7 +138,7 @@ void ReadAnalysisFrames::update(){
// Now create the data collection object and push it back to be stored
unsigned index = my_data_stash.size(); my_data_stash.push_back( DataCollectionObject() );
- my_data_stash[index].setAtomNumbersAndArgumentNames( atom_numbers, argument_names );
+ my_data_stash[index].setAtomNumbersAndArgumentNames( getLabel(), atom_numbers, argument_names );
my_data_stash[index].setAtomPositions( getPositions() );
for(unsigned i=0;i<argument_names.size();++i) my_data_stash[index].setArgument( argument_names[i], getArgument(i) );
diff --git a/src/analysis/ReadAnalysisFrames.h b/src/analysis/ReadAnalysisFrames.h
index 3da79c53d..53d766aaa 100644
--- a/src/analysis/ReadAnalysisFrames.h
+++ b/src/analysis/ReadAnalysisFrames.h
@@ -23,6 +23,7 @@
#define __PLUMED_analysis_ReadAnalysisFrames_h
#include "AnalysisBase.h"
+#include "bias/ReweightWham.h"
namespace PLMD {
namespace analysis {
@@ -38,6 +39,8 @@ private:
std::vector<AtomNumber> atom_numbers;
/// The biases we are using in reweighting and the args we store them separately
std::vector<Value*> weight_vals;
+/// The object that calculates weights using WHAM
+ bias::ReweightWham* wham_pointer;
/// The weights of all the data points
bool weights_calculated;
std::vector<double> logweights, weights;
diff --git a/src/bias/ReweightBase.h b/src/bias/ReweightBase.h
index e3dfd35ad..bc570c4a5 100644
--- a/src/bias/ReweightBase.h
+++ b/src/bias/ReweightBase.h
@@ -40,7 +40,7 @@ public:
explicit ReweightBase(const ActionOptions&ao);
unsigned getNumberOfDerivatives(){ return 0; }
void calculate();
- virtual double getLogWeight() const = 0;
+ virtual double getLogWeight() = 0;
void apply(){}
};
diff --git a/src/bias/ReweightBias.cpp b/src/bias/ReweightBias.cpp
index d5459df96..f3fa47fcd 100644
--- a/src/bias/ReweightBias.cpp
+++ b/src/bias/ReweightBias.cpp
@@ -72,7 +72,7 @@ class ReweightBias : public ReweightBase {
public:
static void registerKeywords(Keywords&);
explicit ReweightBias(const ActionOptions&ao);
- double getLogWeight() const ;
+ double getLogWeight();
};
PLUMED_REGISTER_ACTION(ReweightBias,"REWEIGHT_BIAS")
@@ -88,7 +88,7 @@ ReweightBase(ao)
{
}
-double ReweightBias::getLogWeight() const {
+double ReweightBias::getLogWeight(){
// Retrieve the bias
double bias=0.0; for(unsigned i=0;i<getNumberOfArguments();++i) bias+=getArgument(i);
return bias / simtemp;
diff --git a/src/bias/ReweightMetad.cpp b/src/bias/ReweightMetad.cpp
index e53170be7..2085e65a0 100644
--- a/src/bias/ReweightMetad.cpp
+++ b/src/bias/ReweightMetad.cpp
@@ -70,7 +70,7 @@ class ReweightMetad : public ReweightBase {
public:
static void registerKeywords(Keywords&);
explicit ReweightMetad(const ActionOptions&ao);
- double getLogWeight() const ;
+ double getLogWeight();
};
PLUMED_REGISTER_ACTION(ReweightMetad,"REWEIGHT_METAD")
@@ -86,7 +86,7 @@ ReweightBase(ao)
{
}
-double ReweightMetad::getLogWeight() const {
+double ReweightMetad::getLogWeight(){
// Retrieve the bias
double bias=0.0; for(unsigned i=0;i<getNumberOfArguments();++i) bias+=getArgument(i);
return bias / simtemp;
diff --git a/src/bias/ReweightTemperature.cpp b/src/bias/ReweightTemperature.cpp
index 45f1c6f2b..19fb36ddb 100644
--- a/src/bias/ReweightTemperature.cpp
+++ b/src/bias/ReweightTemperature.cpp
@@ -59,7 +59,7 @@ public:
static void registerKeywords(Keywords&);
explicit ReweightTemperature(const ActionOptions&ao);
void prepare();
- double getLogWeight() const ;
+ double getLogWeight();
};
PLUMED_REGISTER_ACTION(ReweightTemperature,"REWEIGHT_TEMP")
@@ -95,7 +95,7 @@ void ReweightTemperature::prepare(){
plumed.getAtoms().setCollectEnergy(true);
}
-double ReweightTemperature::getLogWeight() const {
+double ReweightTemperature::getLogWeight(){
// Retrieve the bias
double bias=0.0; for(unsigned i=0;i<biases.size();++i) bias+=biases[i]->get();
double energy=plumed.getAtoms().getEnergy()+bias;
diff --git a/src/bias/ReweightWham.cpp b/src/bias/ReweightWham.cpp
new file mode 100644
index 000000000..892c656d7
--- /dev/null
+++ b/src/bias/ReweightWham.cpp
@@ -0,0 +1,122 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2011-2016 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "core/ActionRegister.h"
+#include "ReweightWham.h"
+
+//+PLUMEDOC REWEIGHTING REWEIGHT_WHAM
+/*
+
+\par Examples
+
+*/
+//+ENDPLUMEDOC
+
+namespace PLMD {
+namespace bias {
+
+PLUMED_REGISTER_ACTION(ReweightWham,"REWEIGHT_WHAM")
+
+void ReweightWham::registerKeywords(Keywords& keys ){
+ ReweightBase::registerKeywords( keys ); keys.remove("ARG");
+ keys.add("numbered","ARG","the biases that must be taken into account when reweighting");
+ keys.add("compulsory","MAXITER","1000","maximum number of iterations for WHAM algorithm");
+ keys.add("compulsory","WHAMTOL","1e-10","threshold for convergence of WHAM algorithm");
+ keys.reset_style("ARG","compulsory"); keys.add("hidden","DATA","sneaky trick to remove error for reading numbered args");
+}
+
+ReweightWham::ReweightWham(const ActionOptions&ao):
+Action(ao),
+ReweightBase(ao)
+{
+ std::vector<Value*> targ, fagr;
+ unsigned nbias = 0; wlists.push_back( 0 );
+ for(unsigned i=1;;i++){
+ if( !parseArgumentList("ARG",i,targ ) ) break;
+ log.printf(" bias number %d involves :");
+ for(unsigned j=0;j<targ.size();++j){
+ log.printf("%s ",targ[j]->getName().c_str() );
+ fagr.push_back( targ[j] );
+ }
+ log.printf("\n"); targ.resize(0);
+ wlists.push_back( fagr.size() );
+ nbias++;
+ }
+ plumed_assert( wlists.size()==(nbias+1) ); requestArguments( fagr );
+ parse("MAXITER",maxiter); parse("WHAMTOL",thresh);
+}
+
+double ReweightWham::getLogWeight(){
+ if( getStep()==0 ) return 1.0; // This is here as first step is ignored in all analyses
+ for(unsigned i=0;i<wlists.size()-1;++i){
+ double total_bias=0;
+ for(unsigned j=wlists[i];j<wlists[i+1];++j) total_bias+=getArgument(j);
+ stored_biases.push_back( total_bias );
+ }
+ return 1.0;
+}
+
+void ReweightWham::clearData(){
+ stored_biases.resize(0);
+}
+
+void ReweightWham::calculateWeights( const unsigned& nframes ){
+ if( stored_biases.size()!=(wlists.size()-1)*nframes ) error("wrong number of weights stored");
+ // Get the minimum value of the bias
+ double minv = *min_element(std::begin(stored_biases), std::end(stored_biases));
+ // Resize final weights array
+ plumed_assert( stored_biases.size()%(wlists.size()-1)==0 );
+ final_weights.resize( stored_biases.size() / (wlists.size()-1), 1.0 );
+ // Offset and exponential of the bias
+ std::vector<double> expv( stored_biases.size() );
+ for(unsigned i=0;i<expv.size();++i) expv[i] = exp( (-stored_biases[i]+minv) / simtemp );
+ // Initialize Z
+ std::vector<double> Z( wlists.size()-1, 1.0 ), oldZ( wlists.size()-1 );
+ // Now the iterative loop to calculate the WHAM weights
+ for(unsigned iter=0;iter<maxiter;++iter){
+ // Store Z
+ for(unsigned j=0;j<Z.size();++j) oldZ[j]=Z[j];
+ // Recompute weights
+ double norm=0;
+ for(unsigned j=0;j<final_weights.size();++j){
+ double ew=0;
+ for(unsigned k=0;k<Z.size();++k) ew += expv[j*Z.size()+k] / Z[k];
+ final_weights[j] = 1.0 / ew; norm += final_weights[j];
+ }
+ // Normalize weights
+ for(unsigned j=0;j<final_weights.size();++j) final_weights[j] /= norm;
+ // Recompute Z
+ for(unsigned j=0;j<Z.size();++j) Z[j] = 0.0;
+ for(unsigned j=0;j<final_weights.size();++j){
+ for(unsigned k=0;k<Z.size();++k) Z[k] += final_weights[j]*expv[j*Z.size()+k];
+ }
+ // Normalize Z and compute change in Z
+ double change=0; norm=0; for(unsigned k=0;k<Z.size();++k) norm+=Z[k];
+ for(unsigned k=0;k<Z.size();++k){
+ Z[k] /= norm; double d = std::log( Z[k] / oldZ[k] ); change += d*d;
+ }
+ if( change<thresh ) return;
+ }
+ error("Too many iterations in WHAM" );
+}
+
+}
+}
diff --git a/src/bias/ReweightWham.h b/src/bias/ReweightWham.h
new file mode 100644
index 000000000..dbf16d79e
--- /dev/null
+++ b/src/bias/ReweightWham.h
@@ -0,0 +1,54 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2011-2016 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_bias_ReweightWham_h
+#define __PLUMED_bias_ReweightWham_h
+
+#include "ReweightBase.h"
+
+namespace PLMD {
+namespace bias {
+
+class ReweightWham : public ReweightBase {
+private:
+ double thresh;
+ unsigned maxiter;
+ std::vector<unsigned> wlists;
+ std::vector<double> stored_biases;
+ std::vector<double> final_weights;
+public:
+ static void registerKeywords(Keywords&);
+ explicit ReweightWham(const ActionOptions&ao);
+ void calculateWeights( const unsigned& nframes );
+ void clearData();
+ double getLogWeight();
+ double getWeight( const unsigned& iweight ) const ;
+};
+
+inline
+double ReweightWham::getWeight( const unsigned& iweight ) const {
+ plumed_dbg_assert( calculatedWeights && iweight<final_weights.size() );
+ return final_weights[iweight];
+}
+
+}
+}
+#endif
diff --git a/src/core/ActionWithArguments.cpp b/src/core/ActionWithArguments.cpp
index 9a799a2cf..955013e77 100644
--- a/src/core/ActionWithArguments.cpp
+++ b/src/core/ActionWithArguments.cpp
@@ -47,9 +47,10 @@ void ActionWithArguments::registerKeywords(Keywords& keys){
}
void ActionWithArguments::parseArgumentList(const std::string&key,std::vector<Value*>&arg){
- vector<string> c; arg.clear(); parseVector(key,c);
+ std::string def; vector<string> c; arg.clear(); parseVector(key,c);
if( c.size()==0 && (keywords.style(key,"compulsory") || keywords.style(key,"hidden")) ){
- std::string def; if( keywords.getDefaultValue(key,def) ) c.push_back( def );
+ if( keywords.getDefaultValue(key,def) ) c.push_back( def );
+ else return;
}
interpretArgumentList(c,arg);
}
diff --git a/src/gridtools/HistogramOnGrid.cpp b/src/gridtools/HistogramOnGrid.cpp
index 8adb3b13e..8b4c85fdb 100644
--- a/src/gridtools/HistogramOnGrid.cpp
+++ b/src/gridtools/HistogramOnGrid.cpp
@@ -109,8 +109,8 @@ void HistogramOnGrid::calculate( const unsigned& current, MultiValue& myvals, st
for(unsigned i=0;i<dimension;++i) point[i]=myvals.get( 1+i );
// Get the kernel
- unsigned num_neigh; std::vector<unsigned> neighbors;
- std::vector<double> der( dimension );
+ unsigned num_neigh; std::vector<unsigned> neighbors(1);
+ std::vector<double> der( 0 );
KernelFunctions* kernel=getKernelAndNeighbors( point, num_neigh, neighbors );
if( !kernel && getType()=="flat" ){
--
GitLab