From e29d5504d4f8eb516433f6272752bcd0a23f3802 Mon Sep 17 00:00:00 2001
From: SebastianBruijns <>
Date: Wed, 25 Jan 2023 21:31:55 +0100
Subject: [PATCH] outsourced and combined code, new class for code dealing with
 MCMC convergence analysis

---
 .../dyn_glm_chain_analysis.cpython-37.pyc     | Bin 54287 -> 54049 bytes
 .../mcmc_chain_analysis.cpython-37.pyc        | Bin 0 -> 4481 bytes
 .../process_many_chains.cpython-37.pyc        | Bin 0 -> 3920 bytes
 __pycache__/simplex_plot.cpython-37.pyc       | Bin 2742 -> 2949 bytes
 canonical_infos.json                          |   2 +-
 dyn_glm_chain_analysis.py                     | 256 ++++--------------
 dyn_glm_chain_analysis_unused_funcs.py        |  36 ++-
 mcmc_chain_analysis.py                        | 123 +++++++++
 process_many_chains.py                        | 210 ++++++--------
 simplex_plot.py                               |  38 +--
 10 files changed, 318 insertions(+), 347 deletions(-)
 create mode 100644 __pycache__/mcmc_chain_analysis.cpython-37.pyc
 create mode 100644 __pycache__/process_many_chains.cpython-37.pyc
 create mode 100644 mcmc_chain_analysis.py

diff --git a/__pycache__/dyn_glm_chain_analysis.cpython-37.pyc b/__pycache__/dyn_glm_chain_analysis.cpython-37.pyc
index a82afa9c345c7e3df3fd0494a1e8695fda475388..b4b56a55e7c7e3ff5cceb8da95a87d51c338cbac 100644
GIT binary patch
delta 22891
zcmb_^34B|{@xPuV%kr7SIqW#Zj^iB8;W*BHIu}VG2ZV$iN`fMLpJmyyWUu7JPUMP^
zgsT`FmI4WZ+JOS$F6JubrrbA_wmc|>QYfK7@0PZ|mfv^gJxj6)zaRYm|M;VMd(7_6
z&d$!x?#|}2n;qRZI&zoh<Ye0L@5sT|YKnIp%FWlZ4~>5-E^Hz@II6M8Q`A`ODQ+C?
z8QoaoDbZ}Yvwn=o5xHVOhuu?(*jO=84C>G{n;7gFr;pbs=w%{rzb1x=p<<XAF7icz
z7$FM9NO6i7<<Hu0->-QlqV^(DEJlkGF-DY%v0|JUFD8hxqgkFwy0dGt9-mY{rCalq
zBQjM?OjKn{R8{29LVgAEr-{67?dx!A?Vd`;8P7N-LA|vU`7@9|o!b~CCW*;niYNzZ
z_#0t<xRDus5iYk6jZZx2D6n~Eim9SPO!GTDvpO{~SX7GX2Q<&@l{PU$%sgNdGxys)
zbHpq$`+&_eSIiM}@tr5;iTU`>7YoEfd=~(Fk*ET86|fhgo;<NwEI~bskh4@QL(Vef
zREg!{RNyZbE5yn;0!x58M64ETfVxzy73+|;Oq?Ok#CN$^FV4dERNxF18^lK7tU!2}
zI9qH&c%|4Zwjh6%*eX2uo+h@5?f9N9c8GKET@Bpf!Yg(HcMZb%!Y8T`UMp&ZKz@9k
zeumH+?44<Hm#7tfQ6~bT-k&CH;L4f$nZl1ehn%O6=vuFzB^t!E{dSDjhR_J))pFiO
zD=*lar*G6Zbe-L8gLrK6+dP{w09tTM<5s&(vo(4=+cX<ex8t`%tcrWi(a)jAQwxbk
zKZ<s0twDX$a$BdiquS%uckV}VeWwWRhlJ`obuWK6=-vXG2;(;q1p9O|-ekqAMYE?y
z&tTGJm(acXnR<<0t^2yuM1(sM*RAa6NwaIM+?JE(`mNkuC(W(Xb)gCEpk3@HSZfD?
zx?bBiEsn%m-QPD6X@0%#D`|E5S{)VW$kX-u01c5S9$^!$9ckO`o_bG%zPT%iJ~q^1
zh`^(IiyLhVY<h#((_rg#%75zH35W7Snfi{dxE1MiG#}ecY~9RkHOD<oF2&e3q8b|c
zy<nIf)l$l4q{MntLgkh$q{*fPey}$;f#04;3DHoB?+6Hg+!KaW7Kn4jK1gMgI8W?{
zR5ptP;vl|KoG&iGHv*gy;zDr|aH8U3@m-{~Sb{rPy2T};16t`)ahd3ZZaO3`7s{XM
z*@fA@Ta-XtTU)Qw+w@)gwdJ-QmxEV(JbU%@t{5~|e1a&__YAjnJJ7Q|!q(vEbjeZ3
zZ11*peaCNOMy7X>5-sB1@6}r|FvSP$3v8{+(ZIRl3h_O0Ww!%WS9<p8`;cEu;|->#
z*Pln_NFn;2daFK5Z_!)#YZ$hxI?!IHt9eCyv#m3&`BShyy=#9rNw~?JubYZ4`7#*V
zXDV9pm^J-2hbCIX;gVlidC941mTm9F74fUZVZYsTAa(jd3O11xh-Ps_TqBO657+u#
ztPyI&b)JLb`=U#9i(}$?afA4QxKZ4M*v<G|FODPU7KFNyb0b2O(<N>dw~5=u9pcWT
zBR%JfyF3@@7x<wA{5J8!_)(|lLLR5PF-{k8XaYhPLp^*Kzf17z!0%H0F2k>rF-8O9
zkm&JTZpQS>N#Z?<WA`9*g&64h9<Z<M%xL|bbx6jJ-#~|WyRH)Vc59vKk#}{kxDV<p
zo-Y5>ncn<vuR=kgt2@)%Dm&A6BTd}Tv%yp#516yO+it0j!^r<pO8$fW@{d^B;m7^b
zuhCcPhs8tU;R81Pn)q_4N&N^|JWxD>dF=8W^<3+@&hvdwm#14lq<7D@c@8l#otU>>
zERI(qbRC5_-ISNno@4sq`s>lFNA+Xku@1*}r{@OG4?H(|ZUR?s7LRjtW%@Ci1N$}b
z=s4PZA{ftUyv1`1#O@aTI0(Vdb1Ud>)Q{`8*4i=4pG;8y6LbEHr%bidklvXo@95>5
z)z>68t5;05<DOf^D$i}L59zmcW=8T8YW!*D-0gAwj{USy+zDnrgQo86);MFQcowl+
z^*bA!U3Yb38gi<$>xX7qx5nalB#D8k>+H;s=YSyetiYxpDX?|j-EB*W^>D1HPhPw?
zip>2^@L{KTj+zqJ5n&B>&uyI<&F7G26wfmWNEO#Xa57~z(q1r$)sNH%ICPKIoqI8)
z(Ea+oY+W1yQN0*<T1xpT=m+pQt{(xDlcHlW8Dn`Ls(7(O^W4wFaX+wL`pQt8?s>ps
z{L9FDfY!Hehb5l9;`ho|a_+~-Yr!8iK8OW5ZfpFp=OOg_e*J#&YJu$(8ybA5{$VqQ
zRA>Do`u)~oQzTyNz$%t+6R(RmU^G4~Iy{d8`!W4t@l(+Q5q`77>3JNg><RHRvsv+a
zx5Lt*PqKR32yy(0_&JN?;rgqgSPN|RPtodK|FkKYkcMZVI-i9-IZ*rpwq%j{=T8f0
z;w@oxI3ZTgd7eiv-VWLtOFg9+E}M8KnAZ4$=LO6Pd+QOsP5e^4i?yJ`?<9%uwC^~8
zQZIU5A{~+F(62<iq!#pfUPdoo(I3`d67PxkI~<-@(eKyvSJCjVJM1X<y5|jy;p^5Y
z??fd()!(p$@+I*>gR3)Jehl4lO=|_r=GXN%_S>-@l%V9Bo}cO0V8W8<zN)|307XSz
z=INW=d^+?$#OOCHJT&dD!8j)sMz1vR9A~xhAv5k8OKW_@OnwO}X(0Z7+u=fk&>#9z
zD435^FukHL=7|m$hWv%DpX+g2!Dx;8gvQD93(M^P1ZF?B0WjuE#P6`4m0}n^1=eT`
z%2=#tdEztixqlFE0p1c5J$uCOxj6&9e?xzr(uqWx1H!bhUl)Hk8Adz><&9QP3iH-j
z%=?-6BW8AvJb>l72-<KD)9vhf`>QO~Uo{1_5bV`oHrc(v)|u0M09<}ve^C?-Cuw^}
ze_8xV{8{{^!vXnkE!CgFSdH}jG6m&b{g>je;%^<WcToIKsl~f~WmyG@3cKFxpHlz6
z_&bK?9>~+LhOoL7cZ)B$wyfr#QEe>t@96i(6O`Cmi<yvoB80P3iO6l<(_1V_0AY2#
zr$3YAAx|~4^mkIz#XrpE&6>M@jXcaaVe5{QMQ0ZHK>Snu%R)=9|Bd)B@o(A;Vm5!s
zGnZ!df2>*kWw)d2BN$~c&V|-(nwypZ1A_+Zsne$awjr}~kbF$Dc|N8|xo0P4(8na`
zC}a1W0N*z+x5dHuPxMcyC%PAEb%3lx3`*H^qH{oV9s0Zz=nhNGp5X5Z*s^KByod8J
zzfmQYYuK}dn`SoADq-<Wu#_$K16Bw-1yeY&t5^pKi@yj^jg7;L7D8Kr4YPTmbfWq9
z^b1j;3rr5`FEVQ{6sQWZYAc575c6@2tpwnwgGVs1Ej9c*Xa&q1&u6XYwh?)5%I)2T
zbZT*tfrQVxwQiT^bG@*uu-oSOz5c1G)%DNKN+9*0a^Eru*Tm9)pKt5NG{XRW#`EH{
z&9+cJ*!~BweV}|^gDt7q&a`c{h3s2wodcWC!oX-EE4BIVw61u{D0=?LjbwAB3-wP~
ztnQ)}O#dBXG<E*Olpoaph^lfr(mj9HuSWm>O#LT*VYue#SL+vImwT9(t;0}fNxP|6
zogMD_i(Z5pa`&VEtJtjOByj~)Ml$}_{kGJYDC(12ECwb@jyA)av91Zt^YONixAgWK
zq6ngJ^q^k}^>(GnS2R7q8lKaxWEmb9{9EJSJ%0x?|K5v?`9BbdgP<f(?FQR~Ld|A7
zLFX&$3~iP%;Ec6rde@!3*1K`lS*I>pzQIuE4jwwu7_fHrnm)<HHl4a`&1p-$YnN;}
zZEY(0ql2qwmN|`qk*F`KdqXXa-a0+lq)Q{y?`v%Id23rjHAbc{*i>hQ(}O`y>=~W6
zNy{)azoAtdnpC@oY;)y<ELoy{F=YA>BV%_UTIa0}_#%d*DH1qtlT%dD&<z7711I97
z4rryCoZ7Q*=u~aYG%8q0U`hg%BV@R}-asf2^?EDS8^d1J^3@H)=jBYGoU}c<3`gWN
z^}FGx7tTaVS-PA}8R?Oh>UzB<YB(dlUAm#p%%3@=Z>@4M@aCw)`Mx-XqP8wAW~;a3
z=O7EdU2B~Xby5l~A3H|0+oNd$mLF1FF&G{QX)$c9nn!K6#V{mg_I3Uvhw0z(DiA@`
z$tnV~2+RO5(xfls*X008puUzlOj<`|CiosPv_NX_)uDpCyi<W2aRb=2!P-D=gtn(h
zEnYOd=jDO|yH=<^E?hB}cwT0gYw_W>%It={yh5Hq$r`s})=KRBXXb)nAQG(!H#RNW
zu8xlxr9LmpR%u&@sNarRI)Rvkn$E(<a0Ei4Of#IpKw~6XUT#3%4s~Yf+{uirNoKrR
zMs7rAY7Zz%1<pZ%Vzp>&>6p|OwBX5by}h_&Q*h;`uc>*bx?t?&#U+%BSz(jRi(C?#
zkgrkqdQUnQ89Jm~mk@Uz3b|I@HL74Zb;L+(l7UdvaMgxoqc2)ENX}O8jas3dp~e@@
znoZq<Ohxr3Bdw*$8x1r4VA$EzG>M*Mk<F2Mv_*!zHFdr~D6+ceilSdTrhj#mn9nC4
zCg$3HjL*5_Q@KC{U`V++O6TSn_rSLul8e6;o4@<5urHxO$x0a%^~hMCR-lT<l`lQX
z+>l$SE35$4QfMmy4}on2Xk6rW0y_v)tIly-G((*j_Z_=2M%{UAn2L`d6*tmCUcI)q
zY?wI?=Mbuwz)k`_0@Va+2nYf+Kcr5emVlo?9f1G=7Ug;hH2@eIl!>ehM~v)XxW*T>
zhSqRw@<uVW>%tSf(XSN-Z=Emd^@W62k3^z9%O{j*E}GEjNw93>X-!mpPHi|Sbm4`P
zXwZespfpEQQ!S(|04>)-=&b;K=5ta4nvr#{dadj;oZ<nKMNN|nsX&&g7g!My?~RN)
z%>S@o`Wod0guIu)63QX9*%S<;r-UMxQJOPS7v5c#X6mYpns7@f>TUG3GV=@<$Dy{o
zy%TvT>f$MoMLqFyQ;|p1-10)zIO*lDQkUJ2I{96qrY~Lc5=!eJa3KMrOI`|KktHvq
z1Rl*!3LPTAyh$*~7|i@t4^Ga}o?~i0nEdAdGc_zv{}*ap%ols)r#*Aa@6cuunHX3S
zd#nnqd}%01PaDoaRB!Cp8=ml=s8cFdU^NsK>#}HE&GlxQdblFr#q(*qdb?tzwpM*n
zF~-FcXNej;tu#*S2}T<0qBkH+HfEWl>*TRdhuRO+=n-Q;O*j+*33{le4MQH4hpA*n
zGyt_q#Z8WrTp2E2!XO6`Bdr$Yt@R?w{$yS0s7ywDzfvzu+n^Qn6jeIYM%+r|)2N-e
zg+d8EE(_Gc*#+8YwRv{sw8Vf~>k~JR^~o&)jZzmaH71it_?SPV)q}Iw77n9mF#s^r
zYXg3oO$cf=a)c_GbHA&I@(b0k<}}TiMnP6HW=G9s8MMl)Dd!-8^9itKVwN9K`{zzy
z(T+qY2i+Gm90Ie)wF^aS8d$Q8tVm6usjZwUi^#jEI8QYfXlao;2BO9n#WEx>QU5h}
z#0(auF2Y()fQPM{LOg2DXdBG36%?OAKoXd#PM<e-Ca*m!5eEr@M#J4y<MZy)YoM+K
zhLD0-tRkk~TSRD@Pu{N{m{%Gnnhb4^ksGe94b%jn+FEJFPl>m&1;D}{qZWC7(R4E$
z=!M}5d1ZKa#IU2?^qOj)^g=yYB7G;}vp&c$*$KM7L|UW^z^1u1XV13z0nIfMF*#a2
zzhKVf(FkIMV0}b;NBU;UO^knX&GW1Bg{KZ@WtVDm@j$Ls7cU$!o~9JrU}YJmjJu+~
z2H3bpwn>h!sY&-qL%p!DsN{6QtpPCX5g&e{Uz0|jTE8ey+pdZiot0BS$RxLA1LS5E
zUo<K{hf)?3IFmq-zySiy1TH6_2wXvc<=mjq2LzIR;}FYm1v0275$c~?bGtHg@GHw)
z;?8hq(O+(H<_<|%HFkiWkyYz8=7l+bcmOP`>}GXfaq)m;)nt+FQ$34k$9ZYJ22nPX
z40m-P#424Lp)?X~UYiY9ODNFXqO<6OmAsgIj9g<txH_Wa&P;fjp%M8mrRJL4t(Cr-
zXdoOi7Yj#3Z^m@M8rX|_GtWMYQmRhUuEreScMxh&7cQBiy`%11^2!<>Wli2e2+Vxa
z5)cqd=N@A*;MqwN({MKGz7W%o1*};{%Wz0ah)nnl_1e;@qjpoA<(V0i<POcM18UT=
zBJ5}uEjt|V-!dAaHsl`&FpK_3AzDVU*m5&}qJ+N^*hIi=l|p|e@D~DqBk*-J7(?`3
zn25f{u#9FwW5F&U8$v1H@}b(_)$ZjV@Y;nngZaS}rgXm_X}whSqlJk~62#>G&r}U%
zD!e`^ed?`KvrajRF^CLCM{HWQmW#h~tpYF~Klw$Ew&FBxAgd&;+|(H*Hk*2+bhv6N
zEjZcc;bvk<zri3;y<gq8a;7$1eY$eHkGtwH2jF+eGqSNtkwwzf6bSh(g<)v@=Ps`Z
zj}W;uS6sHT0$yzmCv6gp3`i!)q>Rd^)bUlt+A;O~s^UpJIi?d)JVRJOdAzTthBIlS
zYzh5QHT1N}FiDr5R<W8jG3y;1;@e~zfZ-;o$3%!o9CF#@F{+c)Ae42Bl<Ab#uitZ!
zd%Wl2(_Gs8M2=+}$nR0{D+wgBk~;AzRAzN@&lGjn>H;-<_3$_gY!wyZ{v4$cP7}=8
zZL%VbhPBP?>RQUDuJ&43Y+TJkf{{@hgfUB|i9FeOWuww3h=fca+q5D08)hHabNlL>
zv~M%yFR63au8c!oY;rDvj{qR;zr`o9&*rJWXK+6SRJ^ynY)HzckH-poSeR0-nwG{$
zOQX5H;355l>h7!RBuSOLmVl`pd9~k8@f`$~5V)QIak4Di+-<C(gw+7Orzk8<rkO-+
zlb;g!i~zIpa|->Q05Q=>gUzIyi>5iIyaJ>xM&9#1&#!x5TS0O!i60gh$cAYs8)){0
zC6lF<do_p#mdWN}=T7p7SE)aqIU>#?9HROjATR*fRyT~?sDuF!@xsc$u}#RxVdH^p
zKe7l}k3oa;F%em?5tuKSIhY9M6pDtUuz@mydTlfe`!GVL0jh}@S;=}@i%bsCLfc<{
z??u|LdycOkt0hfySxAW7^Bl9cTPT4$n(QNsHgkZLLMDX5OKoDtO(b+S-guz7yKZyT
zZl^3$s$pa|rC7CN`VLTkLhnPtbk(%sv~&K6(7H0W+)P#dhQNme<`AH{AQw}JIwq-5
zN<}<(d7S@&jGKYWbH=S@7P#H+eD_%S4_nWiO*1u|`l4XI^8KJdU9x%H1ddW^l7aq0
zW7un_729pLKSR~Nv-$LLvutn65S|huirY>QfHIY%QAU=kw>K9%Dq)E2-m<tdfz-Pv
z(b^2H7AkfORWUYE5BFgMnwX$YY?-dzq0Tryc|I20#A<BrOp;TG7FcuJmyM9&5_(XF
z3A_yvTB_86tz~OsgmZ`h%@-rXTzMKHla5F@nCwL#a)=URAZ<Xr{%}}Wd-QzukFC?S
zscMQRJeKFI-5h;3p%a}nH7JeS4E4CDbneLpG>M$xovZ*Krh|E`bKRP(P@}d@S~d--
zAdTnW)6|QT6g37~r$oL;jgCE6$UmpI9-CcOs!O)bpTYGwll2%bY^y~WhIXW`rM5PR
zy^JJtKBfM0^{;KChA*V5T*$H6GqE_!nILF|e5|Kp`v9lA7|AI3L$z|}#vHc6^39Yx
z)Xh7~rqYb_*YdxQs&FuzRuc|}rC|^5HqwH=YCUMA@nNt$1&BE7{m;()v9c6l!@UMm
zf45Hxh|JxTV8_wnsU#_5DW+!oO0|pCcHhF8JgZYAIXR+cS9mCMd2Nq|y=XI_*T$5$
zx>WtiH|HCcHcR1@&P(;gQCgqaLnMy_MMhequRbh6=p95_+U~kQO&xR>^^Z1*hWfnv
zs<~u<OIELIDMX98xdJ>!Q6ApuRGBjZ7b`h~vROUNRFBk5TD6ssSW~fnBE@T@MO&JJ
zy1b6Ec!uz*!VAi^l*aqMX9&|?8<Eda(vxbM*f8k{N^?Y@P1A)QfwiL>t{^<HLPB;E
z7VQvuJ1<wLJ4HTQUUF3;ajW{dC@vtGH*#3ZgDE)6!OoLz6AU#_FKI9aNKD0DJ{VNy
zLSW>O0gsacGinayX@s8^3I`%O>}+48L9U^Gucc%dIyJslZ(C)mKvJ7Ll+uSa!%L?O
z%hk1dgSMw9qqbb*ojt^nxb=-poa<Gcf9jerwU=udRu+L!#LRvkwj}p?oX#O4d5Ivl
zm*vU^O1J@G8BrhjXL5z+AtbF$*x7`FZC)G$dxLt&AFY#n2|r;7_3t3tp)!n%kzjcD
zRZ7&2*Y>GP>bm&AI|VtZIbx{qPE1p%-7^BUGnh!06Q<scTxbv`k`s~|gcWW(@_(pK
z1fGl^AiD6gq4!(~TQDSchGS2lNlv9M784RLh1`2<3&)Ev_lvk`cK#+p@1)WhID;}R
zT%OD+onQ*SgCh6XZE8}(u2ttzaVBIxLZF+MX5x04)2x;m&5U(5KtuTZ9}&qo3jRo)
zXsCe47b>F~Ar$-Q<2Kn4p>(c?SSN`Cl8E4<W-EP}WLQ@#O^vmYJ5)G0RJ*9>VDO*X
zy3s_K-N+33nj$#N3J9HzH_+N|6s{&_P^BP+hobNBJ#7x(OQD_fsq<eg|4l~mSydOl
zH-3;h^7TV0&!_zV>rl!Is5EKhlMUhv=!H7kbm2MNJS_$05kd)tE~b)X*<w(6_Rz?d
zx%v&NWcgu9;!66C0}s%v>b&MTE4aGRR58`bcAdF}pGVQ5Xv}{$|9i6vx?nyfE+Yz7
z02od`uF^3Q%PyCNDO(zz9z0ScEr!D<gsG&LQgO0|jdWjYT@WhMwA@)$vflj<VXOdv
z1<d0cq{hBh#+YXALSzylQp-3gs_r8qymVzHMJw{SdP<gQojt#oH)VZmW%r?a<~ytA
zkaS=Gh7q~NL@riu3X^>xM^I=ar8=OOu~_32Hbh5q@-ekCK6~nHDy&J?+xFdI|4T2@
zzlz1EkkXj9)~@PA6%5VW<N2rX(##d}N=>yB`+aRFv`qave)jUS&00584RZ+0B`}Y`
zd;$vyEF`dqKox<-1eOq3N}x#XJ$Kwv9^da$oSL#M2PR4$q10IfMiBwpL|7yCV8rZ-
ztI;PLuzoGm)Rx($Dt%wh?1b2X%Ck^{Clt)&T%p4%#4$=!7{`kJZ!_ns)%%J{$qYa-
z5;RlfNEpF3wqu+d0*zjnI#t(<RXzL0x@jPgcd>eJ-w17^`q#dZYgPl@NZ(*?^yL{;
z*pACn7(|kc4phB>C|YQiOoUhPX_n!v4oB;NLsZDKRLyw}r*EYy&PZ9%lFI!EVX=mH
z;;cSuDu=a{$WxS8k&h|5gur7efB%#?n;0R)jVwC3u&g5a6k#}-c3<l*nj?nW7Ybn~
z2it&b8|fjyrp-1AZ6~mUGpJ}kl|%}UIs+R$i-K;{7-{m==$VEKtk-3ed|JJ{KfjEO
z9U~JPP1wHC04yRpPV3!VvcV}YrAGdv1{^5WcBmN#O0*ia8DU+;4irxjRQrI49`&|G
z12qk%-Nut$?xd_t%4(&oi1e%H4vdQDQUhFqldTfn)QRl`QVb|W%O;x&fpx(u%5Nlc
zE~OBS3{H^sCLdiPd$GJE?E>6|_!?p18R<>t;SiQl-Y~1$gC#=~i>rBf!W#K!>e_=P
zrF7>A&cm~URx8;|jgcCbQ3B?Ys6IGYc?yZ8k!5;gz*y6R=8Yz+-9@VM{0TWU3m}{~
zP{DHLKfhG_py!hF)3lNABc;ZTT>3{sl7rte{30LWh~|QQX*oxu*)g~|z)jJ$<I<k_
z9n&4B4z52XRvN>Bw<sKLga3m~zzwH9ODsXmOkbLBtF?<Vu}m!QfFvbwv@W}F!pAv7
zp04(F7H&ATz$$^STdW|jFJU`qkLA#VgyA+Z4(^;Q#X@#-scJ8e<<*ab+uy2Kxj0Rn
z-lygfVkA7YGUSEo_nn29e_4kPrQy|zSgX1Z<z=iRpF%jx9#Ib;D$Ky^6~A3L#CrAq
zq3POsHT3dRN3ZYDs5O&*w9I8=Wr9nh-I=QSP=UJQ@`|jpI$-n1hIf_4#w+deyfG+g
zzG?VMFSF#U(Q)Y429>6!Z%*ke>bH7@nip6+8;8`5;%sunOG!=b>)<}PR*t$%4Nu<$
zp5{d6sN2=}s$hbSLNS)cC?8%{Vsrb@*7IVOiM%YlxX6}&iVYQovu*7|+e)KpU9M=l
z*rM{T7@K8zXmw`qmTK7*lL~Q|5gSy$Kx|D=NU?DgQx{z^CQj)KX*{BfVi~dF`{9<_
zKBIl6C~ThvKJ%a7**-h!?wS|N5S|9e05Ede=O86R;E7fHT)(Y-9$wYtLfmGsxM@3V
zz{zT#kCtafC&y-q!q_ZG*!<z}l0(~DqgAmi^v@xx;$j;OaBNm=J~}h2byGB3Y!_J&
zl5Eq(cE?gwHL86<Y<6rE6~Rl29QiT4Ctpr5q%gLCn1u#rfi*eu<yaQSP=ehwU*2qM
z4MrzNGif%k3()rVg&==GEGHIU$gZ;vNMLTHB(^X%4{mz7a#d$;b4hFt5RgyJrZkW%
zOJZ|lS!r-wqY-e}7T99*AX5w3T@B^n5tq|^9{8SDznI=j2`@{=;QA$EC%uyrH8d;0
zn<YZT@%+ln&O=nkTdCgc!Sze=ZpxR)U5YnUwTY^hfu+gm8{h^8)t@Xi52+43uQPE5
z*DnYAnG-xk<MmciB#Q=e%miJhiXa>iSBPcQquAhn@w|Rja1EkS*bmr5h+rzKiQW6P
z_C>KpT`NJxR`RWj4enY6_MaA8)UCzMCb~{X<~UL4cVPOi_S-EAF#(1^nJ~_vE_#dx
zEI&Monao6)=xU2EgSyt3KC`jBu2ob$@bc;tcqgelZc7v%T)!3<zN<x(HCpSs&WN6g
zM|sU1c-KG`b*<?o(yDkpUde&Aen)Q|(Lr5jkyC5T0bj|xqYfdvoH3W_&>bCAzd=M0
z?;SQRmL791?smjnBFassq0kekVk1=%awBv$SB7aBn-QBCn;qR0o6C%wfhz5_F1V_-
zsDi7^>vY%~Y*P%tx_~fSOgCdW=)XzWfV7KL?5bI~;^?_*tm8&4agP(X78~B`{V!D<
za#+V?xPSQuwyZov3)tlZ9s>b>{P-W1Hbse75Qu*d(8~Htm47<_`=Cl*Ntr2zue^%z
z{uj#^)yu1a2kHQ=O|{d!ZtQCd%flwxGCmPHqH+%ptG@<G@+g680W6x{tCH7I?8(1Q
z*nZg)A6Ci#FJ|+lhvzM^82^n=O`b|}e<N6UP*ay3p4&To-?)}{leM6d@1tucbg&9x
zHsiafujU}UQI$ADv^Ez%>K&gat5$%;v~vg;11w_`w`1O#AdV|uK{CuVwwap3a6gfO
zu;lYY`As$nvK!^G5kH2HxxKlbLd(>bNAf3-b&ZyjHa?$-q!{}02ZT&^etA*iq6u3P
z?i}w%xGe3C)MNYd`|JZe8&Nj!>6c%Ai6$-v28>@iO>FXsJsD$*YtRdLVaS`P1==h~
zI`g+S3i3FGzjmLXE<IeJZaBTDf(G5(Y?^*T|MOj-&#Bgz4pYX0p*3(<g^&Djo8NxO
z-t2@H(2@eml4Q+YKHo~TlDt;*ITLeN-U3?DMPjBri^7_`6(7}n-cU9A=)5>@TsE60
z1iZoiF?Eh@Y~FM~MA3%{JVJm~5N}y{^YSrO%cpoGl;-hiJB7|6<Ruh(jKJ*_y@SA=
z1nwemJ!Qm01b;~J#O6lcP0=0#{q{5$Q5x-9%+pvtJ3F7E{Wdlps+6z%Z=>9`h~iWM
z=UBL+mwh%yahtq{P<a0X4>_EjN95a7{C-Nlm%x1lo~E=1C`6kf`6CMPp65}7dZz%n
z{CIIo3{M^oQhkWgYX*3$5O|!x69k?lpwmH`{0RW<xZs=-m4Z@8$;)R5>?H6k0Y8D~
z2*3uk$y_4i1qyXjsPBdd8y;-9qqtVvh4Zqqs{ZFCtEdR=k>bfsl6(<y9C^zYoS(i*
z6^99skkR2aAFGn8l@S6zqZ(eP&>IAPO5jb-AV3==^9bwb6r~N4kwbeW+R8>SyRdN$
z(K&D&H<;Fji4bQ|xf%kmQFJeXn<z?WoP025<e1vVaw0HyT)!Z^8we~Syay@t76F66
z+XP5&$ae@F0Kk<D&ImwIT>g^M-X-uW%5d`VS@)(ih;$*^B6|<ujUjYBL1*LebC1ff
z*F<sduG2%tFiwJZ8~M>hywz=VAnVuZJRLWzNw=GSgGj^870p9U=JxVEBBAeOHU^n}
z6|{2<H#K7OeF3R-I@Ge+Rm1luX?Ath_pi5Wm1<tsnlm*gZjT0QnfMOVCTMQQ2sn2*
z;e|ffF<2YLp@F#J&xL2Y`O5{~FocJp+(7&d(~3D%itrFEQ{B-uGX5_##!gFU(+0(2
zOT%v&ei1qrF@5o|ORysZk9?su*rB`tyY20c7%Zg$u%n^S_*-a;*(n5z+RlFXBu{jT
zL~6yz<QD>PSHM9L=?^RnFWcsTE2YRarHMG#k4f@VpvFn#xsho@?}wnr@>BIf*M@3I
z_<I51=U+|X76PpVA_NW**hQd?KtJ{mB*%MiV@og^FnxLsCY#9yqi|LCZ*^|>p{!;q
z{BO`9_o|V{%JOhqNDj3ja!ZFdgD$2S>W*W>hfMoQMx~ii`IU_6l#&0{nr3i;g)?xu
zNM1(MtZ$|Doo;<+^zTb1J1KFQOVpI>3+=gfb^7(Q&iH~zXAU`-Lst2S$(nu)8N`)f
z>fp@8<j4`?2!4>qY`Yw2?|J3=LmC`ieLooC9*tBiH+08bq>lYyI{TC`U09IVekI)k
zjB?@}H&u=VcK`XY(r#QEr|Xa8l4GP{I>1#kxi}5PVYtjCu$D0ChFVUb%C;edTciI_
z_)7vdLedD7A&PS!?t`32d6NiCCQwZk#0c=xI&umnmJlH7dMPC$jeM1JW4_j{>^IG~
z53%=D-SnB+?ShOGl_oCrdt1LQ$iYNHvI9xlnQZP4(Qf+g4f@&+$ULg&>!>$e(6YFw
zM=!82O6~>on7U$Qnb-5|8Q)(@il|qr>6Ucu3DtVb5_`T~{rHwH^Z1@u>M#*U6D`HX
zHiM!U5g<bf*A2|ZCQ7(~z=Z^o6R)pKvB%{{Bt<lOHpdWBMD~36rEw0u4QU*KhZ$Zk
z4(X&$?~~Z<m)Fp+994h0wM>=VHo<O;R_kvYsdnC$rW$YaX*a4@Zu?LxQ}^C}&vsHS
zau}hrX+r&gt5gvsj37WZu{q$p^7N(8a)&V=Sye#ynHTj?G7~g8#;#tvE5GN!-8ovk
zkZQ>%P(lC|lLT_xY%fix)ZbEj4gy;!%C-wzPHbMVxL-pVOlUFTk<kk$P98a@KL}r6
zCtt|1fbw#AGZmOoN|5RQoJwdiP6^x>l2^GNU{24BdnQ_xGLb~?Hz(3?p=Ekbu$L<Q
zdOFS57^NhkbcJoQHc9HJ1eN`yk7-=0{&jD`HlA({?t_uj7{=SjNRwU@fX_1{^-f3b
zBX!RhYOa)E5I=Ua9>vX~@p2>IhQ%*8)14tF*n9N*wrcSkiMQ1K6|}z^4_fi@A)eSs
z`1cOEW99pMuBZa{`x(JIfG1ZHH~`SQN|Pgk)vh_i2Otxgu7R4Os#ChFm8*&ehH2~6
z$_L8q6YV|nfy=cx?^^pxP}1(FT}+?HYAh@)Kx>HzH;frJM%X^ZEi9W&=TX~#A;40?
zrcV}O_WBdB82y<LrsO3^WPv$>#A<cxgX6UmJ#RdChIR&%mjp~Ga1uA!Eage6uBGa>
z5ST=gKqfgHh7wd6s`;V9c}yA`l|-7%AwVRVPMBn%^<E2;ot47Ca`p5>#agTS?L%Yj
z)9fnu;g|WzmU+EKRLS80N%nmgY01k<jdzuM#;ehDho~`+42v_JOQ_f-1g<B*Q<r<O
zjH2BHP9?y@e<_7l0F-4I8NIxZuv8fCTD*W{f>#k1qcWpuIE+l(Z}r}7aj&P?ph3;P
zm8-`e$<r=TZ$2{GKGUwgcx0DWqQZ|BvR`vzXU>NlkSwM!L^NpTZFFDN#f(2nQy>#K
z<lpJkDmAB_(-W0hf#<6HMF)Sk{((_%R#_=;Rar+c?^O9`neQdv??w8sH}BJ`#60oi
z%mfD~zppwm+9c_2n|kZfIq`o|%x?N09R~vtW+}45tT(OAG|QBgP67{CorXXY3NTVi
zM(^&ME1FBQ*Up!!yfsG{&Ghcj<7Uei5*M&qD^4Ow4yRU>N)e}zG}4!CShc2N=Ct0x
zTq`iEH)U3(;U*k!MaAsi(6soBL~v#zG^aOISz%<Mh}W}ZIZ91SiA+n4R92?sRL)F|
z%rbHU)j=;_@8a=~ciN1;(OJf(KG}1PwAC9bfW@Dg=4Y1qna!WdY38TW{7m2DTDSG=
zDHXHiY_Qo#_pJ3!n=^N|p?>pt>GUF+*`o-|C0=%=j4m_!*GT9~&m1-WiI3KQlLN>s
z7X%ko!vUNvLX0X;B-*8#k2&n}MQZg00QK6_!!G|x!8~3!G<heIfS<52Qxyjl*#K}(
zvVKgTyVQoCT%_Hge)*G5`%q1VpBgrj{O2OkCfw7*2Me#N>GqJ8B13C3)Gbe){cT*N
zJJh77i}fNR%gINjut4j4aF8&v{4yYDaq~9pW;X#4^D*&J<I7v?Snl)QNr^lEX~a}K
zq{509s-gE$88s5`&T#S;LZO4|uBScgwt#Vnd=ox*CHb%?&AYc9F4>Nh>r-MkxI0&T
ziz>U5Rie8qiQ{zEkaTc6F$ifV$lyDXOSwK0p^wnej%O#RQ=WBCO(-(Z!lTF=H8O%u
zHk7xjjnB?q@Rc%*&IJ=$$uqLsk~GSpcJ<t|CF884MLIs*<Jj9i77QeVmQSV)*GBfX
zmK)We=hnFXMd^oB&2yu1%6{;<G4@S%b<1;C#SPb*un*TfFoK;#ET1^@s%EEXE`ehx
zYpK}dReMUz?Ik}xW}oKBR3H4It&KkpkII`c#98Er2e$_7yILXGcT@8;8$crIOng|j
zR;I0h-za60AI7?}w2d6t$Z0G3N9|mQPIsyBiZTtsjVsE&klWNd&*y6s)L)-3iZ4Ny
z;c9BDi!?T}9Yk*AzUr_H$vU(FC&#vke3Y7W6X3l7EjRG?<F|C-CcHs}WhApGF@pp(
z3WtM;xmnAh#60+{Q~!d2>hj*%z$b>rphl%$7?e#iVmRv~c${=t-St92Jg-TH;o}~`
zQ<6}d<vDCV)*lE|Nzx?;Q=v?_T-M5PW3)BOvkD&n@HF68<3t!e^56ru+h|}C#uT4N
z@{>yQqzAhxcAGaJ<^}esq8HQS|3!_^Ih~O~jsSF;6ERN-GbnKY?Nm6`Yx*^eBmD6M
zKB2~@QmBH!I|L+wzZ3X^z&{8~CA<m(l>}xI*h63+KwMr#!A=726QFY%J_|7qB6w%X
z`^h0Rx$+4#Q?2YPJ)NS7Q(bx{?IWLNoN2*_au8}VoDF(g1ZJ4|e=|hn0F*F>nsQLy
z5)G*7FO7^FP7DN|J(B~s;i&cryzr^hYZ}b60n-{{TbQ16Mu~~0g9=XuHlt3qLqDI$
zr2VRq1xFw}t>O(UTaa`-*Vtq@n_B`=IVMiMm_T46ffK}mX%w1CfXq%DT$0ydguU7g
zcT~r@9gH>8E@eYA-Sqtmnuec`s^D{1w)yTQA}kDPT7@+0jO<XDoYv#;QL4dyl%_Qv
zg+n8ouI$vp*+fDmArEW{<Fx``kpN0a$$OM+o$Z;A8hK%nE~+P+9DWi{klrO!%PEDP
zAN`F`2J<;%#2d#6Bbk|AdMHEZoxLtO5Zn5QEGGt}!#l{^rbmt3zG_AO82IL~e{zEo
zq_l*IE_q)-3e?Cu$*HH28@K4N$l#S&8o-yY5O|fqYXoRLmX8p4oxoiL-T*Mt>d^l-
zyhDZn;%R;$)YKA%qmc|(!;zQVW!edqjqq(0+DqVVD(R~A)x_~k)_k2AqR!)F82=rF
zAJZESWUa<<QvuZkTj*Y@ye@$MLqgWnwc*{XnM9B4d(-T7H5lkXr~&7-ayZq7SJ1c&
z@Zv07$E*6FAiEvz+x!;a(9WW*04@(!5J7{9ZCfdH4FPgI<e5WCEy5?|h8^4H?S#gg
zU!6u~kVOWOY+gKI<6x!X;1jrfzrh>(yQxYiUtjR5HjXGvk4qekREKx-A-IvpXO8SD
zMuay7;m8Fu16LmKcK-v__k9AF5fzsbc$R=5@E8HUkf`Ad!s{lWQS?y&L(|QJv9U;o
z8i$_n3YGO*<)&klkS2UlOawt?XiqH5sE66+DO)RhkxPlBl+~Np<Z>$PUI4SIC8#e#
zJGRITKyzzq?`z{%YWa>_94k6;)G|`bwin|+peY8NXwP(Hr{%`M{{`*>_W-=H&CYGg
zF3e4Hr@7O!H@Jtp^V}ubaobq@+T1Sn*=zps`R*a^VeTQ>1Kk;!j(nUF59X8dTt3Df
z0yr5DoLvs3zwVofeNrY`%_WGxT>CKWn(+J?X-;6}Yi>HAO2Y{|VP@dGeW}Z?(%-mE
YD^~Zsv3{_5yq>L<W@Xw3Q}J{D8+G=1yZ`_I

delta 23235
zcmbtc2Vk7VmEM1MC9Rre%T=x`TP`xjMeg8=3$C(_Z45uQ*DK90X;pUR*_ACBnFVY@
zF#?&N;7~%OB(y+)AQDI*^xk41bs-lZgq8v%m%B@HN$&gJzpGuz_~3GeegBj<GjHB|
z^X5%c_pP3*Z}$v1qp&c~g@0N4mFl(~M+b~B<{h2<i>a2keyl1`g=#>D+Z>0OUlpl=
z9fsjjgUnKEyfwipQ^ki2HP~Ke3@ND1Kjc1Sm=l#tjkWg}L*h5(o0BYW*JNu_{giIQ
zEJt{%Dpy0)P&G^q2U>|5sYa>MYHV#j@=rnjG*#SfoD8SNZJsJPlLV)4En{nQkw0D9
z7^9}DQ&fDKIu(ud{wIyh$Uyjvd}ZLj@=H7=E^}t5p$4hxYQ|y1oOPy4%~Z1vyVR^h
zZnHvFsM&{I=4>@b&Bb?)ny2RDJ6A1G3-O%??9<fgaiE_L^!aL$T8!icYKdBk??Saq
zEywpXb%t7j@9Ao#T7~Z-G+3-wt25EyVpOq4twjxMQNt3oPOS&VQnf*C#CI7m2dgdW
zEMP8ITh%tCouRg?v+-S_cBpgkUD+8o!0L0A38YmhFjSqVb|So5?NYmeaHa~VJ@~Fs
zK~;tCTHp**)k*<p9m2zvrD_mfuWD5t@;6u;Rj9$;nWbZ@UNxvj)uh6;S;_@oRa%v*
z0eK$1z#7@L$=a+UsyTkh4XJJkk3?p@WS*7G+^d?MJnJlLOV`$J7u053t;^gFVHu5Q
zH|=n{3|G@R=DCIosV08ssUmZywUdo!mQlJEB|DAQMr+$jSEsS7%G_n`K7`WNZWWIn
zf~s1(tzGinV(lt%sTh9cpg3ST@s?zKkJ@Jjtz4mfe^-^Y%c``3)*dU+ouyiZ*y?0%
z+h^vgWbT1<Zakx~mCQarot=?eV^ve}E^(`N!WuUSthUq%Y09!rNVBY(6Vhs|^;WeN
z2fBr->uV_^HEyI!9qh>3={D=kkhQ(59(@eeKt$kCU6LCYEOJ>P6>o5Ldi8%<CUKY_
z&a=)-g*!cayLOPP+l8x++_UsDh_(T>Q1TaoUv9${m!$1Tx+s+to}Sc&EFJDGus4<4
zi+^!%N+V^d4ynUXz$Pf=2z8125|lEGaEZE9U50Q3;gRZcbp^uB>WI1$xO-KHieH6*
zR#&S|e4{`brH-m=fD}`<x)y0INsTVn6(`l`Y&`^IB{e#!v`+k@q!MorPfd0zTQ5oC
zih{?f{mF!jQ^@Me>bhFDxi3?T!OU%BrrM{DsjsN(F$!*|&6a@`R5zOY)lKTF>Spyd
z)up=C*VQ-FHxau9zb<tva&ANDX5@SWA?AEl-LCFX-%{~A)m=A@He1!V%{Hs87DK7l
zrS87TYaW0NJYU@d^4cXd8KHw1WEbEU$L~V?F2e6(b+36yF#OgbYKnPSV&6gN5)8;M
z0q0WuF6+#19W4Vee^)UFpfloF)#aVOwpNVA`$Sos;nm}efo^wlU|oUy`!n*t+b{n}
za_oGsU;35SYU>L1fO_z-%erzUY>0J4Rbqzvz9@Nz8NbSkS{>F^huk|o=GA7WdDObZ
z>YVQ~FKIp7y2iTtkg?LWYaPaeZC+~~w(Km767=uO_+6*QT36IxP7%}}vkrB8kTy+?
z%}gEc>ULY#c6;W#Jg!Bq&YZn3U{u7cqqH3A2OS>he`n7Ak%+c7<aGM<0*R60hfYfI
zRn}EjtC%S?T;c+0Uoo#wri{f9ei&BNr|-6I5W;Q*Q8!s%!CXkfziNF2Gx69Fu2C)v
z8ibViF7-pw5pK??;t{Fh6I6i;>p!oGR1MRCKkkd;v#ZJ<H*W?Lzh>=Hk6JgY$2vSa
zy=Ir$ZGPSS1}OigdR+KWX5CCd95SF(w?OlsXe?~H)w~r(=~n9&$P_>GHfTVDb&GXd
zjT_eF$zDBp$}u$RX~%ql&-oaxDT5z()GduJb{tH-#k@^bnzy$;Y~9|O7rmf&!gxlM
z@DA%+hd3eK3AKC{P2JgTNcwK|9AdXwcQ$yt?&@|WQz-&TbHrdCNn^y7tHIlutDgcz
z=out(q{P+r?QT~_>~4uc6d8G`81fEz!H3=Ic@m&*Bw-D1^Y+f%y-(2us~3bYq^cVs
z);#?j(q43kwT{%sYbAJ3vOo7?cwmyT?iKIm2#7jT?_k`SwEj=S`jgQ|z~uBGOfngw
z{0?e(slzbu6Ri7y_3{a#*lhMBhwCeGWcF}4bbFE`JvEYlbVAO3GV|Zx^j(?xzh^#x
ze&1)^r;eAn#<<Yn1N9F&F{FCyzi-`_oYcmuS359Uj&P~h)a&q4A5=d!e*o-<tOr&6
zC+c39gP(SI&4;a=t{<wOInAoqx;--M<#avLO;_P8sN<vR=c10OQU4gdg!;!F#e_CI
zVdd06sm7zeU#Q8n1pKcyD`cpIdZWV&t$NDDuvBk0x|+tB;~*|qT>Y{!tLbU;X~Wg&
zZoR>}M*W-m6};Ua*Ltb(o$g(SQS2G>Su}rsZ;#$mr8OYWd=9;M-g?k_R=us>>F}5@
zpx-ZAFQD09ceqjRCG%y7@TH^(ccYG1te2BY`K)@kA-gkPs4HQDuWX$ONBt%1<wI`R
zyYVRcBlEa*CFXPL?hDqB8erL}p5#eq!QOjCwSFVILv6Sca!#Axls3r3EC%?!KIZIw
z;qkLD?nQAh>w}JLG}vKYWi5y0`!EA_m&!l|=0_db5XIA7uY#=Wz{s)cw-l54T5_H7
zTUa-&JJdu~s(z=&sR<Cn?}6oqP|DOKc+P)NAJ-0&CCKZl+_Y8v6KU?p=>5yqOU$4&
zIR=s*{!8l9&%t=P)y&}Ch9vKft3P637wV5JTaCqhVkgP|Nneh}$x`bDM^i^J{G7vX
zOc#ZFKLVR?u%1z4hlBO?Kee7ye^URh{!D}2T4p@~QH?f#mVxne>u2ge)PHusvFaSK
z_ngeqUB9Rm;ia*=5=e<dp6R6eH`ISYHupfG-Wnq2R^6liBGnb_UBv3B10dwS3_|=u
z$OsZ5xqXRvD^)5@0I|9f))VA%++ni}HK+KSPKJ}N{^~UERNeJU<YCAuOsgbYn)^5P
zH}&61tepB^slTg#aP^77{FV${4(fk8gBtFgtJ~A{HogXaa6p~@=XAY8)E%6mM|xe>
zuN(3@i}h5)WxmT{d7uh|=w0eM3VF=mfc4vj_19Y0k_2lP3~P}-ju_0c`J2w7y~ok(
z-N42m1}VEOXQccFUUOFGKs{RWvY03HKU`-bXFE-_YB(MflI9jmk7Njo5=T9;f~ka(
zJ%isKEK#0G=3h|a!gwC28#@Q?U1J?UeK}xrz4eUn_JBf7i0w(nP#a<ar~fhsqJR%<
zf-j#m#P7i_ysyTZAGBWBwg<Vn%uTI$5YI#WgKncc+x*ZP)itWyWqxG6@7Q(gL#Gbt
z{rl3jeBwG-)<+{CKH<v;YHS=RAM9|2M}Y0W1=|Pd3In1rT$^3zy29>rT<FaYAQ(dx
zWH#TOl_^E@chX3qR65^!Uxx3uIfYs8Ax4$`y%4_N`W<Q-(2-;Q!MYs%|AX|O41>HX
ztjnzfSV&(XGuIWcXlW-jWznwa`q<)9q6jvkQj{>2GO$ld7}IS!Kdi>~$z7m|Qe_*g
zg-&D#R<+77-#zc$b5cErGBS0mYszTD<@%7DFO>;no>{bM<K{r+_KksUYtCA=eC1Yi
z$fnJ!D$m%kd~2W7qK)e}^hvM}77sT}`^&{uJ<|rCZP-VKoR?jMN_C0-yCJiOn7R8y
zvARH2C>Tw6nxmmRT-tBX9=f$?DsZA+HgA*}dVJ59ht4o2l(FDM0)H<kMJVA51VZ6Z
zED)G!zd!7_QDi?j{ItT+%*i@n=}1(U+5a5AZqyW{l;>#Hoydu{RMlJ6v4l4o+-Jcm
zkC;2ej|``3)~aU%%(K5eB4{kKi%N#Y%VOSkSG}0)UBeE#V_Ar0w_}z+WV9RY?nSP>
zV|TdP!Mk#!63wuET|rPzFc~0`rGw#GD={D%3&yNKxTPsj(-N-M>{z)+2s86~^P&kO
zl+n}4($$QjDin>aK*^{Nz-1I01C3$EfpPZc#lw5vEh%vu#r8i&tzJlWm%H^+e1wQ{
zcfvh=re4luL)uLD`LB748c}O?q^Ws{X+Jt)tbOC?swsY_hL!jvJfW~EHxk~)P*XHr
zPOn1VdG@Dc{C4PA$uucoXyLb8vz%T{!1R8FcB0VAam)aX(t^==0-jwurp*3rS=q!4
zJiVs3GNx(~8Ohjfe|6$CvxI2|f%Ja>btD{YvI2q66iluKVOSbc5n%|o64^}=Wi>_<
z*}>*!E36WprchW<vELdy)R=1j$Jo`za(nT(dGpyds8P&nPGq$-2VxPSAKf{V%~RWi
zxJ040Vl6rxsICi!!qK%o-x~LhXZDHGarpI_;t@7q1i4%|*^g?I924F`l+qm2)4aI&
zUlGY=6v-E2_2pj}^9l+jKSL<?k0%C=f%d#f)6e*ffuh&38)XEdx9b_&K(LX3ilr+F
zHW6&IZ=G}wtjs?rUEnq+*ulxe?dv9wiJMvBfK^jdKFpE9W}<E(IE!GblqA?ra5e#%
zu6GceLvSvENpK#4sPImPb^#=enCaDZk*HbF7^x06CdF!c&JM&d&g&vm0<n`SLZB`f
z3k1U|U`3<4$DA_Wh<jM^HddUYt!Q1a*-B(NpKKV7221HGW){{&8nFQyNTzw2=4o!O
zfdYnL>uOj_2%yi<_A^xmNHsz|_&GMiGt3y~9%~$!x^ad5xAI@a2NF*-($F&zF3)$^
zF9uBPX>J|NL!?#*n{<@O4Fq!l#Gp1eMqp5hqUSNq8?B4%FV7PGCvvMJE#X+8DcCA}
zPGn0QMlK+10vMX4GLg^T1*ECy!}g}>qwGVcymBIS1^uYg`$<h-x^yek+6ZC<q)Q(F
zNRp+`XMza6ouPvS{$4<)C5naT101F&$$D1E`pdMR{J+Vf%>N~_vV}qJ=xaTjrr&MM
zBhipLhYIGBmH3EEXwMVgP|RxT*DWz7AK0^JuEvCUU}j}`dK_}#B?`g;6^sR&BcX6C
zn()R#u#Yk<EA3Bbjy9Is17=Ogmhm^oo<FNBzQAb_vuPlt92Vz0lJ|;U<Y49qRa;TB
zs5%mkf+{Oq-3E1t=}S=3%w+=sPR5QjsB?*InGB&LQ8TLs<&#s6A9;>VP1lu!%5*gN
zOZ)9vTaAG|b1S@AbFU>HCHmpF7%~W`-kOS@$ggc|3hcK+wRN#5r?g;eD5^)e>@CLz
z+aF&rXvt-m%u8HH-Fv-TGgKmFbg42Z<x`<1ZK;-Ohr84lVeB~j`23Bd3K<>lVCK|>
zYB@$x$f(wX>;(%R$R5W0V*Aevn&*hI6PreMByt@$06o%&iL{qM6NniUS$)aAe&KAR
z$$oj^@How6B0pLkYHpj(qEUT4k!5IRqvaNDK^>}tG0bo>CKB1Ljln9bF_GQoe7&uh
z(;Y>;l11keh<rO367hOtZOxWm$T<BOtq5q1?dsF!&6TO{G{h1fa3$eut_}wFS=F$R
zbX6F{+!%FC?JS}hLEUA)a#~rOv?h!L=72~|O{h8qYf56G8Syr8IAp9e0@)d-ltIFS
z-X*fb0Ug;NO}Il%W=?fgPzPXU{K&?L_zLkUm=se%*WVG1b^*8ypW*G<dwR&o9*meC
zVZXI#!L$(wLQZgYVkvoF!(4wa=oPeOuUx!pxERe$pHIZO)b3h5a<WWbROj+s84qT5
zEZ6|=FHzvgnBGNDoA%p_$BkLU1~?%l+|eL@Dziy@%#sn9=H@Lq%P~4s@A4wOmgsqb
z0KAo!MmRW$f<RzzORzEN?wrCZXA_)3u$!QnfOCPqfFMqAA%RWsD8VxX{=PsuJppOP
z4BNM~)R@pSb?Lc=blw?U#P4O8eZw+;5mVCsNZfvY+5EUn-<L9jYMk&@g~Boy=r1v?
ziGYSCk=+sw?QO9{<H0+bg<fPvQKTwr>3x_69s8i4U~1AeuF=8jSSS*9W*JY^+KX|4
zc{YW6%{VsT2LPq4PQOGTeR~|CkbU#=Q;o;%mzMu%gNW78$A}=4qk#~Li4l=-mNO`&
zBi<$}7#8|5^X=8K@;sT)lm5T`{uwit)UhT}U_V0wq-Bnzi&0~gCc$Oo83XVeYfJ?k
zVffr$#I>@wt++nkpC(LXVhvtLIfr>aW~c!G^S6wcpD^L)1ZxO{kurXN%IMDsenD_D
z>Jme&eHfF$rihN&jVtqpA0?#)DB=QL)A1WY`Vw}>%3q6WV$=!Ggw_le>Apm!=<7$7
zkeen=2>;*d8>$_TO6#Eg@v3~|2K#TTMj3DP3}3z0C=!cgOm|F13GMr=X?~CnBB$_^
zT%nw7fPdG1>CCyt0NcH0XHfd%arpfT$|edhm(izbZVrWOlQkubKBE?BG9x@pp1Bdr
zc9oaNG@bBrh{7sad#Ri{rth{NUE?>dvfo<cKZUgENdz)fL{nvKa6U-naq*bc{>$uB
z)=q<AIA`t5^<sn9v()<p9}w3^^~Z>a>fbQ!DuQ>I;ZZGGf5fQ27l<~Vg6M5MudmHE
zX7uJHoeF&sQ7`VrmwYM2#n+RzZ~daU=wk&@g#ODIdXIp^+Ua^$5VC4>y2iOnG&kjw
zEoTC?)Xc4Egr`g&Mt_!k{s5fE4DdA%0hciZ{|$2w?0ItiZN~pa#4p;{Z9LOof-3Yh
zg69FC9WUThd&BkocgIp)WdD8R$c(Ou7_bh8dzD??(iClJlAbwTdYRSqHD{U%qc0~A
z(QaU<l3)|T9D=I}$g}bSnJV-WCM*Up2a2%QIf~-=NivRqM2zDE!nIc!dX0emN@T%l
zvYeSP!`RG5)^j~?RlaSAHPP&OvcL3B8p;PcT@jrVVnqf2$ZCab5wP@8WLv?evf<vr
zxkwbmbf^XUG2-kw25Aoqb`yLGknEc|0EnZ5jnM#H3Ty#|%|dY#=n>LS)aCZhEhXzr
z86L4nEZAt~Hd-~Y2t357i3uM?qKW);8E0MLEa&`eOh73^Hv_!Y^XQg|M%rE1#mtee
z7dqV)?UbITJISs&@}l?>d6;<uY!ZWMEVINf6G2PYe2&<zW0s@9aJ-u{+}a6!1=N4w
zeGLV(>@RIy8~-IjmE}H}-SxAq?>PdFOii_SrgU~pvrtAwenzF+Zz6;fOJ{iajJy(`
z&o{z1Q8T?7>GF?X{#M`@^#Nd8%T*sYoP8{JB0q-JShCZJtlD5xQ;<s1U#Sv#k|i5j
z!(itAgqmFXp8%iJQmP^S-@r{|1sj{|&_G=OgXR6bpp?0Jk}GYLpbbsJHTLEEw&4M)
z^rip-_NEKloZW24UFj@u-HthGKtA;Tvs#K5HwJaB1-DhkoisVRw@Dm^oTN!k-)s^;
zV-$8$aVYD?u0I?1RAv|xw(N1JpdoS0Uj`|0r|aysUIs9(w{6Dyu(Y?kROnRO`8ct-
z1bdT7DLPhOEK5cs={6EMjg8XB@yT|)F6}jXPTe-w$QH+^$nLzY#J*_zq^T11qp6c{
z&RsL>RHqXAN88uUa7qJR8arv3is`*&#(65>hS#6Htb(mL#gYqJqb*_7V4BUM+;J}F
z<fi@0x%Qi9&o++P19!|=pjlW9m5kxErKPLo>>T7El*q;&A1+A9&JrCGd;gB|4UNR%
z%F*#6a~=Nz`+o^fG}4&Xm_D>nSc!pLrv_>x5tZB#nrQ!S$82nM`Ok?=6aly$8!?%+
z`+I>%nxlKUeg8RS3)eEE->^?3OGru#0uz>jb>dn046G{$oqNiPVx)pJ>HR&d=`)HZ
z1}66-gVAd2Zo+?Am8!?()`RW%xu?&Odc5g+65NYYktQr9>so4R8nKwwv-;LwYX9Nf
zvBTG~s%+%A+<8WUp@)H>)%uQ}31*SkHyp_*c!Pb~?rnt}FsVJKWA=@^%du_!^ltwu
zjgu;u9>QGO=vw^n@2#cD!11LCBa$$xOy35~3HyUuVL06gjD#B-)00p(VPH))Gf-yi
zv$qEp&*k87e0V=XC(*+J;?Bh)0W6F|R?X%1IeW1299ZxLOFN~o0~}yISE95v<R-S8
zLwd5Bg7py%GLPA>?imVy?fpH+7M3wv=2`KAK44T_fjbxy0S;kdZxnaz^ic9drl~3R
z1Hn_)tS8c?1ft_~4-;9jmgYuFU%^xv1L9tYXK^{x<`djYOm|IGe}}QV?5R~-x8BJ#
zPc$4gbCeZ@Lu)0n8?9RChVCHN0s?NG$p%1q3HF3w9}~#~i>QS?2fEt&^j^FL_VFse
zG1>mKYH-$HkZRhw)#Kx4k;aJI7lh-CnH~;+S;!R|_KlpVNKcoetJuMfgrZi$+Z2p8
z=xOZysZ54LTODi-v{hv2bt&pki2QUg8iP}Ty~r)p!YX9Apo3nc*gMc&VA$WP-eGL*
zajWUZEI+cK3S{FKa&PaogVqdVqkYsGHdn+T;!8_G<h7V|PC=Ez{_+8OEfYGKu-Sgo
zn!8~m<5H3<eX}SMZfpzGM0B9h3fIQ!^fu-%Cq7+>{(Y1shmOEpH%9h<)817x**Mp}
zq^4_%zb~>xOPRJE^*sz$$aHyNjQvSyh#jomxJ;-O#SvnD3tSitp;i*4hIBP&qnby^
zQT)KsWJT_}Z^U;ZM>YW9uSlizl|(NSo&%v~J)Y6O5s1%3txT=w$pNRsX|5b)S~~!i
zMlDs&1T15_A1m|@wEMW*WtWBatuaZ15ON+u=!f`l6hpb!VYWzs?Ip4sU@&U?Zzk$<
zQ1BW1ozToiR-suVlmnN2UG&44w3)!EhaqxAh>)UAD}FTK&{@i+lB1)!sV4dmc}JH7
zdroej$Cu^v4l#7F=Ysmb8kMr+<WA%@2AiYU1`R1o+<wsBPfCYV^@B8_$<G9~C31cp
zjcAwsOyi!@$&3@re-SDF#11xnH@=I#JGp4}Zsz~1qSXPGrosL!q5d4bu)iF>Y_~LD
zL@b~?Gl8yRNv???WEp0Zba{3^A#2em8A?)#pXCLVGDKPb(mp@3V6{|7_Dh{I**N<<
zr!wu&0JYoY-%pmpMoed787biWl<?MSs|}JUv;Ptqm9aRK@gt(5`XoF-r5v+9m*xKm
zkjM$P)-}S&I*UFru41X*CdNVltdm5*>)6^!m5z*XShD#PBC=&{99BtRQ#w-0#Qd~w
zMPIQWX)ZVFdOm2rEq|6DxxoCwR`8GZW9@4eP(2`lLKIGJaqAlxrq7@UF*KO<dJ<XK
z-NPwrR+u}q`VM>9!TD+g3manO6Yl+y+LvGAznRNWVL8i_xk+#NLoB^fQNqhx0?0>B
z(%5*N`CQs6PBlo81lxdDX1{T8`^q^^L-SZ^3BgE$Q3Rt2#t@7pz}+|~2tg^qc!CK8
z!|d1vla`Cj>37SzG3jhLP6|oV!wH6xP-AMjH{ou?0v1U)L=Nf(?75fOn~x2_xx(0D
zKS?D3kt-SAX^)E!JBul4Pj(HwdNH)wDc8rl%yE78oU;s^<MtYv!B&MdDk{>&1@dkU
zH3e|R;_VwI+IPh#`Y0ph?XxRy92uI8A_sCO`l};a^8!TvW*w;tR)rcvu{KQ=O5|*H
z_KWlk#Bg1}%VoE|$!T#0KF;Bwp?kxtof8#rRU}pis$}o4!rpmd!@32mVw2;u$e>8i
z53*{e^OEpluOQ~w6Y+v%RLjKlB3pC(A$!n8r^c^ld=Jaya|0~tY3aL}{63R0s`hbI
zC49kf7@G!g?CC%0g#_Y0oyO4V1d9k{n&@Z6Xd|I-5qSE)XV8Zlqs_r;%l_d-xo6zN
zq=K3dR!mh^Y`<lN6W&M|8y0l05_#CofIl4z!3TSt3Hi<f5@P@iuVIsK+rEp-j7|2b
z7mqj2us0*T(r&rfzhMP36Gd>O0&THSb%PwsM8#O>rOe7>Rx7ijdNmWWT6t!b7{qjO
zE(1J@#vzN&wcoyY=(!ze-`Qg7=etdFEQ=Twao5(9F7fRsFq}+S%|YIV2w>t$yCOJ`
z3pT;IH0@1?Du(t>)6UkM*z$+%>kf@CJCFVL1;HTBd%BX9Zz9-2;7qRen};gKP-zqS
z&Y~VJq1EW@-iX1TWuJ0*N}(SZP~CPGEVcI>E;Byqx#Vz`G5SHIRQr(2f3&Jpa#{7)
z74_VD$wgU(^|69>T#?3&*6dw7d*&RO?HP=Pxxy8#dIO7cZ?}O<&?D@-t{Sy<d5NpN
zti6zz+uMEW3|{x+^0~c`PhN(*)FfQwU7;595~NysaC-aT`q8-Bxw3t_TBTNZcnw*0
zOU)zIXxz-p)o0ryuO4NLx944bbk&)BGSxn@eJmH}f|WmtH$K<#YAG&F*WudG@UHUq
z$?c<^cQ2hrcx231e|r09wO(zgbt{kBXcu(OhUdAgb5+U44ucdq6oWda{&Ka>zPWRY
zQE7kDIR=;mj*iN?IW`IX*<{Z=Iy;`xN0dqS20dBi)T1ULb+g*So0l1>nSCtogPW1l
znYH&a^dZ$da`!B7w=kL=^L6Dq>GnHE#~7zwbIs85tp|s;Uff>Mi=K~XJq7x)_MvLj
zeAmIDZDp~nuIyNj+Ga1gW@7$VoV$6Q1^buT;cHGQU517R)i1)UvJ|n5O-je?=8oa<
zSt2sxEEe$$u3w_g<_k3GUUB_Wg$HQK*dWBt!!tA|Z!lsz@f0nExlHZu#T;C}43E>&
zdCS$PbT8B%se4d8QsXBi4@R;FFMS-O;(9!bi{(*0h0tKKC?CsJ@hU00xI$Iq!sJS&
zczv?HxL<tmXVn01uR<%6YUA!?<)ffbZX_>o;=E8z#a1pyvAb4xyAIB2FYY?CeU8&i
z*P7T`sb-QIRhz7ZS0^QRP}iAKJ@&JN#5LeMvDT$H*|pxe@f$}K^ZNAAatIjWQls!f
z#e+r%_cl7XevR_f8lT<pa93~j#q}Fh4Pv!PA#UudjBUc(#<~uNyun@TNghZ#T`tl?
zjW^>JBG_K*NmUXX)U|~dZ`(b%KD#D{gX%7Cd$w~wJvONRELCqWyKd;@25BM-IOm|m
zR+b3+5ZWeXbK8g4xUsX=Xm?&$;kZx9&6ZDLgi_P`H3-W*?|@Q9c}DVO`%l-+!}5RX
zv5Ce7@Cqv8ex$<4$ztVmPNTk>_6yJsZQlRteCZGkXFb6kz;!$br?mzqHWI{}0Fzn&
zzO;6>{QLD%9cJbT0MI60S%i3>@2cr$mib&yP49Kk%HsHHI#Ed)jxSBRe$Ow}QRc-6
zTG9w_E!F!Nd-Bgyc3%C&Z<gx+pO=?2f%Ok_*%w~FaGBrX&KGuukEaE)4Sd>V$FDyv
zCD1Qi$%pAmP{)3d^P%O?A--KJJL`z$^`*5n^6$N+`c@>HMUH#X9KnSe9m3R)JLF4n
zxR1Qg^Hg;kkg!5LKwJ^<`3#lW-8YPwa)N^?OZN;fQy(NUoy_U>_v0h_IGI)boJ=hu
zmi5VnEaf<s`XYuu->=jc1FzSuWZ2M$?C;)mR(t^?+ntr)+hCZsk&tnyb}ge{Cb*6u
zy_6OYQ=CiL#-X5`eY7Tq785zW3wH^lUn01a;4%VnMB{sieL3U39;m*8(IW&y2-1#c
zkkQ|=geH(ZlFwR0i$_Y=Q6FXrdc<&Rv0s7HciqQ(jl1-fMEN~YaGwqyW>i1Ml2<di
zgWxIxx~#gBA^NELC_^;_#|SXR;&KSV^#nH%+(@vJ;3k5v65LF%g5YZee{{tG>6z*-
zg4G1w1Yak}WlP^+sDq(?9_ZZuYunFoL-jXN8vB{L1)E~Gu_js5FGqFUnGDy)kuH9v
z?8`jNX6_)`w+QYexQpQ11oTOr{f4_4rB7-W(mmxOB8tV3W%|OY9Z1JTC0;p;Z6~;$
z(QO31j&EFcG0Z||2qtd<I8N?8L=>-e77Jg)=)DBrA-Io#V@dZgM1R-JiDJ=4cl&-O
zeV5>S1P>7Si7TKKA$a_9fqa;v&m3uP!g}Lcnvgxfbs6>#x=yhlK0L)P>~6P@bZ@X%
zU$WA8#Qy8ow^k|E(*S@A>f0EuAqWxB;5bJfG{d@{z~2j`*@3(m9%yQ5jD-SnBBH%$
zRGe^Zu<}3ssC~~jj^<ah+@mDC!CwE(@*(CxxzQPp=%ye}%q=y`9FBXs!-mc}A*Ujl
zQ*lDh>}1aD6LRJ_Ig96**+SN=<hLUEot^y7>EEF|xz{3P+U%XTjB-C^*loAW+w=zM
z7dCl>O$n#osKcLrTp1*lKzdjiIc&L<Y{3tw9|dmxjM4LtTaFsIW!ZV#NFR5aV83`q
zHO_wZw%OBqYfJeJ%N>j&vQIL}d+{7@w`bnIF)r4l*OcLwa4cBekYV(mVZo;eWHFt?
z(z%RoWW@vYV+=n|@C3n=1bF~h-ANyGK2r(^3JF9{8d=M7g3$w*_zwpC0O*TQT8eOD
zguUqY5xBoH>yFdiFBm=HJAUtUyCnBRF6Wc|o@RgR&Kdo>p`T-|l|XW?L!{{!(xp$N
z-?<Dn=%k&{&$FJB^#E7zA_1KDp(itdCOr}A*&Q?AIpve<-2F9W9D8Pe?XDc-9{aYt
zmb+gv?2qp1a(4XokggyB#$e`EO{#!yA=n3iElOc!6%(QaF@p36^fM}e%#{0iuTzr>
z5OaBk;C;tH_Xu3Ik{mv0XtHzfF3N9aj~zeBo^Uts&>+FyaQ76qS!!Q&_b|Kj?i~EQ
z_OZKz-V2!box4lzG50)f<k`Qv=l-3dO8>zYn+Wy*ppw5b{9gn+2$mCw9rQCKpcJ9x
z1&lhY5&IrvvOukz`r#tZ6hr^ius^<cM4Y=T?MAftQdbQQ;+6Zb@1Q$eATf^<Yo|Bu
zt)$Z;|M*2+uEF?4CjrE;X2+V!t&Gd~s1_^Y4o`9^BPm?YNEPY*TPhn$mV4qd?<aDx
zqY(?q&A|7V?wo_;Ga0dr4h;+AMbe3f{}b6aY)kXxW0qrnd7C%j^nCgv(VFOz=_i(s
z^Nmln6h&!%d%HFjkyO`s+4x`Sp5xX03@nrD1h679Hdmu+r*rTU{D^3ZOaB~rQ_@|L
z9BgKi2Z=(ycL_wBt?CeN;R>pgpX3oMNKHK%eV@SUoJj!J6g{wb2{h8%VUPW{dkg+W
z9pHExjGx?!<eC~kT&uVY57@EC(DzQqNp4es2VoaHl?)y%Lya>)v~wOd#StdA6Z^|D
zopXW+dnp@<%7ikA@8s0Q!DDKjbD5jtte>1fV}BB7rASN(NpOAv-s)v3Y^QgQY!l8E
zF?m+z^yJiD$>ibqa{vxM^&3Xd$M@nIR55EPBH-yd7hjwjdS?~q;@%I~mXIx8vN)^a
zDvLiU^U$R%Eq3E?3~`ZwqdjUMHsTwb@Nhs*+Lkj@e5q>~np{RG=Hyj`QoZmqwUl55
z!2F&E?wgvVREQ+GUz|uh+tRgsZNhmx!@2KIr}M7oO?Fi*K|fN*_EPC+ke{{n(zx8d
z;Qo>+OoM>a@*(+-)I|G|@)>Hsdf(6?b9;qCndJF$^{<ScKi+??F^slZWU`toNeUXD
zt#<HxBhKr$klDogPa|k1;COT9&No;nd4{_+MTA3mDWVv?vVkjSs<_1#s*Y~7-}zp-
z`z@oV@PTWL_#l?)s}{9L!auCJ2jbJ`VWVh(XoW0NeDGca@<QNpV)`6M*#t?cO>kyL
zo$X@rdGlF9?tbHhWH~Go&jMPoi1DL>46>7pTow;7KWNYT{$%6bo^!vyDM_;so~9mc
zl7O^6r%5Gk&UwV!z;gVF{4}X%0co;Yg-^I&=-2?c1E=3%flUOYUVqG*mQMuajRf2<
zN6FOtO<{wd2m(s&4G;N^{dUts6W#9__O%baGF`YV6r(QqA^>Stb0$n5I~v-sJ0G4A
zmr*!}*sTOt6G#`OPxBb<Adt1zCWd57Uk(#<F##sU<J{9o_-gQER_Hvv7j-#9BC$MW
z)dZZ!KTdD`_=^{*-T*5xfHy4IWvL^2#6ItbrS6Xm`+^_t!(SKxD!Y~3wRkdPSBd>j
z?(cj*sZ+_#vpvXx#tU!#2&WZ5Q&cx{;N&5l7k#}(h2eEe`ho1qbGM&5a~{sEcB}}@
zte9PE|KgD{xj>1bFoYFdgyIMbx0B&wrr!jBbP+`|ol1_|HVvU>WC&75Myh~R^n0h$
z!5Iy(%SoKL({d6?3M5+_o<@?rJDpikCPk7y+RRz8b<KvEb7!Rj3zLC)sg!vYrjIy*
znKS37LbKv?dV_O&Lkm)&ikW6UiUiE%D^Y4zMr2lIq@p4tr($kqWS&_Vs%i}2{WczQ
z1!m3Z8=YsK-6wmYnYDiFOkl}puJf7aeCEriVwUr%a6YqBQq{i&Q_UQ6V_?>Th4W4O
z_+w?`UPF?8m0%iq$?dU}bjfIY^5gHV{UV#Gf2Ky!A<K5wSExh=4EG%n76bc%{n*os
z&i<4sUu=hVQg*shQ60m%&_v1VYXMO6cNxClWk3G(2z%9-fh+nvGebLHanM#0Z8?Fz
z7nCxz3Secr9hA7<e)fsWjZS;ZlbwtIjNF8eYjs|1#AV&=CM#AKQN6pWW>zEaNHv#x
z%tCCggmGkO@%50EFq%z!=u_LL`BB3c>~7S)?kT^IxA?G=7_uLEs$|j)PU(jno}G?Q
za?kP)sf}Hn=@;97e`<`e%O3T#Iqm<z-7neQPY=mok2;g0qIjzQ4;5UR!8WfWVyUgy
z27X8t`jBg`4{4M2f3mIT?cL8*%<wZ#fXQBDBR53F-R{kDd;<#TsQu_O3m5ffq{;}o
z6zt7Pm$@cQoqo}<=R7-p5~T@Ql!1%`p7Re*M3}2d*_b!8x5*7UJ=<=7c0=~RG5w(Z
z%Cn_7UijqM3GPn}d+2k=CYsqBB0;<w#tQU(qI^K`F~M&LJ|VaYe#Wpk7W0Up2@z=4
zYTVgm=5RD#S)R2`0^1~Qb^oYa%B+TVvn;<Tgs&n0ySJ?_|3r2#6YfYfkz<ATg>)p0
zTOln?&22J^B?dLwN1rc=bAn4`H@DSAo0>EoA8d_T!K#Q3>oXC<6~(rwzMW*z3Fj7+
zOJgAO0hm00j5OmZDIY_Fm$*X4M+q?;<wtc6F}YAl48|h^PG60os_C^lq`2nM+~x%h
z@pKlC+TobwvWT}niZ^AKFjyozV9h#$R|-+wo5hvWW?X-_!ZkY56l;xXYCgPioXR(A
zJRu}Oe0?KZrH7H1an<2+MCC;+JPcem$M%H%>I*~Szapuh0+_ixM(1_XsI!avJ0=$K
zE{voGoa+E1i9d>9Ji!El#|accg5V8;Hwpep;3gd&f@O?)85&412q3NpGdPT3E5W%0
zL4rDh2MHb_kWId)7?LfT7Z`eppqjO=CHR2R-c9%fp34V$+ZazoaOd%lz({x-thOkg
zDLH?VBdVV^4ErZ9PK}#hFd0wVd7C=nsR}Bypt{bgZg6(llGjZuSpy!|#<=itw%bLk
zH?i%{+UpWmkh?x+KJIGcU7NF&xj`sVP{)LKZ%Zhqe;Oy@;#<?LhMNxX7Qwp&?-6`N
z@HT*15RN#vzVW&=DnrUFtbvQ2OvJwK$HbjTswx1?fz1&-43T|YKn?jm$>ikLE+2))
z@r;#bS%1ik0<KVdUvfRcq+F*sUvb?_8q&^cJ>GgTB?ot40&P~z91z6Am0onbti>JN
zjJR`zvX&HXMa_C0L-C^wiZ|%NU2*Qj<0XiGpLybnE@wV>o19aI4u)Q2+J_9$k4w0*
z0K1w=mlJSO#=HmcD8T~+4-$NzfU};yjo=3a*T)GT0x+}i_Yj)f5?Q$N8U!7oaC1w{
z%!6`N;m$0tjmQdgBk^Um_ZUKn>>6Av#(|j4<P9K(8}%5EwHTFFO{moz?%0d!#ulEP
zO|J{#4}j?Ex;DHZmZZ2g@%Ii=iSD{;ypjlo8-lf#{!i4&VT%)&01hGW2OIEb8BF-h
zxGXJy9wA|z<+$N=zV)97WaWG*LyG}K-D!iO`eCMVnsipge<P`bKtf!O^gyq3vSa#c
zlBDYIvBWC`Q&{Xq*5{S;A(>rrkcdAcfj=`+71=N65$0gofS3DHq`A2fx0$fviz6jm
zD}S9e9U-V?bu|QChHF9cv`5Z|z9tP3oI&s%V#{LVU4{-2979IJuymvSyB{rHaut#i
zS^N_M$tT~<#PRn9r9|>Aa-J?PL3vm7IyhuLw;ew|dA$c)7I^wP#TacAxc&G`Kz_if
zo`S3aaqQ?U@|F0Ce7U}Y0nNBgoBXd})POA9tj#Id>KpDG>>FPYcTL1EVE_Jj?V!8@
zU%?395Z^H0kb;4}+&sfR|JC4fd0d<4_TuYx4{+l>I6^*oc8lHCZ2y3Ph!q>Fy#oyI
q2yCRd{q7-o#l}pe&2&#Uym|JZ*Y3hqpGRKXg5&|%rY$ha_WciK)2YM&

diff --git a/__pycache__/mcmc_chain_analysis.cpython-37.pyc b/__pycache__/mcmc_chain_analysis.cpython-37.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..bdea663ae5d36ec0de32219c100df06016b9244f
GIT binary patch
literal 4481
zcmbVPOLH7o6~4D$^Uz3^Y+05glT3Id83R&M#Rg(*6DXzFm?{F5Dpg2Vt8b6g<LRFE
zy*(;P=~>AX-nfcI7Qo6{W|JL%fL{PR>h4mkyk^G+zH@qJ^sq?+JvI0JIOjX(eD~h>
z78lzB-;>|}E&lUWA^t&?+2^8gAAj>H5+!zolFAv%9S6BPbTe<q%lsW**6qPgpk7nH
z3LeXyhU%(@YNFIsE!BQ3c3Nsdb&$8!qFO?}psuK8<Q=u5t|DJl*VHQVC3Rh`A-|$-
zsCDGa>ZW=H`HI?5w~${|-%zh2zowiIMDO+EqjneH_BfBrbd;NJIqD{drH<mVYs#oh
zO_|1~+t;J48+A=mc1QhgydR~xd2wPq%8GHBbbtBrepFtXGheKGki_LkoAz^nA|CV{
z+tgA1P(@{Ay?msz?bx9ne1zZZZ@!62CC$~}iM}KpA2%mX-r<^C*DBJ;dk{;Si;r?2
z|7-Z0H;|NKAS*Ed{7=NXP>yn+x(bH*?X;!u+!VaT^MB)ePV;-xdeLEOdaiZy!nTS?
zmnk7-7b9cRy*!LZ<GeH&Cl4vfKNJ6YckAKN?pCx{s4cUr$LSz9Tj_q5g}cdqbPy60
ziP_5HEDjmUFv_Fhkx9*Fab(x>az7kKdM`0$Jv8jYIvrAifsx1W@Cu<xkCSje85W7&
zyfYld(a_x8eEIA-q_n~uliTMospL(#;aA|IJ<QNFzbAD-nX4Pr)TE?E4UVRGux^?w
zj9H*|5?wD+vV`Oc{!>M`2TRZ>yO@ns5}G23(|#KB(7J!R3N2t)UeKI&Gqk;so?LP}
zj$T1|uJ6o!DE~iQh^Dd&t)T$-p}*V^PsBs-864+S!eyKXH7YKHU^{t8{#=YwxQ2Dj
zI0F_WDS1eF28dm%^~uK>8<}?aC;sC2t+S5H7C!4CznE)EmYs|G3MwortooJWYVelO
zq~g?inrsA#SE>EXde>21rzR3}W;-6UObygMvlzFs#g@;+trelJGlMxDsc)d=LeNs+
z{d6{r)EE%PlDJ|4)Rdc0ftPaY??qV_O`)v?ol9%m+{v)s!dkV^gv;05N|UPAkA}tm
zf8ByRd5hQPThL|~wjiUhHUu*p1^O&7?J(W2UTqU)#V`qnLt|TcSl7TR>@}BFYYd0H
zDN$mr7i~Y^lpKv>lsEkx1={o@XCfy~DF;r)$IVk`;+FE6Q@N*d;#E!sQe77t0&Q;P
zRv!4#y(UiOs+jng<yNAaait)htow<GpE?!7faXpDUNOLmjk>LJHbm7}6_rn=8)DL|
z8p?a<=szg`1KP{0(6WJjCatPfHTy*2!8Dj^ig`|st(l7*4zjxJTi@(Qpg!Us>t#ul
zvkP0#Jkn*aZG-xhvNEgrrB#?TRHHH+*o9dURV&}K3nm<d=FvFPiL&xT+c}7Ih$Ad4
zhGS!WE?>66ay=r_Iy%A`kQLv79EKu8$hKU<|5bS#>AGw>ASS85kDi#!LLZrk_>)a^
zhyou?e6GPTe&qq|+8j8pDp~o{W9kLyVV;2J!USXr0PoSv%9{coP>Fyq1q?d}yY#Bo
z2}PB60c_ocEIq#_b#^bJ1^q1~*3F{BUR!^g8t63DJxWOJn&#i6(p!|=rsTVne2)^2
z$+N&*O%9@A7-b`zq{oR08OR@CprJ7WSZ(6JPMB?>K6h(6CFYp@1r>H{m)x3M`J7{W
zT=Nv}1&pQ(XLM`lh9i3~pee7x*=a1JSv$KoL!7S%qVKZH)5zN8-Fn4IgZFAo_iDm!
z*r{@XeiLpVPzmU_*wx|S?CR~RUA5TNfqoa7AkIm%`V8^)>s01*^`Z%S3yi*%;W^mZ
z6ld78QNG*&*Dt$nn2m?!kFlvCiW0H|?8ypm&+PCHTIY5kN#}MbP=RORT7?wUathx=
zFv}8zaYErO20rQz>XPe{>uwo5n^A-qG?jA$^nr06x3lk8@?1Mb1Sq6(`yPDOeXxgZ
z@%{bYwD8w^y;jWxZnSG`kR>VwArN)ghE8gxU|3yha*Lu5qJaDL0)}*Fu|FCrIOCKw
z0=kld3}(&K>J`dEx}-+K)S!QDI_*;XWtYFn{b_zMog4Bg*2g+Ohd8z0yQ^|pE<5Y+
zeElw3FW^niz+7<$g*`Emhnto7nfT2cKvq=FiTuc$xD#*UgFECiS%Q*@#7{u(QaoyZ
zC~ExcQ`F=2Z4ebjk}JDR?8PpS*0HEJvB?w|8?dv01sNIrLv-mMQSvq=mqNSEq0KXX
zjc!9UNB6QeR9kcOebnY^ks-)FtQK7gpa@)Mc21GI)Z#jt+^-hOt?B13HM~rK6&Apx
ze3r5tA7iMYYb&|6{t5D!cBYR&r5=`mR+~!zC{cz4kjCdSLF`NdlqG=mtKh_)Nk_z9
zOV_BSn`30#!zfFppwb5R`XxuH`|PckyU}#!MUXWGYK;H%64ilT(2QGWno&Lv(Pf;S
z6rv@L2pj><DG2NY=>&K`^GaX2CsKK*jwC|mIz2+4I?BVHpz+j|q6|*(nEJ9cPQ-;3
zXg=2b+Vk<$y&wbX|FaAYg#+!J6sW&^4MlIf_JiX$iS^Ry3(Q0Zc#_3YKFU)(q(X28
zdc#j9n}v1sE;IA?ufWsGfEKs-4!>6L$VubG?EUGF^nG@3?~~Qf-~ZzH(|f(94T>~=
zI83ZJ97W1{BgC2wDCDEkE>0OZS7BldzgVnaXuKBqGK|NjhZkb0x3QURnDOoaZ|2Lv
zwlcjk(gnMIrVq#~x=4hxhKO62fm{QHuW>%BQ?7r3F}Uo*GsJdJZueGgBgEw>4nynX
zF;X1WWV&GDv^d(N9SoDHCePCEQ$o)ymR8fHTJs&7Z^0Z7baK~mOI)V!Q$oSVe#+aQ
qQ|Ih%(7Xd&$HU|<#U-oa$tF~3y6twm*>w4LrP*qB@M|^wX5&BCxVPf~

literal 0
HcmV?d00001

diff --git a/__pycache__/process_many_chains.cpython-37.pyc b/__pycache__/process_many_chains.cpython-37.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..44eddf35ca7b5b24cc52c61abb188e42c155fc62
GIT binary patch
literal 3920
zcmb^zO-~%jwW_=4gJHmc`L>NQ#@GYi0c<cf#^A&o?;aLwqdhojWz}h>ilL#qr?;xd
zVFp{;2wU09N*t3ttRScCCCVWO=O^TtoT9o#IgN6PBIOW8O5UrU8I0FnbI9~mzxV1@
zy?Wo(-|y<m5%~T5#h*%Z{e=7tiPkR<$bA^{FD4<3FzOSVYG~MohNf+5XxSD(!%qci
zJ57_dnILO-H1V9BV{6P}spr)0WcOK`WdO^wEbDkq>@JpLodCO8o^=5{!@AiS*28*P
zUpd9d2(f#34?DACh+lH+u$T9-v#fvLIHLC1a2Pn<!1-O1GtlC&!7{b`*$^9MBWx73
zJ6BF?Spz)J&hJxu5a0!XL%egJ*u%V&jR7?Rvlju50zD4!9KcHe&jXwQcmd#LfMWnB
z0bT@n1>iWFvM;f#_QZj)`y1Z>J~=Qx{o20FCpoMJFI=;)z?kA^`Q;J93NXfo2%GLF
zY-W_$S2;br#wYmIO|+_Tf|@=-&77dFAJE<3^Y?hMDfjvaX=3k<kbPsHf*$nKgP$1o
zEQ~pOo=>r{B?5X(@;TIto7$hn!yEka5oNR5t2fKUUSPyt+^yGaHitak!w2{rpJ($x
zEgsRMw0(=`4s%DuUg8Uh6n-n&5AqlIB)fqtlQn#4m>hvNaNeRmfAKL1hb_Wy@8k{)
zu|)|!ACn)G(D)H({Mo-A6Au2`F@XCyU$n1;6Z{Ha+%bXvHORiHWliGTz5!D2?8APx
z@T#Uqsl#Qq$Zj1O_FdSq1iZN;s`sV&0;BNEyfDFEH-Um2+g{N&UgCGz?JqO-J=k*x
z_Dp@tf<b-{;=RgOd4W%@5I)qVr!_s?re`!g(x$H?4c5H`>&(M7*x$#NPwLF>LIltv
z<WB>C<vaK@z`yq${OiD9eKQ}ln0pKr{Rt_2uYvXcs6bWc`h)ccj^MKD#}7*ePQwps
zzdhsx5KpP>`vuuV4;s(2&x|myDKnv-8X(hpRHC>Seh5=A&ccvW@QTSceM+_sM!z7#
z&@`C&B{mWAqB*l^j1VoN@n85ga6jl~xx1b$n92wv)e*TO_F}IRsxDVbuN*q1N;Qln
ztfRa53cn%$zB{{9+n9CB5u248qUvpja@N}l0%wD7xlbG!yD^uuk*JioluqD=HK(-Y
zdZ8>vH8mK<TaNFFGM90()Y+_tCG3b&JBHThwFF0cb?$6&KjNae;#W$pFIS6i-raU7
zN@)9J`W2r%rO+7O<rjew5f;1(RG^rGDbhHV!B-Y14AB8VnO?|X4R(0Zi7(V&biwNk
zjHV%)10ryh>Iy21J2K)WZ__L3r3t?&1qG<yZ&{9eC7w*EMZ3+sA$kCB+xve=Arzrg
zXb=!c9~|@?`GV|Nhmbf&n90VY#4DzDg35;ub7MEEc#s;(l+^%yfquaZiH{wkx)XcC
zYM_IRC;Srzsb4#0ydL3W4(3JMQjBa{)lnd{su0z8sw;Y7wdjLaGkT&MaVpU2V`T+g
z(oi-qHPD)e*y%MGQ_XHgOwa0w^&G347=luc17_r@K&wv_Ef(H`)-p6>8+}Y46%5rG
z3pbPw^JE;VRM`yzx9QqM(RQw>oMvjbo(EZp(j%uoN({kTm2&-P>#av<rat0g`-pR|
zAAv#9)}V)VP{5u}-F9@l(RSse8$>>Ld|#@L&`BuhDmrUkH!b7)`cQ0%iCtrWsZ1pJ
zhUZFU>9^_nzREY}oXD?&<az~4EUx*yK2mbSO6Wn;aJ+D{A|2dQ+?LSUsHpe9#!VPK
zAAGd#eBuSp{Jawjue4KsIrif3KYaQ1TKPZU{NvB-8$WF!_t#XW8hV=*5j0t!ua&?0
z75@I%LjJI(Yo@%{|Jx7#T7SOQzzlAnw;I@C8rWYN*j*afU6fhlvVpa&n~>e1JEjAb
z;~hf!H88R&B{^qOWjDPTIkJACI4hyYLci7pIlK2vTMPM9ii?pl#72F&ILm6G8+auL
z%=126Uw2BPBIOfkr??(u#kGh#fIe{r0{!26R*dQ<;PruzDlYTF^7!7f@wigi8Q)tj
zE<Cd{AB)O{&x88;w+jIkXi-t8-tB`7&c;!HkW)yjl#YtBu#c#0(+H{ZXd39@xnYSb
zlU0SvNcRbMJXSxq?(U!-$u)XXo%;I$T&-Sm&p>)TfeQ%qIxfz+w;6i81|9aZjdx0#
zKaHgNhi}y9<h^$)olYg`R60TT!_M=tKpz(Fs}v+qRMXd+N=LOV83X_!>3yy&tP^_>
z^2RCqLDQ*j==wG3Nr8E+dIL9({7US58*r;5hEg;%BR^JF>;+ukRl5HP%wK`ML+D!8
zDU5liEKCh$RU#g$JmL^yoy3q8)s1alilw-XC#LW^m`}hBTT?d^u#kYo1l&r%QUY!#
z;EtWv+vl2aUfe;MG=TI<PEn~yz;$M)Lq}9YsWSS*hl}L)PSrUS4NU6qnE4iRqlMgT
zAq%R5gJ-#bZy;sXJ*a?fWntc{OtWlN8vBMDs}8IIU8|Vn*dL()0xyI##*Z#8l(p$$
zIDFT4RmTr@OFYuwX4Yp_FIE}BBi}7?v8>&Ll_XU<sSafVgzADAztUZgZ;$K>e2!Fo
zzKR(k58;)i$&;sfnl*E|TsCXs(+Bf7dT83{*K5+vm8~6FY8mit0<#u8!+<5voL)Dm
pXANr5QK0*Qp9fjpGzao>FawZ|X1WcWv0w}v$8r(t?ir`}{13l69=QMj

literal 0
HcmV?d00001

diff --git a/__pycache__/simplex_plot.cpython-37.pyc b/__pycache__/simplex_plot.cpython-37.pyc
index 1e8de26b8b6d51af0763e873273c06afd448dab6..1c26491f573cb6355e901c949147f3d921c87a30 100644
GIT binary patch
delta 900
zcmZuv&1(}u6rVTQ-E7h%O`D`?8q*KNkF`QSz>5}As<)O>Q4qY;(AlJ0lQhhxkcLPP
z67Z5DorBhcB)1AGh$6l9RxgTp9rd6mgDCz7`X*`h;#>ANZ{F|y-kUef>sUP=nhOT~
z1ntMO_bYdUFG6o1UNx?7mh91~qjmad6`vvDSGn`}yWe)`om*=`*24ld%CrPwK^AIg
zKx)u5S_GNi4io+l(^uA-I`g)GnMuO>;{+?1|MJSA8mW<HFAIwgbYec7r3sRyEYgjk
zDjFpY)bs~*8-SEUWo@nRSb;2lb;>}gu)4-#4chEyVeCi67do=lsmBIb98ZBd)u#pB
zs%)^KA5YMXumDT6=-++xcwZ*X7#nJ--ehfnCD{p<YN$NzKJZV-##;nP(L`@hc?_-*
zYq<K@KZ>>_Ha!phixzGB@ESQ4m$jL#{o(retNYIna*uYEJvnF(B58rtT&HU1vK{-G
zZUaVhDs{O2F@2Q0o0j0_kU$H}Dov*1xDw7u;Fd=kg$%!KbM8f&_aNlny;b%rEj{bm
zY_*qdaC6yp3Z=YV5l^*bG9h6}8d)c|!EFf((pYec_Vg7=izV?@%b(E-cBUx3+;&$p
zC0m-VU0$x_okG{!M2`oskI3n>V2M3F1ry>>pPxRDLE<wZKtm7(6BI`25DdXtIt&**
zlQ826Lzw!JH^5K%)XTJS0HBNC-hnv{#cczDt0%Z<nLLiG)T=gMcctDvi_SB}oXyqc
zNzV$#X~~fu*C}H9xVT}Sj83Qu9hg+!7{Vp7ZJvTL@ygsWRRfFi9<k_4KwPZ)Mj#^Y
zA*RHhF9k7i=(}l5b(G??-?A>CD?NF;$Y)hhjmRwp|N3;BiHd(bVwowEYFYq?$MGp%
G1pNcJX6K3k

delta 734
zcmZ8e&ubGw7@cpjyP0f$G(R>8u`~^6wM&{xa<G>|MMMuRA*dh}L_=qjrf!omyCpG1
zJc!;D?I5BTvlpR95xn-M7wsP~XD@{wym<5AOsYlvVBfxZ-}m<M%|6M0$|kPF<0g{z
z>-*ux)7Xc^7a#(v0%`*IF%(Qg_EzQNxqJ#gyfu!~D2<LV)e0yG0UE*B2#{1`6E*&y
zs{7hlr{Mr7nMO2`N76}7X8gnmjnH_C#`!V6_$c1O1=PZH>O>2m10(%4eh#M@DVZD#
z@aZ}I_w1u_h9(0!>7y-+rf8aGo~b$Q1+aUvn+*`4E|Dy!7=@e2wuYb0lB|Z?^Go#s
zH?_Um4aeVQU9WVl)81vyLTRJvx(+LC?|Ob|(cA3x=k1B7ZsVhfyc4GL_u6vh?Hg75
z2kmcPRsC2w0@-Z;kyiwGBGhZOhQxO~pLg`>>AZl3Fj|us0c*lo>vWykRe|db{!Vub
z#B=<+1Glp+ChUP?T656xc+I#xJ14&`jQjm2YkLB_)%C~@>uot3zC-sBEAyAe>f*ed
zK_-D1PJjgj6k0e3IXH)la0w?NiA~uvz{ITbF%dohApDTzDq(3nQh%LACDd)gEcq7t
z&|y7K=-y^?#}PVfcH0ii@>Ilv3;cZKtW{Q)2`ng9l(59_MCPE#A4m3xx-T2E7@si<
zVDTk$7Lxq7#525a&On;KGS`jD#L0h}wq2D}gj}b~mX)aXVQTS{h=0QP^=LUWM`j4t
Nw3ysf?uzBs_&2ImxRL+>

diff --git a/canonical_infos.json b/canonical_infos.json
index 6b9451ce..c2ca57ff 100644
--- a/canonical_infos.json
+++ b/canonical_infos.json
@@ -1 +1 @@
-{"SWC_023": {"seeds": ["302", "312", "304", "300", "315", "311", "308", "305", "303", "309", "306", "313", "307", "314", "301", "310"], "fit_nums": ["994", "913", "681", "816", "972", "790", "142", "230", "696", "537", "975", "773", "918", "677", "742", "745"], "chain_num": 4}, "SWC_021": {"seeds": ["415", "403", "412", "407", "409", "408", "405", "404", "410", "414", "401", "413", "402", "400", "406", "411"], "fit_nums": ["773", "615", "107", "583", "564", "354", "142", "184", "549", "185", "924", "907", "105", "531", "9", "812"], "chain_num": 9}, "ibl_witten_15": {"seeds": ["409", "410", "401", "415", "414", "403", "411", "404", "402", "405", "400", "412", "408", "407", "406", "413"], "fit_nums": ["411", "344", "496", "600", "716", "18", "527", "467", "898", "334", "309", "326", "133", "823", "740", "253"], "chain_num": 9}, "ibl_witten_13": {"seeds": ["302", "312", "313", "306", "315", "307", "311", "314", "309", "301", "308", "300", "304", "310", "303", "305"], "fit_nums": ["897", "765", "433", "641", "967", "599", "984", "259", "853", "385", "887", "619", "434", "964", "483", "891"], "chain_num": 4}, "KS016": {"seeds": ["315", "301", "309", "313", "302", "307", "303", "308", "311", "312", "314", "306", "310", "300", "305", "304"], "fit_nums": ["99", "57", "585", "32", "501", "558", "243", "413", "59", "757", "463", "172", "524", "957", "909", "292"], "chain_num": 4}, "KS003": {"seeds": ["404", "407", "413", "403", "414", "405", "400", "401", "402", "410", "415", "408", "411", "409", "406", "412"], "fit_nums": ["846", "256", "845", "945", "293", "406", "420", "109", "690", "421", "54", "866", "784", "81", "997", "665"], "chain_num": 9}, "ibl_witten_19": {"seeds": ["315", "311", "307", "314", "308", "300", "305", "301", "313", "304", "302", "310", "306", "312", "309", "303"], "fit_nums": ["179", "951", "613", "6", "623", "382", "458", "504", "406", "554", "5", "631", "746", "817", "265", "328"], "chain_num": 4}, "SWC_022": {"seeds": ["411", "403", "414", "409", "407", "412", "410", "413", "415", "404", "405", "400", "402", "401", "408", "406"], "fit_nums": ["408", "884", "62", "962", "744", "854", "635", "70", "320", "952", "8", "67", "231", "381", "536", "962"], "chain_num": 9}, "KS022": {"seeds": ["315", "300", "314", "301", "303", "302", "306", "308", "305", "310", "313", "312", "304", "307", "311", "309"], "fit_nums": ["899", "681", "37", "957", "629", "637", "375", "980", "810", "51", "759", "664", "420", "127", "259", "555"], "chain_num": 4}, "CSH_ZAD_017": {"seeds": ["401", "409", "405", "403", "415", "404", "402", "411", "410", "414", "408", "406", "413", "412", "400", "407"], "fit_nums": ["883", "803", "637", "806", "356", "804", "662", "654", "684", "350", "947", "460", "569", "976", "103", "713"], "chain_num": 9}, "CSH_ZAD_025": {"seeds": ["303", "311", "307", "312", "313", "314", "308", "315", "305", "306", "304", "302", "309", "310", "301", "300"], "fit_nums": ["581", "148", "252", "236", "581", "838", "206", "756", "449", "288", "756", "593", "733", "633", "418", "563"], "chain_num": 4}, "ibl_witten_17": {"seeds": ["406", "415", "408", "413", "402", "405", "409", "400", "414", "401", "412", "407", "404", "410", "403", "411"], "fit_nums": ["827", "797", "496", "6", "444", "823", "384", "873", "634", "27", "811", "142", "207", "322", "756", "275"], "chain_num": 9}, "ibl_witten_18": {"seeds": ["311", "310", "303", "314", "302", "309", "305", "307", "312", "300", "308", "306", "315", "313", "304", "301"], "fit_nums": ["236", "26", "838", "762", "826", "409", "496", "944", "280", "704", "930", "419", "637", "896", "876", "297"], "chain_num": 4}, "CSHL_018": {"seeds": ["302", "310", "306", "300", "314", "307", "309", "313", "311", "308", "304", "301", "312", "303", "305", "315"], "fit_nums": ["843", "817", "920", "900", "226", "36", "472", "676", "933", "453", "116", "263", "269", "897", "568", "438"], "chain_num": 4}, "GLM_Sim_06": {"seeds": ["313", "309", "302", "303", "305", "314", "300", "315", "311", "306", "304", "310", "301", "312", "308", "307"], "fit_nums": ["9", "786", "286", "280", "72", "587", "619", "708", "360", "619", "311", "189", "60", "708", "939", "733"], "chain_num": 2}, "ZM_1897": {"seeds": ["304", "308", "305", "311", "315", "314", "307", "306", "300", "303", "313", "310", "301", "312", "302", "309"], "fit_nums": ["549", "96", "368", "509", "424", "897", "287", "426", "968", "93", "725", "513", "837", "581", "989", "374"], "chain_num": 4}, "CSHL_020": {"seeds": ["305", "309", "313", "302", "314", "310", "300", "307", "315", "306", "312", "304", "311", "301", "303", "308"], "fit_nums": ["222", "306", "243", "229", "584", "471", "894", "238", "986", "660", "494", "657", "896", "459", "100", "283"], "chain_num": 4}, "CSHL054": {"seeds": ["401", "415", "409", "410", "414", "413", "407", "405", "406", "408", "411", "400", "412", "402", "403", "404"], "fit_nums": ["901", "734", "609", "459", "574", "793", "978", "66", "954", "906", "954", "111", "292", "850", "266", "967"], "chain_num": 9}, "CSHL_014": {"seeds": ["305", "311", "309", "300", "313", "310", "307", "306", "304", "312", "308", "302", "314", "303", "301", "315"], "fit_nums": ["371", "550", "166", "24", "705", "385", "870", "884", "831", "546", "404", "722", "287", "564", "613", "783"], "chain_num": 4}, "CSHL062": {"seeds": ["307", "313", "310", "303", "306", "312", "308", "305", "311", "314", "304", "302", "300", "301", "315", "309"], "fit_nums": ["846", "371", "94", "888", "499", "229", "546", "432", "71", "989", "986", "91", "935", "314", "975", "481"], "chain_num": 4}, "CSH_ZAD_001": {"seeds": ["313", "309", "311", "312", "305", "310", "315", "300", "314", "304", "301", "302", "308", "303", "306", "307"], "fit_nums": ["468", "343", "314", "544", "38", "120", "916", "170", "305", "569", "502", "496", "452", "336", "559", "572"], "chain_num": 4}, "NYU-06": {"seeds": ["314", "309", "306", "305", "312", "303", "307", "304", "300", "302", "310", "301", "315", "308", "313", "311"], "fit_nums": ["950", "862", "782", "718", "427", "645", "827", "612", "821", "834", "595", "929", "679", "668", "648", "869"], "chain_num": 4}, "KS019": {"seeds": ["404", "401", "411", "408", "400", "403", "410", "413", "402", "407", "415", "409", "406", "414", "412", "405"], "fit_nums": ["682", "4", "264", "200", "250", "267", "737", "703", "132", "855", "922", "686", "85", "176", "54", "366"], "chain_num": 9}, "CSHL049": {"seeds": ["411", "402", "414", "408", "409", "410", "413", "407", "406", "401", "404", "405", "403", "415", "400", "412"], "fit_nums": ["104", "553", "360", "824", "749", "519", "347", "228", "863", "671", "140", "883", "701", "445", "627", "898"], "chain_num": 9}, "ibl_witten_14": {"seeds": ["310", "311", "304", "306", "300", "302", "314", "313", "303", "308", "301", "309", "305", "315", "312", "307"], "fit_nums": ["563", "120", "85", "712", "277", "871", "183", "661", "505", "598", "210", "89", "310", "638", "564", "998"], "chain_num": 4}, "KS014": {"seeds": ["301", "310", "302", "312", "313", "308", "307", "303", "305", "300", "314", "306", "311", "309", "304", "315"], "fit_nums": ["668", "32", "801", "193", "269", "296", "74", "24", "270", "916", "21", "250", "342", "451", "517", "293"], "chain_num": 4}, "CSHL059": {"seeds": ["306", "309", "300", "304", "314", "303", "315", "311", "313", "305", "301", "307", "302", "312", "310", "308"], "fit_nums": ["821", "963", "481", "999", "986", "45", "551", "605", "701", "201", "629", "261", "972", "407", "165", "9"], "chain_num": 4}, "GLM_Sim_13": {"seeds": ["310", "303", "308", "306", "300", "312", "301", "313", "305", "311", "315", "304", "314", "309", "307", "302"], "fit_nums": ["982", "103", "742", "524", "614", "370", "926", "456", "133", "143", "302", "80", "395", "549", "579", "944"], "chain_num": 2}, "CSHL_007": {"seeds": ["314", "303", "308", "313", "301", "300", "302", "305", "315", "306", "310", "309", "311", "304", "307", "312"], "fit_nums": ["462", "703", "345", "286", "480", "313", "986", "165", "201", "102", "322", "894", "960", "438", "330", "169"], "chain_num": 4}, "CSH_ZAD_011": {"seeds": ["314", "311", "303", "300", "305", "310", "306", "301", "302", "315", "304", "309", "308", "312", "313", "307"], "fit_nums": ["320", "385", "984", "897", "315", "120", "320", "945", "475", "403", "210", "412", "695", "564", "664", "411"], "chain_num": 4}, "KS021": {"seeds": ["309", "312", "304", "310", "303", "311", "314", "302", "305", "301", "306", "300", "308", "315", "313", "307"], "fit_nums": ["874", "943", "925", "587", "55", "136", "549", "528", "349", "211", "401", "84", "225", "545", "153", "382"], "chain_num": 4}, "GLM_Sim_15": {"seeds": ["303", "312", "305", "308", "309", "302", "301", "310", "313", "315", "311", "314", "307", "306", "304", "300"], "fit_nums": ["769", "930", "328", "847", "899", "714", "144", "518", "521", "873", "914", "359", "242", "343", "45", "364"], "chain_num": 2}, "CSHL_015": {"seeds": ["301", "302", "307", "310", "309", "311", "304", "312", "300", "308", "313", "305", "314", "315", "306", "303"], "fit_nums": ["717", "705", "357", "539", "604", "971", "669", "76", "45", "413", "510", "122", "190", "821", "368", "472"], "chain_num": 4}, "ibl_witten_16": {"seeds": ["304", "313", "309", "314", "312", "307", "305", "301", "306", "310", "300", "315", "308", "311", "303", "302"], "fit_nums": ["392", "515", "696", "270", "7", "583", "880", "674", "23", "576", "579", "695", "149", "854", "184", "875"], "chain_num": 4}, "KS015": {"seeds": ["315", "305", "309", "303", "314", "310", "311", "312", "313", "300", "307", "308", "304", "301", "302", "306"], "fit_nums": ["257", "396", "387", "435", "133", "164", "403", "8", "891", "650", "111", "557", "473", "229", "842", "196"], "chain_num": 4}, "GLM_Sim_12": {"seeds": ["304", "312", "306", "303", "310", "302", "300", "305", "308", "313", "307", "311", "315", "301", "314", "309"], "fit_nums": ["971", "550", "255", "195", "952", "486", "841", "535", "559", "37", "654", "213", "864", "506", "732", "550"], "chain_num": 2}, "GLM_Sim_11": {"seeds": ["300", "312", "310", "315", "302", "313", "314", "311", "308", "303", "309", "307", "306", "304", "301", "305"], "fit_nums": ["477", "411", "34", "893", "195", "293", "603", "5", "887", "281", "956", "73", "346", "640", "532", "688"], "chain_num": 2}, "GLM_Sim_10": {"seeds": ["301", "300", "306", "305", "307", "309", "312", "314", "311", "315", "304", "313", "303", "308", "302", "310"], "fit_nums": ["391", "97", "897", "631", "239", "652", "19", "448", "807", "35", "972", "469", "280", "562", "42", "706"], "chain_num": 2}, "CSH_ZAD_026": {"seeds": ["312", "313", "308", "310", "303", "307", "302", "305", "300", "315", "306", "301", "311", "304", "314", "309"], "fit_nums": ["699", "87", "537", "628", "797", "511", "459", "770", "969", "240", "504", "948", "295", "506", "25", "378"], "chain_num": 4}, "KS023": {"seeds": ["304", "313", "306", "309", "300", "314", "302", "310", "303", "315", "307", "308", "301", "311", "305", "312"], "fit_nums": ["698", "845", "319", "734", "908", "507", "45", "499", "175", "108", "419", "443", "116", "779", "159", "231"], "chain_num": 4}, "GLM_Sim_05": {"seeds": ["301", "315", "300", "302", "305", "304", "313", "314", "311", "309", "306", "307", "308", "310", "303", "312"], "fit_nums": ["425", "231", "701", "375", "343", "902", "623", "125", "921", "637", "393", "964", "678", "930", "796", "42"], "chain_num": 2}, "CSHL061": {"seeds": ["305", "315", "304", "303", "309", "310", "302", "300", "314", "306", "311", "313", "301", "308", "307", "312"], "fit_nums": ["396", "397", "594", "911", "308", "453", "686", "552", "103", "209", "128", "892", "345", "925", "777", "396"], "chain_num": 4}, "CSHL051": {"seeds": ["303", "310", "306", "302", "309", "305", "313", "308", "300", "314", "311", "307", "312", "304", "315", "301"], "fit_nums": ["69", "186", "49", "435", "103", "910", "705", "367", "303", "474", "596", "334", "929", "796", "616", "790"], "chain_num": 4}, "GLM_Sim_14": {"seeds": ["310", "311", "309", "313", "314", "300", "302", "304", "305", "306", "307", "312", "303", "301", "315", "308"], "fit_nums": ["616", "872", "419", "106", "940", "986", "599", "704", "218", "808", "244", "825", "448", "397", "552", "316"], "chain_num": 2}, "GLM_Sim_11_trick": {"seeds": ["411", "400", "408", "409", "415", "413", "410", "412", "406", "414", "403", "404", "401", "405", "407", "402"], "fit_nums": ["95", "508", "886", "384", "822", "969", "525", "382", "489", "436", "344", "537", "251", "223", "458", "401"], "chain_num": 2, "ignore": [10, 12, 4, 1, 0, 3, 2, 6]}, "GLM_Sim_16": {"seeds": ["302", "311", "303", "307", "313", "308", "309", "300", "305", "315", "304", "310", "312", "301", "314", "306"], "fit_nums": ["914", "377", "173", "583", "870", "456", "611", "697", "13", "713", "159", "248", "617", "37", "770", "780"], "chain_num": 2}, "ZM_3003": {"seeds": ["300", "304", "307", "312", "305", "310", "311", "314", "303", "308", "313", "301", "315", "309", "306", "302"], "fit_nums": ["603", "620", "657", "735", "357", "390", "119", "33", "62", "617", "209", "810", "688", "21", "744", "426"], "chain_num": 4}, "CSH_ZAD_022": {"seeds": ["305", "310", "311", "315", "303", "312", "314", "313", "307", "302", "300", "304", "301", "308", "306", "309"], "fit_nums": ["143", "946", "596", "203", "576", "403", "900", "65", "478", "325", "282", "513", "460", "42", "161", "970"], "chain_num": 4}, "GLM_Sim_07": {"seeds": ["300", "309", "302", "304", "305", "312", "301", "311", "315", "314", "308", "307", "303", "310", "306", "313"], "fit_nums": ["724", "701", "118", "230", "648", "426", "689", "114", "832", "731", "592", "519", "559", "938", "672", "144"], "chain_num": 1}, "KS017": {"seeds": ["311", "310", "306", "309", "303", "302", "308", "300", "313", "301", "314", "307", "315", "304", "312", "305"], "fit_nums": ["97", "281", "808", "443", "352", "890", "703", "468", "780", "708", "674", "27", "345", "23", "939", "457"], "chain_num": 4}, "GLM_Sim_11_sub": {"seeds": ["410", "414", "413", "404", "409", "415", "406", "408", "402", "411", "400", "405", "403", "407", "412", "401"], "fit_nums": ["830", "577", "701", "468", "929", "374", "954", "749", "937", "488", "873", "416", "612", "792", "461", "488"], "chain_num": 2}}
\ No newline at end of file
+{"SWC_023": {"seeds": ["302", "312", "304", "300", "315", "311", "308", "305", "303", "309", "306", "313", "307", "314", "301", "310"], "fit_nums": ["994", "913", "681", "816", "972", "790", "142", "230", "696", "537", "975", "773", "918", "677", "742", "745"], "chain_num": 4}, "SWC_021": {"seeds": ["415", "403", "412", "407", "409", "408", "405", "404", "410", "414", "401", "413", "402", "400", "406", "411"], "fit_nums": ["773", "615", "107", "583", "564", "354", "142", "184", "549", "185", "924", "907", "105", "531", "9", "812"], "chain_num": 9}, "ibl_witten_15": {"seeds": ["409", "410", "401", "415", "414", "403", "411", "404", "402", "405", "400", "412", "408", "407", "406", "413"], "fit_nums": ["411", "344", "496", "600", "716", "18", "527", "467", "898", "334", "309", "326", "133", "823", "740", "253"], "chain_num": 9}, "ibl_witten_13": {"seeds": ["302", "312", "313", "306", "315", "307", "311", "314", "309", "301", "308", "300", "304", "310", "303", "305"], "fit_nums": ["897", "765", "433", "641", "967", "599", "984", "259", "853", "385", "887", "619", "434", "964", "483", "891"], "chain_num": 4}, "KS016": {"seeds": ["315", "301", "309", "313", "302", "307", "303", "308", "311", "312", "314", "306", "310", "300", "305", "304"], "fit_nums": ["99", "57", "585", "32", "501", "558", "243", "413", "59", "757", "463", "172", "524", "957", "909", "292"], "chain_num": 4}, "KS003": {"seeds": ["404", "407", "413", "403", "414", "405", "400", "401", "402", "410", "415", "408", "411", "409", "406", "412"], "fit_nums": ["846", "256", "845", "945", "293", "406", "420", "109", "690", "421", "54", "866", "784", "81", "997", "665"], "chain_num": 9}, "ibl_witten_19": {"seeds": ["315", "311", "307", "314", "308", "300", "305", "301", "313", "304", "302", "310", "306", "312", "309", "303"], "fit_nums": ["179", "951", "613", "6", "623", "382", "458", "504", "406", "554", "5", "631", "746", "817", "265", "328"], "chain_num": 4}, "SWC_022": {"seeds": ["411", "403", "414", "409", "407", "412", "410", "413", "415", "404", "405", "400", "402", "401", "408", "406"], "fit_nums": ["408", "884", "62", "962", "744", "854", "635", "70", "320", "952", "8", "67", "231", "381", "536", "962"], "chain_num": 9}, "KS022": {"seeds": ["315", "300", "314", "301", "303", "302", "306", "308", "305", "310", "313", "312", "304", "307", "311", "309"], "fit_nums": ["899", "681", "37", "957", "629", "637", "375", "980", "810", "51", "759", "664", "420", "127", "259", "555"], "chain_num": 4}, "CSH_ZAD_017": {"seeds": ["401", "409", "405", "403", "415", "404", "402", "411", "410", "414", "408", "406", "413", "412", "400", "407"], "fit_nums": ["883", "803", "637", "806", "356", "804", "662", "654", "684", "350", "947", "460", "569", "976", "103", "713"], "chain_num": 9}, "CSH_ZAD_025": {"seeds": ["303", "311", "307", "312", "313", "314", "308", "315", "305", "306", "304", "302", "309", "310", "301", "300"], "fit_nums": ["581", "148", "252", "236", "581", "838", "206", "756", "449", "288", "756", "593", "733", "633", "418", "563"], "chain_num": 4}, "ibl_witten_17": {"seeds": ["406", "415", "408", "413", "402", "405", "409", "400", "414", "401", "412", "407", "404", "410", "403", "411"], "fit_nums": ["827", "797", "496", "6", "444", "823", "384", "873", "634", "27", "811", "142", "207", "322", "756", "275"], "chain_num": 9}, "ibl_witten_18": {"seeds": ["311", "310", "303", "314", "302", "309", "305", "307", "312", "300", "308", "306", "315", "313", "304", "301"], "fit_nums": ["236", "26", "838", "762", "826", "409", "496", "944", "280", "704", "930", "419", "637", "896", "876", "297"], "chain_num": 4}, "CSHL_018": {"seeds": ["302", "310", "306", "300", "314", "307", "309", "313", "311", "308", "304", "301", "312", "303", "305", "315"], "fit_nums": ["843", "817", "920", "900", "226", "36", "472", "676", "933", "453", "116", "263", "269", "897", "568", "438"], "chain_num": 4}, "GLM_Sim_06": {"seeds": ["313", "309", "302", "303", "305", "314", "300", "315", "311", "306", "304", "310", "301", "312", "308", "307"], "fit_nums": ["9", "786", "286", "280", "72", "587", "619", "708", "360", "619", "311", "189", "60", "708", "939", "733"], "chain_num": 2}, "ZM_1897": {"seeds": ["304", "308", "305", "311", "315", "314", "307", "306", "300", "303", "313", "310", "301", "312", "302", "309"], "fit_nums": ["549", "96", "368", "509", "424", "897", "287", "426", "968", "93", "725", "513", "837", "581", "989", "374"], "chain_num": 4}, "CSHL_020": {"seeds": ["305", "309", "313", "302", "314", "310", "300", "307", "315", "306", "312", "304", "311", "301", "303", "308"], "fit_nums": ["222", "306", "243", "229", "584", "471", "894", "238", "986", "660", "494", "657", "896", "459", "100", "283"], "chain_num": 4}, "CSHL054": {"seeds": ["401", "415", "409", "410", "414", "413", "407", "405", "406", "408", "411", "400", "412", "402", "403", "404"], "fit_nums": ["901", "734", "609", "459", "574", "793", "978", "66", "954", "906", "954", "111", "292", "850", "266", "967"], "chain_num": 9}, "CSHL_014": {"seeds": ["305", "311", "309", "300", "313", "310", "307", "306", "304", "312", "308", "302", "314", "303", "301", "315"], "fit_nums": ["371", "550", "166", "24", "705", "385", "870", "884", "831", "546", "404", "722", "287", "564", "613", "783"], "chain_num": 4}, "CSHL062": {"seeds": ["307", "313", "310", "303", "306", "312", "308", "305", "311", "314", "304", "302", "300", "301", "315", "309"], "fit_nums": ["846", "371", "94", "888", "499", "229", "546", "432", "71", "989", "986", "91", "935", "314", "975", "481"], "chain_num": 4}, "CSH_ZAD_001": {"seeds": ["313", "309", "311", "312", "305", "310", "315", "300", "314", "304", "301", "302", "308", "303", "306", "307"], "fit_nums": ["468", "343", "314", "544", "38", "120", "916", "170", "305", "569", "502", "496", "452", "336", "559", "572"], "chain_num": 4}, "NYU-06": {"seeds": ["314", "309", "306", "305", "312", "303", "307", "304", "300", "302", "310", "301", "315", "308", "313", "311"], "fit_nums": ["950", "862", "782", "718", "427", "645", "827", "612", "821", "834", "595", "929", "679", "668", "648", "869"], "chain_num": 4}, "KS019": {"seeds": ["404", "401", "411", "408", "400", "403", "410", "413", "402", "407", "415", "409", "406", "414", "412", "405"], "fit_nums": ["682", "4", "264", "200", "250", "267", "737", "703", "132", "855", "922", "686", "85", "176", "54", "366"], "chain_num": 9}, "CSHL049": {"seeds": ["411", "402", "414", "408", "409", "410", "413", "407", "406", "401", "404", "405", "403", "415", "400", "412"], "fit_nums": ["104", "553", "360", "824", "749", "519", "347", "228", "863", "671", "140", "883", "701", "445", "627", "898"], "chain_num": 9}, "ibl_witten_14": {"seeds": ["310", "311", "304", "306", "300", "302", "314", "313", "303", "308", "301", "309", "305", "315", "312", "307"], "fit_nums": ["563", "120", "85", "712", "277", "871", "183", "661", "505", "598", "210", "89", "310", "638", "564", "998"], "chain_num": 4}, "KS014": {"seeds": ["301", "310", "302", "312", "313", "308", "307", "303", "305", "300", "314", "306", "311", "309", "304", "315"], "fit_nums": ["668", "32", "801", "193", "269", "296", "74", "24", "270", "916", "21", "250", "342", "451", "517", "293"], "chain_num": 4, "ignore": [9, 11, 0, 1, 14, 2, 12, 13]}, "CSHL059": {"seeds": ["306", "309", "300", "304", "314", "303", "315", "311", "313", "305", "301", "307", "302", "312", "310", "308"], "fit_nums": ["821", "963", "481", "999", "986", "45", "551", "605", "701", "201", "629", "261", "972", "407", "165", "9"], "chain_num": 4}, "GLM_Sim_13": {"seeds": ["310", "303", "308", "306", "300", "312", "301", "313", "305", "311", "315", "304", "314", "309", "307", "302"], "fit_nums": ["982", "103", "742", "524", "614", "370", "926", "456", "133", "143", "302", "80", "395", "549", "579", "944"], "chain_num": 2}, "CSHL_007": {"seeds": ["314", "303", "308", "313", "301", "300", "302", "305", "315", "306", "310", "309", "311", "304", "307", "312"], "fit_nums": ["462", "703", "345", "286", "480", "313", "986", "165", "201", "102", "322", "894", "960", "438", "330", "169"], "chain_num": 4}, "CSH_ZAD_011": {"seeds": ["314", "311", "303", "300", "305", "310", "306", "301", "302", "315", "304", "309", "308", "312", "313", "307"], "fit_nums": ["320", "385", "984", "897", "315", "120", "320", "945", "475", "403", "210", "412", "695", "564", "664", "411"], "chain_num": 4}, "KS021": {"seeds": ["309", "312", "304", "310", "303", "311", "314", "302", "305", "301", "306", "300", "308", "315", "313", "307"], "fit_nums": ["874", "943", "925", "587", "55", "136", "549", "528", "349", "211", "401", "84", "225", "545", "153", "382"], "chain_num": 4}, "GLM_Sim_15": {"seeds": ["303", "312", "305", "308", "309", "302", "301", "310", "313", "315", "311", "314", "307", "306", "304", "300"], "fit_nums": ["769", "930", "328", "847", "899", "714", "144", "518", "521", "873", "914", "359", "242", "343", "45", "364"], "chain_num": 2}, "CSHL_015": {"seeds": ["301", "302", "307", "310", "309", "311", "304", "312", "300", "308", "313", "305", "314", "315", "306", "303"], "fit_nums": ["717", "705", "357", "539", "604", "971", "669", "76", "45", "413", "510", "122", "190", "821", "368", "472"], "chain_num": 4}, "ibl_witten_16": {"seeds": ["304", "313", "309", "314", "312", "307", "305", "301", "306", "310", "300", "315", "308", "311", "303", "302"], "fit_nums": ["392", "515", "696", "270", "7", "583", "880", "674", "23", "576", "579", "695", "149", "854", "184", "875"], "chain_num": 4}, "KS015": {"seeds": ["315", "305", "309", "303", "314", "310", "311", "312", "313", "300", "307", "308", "304", "301", "302", "306"], "fit_nums": ["257", "396", "387", "435", "133", "164", "403", "8", "891", "650", "111", "557", "473", "229", "842", "196"], "chain_num": 4}, "GLM_Sim_12": {"seeds": ["304", "312", "306", "303", "310", "302", "300", "305", "308", "313", "307", "311", "315", "301", "314", "309"], "fit_nums": ["971", "550", "255", "195", "952", "486", "841", "535", "559", "37", "654", "213", "864", "506", "732", "550"], "chain_num": 2}, "GLM_Sim_11": {"seeds": ["300", "312", "310", "315", "302", "313", "314", "311", "308", "303", "309", "307", "306", "304", "301", "305"], "fit_nums": ["477", "411", "34", "893", "195", "293", "603", "5", "887", "281", "956", "73", "346", "640", "532", "688"], "chain_num": 2}, "GLM_Sim_10": {"seeds": ["301", "300", "306", "305", "307", "309", "312", "314", "311", "315", "304", "313", "303", "308", "302", "310"], "fit_nums": ["391", "97", "897", "631", "239", "652", "19", "448", "807", "35", "972", "469", "280", "562", "42", "706"], "chain_num": 2}, "CSH_ZAD_026": {"seeds": ["312", "313", "308", "310", "303", "307", "302", "305", "300", "315", "306", "301", "311", "304", "314", "309"], "fit_nums": ["699", "87", "537", "628", "797", "511", "459", "770", "969", "240", "504", "948", "295", "506", "25", "378"], "chain_num": 4}, "KS023": {"seeds": ["304", "313", "306", "309", "300", "314", "302", "310", "303", "315", "307", "308", "301", "311", "305", "312"], "fit_nums": ["698", "845", "319", "734", "908", "507", "45", "499", "175", "108", "419", "443", "116", "779", "159", "231"], "chain_num": 4}, "GLM_Sim_05": {"seeds": ["301", "315", "300", "302", "305", "304", "313", "314", "311", "309", "306", "307", "308", "310", "303", "312"], "fit_nums": ["425", "231", "701", "375", "343", "902", "623", "125", "921", "637", "393", "964", "678", "930", "796", "42"], "chain_num": 2}, "CSHL061": {"seeds": ["305", "315", "304", "303", "309", "310", "302", "300", "314", "306", "311", "313", "301", "308", "307", "312"], "fit_nums": ["396", "397", "594", "911", "308", "453", "686", "552", "103", "209", "128", "892", "345", "925", "777", "396"], "chain_num": 4}, "CSHL051": {"seeds": ["303", "310", "306", "302", "309", "305", "313", "308", "300", "314", "311", "307", "312", "304", "315", "301"], "fit_nums": ["69", "186", "49", "435", "103", "910", "705", "367", "303", "474", "596", "334", "929", "796", "616", "790"], "chain_num": 4}, "GLM_Sim_14": {"seeds": ["310", "311", "309", "313", "314", "300", "302", "304", "305", "306", "307", "312", "303", "301", "315", "308"], "fit_nums": ["616", "872", "419", "106", "940", "986", "599", "704", "218", "808", "244", "825", "448", "397", "552", "316"], "chain_num": 2}, "GLM_Sim_11_trick": {"seeds": ["411", "400", "408", "409", "415", "413", "410", "412", "406", "414", "403", "404", "401", "405", "407", "402"], "fit_nums": ["95", "508", "886", "384", "822", "969", "525", "382", "489", "436", "344", "537", "251", "223", "458", "401"], "chain_num": 2, "ignore": [10, 12, 4, 1, 0, 3, 2, 6]}, "GLM_Sim_16": {"seeds": ["302", "311", "303", "307", "313", "308", "309", "300", "305", "315", "304", "310", "312", "301", "314", "306"], "fit_nums": ["914", "377", "173", "583", "870", "456", "611", "697", "13", "713", "159", "248", "617", "37", "770", "780"], "chain_num": 2}, "ZM_3003": {"seeds": ["300", "304", "307", "312", "305", "310", "311", "314", "303", "308", "313", "301", "315", "309", "306", "302"], "fit_nums": ["603", "620", "657", "735", "357", "390", "119", "33", "62", "617", "209", "810", "688", "21", "744", "426"], "chain_num": 4}, "CSH_ZAD_022": {"seeds": ["305", "310", "311", "315", "303", "312", "314", "313", "307", "302", "300", "304", "301", "308", "306", "309"], "fit_nums": ["143", "946", "596", "203", "576", "403", "900", "65", "478", "325", "282", "513", "460", "42", "161", "970"], "chain_num": 4}, "GLM_Sim_07": {"seeds": ["300", "309", "302", "304", "305", "312", "301", "311", "315", "314", "308", "307", "303", "310", "306", "313"], "fit_nums": ["724", "701", "118", "230", "648", "426", "689", "114", "832", "731", "592", "519", "559", "938", "672", "144"], "chain_num": 1}, "KS017": {"seeds": ["311", "310", "306", "309", "303", "302", "308", "300", "313", "301", "314", "307", "315", "304", "312", "305"], "fit_nums": ["97", "281", "808", "443", "352", "890", "703", "468", "780", "708", "674", "27", "345", "23", "939", "457"], "chain_num": 4}, "GLM_Sim_11_sub": {"seeds": ["410", "414", "413", "404", "409", "415", "406", "408", "402", "411", "400", "405", "403", "407", "412", "401"], "fit_nums": ["830", "577", "701", "468", "929", "374", "954", "749", "937", "488", "873", "416", "612", "792", "461", "488"], "chain_num": 2}}
\ No newline at end of file
diff --git a/dyn_glm_chain_analysis.py b/dyn_glm_chain_analysis.py
index 3fefa03f..7d57a217 100644
--- a/dyn_glm_chain_analysis.py
+++ b/dyn_glm_chain_analysis.py
@@ -1,5 +1,4 @@
 import matplotlib
-# matplotlib.use('Agg')
 import os
 os.environ["OMP_NUM_THREADS"] = "6" # export OMP_NUM_THREADS=4
 os.environ["OPENBLAS_NUM_THREADS"] = "6" # export OPENBLAS_NUM_THREADS=4
@@ -12,7 +11,7 @@ import pyhsmm
 import pickle
 import seaborn as sns
 import sys
-from scipy.stats import rankdata, norm, zscore
+from scipy.stats import zscore
 from scipy.optimize import minimize
 from itertools import combinations, product
 import matplotlib.gridspec as gridspec
@@ -20,6 +19,7 @@ from matplotlib.ticker import MaxNLocator
 import json
 import time
 import multiprocessing as mp
+from mcmc_chain_analysis import state_size_helper, state_num_helper, gamma_func, alpha_func, ll_func, r_hat_array_comp, rank_inv_normal_transform, eval_r_hat, eval_simple_r_hat
 
 
 colors = np.genfromtxt('colors.csv', delimiter=',')
@@ -69,7 +69,7 @@ class MCMC_result_list:
     def state_num_dist(self):
         state_nums = np.zeros((self.m, self.n))
         for i in range(self.m):
-            state_nums[i] = state_num_func(self.results[i])
+            state_nums[i] = state_num_helper(0.05)(self.results[i])
         return state_nums
 
     def return_chains(self, func, model_for_loop, rank_norm=True, mode_indices=None):
@@ -149,14 +149,11 @@ class MCMC_result_list:
                 counter += 1
                 step = len(ind) // min_len
                 up_to = min_len * step
-                # print([j + i * 640 for j in ind[:up_to:step]])
                 chains[counter] = func(self.results[i], ind[:up_to:step])
-            # print(chains)
-        # print(np.unique(chains.flatten(), return_counts=True))
+
         self.chains = chains
-        self.folded_chains = np.abs(self.chains - np.median(chains))
 
-        self.rank_inv_normal_transform()
+        self.rank_normalised, self.folded_rank_normalised, self.ranked, self.folded_ranked = rank_inv_normal_transform(self.chains)
 
         self.lame_r_hat, self.lame_var_hat_plus = r_hat_array_comp(self.chains)
         self.rank_normalised_r_hat, self.var_hat_plus = r_hat_array_comp(self.rank_normalised)
@@ -178,15 +175,6 @@ class MCMC_result_list:
         print("Effective number of samples is {}".format(self.n_eff))
         return chains
 
-    def rank_inv_normal_transform(self):
-        # Gelman paper Rank-normalization, folding, and localization: An improved R_hat for assessing convergence of MCMC
-        # ranking with average rank for ties
-        self.ranked = rankdata(self.chains).reshape(self.chains.shape)
-        self.folded_ranked = rankdata(self.folded_chains).reshape(self.folded_chains.shape)
-        # inverse normal with fractional offset
-        self.rank_normalised = norm.ppf((self.ranked - 3/8) / (self.chains.size + 1/4))
-        self.folded_rank_normalised = norm.ppf((self.folded_ranked - 3/8) / (self.folded_chains.size + 1/4))
-
     def rank_histos(self):
         count_max = 0
         for i in range(self.m):
@@ -385,42 +373,6 @@ def state_appear_and_dur(m, n_sessions):
     return state_appear, [(1 + session_dict[s][1] - session_dict[s][0]) / n_sessions for s in observed_states]
 
 
-def r_hat_array_comp(chains):
-    m, n = chains.shape
-    psi_dot_j = np.mean(chains, axis=1)
-    psi_dot_dot = np.mean(psi_dot_j)
-    B = n / (m - 1) * np.sum((psi_dot_j - psi_dot_dot) ** 2)
-    s_j_squared = np.sum((chains - psi_dot_j[:, None]) ** 2, axis=1) / (n - 1)
-    W = np.mean(s_j_squared)
-    var_hat_plus = (n - 1) / n * W + B / n
-    if W == 0:
-        # print("all the same value")
-        return 1, 0
-    r_hat = np.sqrt(var_hat_plus / W)
-    return r_hat, var_hat_plus
-
-
-def eval_amortized_r_hat(chains, psi_dot_j, s_j_squared, l, m, n):
-    psi_dot_dot = np.mean(psi_dot_j, axis=1)
-    B = n / (m - 1) * np.sum((psi_dot_j - psi_dot_dot[:, None]) ** 2, axis=1)
-    W = np.mean(s_j_squared, axis=1)
-    var_hat_plus = (n - 1) / n * W + B / n
-    r_hat = np.sqrt(var_hat_plus / W)
-    return max(r_hat)
-
-
-def r_hat_array_comp_mult(chains):
-    l, m, n = chains.shape
-    psi_dot_j = np.mean(chains, axis=2)
-    psi_dot_dot = np.mean(psi_dot_j, axis=1)
-    B = n / (m - 1) * np.sum((psi_dot_j - psi_dot_dot[:, None]) ** 2, axis=1)
-    s_j_squared = np.sum((chains - psi_dot_j[:, :, None]) ** 2, axis=2) / (n - 1)
-    W = np.mean(s_j_squared, axis=1)
-    var_hat_plus = (n - 1) / n * W + B / n
-    r_hat = np.sqrt(var_hat_plus / W)
-    return r_hat, var_hat_plus
-
-
 class MCMC_result:
     def __init__(self, models, infos, data, sessions, fit_variance, save_id, seq_start=0, dur='yes', sample_lls=None):
 
@@ -533,21 +485,6 @@ class MCMC_result:
         return np.array(glm_weights)
 
 
-def gamma_func(x): return x.trans_distn.gamma
-
-
-def alpha_func(x): return x.trans_distn.alpha
-
-
-def largest_state_func(x): return x.assign_counts.max(1)
-
-
-def state_num_func(x): return ((x.assign_counts / x.n_datapoints) > 0.05).sum(1)
-
-
-def ll_func(x): return x.sample_lls[-x.n_samples:]
-
-
 def find_good_chains(chains, reduce_to=8):
     delete_n = chains.shape[0] // 2 - reduce_to
     mins = np.zeros(1 + delete_n)
@@ -943,6 +880,7 @@ def state_development(test, state_sets, indices, save=True, save_append='', show
     ax0.plot(1 + 0.25, 0.6 + 0.2, 'ko', ms=18)
     ax0.plot(1 + 0.25, 0.6 + 0.2, 'wo', ms=16.8)
     ax0.plot(1 + 0.25, 0.6 + 0.2, 'ko', ms=16.8, alpha=abs(num_to_cont[1]))
+
     if test.results[0].type != 'bias':
         current, counter = 0, 0
         for c in [2, 3, 4, 5]:
@@ -1095,12 +1033,30 @@ def state_development(test, state_sets, indices, save=True, save_append='', show
 
     durs, state_types = state_type_durs(states_by_session, all_pmfs)
     dur_counter = 1
+    contrast_intro_types = [0, 0, 0, 0]
+    state, when = np.where(states_by_session > 0.05)
+    introductions_by_stage = np.zeros(3)
+    covered_states = []
     for i, d in enumerate(durs):
         ax0.fill_between(range(dur_counter, 1 + dur_counter + d), 0.5, -0.5, color=type_colours[i], zorder=0, alpha=0.3)
         dur_counter += d
 
-    ax2.set_title('Psychometric\nfunction', size=16)
+        # find out during which state type which contrast was introduced
+        for j, contrast in enumerate([2, 3, 4, 5]):
+            if contrast_intro_types[j] != 0:
+                continue
+            if test.results[0].infos[contrast] + 1 < dur_counter:
+                contrast_intro_types[j] = i+1
+
+        # find out during which stage which state was introduced
+        for s in range(len(state_sets)):
+            if np.sum(state == s) == 0 or s in covered_states:
+                continue
+            if when[state == s][0] + 1 < dur_counter:
+                introductions_by_stage[i] += 1
+                covered_states.append(s)
 
+    ax2.set_title('Psychometric\nfunction', size=16)
     ax1.set_ylabel('Proportion of trials', size=28, labelpad=-20)
     ax0.set_ylabel('% correct', size=18)
     ax2.set_ylabel('Probability', size=26, labelpad=-20)
@@ -1142,68 +1098,13 @@ def state_development(test, state_sets, indices, save=True, save_append='', show
     else:
         plt.close()
 
-    return states_by_session, all_pmfs, durs, state_types
-
-
-def comp_multi_r_hat(chains, rank_normalised, folded_rank_normalised):
-    lame_r_hat, _ = r_hat_array_comp(chains)
-    rank_normalised_r_hat, _ = r_hat_array_comp(rank_normalised)
-    folded_rank_normalised_r_hat, _ = r_hat_array_comp(folded_rank_normalised)
-    return max(lame_r_hat, rank_normalised_r_hat, folded_rank_normalised_r_hat)
-
-
-def eval_r_hat(chains1, chains2, chains3, chains4):
-    rank_normalised1, folded_rank_normalised1 = rank_inv_normal_transform(chains1)
-    rank_normalised2, folded_rank_normalised2 = rank_inv_normal_transform(chains2)
-    rank_normalised3, folded_rank_normalised3 = rank_inv_normal_transform(chains3)
-    rank_normalised4, folded_rank_normalised4 = rank_inv_normal_transform(chains4)
-
-    r_hat1 = comp_multi_r_hat(chains1, rank_normalised1, folded_rank_normalised1)
-    r_hat2 = comp_multi_r_hat(chains2, rank_normalised2, folded_rank_normalised2)
-    r_hat3 = comp_multi_r_hat(chains3, rank_normalised3, folded_rank_normalised3)
-    r_hat4 = comp_multi_r_hat(chains4, rank_normalised4, folded_rank_normalised4)
-
-    return max(r_hat1, r_hat2, r_hat3, r_hat4)
-
-
-def eval_simple_r_hat(chains):
+    return states_by_session, all_pmfs, durs, state_types, contrast_intro_types, smart_divide(introductions_by_stage, np.array(durs))
 
-    r_hats, _ = r_hat_array_comp_mult(chains)
-
-    return max(r_hats)
-
-
-def find_good_chains_unsplit(chains1, chains2, chains3, chains4, reduce_to=8, simple=False):
-    delete_n = - reduce_to + chains1.shape[0]
-    mins = np.zeros(delete_n + 1)
-    n_chains = chains1.shape[0]
-    chains = np.stack([chains1, chains2, chains3, chains4])
-
-    print("Without removals: {}".format(eval_simple_r_hat(chains)))
-    if simple:
-        r_hat = eval_simple_r_hat(chains)
-    else:
-        r_hat = eval_r_hat(chains1, chains2, chains3, chains4)
-    mins[0] = r_hat
-
-    for i in range(delete_n):
-        print()
-        r_hat_min = 10
-        sol = 0
-        for x in combinations(range(n_chains), n_chains - 1 - i):
-            if simple:
-                r_hat = eval_simple_r_hat(np.delete(chains, x, axis=1))
-            else:
-                r_hat = eval_r_hat(np.delete(chains1, x, axis=0), np.delete(chains2, x, axis=0), np.delete(chains3, x, axis=0), np.delete(chains4, x, axis=0))
-            if r_hat < r_hat_min:
-                sol = x
-            r_hat_min = min(r_hat, r_hat_min)
-        print("Minimum is {} (removed {})".format(r_hat_min, i + 1))
-        sol = [i for i in range(32) if i not in sol]
-        print("Removed: {}".format(sol))
-        mins[i + 1] = r_hat_min
-
-    return sol, r_hat_min
+def smart_divide(a, b):
+    c = np.zeros_like(a)
+    d = np.logical_and(a == 0, b == 0)
+    c[~d] = a[~d] / b[~d]
+    return c
 
 
 def find_good_chains_unsplit_greedy(chains1, chains2, chains3, chains4, reduce_to=8, simple=False):
@@ -1212,8 +1113,8 @@ def find_good_chains_unsplit_greedy(chains1, chains2, chains3, chains4, reduce_t
     n_chains = chains1.shape[0]
     chains = np.stack([chains1, chains2, chains3, chains4])
 
-    print("Without removals: {}".format(eval_r_hat(chains1, chains2, chains3, chains4)))
-    r_hat = eval_r_hat(chains1, chains2, chains3, chains4)
+    r_hat = eval_r_hat([chains1, chains2, chains3, chains4])
+    print("Without removals: {}".format(r_hat))
     mins[0] = r_hat
 
     to_del = []
@@ -1224,7 +1125,7 @@ def find_good_chains_unsplit_greedy(chains1, chains2, chains3, chains4, reduce_t
             if x in to_del:
                 continue
             if not simple:
-                r_hat = eval_r_hat(np.delete(chains1, to_del + [x], axis=0), np.delete(chains2, to_del + [x], axis=0), np.delete(chains3, to_del + [x], axis=0), np.delete(chains4, to_del + [x], axis=0))
+                r_hat = eval_r_hat([np.delete(chains1, to_del + [x], axis=0), np.delete(chains2, to_del + [x], axis=0), np.delete(chains3, to_del + [x], axis=0), np.delete(chains4, to_del + [x], axis=0)])
             else:
                 r_hat = eval_simple_r_hat(np.delete(chains, to_del + [x], axis=1))
             if r_hat < r_hat_min:
@@ -1237,23 +1138,11 @@ def find_good_chains_unsplit_greedy(chains1, chains2, chains3, chains4, reduce_t
         mins[i + 1] = r_hat_min
 
     if simple:
-        r_hat_local = eval_r_hat(np.delete(chains1, to_del, axis=0), np.delete(chains2, to_del, axis=0), np.delete(chains3, to_del, axis=0), np.delete(chains4, to_del, axis=0))
+        r_hat_local = eval_r_hat([np.delete(chains1, to_del, axis=0), np.delete(chains2, to_del, axis=0), np.delete(chains3, to_del, axis=0), np.delete(chains4, to_del, axis=0)])
         print("Minimum over everything is {} (removed {})".format(r_hat_local, i + 1))
     return to_del, r_hat_min
 
 
-def rank_inv_normal_transform(chains):
-    # Gelman paper Rank-normalization, folding, and localization: An improved R_hat for assessing convergence of MCMC
-    # ranking with average rank for ties
-    folded_chains = np.abs(chains - np.median(chains))
-    ranked = rankdata(chains).reshape(chains.shape)
-    folded_ranked = rankdata(folded_chains).reshape(folded_chains.shape)
-    # inverse normal with fractional offset
-    rank_normalised = norm.ppf((ranked - 3/8) / (chains.size + 1/4))
-    folded_rank_normalised = norm.ppf((folded_ranked - 3/8) / (folded_chains.size + 1/4))
-    return rank_normalised, folded_rank_normalised
-
-
 if __name__ == "__main__":
     fit_type = ['prebias', 'bias', 'all', 'prebias_plus', 'zoe_style'][0]
     if fit_type == 'bias':
@@ -1268,16 +1157,13 @@ if __name__ == "__main__":
     # test.r_hat_and_ess(return_ascending, False)
     # test.r_hat_and_ess(return_ascending_shuffled, False)
     # quit()
-
-    import pyhsmm.util.profiling as prof
-
     good = []
     bad = []
 
-    check_r_hats = False
+    check_r_hats = True
     if check_r_hats:
         subjects = list(loading_info.keys())
-        subjects = ['GLM_Sim_11_trick']
+        subjects = ['KS014']
         for subject in subjects:
             # if subject.startswith('GLM'):
             #     continue
@@ -1322,24 +1208,6 @@ def dist_helper(dist_matrix, state_hists, inds):
     return dist_matrix
 
 
-def state_size_helper(n=0, mode_specific=False):
-    if not mode_specific:
-        def nth_largest_state_func(x):
-            return np.partition(x.assign_counts, -1 - n, axis=1)[:, -1 - n]
-    else:
-        def nth_largest_state_func(x, ind):
-            return np.partition(x.assign_counts[ind], -1 - n, axis=1)[:, -1 - n]
-    return nth_largest_state_func
-
-
-def state_num_helper(t, mode_specific=False):
-    if not mode_specific:
-        def state_num_func(x): return ((x.assign_counts / x.n_datapoints) > t).sum(1)
-    else:
-        def state_num_func(x, ind): return ((x.assign_counts[ind] / x.n_datapoints) > t).sum(1)
-    return state_num_func
-
-
 def state_glm_func_helper(t, mode_specific=False):
     if not mode_specific:
         def temp_state_glm_func(x):
@@ -1364,26 +1232,6 @@ def state_glm_func_helper(t, mode_specific=False):
     return temp_state_glm_func
 
 
-def sample_statistics(test, mode_indices):
-    # prints out r_hats and sample sizes for given sample
-    test = pickle.load(open("multi_chain_saves/canonical_result_{}.p".format(subject), 'rb'))
-    test.r_hat_and_ess(state_size_helper(1), False)
-    test.r_hat_and_ess(state_size_helper(1, mode_specific=True), False, mode_indices=mode_indices)
-    print()
-    test = pickle.load(open("multi_chain_saves/canonical_result_{}.p".format(subject), 'rb'))
-    test.r_hat_and_ess(state_size_helper(), False)
-    test.r_hat_and_ess(state_size_helper(mode_specific=True), False, mode_indices=mode_indices)
-    print()
-    test = pickle.load(open("multi_chain_saves/canonical_result_{}.p".format(subject), 'rb'))
-    test.r_hat_and_ess(state_num_helper(0.05), False)
-    test.r_hat_and_ess(state_num_helper(0.05, mode_specific=True), False, mode_indices=mode_indices)
-    print()
-    test = pickle.load(open("multi_chain_saves/canonical_result_{}.p".format(subject), 'rb'))
-    test.r_hat_and_ess(state_num_helper(0.02), False)
-    test.r_hat_and_ess(state_num_helper(0.02, mode_specific=True), False, mode_indices=mode_indices)
-    print()
-
-
 def state_type_durs(states, pmfs):
     # Takes states and pmfs, first creates an array of when which type is how active, then computes the number of sessions each type lasts.
     # A type lasts until a more advanced type takes up more than 50% of a session (and cannot return)
@@ -1452,6 +1300,22 @@ def pmf_type(pmf):
 
 
 if __name__ == "__main__":
+
+    # visualise pmf types
+    lapses = [0.1, 0.2, 0.25, 0.33, 0.4, 0.45, 0.5, 0.55, 0.66, 0.9]
+    test_pmf = np.zeros(4)
+    for i, lapse_l in enumerate(lapses):
+        plt.subplot(1, 10, 1+i)
+        if i != 0:
+            plt.gca().set_yticklabels([])
+        plt.ylim(0, 1)
+        test_pmf[:2] = lapse_l
+        for lapse_r in np.linspace(0.02, 0.98, 33):
+            test_pmf[2:] = lapse_r
+            plt.plot([0, 1, 9, 10], test_pmf, c=type_colours[pmf_type(test_pmf)])
+    plt.show()
+    quit()
+
     fit_type = ['prebias', 'bias', 'all', 'prebias_plus', 'zoe_style'][0]
     if fit_type == 'bias':
         loading_info = json.load(open("canonical_infos_bias.json", 'r'))
@@ -1505,6 +1369,8 @@ if __name__ == "__main__":
     num_sessions = []
     state_appear = []
     state_dur = []
+    contrast_intro_types = []  # list to agglomorate in which state type which contrast is introduced
+    intros_by_type_sum = np.zeros(3)  # array to agglomorate how many states where introduced during which type, normalised by length of phase
 
     n_points = 150
     state_trajs = np.zeros((3, n_points))
@@ -1530,12 +1396,14 @@ if __name__ == "__main__":
 
             mode_indices = pickle.load(open("multi_chain_saves/mode_indices_{}_{}.p".format(subject, fit_type), 'rb'))
             state_sets = pickle.load(open("multi_chain_saves/state_sets_{}_{}.p".format(subject, fit_type), 'rb'))
-
             # lapse differential
             # lapse_sides(test, [s for s in state_sets if len(s) > 40], mode_indices)
 
             # training overview
-            states, pmfs, durs, _ = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, show=0, separate_pmf=True)
+            states, pmfs, durs, _, contrast_intro_type, intros_by_type = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, show=0, separate_pmf=True)
+            intros_by_type_sum += intros_by_type
+            continue
+            contrast_intro_types.append(contrast_intro_type)
             # state_development_single_sample(test, [mode_indices[0]], show=True, separate_pmf=True, save=False)
 
             # session overview
@@ -1794,12 +1662,6 @@ if __name__ == "__main__":
 
         plotSimplex(np.array(abs_state_durs), c='k', show=True)
 
-        plt.hist(abs_state_durs, label=['1st type', '2nd type', '3rd type'])
-        plt.legend()
-        plt.tight_layout()
-        plt.savefig("rel state type hists")
-        plt.show()
-
     if False:
         ax[0].set_ylim(0, 1)
         ax[1].set_ylim(0, 1)
diff --git a/dyn_glm_chain_analysis_unused_funcs.py b/dyn_glm_chain_analysis_unused_funcs.py
index 17daace5..c2e51ba4 100644
--- a/dyn_glm_chain_analysis_unused_funcs.py
+++ b/dyn_glm_chain_analysis_unused_funcs.py
@@ -99,6 +99,40 @@ def find_good_chains_unsplit_fast(chains1, chains2, chains3, chains4, reduce_to=
     return sol, r_hat_min
 
 
+def find_good_chains_unsplit(chains1, chains2, chains3, chains4, reduce_to=8, simple=False):
+    delete_n = - reduce_to + chains1.shape[0]
+    mins = np.zeros(delete_n + 1)
+    n_chains = chains1.shape[0]
+    chains = np.stack([chains1, chains2, chains3, chains4])
+
+    print("Without removals: {}".format(eval_simple_r_hat(chains)))
+    if simple:
+        r_hat = eval_simple_r_hat(chains)
+    else:
+        r_hat = eval_r_hat(chains1, chains2, chains3, chains4)
+    mins[0] = r_hat
+
+    for i in range(delete_n):
+        print()
+        r_hat_min = 10
+        sol = 0
+        for x in combinations(range(n_chains), n_chains - 1 - i):
+            if simple:
+                r_hat = eval_simple_r_hat(np.delete(chains, x, axis=1))
+            else:
+                r_hat = eval_r_hat(np.delete(chains1, x, axis=0), np.delete(chains2, x, axis=0), np.delete(chains3, x, axis=0), np.delete(chains4, x, axis=0))
+            if r_hat < r_hat_min:
+                sol = x
+            r_hat_min = min(r_hat, r_hat_min)
+        print("Minimum is {} (removed {})".format(r_hat_min, i + 1))
+        sol = [i for i in range(32) if i not in sol]
+        print("Removed: {}".format(sol))
+        mins[i + 1] = r_hat_min
+
+    return sol, r_hat_min
+
+
+
 def params_to_pmf(params):
     return params[2] + (1 - params[2] - params[3]) / (1 + np.exp(- params[0] * (all_conts - params[1])))
 
@@ -214,7 +248,7 @@ if __name__ == "__main__":
         plt.savefig("pmf fit scatter")
         plt.show()
 
-
+        # New things
         xy = np.vstack([short_pmfs[:, 0], function_range, short_pmfs[:, 1]])
         z = gaussian_kde(xy)(xy)
         plt.figure(figsize=(24, 24 / 3))
diff --git a/mcmc_chain_analysis.py b/mcmc_chain_analysis.py
new file mode 100644
index 00000000..2a99a676
--- /dev/null
+++ b/mcmc_chain_analysis.py
@@ -0,0 +1,123 @@
+"""
+    Functions to extract statistics from a set of chains
+    Functions to compute R^hat from a set of statistic vectors
+"""
+import numpy as np
+from scipy.stats import rankdata, norm
+import pickle
+
+
+def state_size_helper(n=0, mode_specific=False):
+    if not mode_specific:
+        def nth_largest_state_func(x):
+            return np.partition(x.assign_counts, -1 - n, axis=1)[:, -1 - n]
+    else:
+        def nth_largest_state_func(x, ind):
+            return np.partition(x.assign_counts[ind], -1 - n, axis=1)[:, -1 - n]
+    return nth_largest_state_func
+
+
+def state_num_helper(t, mode_specific=False):
+    if not mode_specific:
+        def state_num_func(x): return ((x.assign_counts / x.n_datapoints) > t).sum(1)
+    else:
+        def state_num_func(x, ind): return ((x.assign_counts[ind] / x.n_datapoints) > t).sum(1)
+    return state_num_func
+
+
+def gamma_func(x): return x.trans_distn.gamma
+
+
+def alpha_func(x): return x.trans_distn.alpha
+
+
+def ll_func(x): return x.sample_lls[-x.n_samples:]
+
+
+def r_hat_array_comp(chains):
+    m, n = chains.shape
+    psi_dot_j = np.mean(chains, axis=1)
+    psi_dot_dot = np.mean(psi_dot_j)
+    B = n / (m - 1) * np.sum((psi_dot_j - psi_dot_dot) ** 2)
+    s_j_squared = np.sum((chains - psi_dot_j[:, None]) ** 2, axis=1) / (n - 1)
+    W = np.mean(s_j_squared)
+    var_hat_plus = (n - 1) / n * W + B / n
+    if W == 0:
+        # print("all the same value")
+        return 1, 0
+    r_hat = np.sqrt(var_hat_plus / W)
+    return r_hat, var_hat_plus
+
+
+def eval_amortized_r_hat(chains, psi_dot_j, s_j_squared, m, n):
+    psi_dot_dot = np.mean(psi_dot_j, axis=1)
+    B = n / (m - 1) * np.sum((psi_dot_j - psi_dot_dot[:, None]) ** 2, axis=1)
+    W = np.mean(s_j_squared, axis=1)
+    var_hat_plus = (n - 1) / n * W + B / n
+    r_hat = np.sqrt(var_hat_plus / W)
+    return max(r_hat)
+
+
+def r_hat_array_comp_mult(chains):
+    _, m, n = chains.shape
+    psi_dot_j = np.mean(chains, axis=2)
+    psi_dot_dot = np.mean(psi_dot_j, axis=1)
+    B = n / (m - 1) * np.sum((psi_dot_j - psi_dot_dot[:, None]) ** 2, axis=1)
+    s_j_squared = np.sum((chains - psi_dot_j[:, :, None]) ** 2, axis=2) / (n - 1)
+    W = np.mean(s_j_squared, axis=1)
+    var_hat_plus = (n - 1) / n * W + B / n
+    r_hat = np.sqrt(var_hat_plus / W)
+    return r_hat, var_hat_plus
+
+
+def rank_inv_normal_transform(chains):
+    # Gelman paper Rank-normalization, folding, and localization: An improved R_hat for assessing convergence of MCMC
+    # ranking with average rank for ties
+    folded_chains = np.abs(chains - np.median(chains))
+    ranked = rankdata(chains).reshape(chains.shape)
+    folded_ranked = rankdata(folded_chains).reshape(folded_chains.shape)
+    # inverse normal with fractional offset
+    rank_normalised = norm.ppf((ranked - 3/8) / (chains.size + 1/4))
+    folded_rank_normalised = norm.ppf((folded_ranked - 3/8) / (folded_chains.size + 1/4))
+    return rank_normalised, folded_rank_normalised, ranked, folded_ranked
+
+
+def eval_r_hat(chains):
+    r_hats = []
+    for chain in chains:
+        rank_normalised, folded_rank_normalised, _, _ = rank_inv_normal_transform(chain)
+        r_hats.append(comp_multi_r_hat(chain, rank_normalised, folded_rank_normalised))
+
+    return max(r_hats)
+
+
+def eval_simple_r_hat(chains):
+    r_hats, _ = r_hat_array_comp_mult(chains)
+    return max(r_hats)
+
+
+def comp_multi_r_hat(chains, rank_normalised, folded_rank_normalised):
+    lame_r_hat, _ = r_hat_array_comp(chains)
+    rank_normalised_r_hat, _ = r_hat_array_comp(rank_normalised)
+    folded_rank_normalised_r_hat, _ = r_hat_array_comp(folded_rank_normalised)
+    return max(lame_r_hat, rank_normalised_r_hat, folded_rank_normalised_r_hat)
+
+
+def sample_statistics(test, mode_indices, subject):
+    # prints out r_hats and sample sizes for given sample
+    test = pickle.load(open("multi_chain_saves/canonical_result_{}.p".format(subject), 'rb'))
+    test.r_hat_and_ess(state_size_helper(1), False)
+    test.r_hat_and_ess(state_size_helper(1, mode_specific=True), False, mode_indices=mode_indices)
+    print()
+    test = pickle.load(open("multi_chain_saves/canonical_result_{}.p".format(subject), 'rb'))
+    test.r_hat_and_ess(state_size_helper(), False)
+    test.r_hat_and_ess(state_size_helper(mode_specific=True), False, mode_indices=mode_indices)
+    print()
+    test = pickle.load(open("multi_chain_saves/canonical_result_{}.p".format(subject), 'rb'))
+    test.r_hat_and_ess(state_num_helper(0.05), False)
+    test.r_hat_and_ess(state_num_helper(0.05, mode_specific=True), False, mode_indices=mode_indices)
+    print()
+    test = pickle.load(open("multi_chain_saves/canonical_result_{}.p".format(subject), 'rb'))
+    test.r_hat_and_ess(state_num_helper(0.02), False)
+    test.r_hat_and_ess(state_num_helper(0.02, mode_specific=True), False, mode_indices=mode_indices)
+    print()
diff --git a/process_many_chains.py b/process_many_chains.py
index 7aba22fb..88840671 100644
--- a/process_many_chains.py
+++ b/process_many_chains.py
@@ -1,3 +1,7 @@
+"""
+    Script for taking a list of subects and extracting statistics from the chains
+    which can be used to assess which chains have converged to the same regions
+"""
 import numpy as np
 import pyhsmm
 import pickle
@@ -5,25 +9,7 @@ import json
 from dyn_glm_chain_analysis import MCMC_result
 import matplotlib.pyplot as plt
 import time
-
-def gamma_func(x): return x.trans_distn.gamma
-
-
-def alpha_func(x): return x.trans_distn.alpha
-
-
-def state_size_helper(n=0):
-    def nth_largest_state_func(x):
-        return np.partition(x.assign_counts, -1 - n, axis=1)[:, -1 - n]
-    return nth_largest_state_func
-
-
-def state_num_helper(t):
-    def state_num_func(x): return ((x.assign_counts / x.n_datapoints) > t).sum(1)
-    return state_num_func
-
-
-def ll_func(x): return x.sample_lls[-x.n_samples:]
+from mcmc_chain_analysis import state_size_helper, state_num_helper, ll_helper
 
 
 fit_type = ['prebias', 'bias', 'all', 'prebias_plus', 'zoe_style'][0]
@@ -45,117 +31,79 @@ func8 = state_size_helper(4)
 func9 = state_size_helper(5)
 func10 = state_size_helper(6)
 func11 = state_size_helper(7)
-model_for_loop = False
 dur = 'yes'
 
-def temp():
-    m = 16
-    for subject in subjects:
-        print(subject)
-        n_runs = -1
-        counter = -1
-        n = (loading_info[subject]['chain_num'] + 1) * 4000 // 25
-        chains1 = np.zeros((m, n))
-        chains2 = np.zeros((m, n))
-        chains3 = np.zeros((m, n))
-        chains4 = np.zeros((m, n))
-        for j, (seed, fit_num) in enumerate(zip(loading_info[subject]['seeds'], loading_info[subject]['fit_nums'])):
-            counter += 1
-            print(seed)
-            info_dict = pickle.load(open("./session_data/{}_info_dict.p".format(subject), "rb"))
-            samples = []
-            mini_counter = 0
-            while True:
-                try:
-                    file = "./dynamic_GLMiHMM_crossvals/{}_fittype_{}_var_{}_{}_{}{}.p".format(subject, fit_type, fit_variance, seed, fit_num, '_{}'.format(mini_counter))
-                    lala = time.time()
-                    samples += pickle.load(open(file, "rb"))
-                    print("Loading {} took {:.4}".format(mini_counter, time.time() - lala))
-                    mini_counter += 1
-                except Exception:
-                    break
-            if n_runs == -1:
-                n_runs = mini_counter
-            else:
-                if n_runs != mini_counter:
-                    print("Problem")
-                    print(n_runs, mini_counter)
-                    quit()
-            save_id = "{}_fittype_{}_var_{}_{}_{}.p".format(subject, fit_type, fit_variance, seed, fit_num).replace('.', '_')
-            # for file in os.listdir("./dynamic_GLMiHMM_crossvals/infos_new/"):
-            #     if file.startswith("{}_".format(subject)) and file.endswith("_{}_{}_{}_{}.json".format(fit_type, fit_variance, seed, fit_num)):
-            #         print("Taking fit infos from {}".format(file))
-            #         fit_infos = json.load(open("./dynamic_GLMiHMM_crossvals/infos_new/" + file, 'r'))
-            #         sample_lls = fit_infos['ll']
-
-            print("loaded seed {}".format(seed))
-
-            lala = time.time()
-            result = MCMC_result(samples[::25],
-                                 infos=info_dict, data=samples[0].datas,
-                                 sessions=fit_type, fit_variance=fit_variance,
-                                 dur=dur, save_id=save_id) # , sample_lls=sample_lls)
-            print("Making result {} took {:.4}".format(counter, time.time() - lala))
-
-            # if model_for_loop:
-            #     for i in range(n):
-            #         chains[j, i] = func(result.models[i])
-            res = func1(result)
-            chains1[counter] = res
-            res = func2(result)
-            chains2[counter] = res
-            res = func3(result)
-            chains3[counter] = res
-            res = func4(result)
-            chains4[counter] = res
-
-            # res = func5(result)
-            # chains5[j] = res
-            # res = func6(result)
-            # chains6[j] = res
-            # res = func7(result)
-            # chains7[j] = res
-            # res = func8(result)
-            # chains8[j] = res
-            # res = func9(result)
-            # chains9[j] = res
-            # res = func10(result)
-            # chains10[j] = res
-            # res = func11(result)
-            # chains11[j] = res
-
-        # func2 = state_size_helper()
-        # func5 = state_size_helper(1)
-        # func6 = state_size_helper(2)
-        # func7 = state_size_helper(3)
-        # func8 = state_size_helper(4)
-        # func9 = state_size_helper(5)
-        # func10 = state_size_helper(6)
-        # func11 = state_size_helper(7)
-        # plt.plot(chains2.flatten()[::10])
-        # plt.plot(chains5.flatten()[::10])
-        # plt.plot(chains6.flatten()[::10])
-        # plt.plot(chains7.flatten()[::10])
-        # plt.plot(chains8.flatten()[::10])
-        # plt.plot(chains9.flatten()[::10])
-        # plt.plot(chains10.flatten()[::10])
-        # plt.plot(chains11.flatten()[::10])
-        # plt.axhline(1643, color='r')
-        # plt.axhline(3438, color='r')
-        # plt.axhline(2216, color='r')
-        # plt.axhline(811, color='r')
-        # plt.axhline(2721, color='r')
-        # plt.axhline(743, color='r')
-        # plt.show()
-
-        pickle.dump(chains1, open("multi_chain_saves/{}_state_num_0_fittype_{}_var_{}_{}_{}_state_num.p".format(subject, fit_type, fit_variance, seed, fit_num), 'wb'))
-        pickle.dump(chains2, open("multi_chain_saves/{}_state_num_1_fittype_{}_var_{}_{}_{}_state_num.p".format(subject, fit_type, fit_variance, seed, fit_num), 'wb'))
-        pickle.dump(chains3, open("multi_chain_saves/{}_largest_state_0_fittype_{}_var_{}_{}_{}_state_num.p".format(subject, fit_type, fit_variance, seed, fit_num), 'wb'))
-        pickle.dump(chains4, open("multi_chain_saves/{}_largest_state_1_fittype_{}_var_{}_{}_{}_state_num.p".format(subject, fit_type, fit_variance, seed, fit_num), 'wb'))
-
-
-temp()
-# import pyhsmm.util.profiling as prof
-# prof_func = prof._prof(temp)
-# prof_func()
-# prof._prof.print_stats()
+m = 16
+for subject in subjects:
+    print(subject)
+    n_runs = -1
+    counter = -1
+    n = (loading_info[subject]['chain_num'] + 1) * 4000 // 25
+    chains1 = np.zeros((m, n))
+    chains2 = np.zeros((m, n))
+    chains3 = np.zeros((m, n))
+    chains4 = np.zeros((m, n))
+    for j, (seed, fit_num) in enumerate(zip(loading_info[subject]['seeds'], loading_info[subject]['fit_nums'])):
+        counter += 1
+        print(seed)
+        info_dict = pickle.load(open("./session_data/{}_info_dict.p".format(subject), "rb"))
+        samples = []
+        mini_counter = 0
+        while True:
+            try:
+                file = "./dynamic_GLMiHMM_crossvals/{}_fittype_{}_var_{}_{}_{}{}.p".format(subject, fit_type, fit_variance, seed, fit_num, '_{}'.format(mini_counter))
+                lala = time.time()
+                samples += pickle.load(open(file, "rb"))
+                print("Loading {} took {:.4}".format(mini_counter, time.time() - lala))
+                mini_counter += 1
+            except Exception:
+                break
+        if n_runs == -1:
+            n_runs = mini_counter
+        else:
+            if n_runs != mini_counter:
+                print("Problem")
+                print(n_runs, mini_counter)
+                quit()
+        save_id = "{}_fittype_{}_var_{}_{}_{}.p".format(subject, fit_type, fit_variance, seed, fit_num).replace('.', '_')
+        # for file in os.listdir("./dynamic_GLMiHMM_crossvals/infos_new/"):
+        #     if file.startswith("{}_".format(subject)) and file.endswith("_{}_{}_{}_{}.json".format(fit_type, fit_variance, seed, fit_num)):
+        #         print("Taking fit infos from {}".format(file))
+        #         fit_infos = json.load(open("./dynamic_GLMiHMM_crossvals/infos_new/" + file, 'r'))
+        #         sample_lls = fit_infos['ll']
+
+        print("loaded seed {}".format(seed))
+
+        result = MCMC_result(samples[::25],
+                             infos=info_dict, data=samples[0].datas,
+                             sessions=fit_type, fit_variance=fit_variance,
+                             dur=dur, save_id=save_id) # , sample_lls=sample_lls)
+        print("Making result {} took {:.4}".format(counter, time.time() - lala))
+
+        res = func1(result)
+        chains1[counter] = res
+        res = func2(result)
+        chains2[counter] = res
+        res = func3(result)
+        chains3[counter] = res
+        res = func4(result)
+        chains4[counter] = res
+        # res = func5(result)
+        # chains5[j] = res
+        # res = func6(result)
+        # chains6[j] = res
+        # res = func7(result)
+        # chains7[j] = res
+        # res = func8(result)
+        # chains8[j] = res
+        # res = func9(result)
+        # chains9[j] = res
+        # res = func10(result)
+        # chains10[j] = res
+        # res = func11(result)
+        # chains11[j] = res
+
+    pickle.dump(chains1, open("multi_chain_saves/{}_state_num_0_fittype_{}_var_{}_{}_{}_state_num.p".format(subject, fit_type, fit_variance, seed, fit_num), 'wb'))
+    pickle.dump(chains2, open("multi_chain_saves/{}_state_num_1_fittype_{}_var_{}_{}_{}_state_num.p".format(subject, fit_type, fit_variance, seed, fit_num), 'wb'))
+    pickle.dump(chains3, open("multi_chain_saves/{}_largest_state_0_fittype_{}_var_{}_{}_{}_state_num.p".format(subject, fit_type, fit_variance, seed, fit_num), 'wb'))
+    pickle.dump(chains4, open("multi_chain_saves/{}_largest_state_1_fittype_{}_var_{}_{}_{}_state_num.p".format(subject, fit_type, fit_variance, seed, fit_num), 'wb'))
diff --git a/simplex_plot.py b/simplex_plot.py
index 41db6f5e..bdaad788 100644
--- a/simplex_plot.py
+++ b/simplex_plot.py
@@ -4,7 +4,7 @@ Visualize points on the 3-simplex (eg, the parameters of a
 contained within a 2D triangle.
 Adapted from David Andrzejewski (david.andrzej@gmail.com)
 """
-import numpy as NP
+import numpy as np
 import matplotlib.pyplot as P
 import matplotlib.ticker as MT
 import matplotlib.lines as L
@@ -14,38 +14,42 @@ import matplotlib.patches as PA
 
 
 def plotSimplex(points, fig=None,
-                vertexlabels=['1', '2', '3'], save_title="test.png",
+                vertexlabels=['Type 1', 'Type 2', 'Type 3'], save_title="test.png",
                 show=False, **kwargs):
     """
     Plot Nx3 points array on the 3-simplex
     (with optionally labeled vertices)
 
     kwargs will be passed along directly to matplotlib.pyplot.scatter
-    Returns Figure, caller must .show()
     """
     if fig is None:
         fig = P.figure(figsize=(9, 9))
     # Draw the triangle
     l1 = L.Line2D([0, 0.5, 1.0, 0], # xcoords
-                  [0, NP.sqrt(3) / 2, 0, 0], # ycoords
+                  [0, np.sqrt(3) / 2, 0, 0], # ycoords
                   color='k')
     fig.gca().add_line(l1)
     fig.gca().xaxis.set_major_locator(MT.NullLocator())
     fig.gca().yaxis.set_major_locator(MT.NullLocator())
     # Draw vertex labels
-    fig.gca().text(-0.05, -0.05, vertexlabels[0])
-    fig.gca().text(1.05, -0.05, vertexlabels[1])
-    fig.gca().text(0.5, NP.sqrt(3) / 2 + 0.05, vertexlabels[2])
+    fig.gca().text(-0.06, -0.05, vertexlabels[0], size=24)
+    fig.gca().text(0.95, -0.05, vertexlabels[1], size=24)
+    fig.gca().text(0.43, np.sqrt(3) / 2 + 0.025, vertexlabels[2], size=24)
     # Project and draw the actual points
     projected = projectSimplex(points / points.sum(1)[:, None])
-    P.scatter(projected[:, 0], projected[:, 1], s=points.sum(1), **kwargs)
+    P.scatter(projected[:, 0], projected[:, 1], s=points.sum(1) * 3.5, **kwargs)
+
+    # plot center with average size
+    projected = projectSimplex(np.mean(points / points.sum(1)[:, None], axis=0).reshape(1, 3))
+    P.scatter(projected[:, 0], projected[:, 1], marker='*', color='r', s=np.mean(points.sum(1)) * 3.5)
+
     # Leave some buffer around the triangle for vertex labels
     fig.gca().set_xlim(-0.05, 1.05)
     fig.gca().set_ylim(-0.05, 1.05)
 
     P.axis('off')
 
-    P.savefig("test.png", bbox_inches='tight')
+    P.savefig("dur_simplex.png", bbox_inches='tight')
     if show:
         P.show()
     else:
@@ -59,22 +63,22 @@ def projectSimplex(points):
     N points are given as N x 3 array
     """
     # Convert points one at a time
-    tripts = NP.zeros((points.shape[0], 2))
+    tripts = np.zeros((points.shape[0], 2))
     for idx in range(points.shape[0]):
         # Init to triangle centroid
         x = 1.0 / 2
-        y = 1.0 / (2 * NP.sqrt(3))
+        y = 1.0 / (2 * np.sqrt(3))
         # Vector 1 - bisect out of lower left vertex
         p1 = points[idx, 0]
-        x = x - (1.0 / NP.sqrt(3)) * p1 * NP.cos(NP.pi / 6)
-        y = y - (1.0 / NP.sqrt(3)) * p1 * NP.sin(NP.pi / 6)
+        x = x - (1.0 / np.sqrt(3)) * p1 * np.cos(np.pi / 6)
+        y = y - (1.0 / np.sqrt(3)) * p1 * np.sin(np.pi / 6)
         # Vector 2 - bisect out of lower right vertex
         p2 = points[idx, 1]
-        x = x + (1.0 / NP.sqrt(3)) * p2 * NP.cos(NP.pi / 6)
-        y = y - (1.0 / NP.sqrt(3)) * p2 * NP.sin(NP.pi / 6)
+        x = x + (1.0 / np.sqrt(3)) * p2 * np.cos(np.pi / 6)
+        y = y - (1.0 / np.sqrt(3)) * p2 * np.sin(np.pi / 6)
         # Vector 3 - bisect out of top vertex
         p3 = points[idx, 2]
-        y = y + (1.0 / NP.sqrt(3) * p3)
+        y = y + (1.0 / np.sqrt(3) * p3)
 
         tripts[idx, :] = (x, y)
 
@@ -87,7 +91,7 @@ if __name__ == '__main__':
               '[0.8  0.1  0.1]',
               '[0.5  0.4  0.1]',
               '[0.33  0.34  0.33]')
-    testpoints = NP.array([[0.1, 0.1, 0.8],
+    testpoints = np.array([[0.1, 0.1, 0.8],
                            [0.8, 0.1, 0.1],
                            [0.5, 0.4, 0.1],
                            [0.33, 0.34, 0.33]])
-- 
GitLab