From edf0f7a39190ecdc1eedac47b96c98899695d8d2 Mon Sep 17 00:00:00 2001
From: SebastianBruijns <>
Date: Fri, 7 Jul 2023 16:32:06 +0200
Subject: [PATCH] code update
---
.../dyn_glm_chain_analysis.cpython-37.pyc | Bin 58714 -> 53646 bytes
__pycache__/index_mice.cpython-37.pyc | Bin 0 -> 1756 bytes
.../mcmc_chain_analysis.cpython-37.pyc | Bin 5771 -> 6838 bytes
behavioral_state_data.py | 7 +-
behavioral_state_data_easier.py | 9 +-
canonical_infos.json | 2 +-
combine_canon_infos.py | 33 +
dyn_glm_chain_analysis.py | 638 +++++++-----------
dynamic_GLMiHMM_consistency.py | 1 +
dynamic_GLMiHMM_fit.py | 27 -
index_mice.py | 8 +-
multi_chain_saves/.gitignore | 23 +
multi_chain_saves/work.py | 15 -
...ocessing.py => raw_fit_processing_part1.py | 81 ++-
raw_fit_processing_part2.py | 238 +++++++
remove_duplicates.py | 86 +++
test_codes/bias_shift_test.py | 38 ++
17 files changed, 730 insertions(+), 476 deletions(-)
create mode 100644 __pycache__/index_mice.cpython-37.pyc
create mode 100644 combine_canon_infos.py
create mode 100644 multi_chain_saves/.gitignore
delete mode 100644 multi_chain_saves/work.py
rename raw_fit_processing.py => raw_fit_processing_part1.py (60%)
create mode 100644 raw_fit_processing_part2.py
create mode 100644 remove_duplicates.py
create mode 100644 test_codes/bias_shift_test.py
diff --git a/__pycache__/dyn_glm_chain_analysis.cpython-37.pyc b/__pycache__/dyn_glm_chain_analysis.cpython-37.pyc
index 4a89c93f30137d03c0b66cc48856b5c8f6b5a111..797f6ea4ca71054cc38cb50347dc59cdea7fef0b 100644
GIT binary patch
delta 22569
zcmbt+34C1DdGFjii$<fxvMt-P<wf2tud=+$`@S3RHo|5+*13{L%Zzm9%9g~4F}6ul
z8iULMLWltw0|5dUB9gFWX+qdS2;rrW+)zlGCM|(9O`4=hn!LRK|J>1xBut^N((itA
z_Vqj8`L=tVm%TS$^cJouC<uA*cckcEvH8%s!V+_S@Z7XttS~&nBaFrg)f4es&^W1S
za`oh<Db-U9kL9bMDhkCQQFPX<o`zVt7%YaIHBe({^>k~7HPf1<gfTuf8#T+scrigt
z6qCecF-1%j(?q$Lem<{yj^*pFu;$dyy<k*VB0Nvb5W~d?Q6fg7R;d^*#)z?^tS%4v
z^O3(m6kjkttW2$0y-=5#rpu&TsZ2LhR+ozci*zUB#7r?u%ocOdN#<X7vY7d5p*eZb
z5^+gXh`DEs>MEU9Z<eN(ib^r>jHi0pMvs^;7M$^j1*gsG<zk^&gwz#CT`ZO$bqP{e
zipNBiSc>|qkiJYTNBVN4uNEuBN|atBR*BX4T`L|J>%|6?TPHS(O-NZUHj6Fz-5|D#
zZTQ_Nwu>G3-INknh@A**7Q4i5{B99n5PQU4RNacTUx#O{SR)44dBncjpx7@CAZMF6
zD6T}=?V?9ii>r{bLp&i4iNgq8A-*WC7GZ>TiX$R|QoF=aQG?&zLI?}Ldql0M!|z^k
zOhnJ5Jk|R|y=Xwfe(@#ID4J09fQX4WQVxnQixwe~a;11uB!rDnwRM$fZ7{onvO*je
zCq$cQ7aete;Th?vK4cvdCy?irk6I(U4_jA@Yee$2iH#1&Mj`LG&U-PP*U7vTGsU$c
zRcBTou_9Ku`{)JFX`{MEoU8-1gjI8z_;;N+B~FVo*xKvs0=hFxe4@G*+gK<525lVE
zp>l+xb)M>aJPmjn@igIy;fd=K6HuaAe6qU5jajLdT8>ZZ=-(og5QD326mN}ss*mG2
zfp8n1c03(;uECSU(~0L=JSjXU@mz=J6rR)KhC^ocnd<9xi+_g}KY`T0!ShKxf7@lW
zIabX&PnU7{wdxxzkvi?cxf9k6oTZhk7wuEse`nob9lBtOvwHD9C2p)f+rBuRdlPcc
z>U?iH|HgDq*9UTLI^c=zHqe9NX`i1ib*^9QDWkf#`px|^jA~Xcaw~WFPCR|3CuSb>
zxX{{WeM%^C%NdV#(>iRIbt5*W49I%x*}$RP>aOZ@)i+nG>RYOBt-h`L)77_IpV8BF
zr*~W5ez^a4Sojw=jP#80AiNGue>&T=xZSz|4G%_%HR3a8gVmqy^0hpjZg7v)Wu3zT
zJa3&B<AGcOIn5H{4!xGWYuXz>FUoX@D(j3WTjZ%egF2pu6f14BuGdM|V<XO3*Pk}A
z71bEVXS2hwZanR+z7tD$mvtxb`aI_Vl5{6Fh^6)XIIsh~n|N1xb6PQ=yEO`XM0aCB
z)7{l~>kYr#Y88esQrDZpQ)`a&xKYH)7&~pI=a4S;fl}i?q?E3mY3l=}y1PZ$1tXL1
zsZ9;&Y(Rkz3}!%npALWl7aA`#-BW!J=64UUf1cHf6}ktTcTun7EHUA<L2Ln8!+)Ju
z+<VraLgAeHi&m@Ku((g$FFtqHi`2aCdoS>8<ka6M9uN=m7DD*`KH<-ahs5WR!wtx(
ze?UB(jdeemuKSR9MBI)cJgVEw>Ha*^g|Wv|{?xm?9Eb04jD^miE4J?p@7T9Jym!;?
z^{dzIan!kjVbh&K+qZ7(mppv``nB6OuL*BoeQ@*kzQtcFxN@k^DN5K8+X}~8o5IJe
z#%4>Zy26DQpDi?to!k>q`&jsBG?H+<&57uP9yvq3RunO2smj5BGsYPttw!s3BxZ-3
z<I$L%kTWng7tB;|3@$U~_x$JJ5yr5EOkR}9)1_vri9^3yy#$F);p#-fl6Ev6TQ6l?
zI$pFTS;q-9#f8<Fkc$~9sEN0tt4Jadt&1hfy>gi@QX8qUSxH*``tMm<JlQZ7s>rY{
zv#U_BJjX4clW09!Z`Ih2FA+IzIYp7i#<05pfD&_9tA03aLQ$WYpt9Vi@`it9^;QO#
z66kUJWh7Q-$wDS@xV=LPWfnBy7}36q(vv%4w{OD1-WJMD$DCaut6A63!|PqV&eE<)
z974G%58;PPBFYWN*BEU|&vGUkIjrW6UOr0~G^C!TTV7tx+yQeqjI_8%&3<gS`nSco
z>gmyw%KCO;Huf8;+x-6>y=i=(Lh{Jv+EsO9UR$)4IoY+zwSZ5S1Y+r>k$TRnP{s*F
znwzbdaJ)^?n5<H<v8BTYGkp+NvmqRd%ce+UG+_xx-92{B{vu|D8Y4|sScZ>9Y-g|z
zACKr*xVf=4;S3v4GhMPIlWd7_zY3#O_t@b%>0)a2xCN6Pe=Ka()|QWOk;GTq8mn;%
zHTH#T<1*YBk2gE&v*Siiuj7D>dM5yUYi+dk<~crWPXfDR#S%^q<`qWBQE!hMa}_rN
zsAgNuj=!}zY{wk~s3(ub6Hb0(ye85Zt~nNoqS_EQ)^D8N{bAejqj1>Pn+a?PqdgZ>
zX7=nad&eAl13Sp8jW-I+k~5a|YSY9qsm(~=QC=W9UfDoUPoVcfZ_^Y+o!n#5gdMMw
zktQjby`70XBF&9)jFY+Y3a0rI$Kof-gX!h5+nX&Xuh00MT=#G`McUjIDGlfZ9DuO=
z6k8@1rLnu`uM=N1Dm7Xeauw2Xa5R>)W9&UTj@njJ-zuzpSp9MGT>*_}3)P)dN&-ij
zYN=<Yj5ZFcA557TsAa--WlWuv+Jzt{f+HS|3fYBNd0u+(dTKe?v1pBza0b=HV+r(W
z#cJBaO%YpumNj$jC@_xI^>j31I{|&JF{Om#ujP&p7-65bdhVSX_6Hu}cy3j>vySL9
zs*9&jbo$UcK@n||mT0YU+j38J&uz0lV^nFXX|V}+NG?^8iHN*|L(BmJ9*bZf{Q5`M
zup~ECo>$R1^9)<vIcG$Q`|ac<YNE~Uv)MpGet|{wJ_l-}b*<9EGSx(E98`HD3kTX7
zBS)=9C(!Qx`r1HrG!n>rSam0XUg?VraiQfog53lpjFN{|lOspXt5~p-bXx8KaJ-nU
z6Kt-Dgf#}`3oLy2IH+-4Cge3NwuQxvh+7Cb3HglrT*agsj@mIgoWgi*ZL}tebHiaJ
z`o!BcZPe(VU`w|WkhD5p%+m?P!ZLm$;g}dfPR&u!51?uq2OnqAnLV@TzHE$NiE!^8
zJ;d}rh|B0<HFDm*A{Pr8*9%;%E-Wmbrg2%LN_nn}YXLjbfSc9HclVBzFCrcF^t^Fn
z*0OR3fMX^ictqch)bHoj74T$ZB%@gH5mi5b`kLhot|Qn+a23Jd65K}cFu~gdKPSkJ
zTZc4eEl_XFA7?D+`ThJFqkzeonX1}_WqH|xht;`-OH);B=S~9Rl@mM~jp<V@&oj*?
z;2CrRt+8lJtEEv8Q(wxA_n9#$eiS$RaZo1r=*XWiwa{IST8Wz*mwtMa5>^Wk)k-9K
z%RW@fcG9C{bMlu2dLr+r+C_7XcY1DLwAdJ=QOlU^?)XHHdi*Ps)U3}JOw(2S-;0;u
z(Zx&d+43unEFZ+c1{G*-j>hWJYvUOGwuVFCP1<01JP8wlfCsn6Dg+!~4A+!ArzTaE
zEx(ru{%EWggx)72U@i1wE=IDuqIyDl3GZOh5ZSNvBHg8~tD0#z>d~rkn^z+<jZH3k
z{e><3Usefnv9Nv#`Ddo-9=)O!gd9Wbv;BI(miG)_8Zc()9oFmOl3e*|R=u4dljW}P
zVYPAD`qWnD=t19wkW=6?UtyQF>&<tk?vF?rX?JJF6QPl#-0$90)l48#IJvcr5j#v0
zFE6(foNS!y8bZKh6dOHx%O5sA{yx65VrRj6)Reyf07vznTDEeQI|seUV4rZg1ZrBF
z60J?{a%$}UH`eKo(~|fg?<UY2r>UzZ#+vkze#-j@xH#qc?nDnTVLw1GPmK+gt8Mol
z(Wm0yS>o3OdQpDE(0>qcIfyU9a;cKL9eV4w^^~vrwUK%cCBOn24dq4d&@~3>!Rbka
zxa-+VK+lPrCb@?oe*$-^bn<N(Z4IX{KlcdFVWDRUf@mT=7pD*fBaMkL<N*i(F(*}^
zi8XIX-ka`z=D8=CYyub`oM1{3PA&s>+>SIlxs6t>9mg%2AOS`@2`4Yxj>{`(DmW1R
zzKrqQL{ylLfXC-EdP>&LF;c4#anVlCIXhjAsCqH<K@NfNgxL{J$28V?l0a{<o}|mj
z$V*IhF%t|&bDs;p|IYM|E*Hi1o;<C-yKd9r-=X>rFh&2uX8x6c(wSUNu$o{S0jDCL
zCYaG5=$3wqv;&+daXx5-Mg)VwlHg?d+a7!UNn;qRaM$}`h05J{bh4WxOHsE^!)((r
zYb#%6@)$#PZ9F<V(`_%GWwbj+Ey&ww78;vrK#%9JD%-TkcuXDKG-vTsXv}3Pwlb7O
z{~|{9InQH%0SIXpuuC8UHmgp3ZBzNe0he&D%dmqvcD#vrV|LvA#xYZk-MrbHW~k`q
z*m!-kO!sJMz9KUbcN3PYKWv^>b-BgQ7Stg<v2x7w$9Oc+4jQuDP}N(ekDrbt41pJ%
ze2HLSU3K4<@#YLey||^SKo|C<3#&m}M-*-8UuC9RxOH~UY_z#SKBpSCmQMS=&*MpY
z>LI<4^wb-jxx(lS)p<Jex_oU<+GaAB%7rf9;fH!2-8$bKZ;!UeBny%wgh5S7o=1m^
z7#;x?#7Iwx`s4PsRUcW(Q;hP^SWid}gRUDo@3T3E$M)5aRgF8wj~Qo|iJ}sZm_qFY
zQpY3S8=jyZ+A$_I@q(v*l9-t+Nsg>F#jLYt=cweU`pHmKLu*dWMRHX46fp{F>k25f
zi|YIr&=N{a6>~>=x~Czhhzcg8lv7cS_UBWRfxf1T1yEA!HtZSQGpW}=>MXH{3TvvV
z&-x5cz3uCs&8(ysTG`*Dz#OrxdqmPq<|hLeyh*=Ut_R`k^2q?oR<LX=h|pYJD=5n{
zOI`_8{)!7mqRgJ(y}({5R=Vwoq7fdbDp#pzub7&uwHKW-?8Ws<lBHrb6n1N&<ynVG
zk4}#6t`e&+7~M-Rc*L5Kp3X6ulx69ZvFLSqa!maSJS&r9>sMZICptEn&pt3Mce3eG
zwv-+4B=exiS-(R9OMrUs4#_7Djc3G=Ju0^I2#()(cP8dKxx1rv$LtdkDeC6wzn`uw
z_sJ(v(lMb2fK9|BR{}VObb{NgT9w#!Yy|Na4Gp9V>NmSOJ~Z$2?)Ks@aunm(n!hPh
zAD18&kEx<P!?rvAy}SzE9gPAT-W|?jqu7U8SF8VizOwF)&+zX-jLyjcTi0a8V@|LZ
zyi_9EVSz8Ii8nf4h>IPol`*GS-MOd9@0h#yI%?3~;Y(M!8iM!Pu%;>GZnT6VpJPZL
z&6ymSFJW0i&f>JS>{q3Z?42=X6$|O>U@Jr9sT{xE+T3W#ZYF3dR;wQ0J9M7DvGt|F
zo88H!@`o&6LjEU9?j(4bIc9A_ew(ou)W7WAQ}I00ya`J8!b*TSw46YrRR=YOyvR~o
z2<l=;3qgNW@qLrVu3#mzSva}t!2fOFvN$<)GTz#Jv|Y|tkMA3578vUJedCJAz{<S@
zx{dwncl*ZWINr6}H^?QbbpJ?yF=wxq?4OY8qjPgFXz*m<q9NwGQ5PxJ>TQh0p+|9i
zO_4-{+{8uO%w({QHIcS(dnF`DJ51KpF-0vl-HKp-P*r?Tt}Wym$p|ia<?Ig`N%gb+
zODFC$i;bP;1hdG<Hww+57ym{f{)@dsj3WF7y+hRc0~3rJdyXC0V@%aMo4ru<`UfJi
zjv9m6oVt6vsDC;*Z;g}NEUn|N-Z4b&W4N|4+RVX$<CPDx^g{rSk9mCvwx;v9Vv+~c
zfh!jnxAvU9vcXJgq8)NktN#Gxd^SZkqtD$YnMNEbFO*%(I!ACb0K^lmDFsUe^iWJ0
zQUiY5g3NH#f+VZ6<w?{90zyRsJ%NZWuSMb8*hXgH@_X!KHPf=U4F~TGYHq3xOREmD
zz02CGCk`!Na~>sYa*)rzL3ncUti_Y~F#v9qoTT@(r!(00n$5cg3VPf%riNbLFNeRR
z<{jRa5@u(9+dA7Lyux=z7o=WW7npg->@4Uk#C2md<jC_aH-LKi7kN7Klm293G8c3=
zFPW1Jo(>qEmK#990?A;~KzdFxPXx~TL8)`>Kr(dNK%JbH2`CW+jmvRqoScYJ));E-
zBMauHsnaT?q`8s*Mg&2jd(^J0XXeV0EN{wd)pb`7U!c$aLZrx51p4S}B5{N(s+mCU
zbNL|I&V>6}=(Kv}>IJE58P_#CkO=WK*4Xaq_pDf*eN5iT{F5wSO`uu!0Xw8cL>X@c
zO2$vz3IQ-N0YyEg&j^@hMi~%4-!P33;$w{<LP2A=i34HEi-t#)*|)4AIys@V>=v#E
zL(DxtVUg7LMdtqf%DYhQ61pQ(_j|i>`&<r?B-sgN>*IFN3$c}%dY-tKxo~<Rbi#mU
zVb<|HZmLVIGg1Lm;vX4NA{(hv;uZiKC<U}KlR&GaH|Z;Z)f84B^Q;N133TSP#cVGy
zBImRTR6BySU}rAwXs;jTY#(a+lR+S4fb~xUYl_iIeU9*r^ptuQdAy!w9+dW%uz3H%
zZc){8e00hwke1ylW4)rMo<7Eo?d1%E?!3wSYkY9L9noevi_wwDcVz&L6d`w4x&Q7?
z>pjmRhCZa6HN=?Kqv<341o3acS(V*r_=~0*enRchBtj4AKK8~T$$J5?@}wjhC$Tqr
z12j5n^b0iLYOm|N=bC_BK*eWF_2(0FQ*6VfKWrDbI*^Ey)!TWRnYuONhU8A_F-bPy
z_SJv0fpsqZYN|~<O3cm2*gQti4ElnBATAmCdDN`&G8KP8k9xT6hH^cIY(v?Q9xS{5
zRvInC?PoS^MMbyoRSXSS{EtAy7dR+97?g^%j~Wt{AgWNm*JMyr9q1Ud{&Fa|st*cQ
z_D4bBlusj|Un8)KegQT8WNGSiTq!LdB0v+-6b?FrMe+gW517NpMj-hhYd^sj2P_&f
z=Vi=O?YO3G!^144$MBIDBOhYf0gZ9DziX=eWTlHPZ0~=KDFrU3q^Z(#+?Ncpc!U0+
z56JQeI?5mmb9=s${EM+_KUeYtbm)%WWx&uh_Uc(=Q0Oi60Ltd^bc|KsJ2^}ZN^Mxi
z^Yp=nJ`O>9+f?UL)0;9_gN8B$tsaq&BxvjUVAI0>*z{op4OD7`T6r?y;#7t_eSk8#
zIQ74xO#dxnc|AWd)i+Mw5Y|J`bPLPs(RZ+X6hW<^0%59RDydU>U;tr~ni^`ls2jNq
z!u{M-pTBO=2HoBSHpE?b`_Nr(V6+gu)$P-N|F9Cb{YEA(;}Epa;sck{4s=+irkyHt
zwRD<5YE&VibG#8DTtQ$5tA?1C6KOlv2nO8MzU^XqEx{iN$bCbBr4g)y9o^5=%aKPt
zUOjwGMi)*+L1^KR>eJPWqeM@bxgyO8*n*<q{NFUyd#9!uPxKU@{#;%q*9Apa@#iDe
zsVMYcruyZL8yA!4VEsoRB_s7QOg>eXGBgID+zXo|<R<WtF&ZwQPN=_W(fkTlHso6T
zIOd6X-FIH$U&Sgk$TYpOX>sh`m6~mmpH`uVCx;v?ccdM~mMLG?zLmS&cJ{E&T!Ko1
zc?9zb77#2XSVXXxU<pAL!BTa5*NioKD=#rlJ_JjkX&v`geuZET!8l|&#&k%ACJgy{
z+<~TuY=F6GnE_S(P_=c|U{!vuU{MCmFvvq(1|lDHUxmVgwGB3ezD@?3avVm1U8?Ea
zxUrA1sHtb~Vi<7`>!e3?Kz;Yo>83oW-Z(eI*r--t8kO3K?2g>YWi@fI?{Egl+bQEm
zBS)i+QM+AkK&F$k$7R9gR<>l)xLa#W3d!ha6Mp1ohK+1^gBFV+l*H|0D8nJi9eNlJ
ztc^CrMcFX~CPkPQ)EzhP+p?c|SGba_M*JzPp)a5F43V#e?I2oXyEH>fQ2NwvW%m07
zS`Yo4TB2s9v}`+vb@FISu;W^?bxm`BVKPq6ah@C}7>UK8orS`VQnfrppoQwg3|&nS
zCeW<Q07^yj%PB&G=VLI478A{p8Y|=kI7!(oUsAtPB`Gc7g&-=`Ks<^<7N?NcD?tIT
zvG8ADC;HsR8EPhIA)sbNpK53D22|Y<w>QA`L%7}{&@S6yJ8GlJI+P<TZNnI-ZI+O2
z8HR~3PX}oelu`AGTgqk~W7$E>Z?~g04Z3_nlN2ek%`ht=8`Rft8EX`%cWxP$y2xeL
zvZ<yuU5sinGt(tAn!@a5>NOm-&gag;=4UlW=l~H0f1Z^mU`38JK{w&#H0wJy0X|c6
zWw1JFd9ybTB|0R=p=k3KwoZZ9MQ2BK-a28@NyM>}P+*|>guITOpCULzKvn{lsBhj{
zIff^}$%`~LAB(^(;kec4>f1GC-=c=xHqA8z8uA6!nW;A3Hp%>>+0%5J&lvv~CXp`H
z<e@wM5$Ox?ti_XPEAF}Z({=u(!S0!MuzMC{MsI_+el}EsRcDP0__ycSp(G?s{KGXY
z=?PbK<+OF#`Amb_Gw1NBp80p)X&Q6X>vyLBcX$7==Jn*T`VsaVG*SxA4Z<o!hdgLM
z&14aT%?%Je@P9pGu2xh`fhKdKSV8@wSk*Z%IlO)}JbO1K=h1nU+QdvNqr_<V(uHJ^
zI)7oLHzfb8p1d&Cn52Gu;kxl~a24B=<;jVe(s`wH`r2_8B5zlyWs*8|&%`nCdaW}d
zK3ARfcFsb{iewRFmfVJ%u6&t~cE+f8?r9wCwl&+Wi}5U0g%?Mv){FCvUFxBWOLETG
zQ!ta=>K7O1rl=xJkA6UVucPJ|v8UT7O6`@&`N^3So}jce#NOVLVqb{|`dRqb+N)TG
z8e%x99pso%dv&JV6!dbXsMbo;zUiA*^;xE211j`c!Y+TydC2s=+hKH@;;IWqS7Bn6
z`qz7>!quSYzG3BuI!oG4B^PH(<u!P^3gn-XC888sfs*z~w!b@I=ZM2<$9<FYUN@k#
z#nKj>=up?*H)C=K+ApqO4L9K4figCwPi(P0N4<LANH0pLU))!^Al%^7TVZxCOI8t0
zfPKEMLiuCFO8do*AXe(;sjB;1jIet6{&}mfMt^8Cxx!utR9i8^lU$1@1WYy2@74CI
zWFcFDcXUKA=!)c8%wt8{emfsd%6Ve5D*N1+)XiOkTG}#mE={guUub6qW<N-_B@1<I
zWpd&	P`9c-k7Xrvmcwq0aR<JA;7S>l@(Hi#eEGgA%VK*C*FB1kl27Q0tZCDmS0w
z>DT@%>RF88l>=;>n&sj3>*+Wxgx>05^&8+71jlJ{jK>4>-5{c3SUPt&V(_8P<PNLf
z2v6$TOzuX}gxubiHuXJ@BBt9fuGgul-cy=R9*$%$&aYc&Sp8-!rX~aWtj5z-kt&*L
zfxxh~h!)rlwu*JyjxcOMeE8*B5K@?DJaF$e&y4Z(T~Sh&(?;ip<c98T7{=A=w-1b8
zv{h`aGt+}8Ez#x?-Boe-b_~NSrqp@ec>NWOiPAc+y;BIa?ZL^Z%hY9FdWClj8?n}O
zY)|)IdmnWE$4UQq2)lRnw&%v{_rv8Ki{rM>TK3@X1KMig?HVEv*<Nv?+XpKLC?xah
z52i#L5;8{*DmNcYAOX>?yY!<{BbvC9O~3>)xVxGm`zjn!AC73T+>b_vh@i;jTDLvj
zOMDMM3c04JbFh6V=}8Vwny0;;8&Uf(sODzu-4XPtsatUhrX+3C7~Fkza-&Opg_p#a
zGeYB|LCiWp20;onc3uNw?9R(=>`0%+F4tD+hqT2rjJ8s2tZPWiZ3i+-g>fDfNw<4<
z=2+6s#&FY-PUyLDoR*4fFBqK_{bv59%*-G9keMIh%zNjZnRIGIa#MFCS%GG!51fC$
zwDdeKXiP%a>3aQFr%w#?yWr&v)2nq<z^$As+Fc_~GSoG+<;DY^6-dwsqGx?w-&sQ~
zJ1SY4tV%9TuCzsRHHNbzSq0QFldF@<YjdH3J*751ROw18X(MM4*Lg38t(l4p-|bpe
z$yDF%RS!ROm3fd1y7r|>u@IPDvY#JS`9W`5iv!P?k0IvBFOtdl5&_LUb$Ik2diydH
zpCs7z6ks~*w<|N&F!?ky`^ZG{8J7E)LeZ4evnY)*YC(n~OS*^GdM#Qnt(0G74__nr
zdbYriSIWO<?3=%vcIc*W{A{KC{}ZPQ(W)y?vHdc&;*qM2>4jc}HTYOU<=ffS-isyE
z!lePr9opLYKW*!e*|yNsH!!Go$@wD86^eGL7ay5%j3T2mC@n(3{U}@m+YC5sSG5?2
zC#8&nxnIXAeH`IXK8KRf^L&dHv<&+^L(9~WM@s@)m|Ui=eRSIFcgc1_j0{T3F7zAQ
zsLK&}s6HAf@~z9u?5kX<#$KA>%6<iu{i5&x^{;=$&Tmx5F3p)kflt0f;0k;4Wrjag
z*2CJ4RPS6GSCKh6A8V)1`rv!Lah=|jw1iuYDGk);@%){I(?Hs}Onw_<!ASdvlkyey
z)ML9-lrmhq6nx)BM=o5y%CsjAWadB&zBciu=o##2N8ZG;cx_qQTxjfcOVT+eQJ(hB
z&I)0!6^_FiAdYF_bSpEo2>QM!h2T-fR}oNHb=`@Ma9FNf%CuJ*{Vu_41m7dbh^i7?
zVA^|>BC7m9jDDYBfMiO=v;;{^TXQKSr!K2Wt6?+Rp7<S>{u-hXe4%BJOWU>dX)*dJ
zGY5#f|G+fJ3-D*r@@Yb@WMe;I@*4zi61>E;A0o8W^$v5GIE0B1)jY0D`xerAPXz6N
zVQfu8d&kIsWbr!H(lay_W$-5~3dxx-KuF4{+6hbkGed&lX9NudKPUJX0-gl<uMAye
zsJ~$9dfE~j!WtX2SFgN@1;0s<!?ObEZ?{K2&-e=f5LTg+fi*aK9BQue%>y-PmyvJf
zyXfC3(9_e)MYYC^o^h!zj!SD0iQTl5Ytz~(CZ8s7*z3PB^ecjYC-^nNZwO8@_X&dk
zKoo?=m2fi#vg?>yK$({+jsz4qP-w(@t>5<=M9Xt*=vKC&EqEOa{gP>S5@akqt{(La
z)3n*<w+v}1cnfR%h|%8>{GQ-H2}on)9|*Mkd_O|)dy6({8UK$=AN?o7|01}VxoMlR
zYYtZ8iScmk{KON}053f;BIUre4)3_EgX|yBvlHyKx4TmRpV>}`8%l)&96$6HW{8iz
zww=K)0Tz_-<IPP_y4;DgpCBrF%tDwThZ|#kqYR%h-56plgpo4e7;UI2Uwq!EQosD-
zRHIT2`_kPznLNbEH%FuRC=JU@IHkcxln?0Ul^BbS3dD-iIdTFpDq{!;qCY<9+>_X=
zUPzXMdHtkEMeO2~ZF~}}c`1y=A{RzKz>%;rmWOuKpm65Tj}XsNVmN>>Y)`jgo2)2I
zzV#U8Yf#=Day;>P6pDcO4j@jN#2T*hS1*%mXBdXH_zozmYiJmqJ|QaBJ|V?b()Lbd
zL%o0=j2*5j_#;eK-|eN;9;2hTMLRKjnSM7Sv0`@Z2W(h*=poeG-2(`~8j<?%?W?Ey
zm?Lc09MN6`Mm25L>I~tX#}Z1@kT`l%1DzrFtw&m7!JV7h8tduK$&?7DEwQ=6S3*;4
zGL%hVybg=4Doze_EgE{cxlHAu%UQ`>{B%7ooq+2}=@`u}HCu;z+6@QJf!Pyz5{<bh
zQd0}90v&KRH62(onLFJd-L;|fb*?>lD3}Y&;4Eyy?IZx~?HJQb`&7jVcWbXgnZyUS
zcA!D_L(OYH>a)pZy^Z(TWcUGvGaIcD6BqVHPP&TM<p0oM3wBE@6*bZpVD7JE=b2AO
zo7|wc81t3PtKn(kUp}7fxkKNy&_gU5>=XVvf2ZHZAv}%4l=Q2XXBMS!pdfJvL@w2I
zJcpeDc#q_;43rYj0m*VYAv1R7wtdOYNd|N(<)%<)9;A$15kgByO@>6sm6Ujz8%+5H
z<c5-XOCG3b!bv2L!X0vaU4{(V!TLPg*X>Vwx6mNVZr#q|Rls<N$}~J~3?YR5$sFY5
zY5fJdg$5=UUE~XdwTod$z4zrYgM1>t!HX_7A#YI9*E8bDfMLpjdh9Dprf6!AX1G2~
z;Q6IbUHg<6q>CDo+yQp*koxsgrNi2}j~e^*hLAn`l~0X-dU#4>=Nx43C?D){MRfVn
zwH;rB)t+z<%H2$~bq;}Mjc#OUCxPaLu3*SwN16r5T<E>FtQkB3`eNm!=vEAEWul*P
zUbL?A_asBS@92%KU5I+Am~KMPh11}h4o$dCv~oJ}Q=iIvregGL;!^;*{3AXQw=^BU
z*`p3Tvu#^jC9I!N<N`ycxEPWI#=v2ywq>vF|D=Ka_MY}0@^<DDPw3_dnLW^*qcswN
z?jU~!rP5HxAJnC1r>>4@hE>0RrdiZuR#cW{-Ji<zeMisYL>b%w$|7ZnW~(NO7l6vr
zvr~CrnW-LoG8ozygN-PfIA#fT>B-&0<rFT?RDx*)<pdfXR(u?E7>t>x`R|SinU&(@
ztByp5spVf?(kl@6q6$_@PDdv|o9s&W$C{zpMt9(qjw)8kAdH;J5gujQBivGrB^Mc`
zMWGiDX0cEQ^HZ8c=&m#xg_x7&vF>FV%?4@CCOzlu$TSXs_ekS_W-I(0AMadGF1+}_
zGmHY)Ko_PVAggego`Wtj=+U$3nLXyP4_#%~JN+XLUAZd!>alrxY}qr%DSyO5&QK42
zb+kTx(LX_qlVe%OJ%>XsyR4c$7;p@;tNd}#Vo^V?xB-3Y_0-2Qy9fF}DlUI^H#93c
zjiWqC!1F2fR$hKz@Jgb)FKPFSvFg&-mKqi6y|1;dqjUvf`Fm_)HG9!Sh+N~q8?q4F
zwSim8TfqZ<t=Jf4NT#aMFIcHTtgNXLZH9PN2Mh#c7=k2Q(=<hol!MF?!1VhY?wjJ<
z*}kS2XnmR?D~Td_^dI|W-da3~=K-RgPB<lYdg1u)Nr4T)1x1?%sLB1EzP6LLne|@8
z%>peJa*lz{_&Wp8QTmfU-X*dUArmxbzwqLQ@q>>9By%H7*O`pW5?f}vzCr@@gk(>`
z#0^!M^rB?S)2XYYn$dvEuB8&%b7VC5aE1AVkJ?w^-{Yp4vLu<a$J0`Bz~v1CX}?D8
z>M#*Ab8=`uZp1sQqY%SlY9aF}gk2$3%ZC=Tg0A>+jvS}!LYiMo+h+RRFIp%y*>?Th
zU?kUAL}A16#;g+qHrW}h<*&-H`(|MzR%7kvmXx79T@roRma$;lgZM83lSIfcoD*@N
z7Y*F49Msag^`iRq3!fYT{vc^`C56Ai2VJruPl6fGQDrX``*h&37e^R6s@{6B0HGr<
zjq0B<(*Jc`@Z?M5b@n*G-tZ{(>`Sv#X~(^-=-GADet-+3mye4E7m!%H7PK=f$NUnE
z+?9>=%4iQt{Z1MdOd~EuuK{Z%^H9m+>X00G;K5XE_$B~wLFU7mxx}b6ra&|*0m++Y
z5>7EtspsaGW5%2T1JuOiC<m>#wgw^awwik-)pgXIS4yk&C^P$)y{g?e_b4*<dgO5z
z`wmx~d8N!e=~wr^GRauc^X*sajg(#&j@(r{Wcu#c?O?oK+8EWY$Y}>+C*(4w6qxY^
zg=6?)0R(l|W_+2U$z@#G+djB@>WNp2%v=2G`Bx_yYt(zMzCZjD3T3@}4DEZ-F;)Bb
z8?Oyd=>cYWo%=ZQ?C>?ibCGG9*QsTQ7oBTb&xHdxd}6iC?5Z<@wNV)Txk`I8MfFDO
zJFb#x#Dwe`j&P;oLcZ+H*}9G0?^h4JRu))~N!FO^-(M|O|Mc1q_i1$zm5FdLZ-E*n
z7cDv78c{1XY0QOR&|Ex~aOTHvGk#<Eja5o5n7ql-{$_k(pfTG$?C6EBk1G8#(X9vt
zJdmdgfmua3deiYMujpC(x-yL(_2&2Ir`})_KHF-uGjCGJ05^j7Yg+r-k-s+WgkN?|
zXPH(`YBgou`@5bTd;NQlf30rZ_3t70>y`Q{EFaoO2d9ry6*_T`nBE3G=_}aD_cblj
zTcCHTk-OhSpl`K5vEG3F;#NTpWJSLDpME*OEzfPidca41uKSZ-r!{Pc8d)cgH9+tr
zFp*0Z?7d#p;5SOlZ}`=WH)f{xV`w!L+4wgIPaYmkoNh-OxII9);Td3b`f#QADe9c@
z;0E!!?grWf7-(AoepR~~9HvhwBpTa~n<PkUSB6W?rRy?ioDVliC>{1MOL<x<U{y8R
zrsE|o+y?sFN2B3C|A$PI@26SdAX!{oO2Gy)z<vy&cW^;HM0n4;ruzoGaRfoD0Odm;
zRtjdTLfk{1R6|hyMg-t|Uy$^dz(qQC$PU!!CB1e&suv(NcmbR^uP>2v*36dc4e1*(
zu-J94PPc+PF<4*7L5V;KCH_q7XfrQa2w}{?1$YXpkm8a+CWADEWO)6b!u|N6bAwdo
zZ;>bw0CS(~-i|fcCjEzv)O>)}uZ%(@itRf*q|>N_y`>H}>UnIMGFyDqx<HwuIjEFY
z73$Ls4AIgUe%#xT3GPuN8~dlwGWtl7XIk+V&zsv{$M_QxZ@l7V3Vx=ekFcIC?e?Vm
z;7yK$_841c&*iBi*D3OSrex30QAV{gS|1{EV{#SKwElZ5Lz<_GFr%AM%>?NCA^kdn
zX5a5;R+5F5Fhm7N?_HoBzE0^gsn=itD__h?*Ask#p#NGp1^p)dLoNsTP7inhA74ZO
zN1tyL;SHT3xa_FfD>SFz3Y=gTB6orzuf*ax!Em@~6h2B44#P=FT1WYk&8OL!YYS6f
zes20Is&qPvYa+3DEXuE#=r`SG*C*n!j=+v9_syzU2oJggYr_?l^PJGeZQH|ppd_8U
zsE)L1*OPaxREakxSN@K3Ur+E&GzZa1dp{B&j1;K+-LOt)8MW)pVwLx1NufN57(c@S
z7x9jARa9FHua0wX4l|}Qtm@zVjCwU543p?ll4n~VCi`Cz$RLe=vmYo(`lc1!au61z
zqp}Go6m`GUqTdlt%yDm3wW(=_`qdAn9gd(F@aN-ncQD&slvWP7-)c(^k4K-|O(RDq
zCplU@{;gr^>)#4shTr?)lFk2vzK&)4>PWkE6TQOMbi#S`sJiP%lc0zG=8wv;{n{1t
zS3mHpz*{4Gg<EQ2dsmNSeE_!U)%t6;3&<r}8sX|W`A}=)MsqDGy|J{{9a)EWsbE0k
zyb;_l#Lx<dHN-_c2)>!xX4QmO&1J+?kGwV2hqC+CKfX0F^&Cf{#k~#Xcr}rsz1noV
z0ks($@MG@XfVw-eL+RWd<qvDyk7MHF3^`WpcvRwjML)bIo7**y;ta-X?&KUAqeo{q
zx9bna1l(V_9SS{IHJaOxC7PPFk&LRA$kDirN#2}zO%RHoZ*q(jfLvyx(ER{QycwUB
z!5fj1`d=rPFXY?sq)f;>mgIssL;LhxC!>N24{bxsjA;Z1ML$LGx-<kiKAt-HS>=Cw
z_V_|J#M{#S2FqaP2bv=>fg_yY=v33&lPahzz$Xge3)zf2*ygbK83};`@;d}sCj50>
z>F?hjKT^LT8Dwq^!AtDK)VanxrKwqrRuIf5SU~U&!EXtsu|_$;41(DNHxk@RaGu~I
zfu^ilGJk^6zYx&6pjDlc`YEl2kQEG3l+_X^1Q~3K=I`YSCaxhMlP5KULsmjlHupvS
zb&R&KX_6OyhAAS*6O)8F`LTGoE*^)oaSfJ?9_Ly!1D4(WY*@ne1^qK;6<{uDzEhNH
zjzj0IP2hn1jMz_^o;J$3(pi|FlahaBdcOYJN18>KKVeeZ3l!RAIT{0WhVUb3@kID&
zyJlU({9*--Q5l1hjMU>R3@n-xi>Kh$YS~U<1fM-f8-ieRL>$($J_5jq{&C%MBZh&O
zNYiFE5yk6@Zas3%kwkkF3{KQKI)hoJp9J9y?@gM~D^-N~4U~ivVum{dx3TnWU|@vd
z>nGz=PEl6P2eV_QQ(iXi+n%zRbIXI!+bau0+QkR163_%TWs*j^v(=&kJA?yr9N!oL
zKD(yFKH0rr;lbJ`_J$WYX@baTcy5K3Gu&4Rs3mnj{uNH^8*$Xyu}QA#EhFW*R>PZ7
zS_`C(+VxMxu;}&QI42*{ee`&GG2|#t%Gc5wwcYD~3!7btM0rr>2oND@F-|@LLypo}
zxD6jVi`%l6={K>4X3t+{=-ULe;kwq`%8LwYp-Nw3S^hwW$S^s<Wx&MNMr#E|?nzWD
z>!-7f5~Y4Re`{^9G&nIhGB`XqG&muY@=P9;f&pW0a71WAXk=(wFbA*x3`LpZ;E>SZ
zU|wiMXj}*^MF{K!?17M@^UVO53JBIdRsYY6c4-TtkBSJhFaRZ@S%R<$J=9pJ6ivuM
pqw!l5sxSh95-=tcz@!B5?j-Du`k%l0@eW7X{FKM9Dt|WZ{{bE{09F71
delta 27705
zcmb`v31D1Dl{VV<_M%?2Tb5+YvfI1meRtwb-tBmc6Wg(U?WEmOT}iF2<?5C#Dcxy;
zvxA7k4Ot*@Bm!9=BoGA%frPMR1u_W?TQkEDm>Cja!oZuEpT{upzOSm)TI?|4y(jyg
zuJznnPMtdE)T#TS2i@6M-Njpqih?fu-QTlC*nM{tSDU+w?wI$hD?4UbGdgBkGdpHk
zvpQy5vpeQka}1Yo38Q_kH4nc<?R6dVt@#}btObV4_Ovb(#iB%%_M6rs#1@M(QQmK$
z#sq7Lz0_W2FLw?Y(-*Bk&9Im*W{8<$mY6N(h`C~(s1x(=3t1~|&%i2sW$Wre!&-yz
zCC+)H#(OD3Yn|teY1`K^B<2koVu7d<)nXEG)rcvgR)m`aO{TS;b;LBWP%ILQ#S-M^
z|22nzbA!S;!<?PjD1I)Mie)z#R=ujT02uYZd=24JbvB9RV#N(EYx7Q*SSeQB;1a9O
zn${MvTC73pR-|4cE=B64NZlrWA=Zj@z`7mj>%|77Z$SDEu~F2c^iHu!Y{u`djQD4<
zUF<-?%fwEx3rV}hWnwpe_lP}WFMjumePTa;_le8J0sQV4SBNX|d$~9$4&nEJI4rKh
z?-k;RIEvpZ#W7*w_n`O}aa^1L=0li(ZFshdEuyU1C9ZDri)%y#Ifq5&U&TogMer)|
zOVJ=25jrA-uz}#HcuO>iW~3amEz#0up7Y8lMNG7cHqkCRn!Um`$z>h4kBe61x#fNK
zq=6Im)uL0x&zczHYjW|fTs+b_A7x{z?2Ea)Q<-sv<#1t<5J}N3PKz_5NAxxqSSMpH
zD<Za7Q9H_poHKe(*eAES&KXx<X*Jjl1C63@&^YH#yraU`BJ7c+=iH~>v4u#T1^L;+
zZcxA1*bS3hB8_Lkq>Rh5ZJltPN^n_CB4ah%0o9lH4Yb$|_Hnz}ZnEt`uee_IESAfC
z|3!0KbGc_Pn%kDky+P+@#-VPvW5N)|ji$Jfu*t*_$L!V%(pv4d3)0%`_6yS5?PE6T
z+sLzHtsR^)F@3U2+|=(qVOpJ5CT_>=h~2Q>Wq0Pr;pX)&4%-n9SVRq&o?)-N4zt#U
zQk>&kFcKyPl_?RVbPcC;F3!!eS3Z_6e`+{4Uta2z3{adp>jHDQRon)mli~y7gZNE|
ze(@pvCdKXI9DcjS+u{!KVG!eK;fOntaz^}0+$HWtN>7eFi{y@XB~PDRjaOeto*A}W
z)9VO*`sB%6;ZyHQp7)4*n@y|N?z68MxOUKW*04^Bk2V8h)J}5LL`Hl}+$Zh_=Q-c(
zQ){JJJYc27fR(n>%`W~EgPHq0)^%V!8S!x}nfIws9YWVPyR7%)Ig958JU8OG3D3={
z#0-?UMLcNTs$=%8O<Zudsn{nF`hY02K8Vu&STi5Ob34N4@Z5pt!+0D#cjEa-GwA;=
z1n$Oj51xDRd=$^e6sk`G)qRNFf6mi$m0E0`tM}#>+e27j1&Kp;gLoJ#x<Dq+6`VR`
z3r^+mLOW;lzMiA@M>Mq$nz<EpzIUb0$%seCWj&Cq@~Ls@14mq)R~kTXxP~i>Pme1z
zT=%hY>8x9->z*q(bLg0>(>&@@Fypch+vmk+#N#)(?162VNBaTs1oow{*n+uTF6N#!
z&T^OeY=6KSRAoOmrtFiW%Fa<`r;jQ-G^XrRqsq=!W#?o@75@B~!e1Cwc%~}+#ZhIS
z9#fXhEfZhr_nin>AGaQ~K4E>*ddPa%dc=Cv`jp*2@M#-UX}H*Td!K!tEbXit6LrAO
zoEOiisE6|=zO2^!V`p9BD{a1W0eJvL&WmS<OB?};Q-?e*`!j99b0K-&a9NK7t=TsZ
zbN)E%s#gf{oWcYCK~}56DkQ$z<^|S=@TTW=>CBl!oqj~W#wbYb3Gnr6?9W0-$TQ_!
z2>kl2dDc)Ie=T^|=kPRvd4F!0ZB25iCja4F;MCvS*NLxBa-rX!W4{ptb4RQx7XvQD
zY{W4Cn$v%SBm;RubNDCGlpztep6coAwZu0yIcZicUO?=rL1WNoecqltFnQ2reZd~m
zSnSX1`fE@doqCZabs75$)viG_AB_POApa?lWVOeYvA>9+49cf8bso@E={grYwdIIQ
zyfpg7gWiFs#mnNALDTvYkiV*sZL^2S^ToGrbWd{GPa*~m-};Q&vQqY!#JBqktS{S7
zptoOE6o9b-W8Y#wVP6Mv<h-2;fPbF{{~lSP-^p?I^8;V8r-|1<{aX~L*YUp_6Q2r7
z8Ofb0!eeqll=+&|br_vObFIgZKdu*G=BOemv`J6H>$wT&M-Sg%55Yda+g{Z1to1Aj
z)PC0P&+Tz8?v43am@(4FUoAELf>Ig1*{bS!WET`0cus@|jeNd(_2`$4%fFy8<MK!0
z&!B(bYyWD;^Vaid`SW1Am)QLv!@tLR{=Qn{%f$>*JqC6g|24bC5Bj|*ThRJ7tZdaV
z@k8+=@#B83=+MB|2g!m9TE8K(;wKcl5&q_=@C)J}#6Kd3n|DEL=0)+-k?6onxymn#
ze-b}K6W&xv3kF_cx-bs9>VEw;RfPRnW3^Rs`Mv{@{fGBO4(+;f`{r#2t%}PJY~R0S
z@8*MJQcL#j**hko@~Z7y_g=OovTyUz%L@06s_#7GKejQDElVV$Njs8=_1TdYyS>Yn
z*;1X<+1)Xm5^RokbVMUf-JOluV6?rfMTea^fz{KpLHl&HJt8A5(PXwvg%h!ku6A1|
zJ$NEuR#<^Ev1Ci6Ar?(!-Cc>;qb@n$=?X=Sh0b?FH`iJv5)#wtXlF9g6_0f$6LLPf
zt-%84bYa+7K6H0sl`(N8lUL>QRH+5dH;SIK)*vxkyg8AuWil4;+%9EYX5GLgStsl3
zhzq+tAy+d})EMu^tD=cSthqB$=ay?#k)~*4l9i-A=KVuIDV}W@D;=YB_u@-Xu&zLt
zFGzGZwAziytS1pYZCj<$_V!3_ER#`*K5TY+OJ|gh>Io{#-Oi^=pWS>JgO?Jhc6(*C
zv)Pu^2UT-(O|pXdiAZFP*yuqSdbI3H&y2&v80yTdxp=9(3h_EaHE($EYAespViW32
zc@#e!5>aPlJ?*iMTrU?8$uZ}R$_?{WK|`uu>hjWJ?zlc2b5>XF3XjHWwvTC@#(2JJ
z*Yr_^<cW)-DyV*a&A8zR*dRSyf-Kc&NY%}asGapiySnU7k#%>(I^`u!@uZr{5~PEE
zqMdD#&baJ|w#O2-u$=lyORg$rR<J$Vkt0Q`OodNJ6|F|P+Pf3h#Bnt>(OT8{WLreW
zRG8vynp9bkE9Sg8Y1M4Y+ZnN&n(C@FU3j=ktYSrQk*2tew8!IJma}g1r1_2PgVCDx
zfyBDolWM-Q9?VDr6JmEJtOE2af{^80H@WsW#~##`w7V>CcUL4C&l*X~BU|DLtFS%Z
zh(*`f5{;o+xsHu#r~9f%GV4X*NK(xo2;eyIY2sKg^k0+z$()en0EL?3?E<}IkEL+l
zs;$lJMEd@^BFXkjn?Mk#c~De42hnVxC6-9Wn`N{^HZgk-6S*q7+T&Oi%$56?=1H`~
z&(!&I!;|dovaQgl_FDnH9y_8vdPHh`DgiVfk=?{RivETknEJA@T#=O_NfTHrim>@&
z%spNlOWGZys}RF6=icynpQ5vs&W7pLzH68ob&gM;VjOmQrqA>>Fkz1~IDJ;;as<&4
ztm;Ti$X>+iLb=APt`%U$VvTmfDrt;&Ch$(Xv#~eQ5lzbXqpTH3#z1kbuDYYvdDf?v
zGP;zoyiJ_(aV^Y_;2PRCBjWYl!S<w_N9Uegl<$6S5$D!<)-Pfm(iYu~8kb(kL#yZA
zYphj_GD;L|QZ}F&iA+T0^^6yQ09&G6cGjzYWP~L-sq(B7n7_hkaW>4a%5c83p+sY>
zt9LOGB;;KzqUPDx6l?C5HioG&n#4ktXIR+R(;jWG+q1r2{p;xgpHW00Z)Mf}1Zt#j
zV~7JS`v?vYkQYiWTE&5~&iCf8swZofR{&()=xx^D)fkN^3Mv#B_{eE6;J8f4b{5;k
zVn)<bg^-hw4>&s)%xYw-vqqm)9B*ohHO8=R*sR2;c&}oIirkxtDMdhTnsuX}Szl*F
z#?K_OW~{?1Xlwxg0IB9^@FOg`VCc04uNu=gAUr%rw=+FIb27(ezqIIZsis0^s{yWX
zHZQA~r)XJ`N?kzFn&nGI+w4eZwouO<IbT#-j<tAN?PgYP1IU_*C>}BTMb2G|n~S)z
z(UQq5c!v{QGJnfD2DcFGCOATnAV?FOBY2kJ8w4ZmRw32+70&4;(~K2EcP(i&ikO`5
zspDG~4viE%<}6#bF2hYk-aw$3g}))zsaCZ-%d{2(uEDIYyEAsG+g2oqt}FKP3NuRL
z4cN?2gEQ%+BfrYjVm%s75}P&l`P?KW>{B3AJCPVJdwUIEkZUE;$yW(fN1k(h%a<9?
z4PCPQ5~D<sma$mR_)M4c#-p>GJx<ZZw!(|g=#N%>eD}B6vO+L}B)C9VSFE!+H#S*g
z%+#<6#{k+4mn&f+;B#TuP&=Ea6Wf&RbH2GUyx{{x=Z$qX#j!7lD5MEBn5$V`ji~C7
z8bY<2gOt>AgLI=aebqw4a`vy9c3C|t<fuuLSC$xlLhODH7RE0jf6O%XMz`pe@*f#h
zb2&1h>xRCt%4aOdr{-9$ypClv7pZW}`Q_^E8Aa!+(Kn)uRix#ui018T@^#m}Q7NOn
zx@TMwia6@LdZ*$dR77C~n%bku2wA+m*i5jqagJ*c0hdu>41IddW5#=*$IPV%inj7)
zF9AS^dfEAhOBd-rs6mFrgw4g**xixn?$E=j=HgqdGnS?$=|SE^pe8QDP?SKizC8?m
zkbr|zSExIBhzSP)tep6b36-0ZdLL1%;yWzy8i5*=zhmfi0uBf1McA4v>G@m0tldMe
zto@ylc^P3yeMVbdscyQWpphl|&n&B&sWt)CCr+B=9D@A`?5Q$an3S>ZNN)4fOL#F0
z4G{bSAlDbG7zLy4i3pSl$o#N-7Aa=U9g_Q|-VsB3B~wgDU>mVdML@<@fPrK@8Ev-$
z?RHZ#j$JfifdK<gA{!dPqh$rf1;^9t$Iz79NeZM`m&aoaJ-uOxk>NViWTz%>q`Qhx
z)nKTF90cQu@I_c1^H}FT0yV{|lUk6G4pTWX*zoO+u2CEQ*qIq_T>!<^oSb*!^}DYA
zE^6+FQ1nfrd4b?X0>!3D%bHEIEAl)+e)Wx--8YbNWSHjtMzG57_gDL8%QuE9Hoec7
z$SR!mF;w9^y18Mt&XMnoMa<S^&8FqA0j2zRm$P(B!{YpFhvh7zH)}M3z5SMTf0sv~
zdUnite#;u;PUjC>mRxc_vlKO{;T<7=wZW@(?&JG>(93L4*GWVmW;HsIt#zx%ZNdR9
zutPY`x)bsCk#>)1#{%cMt(Teq?Q#O!I;X3pW$LA+<cfSp^dziv?%6hP?ZpOvq@W6^
zj{OOZ`xytCOxthBKe?PgZks>-f7oVI_wq9cUoc6|_U+TnKf9ck?Q4rvVNb5G^ThV5
z(p_V#EO36XeQ^P72(BIS5ht{xX5P0wE?3Ic3R?ioOh!5&jC8Qsl@7r|{Mn?L3Md1e
z=jx9S?ccG|tV$K7s)RwKOSMadOBt?m9^AEc>AU2)DqK}Aoa4<%O~kn%&IkJo3|G?A
zI>iZJHodksIaQQayTlwCHjp|E@!@dTx#qIk%=AH5>kP3lRh^pDWQs-oW_ogJa_dZ-
z0^)>`Mvv6wfmvcQ&IFg@#IUs4OS2eC%ofXF6r6*cQrfJLGMDFlz`v3<5WH)iScS7b
zg(F!vFrO9?q%II^c<x8ch%L<U*pCN(y3)WR7E8HtJoXN<7mJMpRVg!7nDPy}Q(jT8
z8stIlpOc<}B`n+NM`)?4<(J|7P%j%;E;bDsiEwhoz{=z*v0390rByI1iY+{rO4p=n
zgwwHmYUY~c>h~MTHLaH*VJnPA+hH`?fsvS!nlf;y*g9wotQ~ZTZIfK-+I-5oT*_3u
zetoL8bpxJ_sj00S2X*hKqzd^mbWV>zuDPeeM_j28>_R*DOOOd@@9lCYxF=-X?NshL
ziG_D&Phy$n-wVEZCMreq3ibDk`Z|xi2gS1{bpA3_l;r?fLt23=W6dp2&EA$O(mHm%
z@f^YV=-$3}&0DaqH*+s(`1eHa?TEI<CD_EBDt4wN*4UB_C_RA2CCht=o5I^uko)!@
z0b*ZABVf3I=G&B%#Ixh=O29g^`J?yp@9lcB0sQ@}zRsO3fXLTj$2+tBCdgUDm@RH-
zD314Fw}?bL8(ZSCZ?iL4SK%z)-{G~)D-T)D!2Zg0_1cE;3h-zJmhz9%VdNtWDHin?
zgrKI{wveET+%9p6^Sk{EPN-c&Z5WD=wXmc&*$w-ayqWbCpX0usEz%oWP6E_ByxL@L
zV!8bUkF%WFl#ov__G#yq%MY&l6w}-Z>~00ZPC$6H<zr0pwZjD7DdcS|dO1Nafzp(I
z;Dip$!Z9MTbqOo54MO4$j!(9rS;o7&8hYg-X1WtLbS7wLaS5$O{wvdbu-0|kiMnZ7
zPtt}+lQS|UZ0L8ChUJw6>Rof4I}S{nZMnDZ+o26rMQUT$A{+z)Ph)2^A%BhXmh=4s
zGcu#>cv=F0&<8OcdbEye;+u#OhRq}nLrvDx5lytoNgVUZOom+A810GlF1K*V6`|~#
zHANF>%Z_4zU<<iWaW1CZEu)YkLNMvXlfMI=#C6W~SFH2yFe{85&eyJ}C^X8AQlr@L
zo5e<{^G{dIFb)lwR~|Ixs(ol=tKv+-$RwL-i{J?6#5N~&a0P_w&mOFu^&zI-&T<}_
zGe@y}wNqGl`(fvI2Ui)^(A-09W=667peDVsFOtiN>M;8rWk?B!O$?E0>Wbx6%;Nqo
zj{;aZfeAtt2<RA7${QiE+t4~1Fp*jVka5%nfx^g=$a+O=LlX)gCyso><x_k|1=G|H
zrvMwDE$QlLib%T|)&Zp)bUFKvY}gV($;JZY^RE~WBpKIMJPFD$*pmuU?z676zvt^o
z?pCnzW4AIDbaRKg`Vr?lNA~XSS(|hTQ@C$Xg=k|}nKMq9>B4joxU|t$An!SKQ_73Y
zDv%1L3R3>Fu)UwU2~5wI@}~@>7NkPL<D5D=$?!V!!yA^d?O<D6r2o5Cf<g~Ce?GeK
za_rf=(J3gCx}!Prgfw;y%4mi$vhL?FZ3%}<$%vdpOwJA_^fPO}bKS9386wm*+mL9L
z#5>!2wfWEPY)-bwlgw9w`=zWjZse8nDC6ysruO(52U132DoVLbkKr@JMi``9XqZOO
zoND;-TVaBnO=$yVqnqViu8b=P7{bGZL^ZN7Gw$AS(luZtO<@kWF@W$Q5N>W?7<5nC
zi&L@dJlILwg7V`A7#}uISF1N=2v5o_y#4MIp0qDrkoMo|O8N#0Zgw>p!Y2xD@F)G4
zadk|BjnE_fBG8@*bOh4@APA&9XAM4SFJ2u=c~agcb4&w*`p8S43nvRjC}mWE8nC4o
zv2@BMKnp4JtPd%LqDbWwpLL01(8rI2>)=Sal96)FV+)=$q&cQ&>2Oh|)|Ks-H5jJ8
z+76guW7>?W1o*^~p5DI5wQ2lc+$HNcXLq0=2xbUB@<|#9COE&IJ9T${0+Dx=d7N92
z;zt-(nTj7AL-cc;Pdrff`{p)Eo98Vw59d@V=hQ<JGkQNal!Yy8o@wYS-G>t4oU#>;
zH-)n>xk)VZzVaPN4M#iU$rf9NiKrfJ;9ArhF?1=gL=k}nsU;>wq^?A6WJl@=HW851
zWCPoK8tpD+M0pnEliLX&A{xaNA7<z-mh_13jxNQcvZbJ8m=0kvjwE1?!VV9;v#UE<
z7t{-a8whsjSW+43FkZm%G{Kh$K2Pukf-e#fxf-fEBcsPnaiIwp4pHJ$96~B%F2iG*
zYC%*RQ?Musjag=)S#Co6K@NWl`D@5dbhOdO8u%;6ljmB?feUNH2WG?rPyB24X3_<M
z<nB)g(m@E=Mo%*3J7*;BOuFb~0nP=A#;*ZJR|w(E+h)q&Ak~*oElL%M0#IlJ1ar`A
z+K*iS1-b9b<pz)o@jZvJC<P6fmCxt0gUAk!%GSAG$>oNS8&bI@8-$mkNnntHR4^54
z_NEJvUx@WyFbEb4A-5Fu3Uj3jgt?yGDMEgs&QE$$MkeJ-6`Xaai;-1~tl|OhpevV+
zrJpLsDJhhRbcxOpW}8ROAiHO3sszHke{ggxpkhcLl{Xwi-dQ+6;7z4$o+v{nGOji=
z?K?H&2$Rc)doIf0?PAEEB5i_71G-i{E}>k{OH(DB!Bi=a<kKZ+Qm#a*Fj<)LzyYfQ
z{Vf_OPL(7}#6)z>O(rww7L}R}rOZ~~$cU<bBVDEnmZ4zvg#{Jkn6!U~Ua6k_`0EQE
zh&8r_adZirw}*?a&)XB;80LDU(wQ|o<4jo}uJ7~qGJ`*i1rNx0gRa<Dx;5UJ5nPCg
za8%miYtwz9y+j#~B^W84$>pDj@=P_c7AbNPaQFFk!PMAWuNHn^*}-mDOOxS--mpH#
ztnVw>q0Tew`%1_n!d)XToOdDDGp|iM&%aoakxf`O;9}g^`W9bU?%K3sT<B`<+#_Zu
zhS_CZnU-VPfg<b0Mh0__;xrW;nV9Mzj*#1%V(sk^v?)nz>H8*zP7-9A2$h($0ii68
zxuOE*RB%e1q$w`eN(QF57;V;d<)hdOx8!U{HLy$CO|c#=&vqc$^5VFot2gUqQy>`i
zL}*OlPqZiN#o-NOy^#!}R**x5Q#duL;C5Y*QZa(Em2KNZa2eZ>-vG2hfi~Tj$fek2
za4EHA2Mr%KY(6;%{-k0T`3<5*rW6?vQcVi1-19f0VzFXsWN<7JGzv{JKFJpj@BS_r
z3b4E6#4-*H*0U19+B|7DLRh78Aw>a@zY@@NTAOFkRQXsua&OX=8B+j>ZBTes0jwn1
zo=3n%ZwjR45qOdR!EyOW%hh*d&-H1YRS~Z~nBW#HsSjeg%wc#*27Q9Esx}L?ZFD2S
z@}7yywhJU}u0s!mokRtyx-cSsqk8mXN|tEk$ndWUPY{n1T;U4@Mxz8H11lZUM#>GA
zU}ltR@5yDkM`x)VZ_<<W4)}yec$FAiAbgN@{h|OcK#`U#6FUn$5Zi+tAro6iM`0Q^
z6^J!n;h*HHb)~(jg4QA(OL>#t*5Z`6wS;&?KxBgb7+RQILZT2R1b?!0piHr9v`~Zw
z-87#dt(+E%TvLkwdQ+06trO@T4)>CByg>bH_Qnh`K~(g6AO>Q;Lvts#S9DBFhcN!<
z%$_;PN>Qmjy_>P^!3oT~x>Sql=zUeby^}>v6PlDRMvtno{}j<Dyx$Ej>}{O{z8nx!
z`%Tm<1!J9@Ds3%MJbW^6)TGMv=7KuY;7Df5V%W9HdooLtQ&Od=va=@C25JTq(iN%l
zfm&3Vnkq?6XoKj=p;F#e#;N<TT}-AM6nB$>X{mC^nbSv2+6>L&lVR*rGu1eibCOaX
zF-w&yf;e2N{>`T0OB7_%6X6+;xlK+-L&{R6+6fNsK-gnMf!bqA;W_`;OO*9kqq?y-
z*NqA>ug!y=Ho(>vSmR1hOoaw!pkp&fbq)?vs(0AL6iL*erzvb*>^U5T%~K7~g@TYj
z0#)=po0Tf23%^)E`&ZAbln<jeC0#kH#Ozch-STPqLhaW_*B+SDtj2Mq^uXNA_$1ij
z77w_Sb(p8g6^g!T%Twr#Q}c*Ut$<me=aGhlRYR)T5S2T+5$;n*hii>AfTafJrhFok
zA7*7J)Ya!BBcCpo>DO!20}GIc#V35qjHt^mq$yFC0`)KA^#NtTEFD;^)(zLxYE30#
z&7gZ=2~CQFpe(~Rdw6lF&3^6vzKd*?dZ(;W%E86ANax;{7C8T2sdWTwrc@h;wFFuY
z#%~?e5hjw8Vf&ztFh?7bu%E-QC6rOaIXSQ5N+o&7nHn2-n@cB98o*nAiJz>$F-|cd
zX~XXC>#&n8anU!qi}v-fKCn+Lxi+1Ywyig9=lT^(GTg}CV=JePv{Kfh+M#UWnV5iQ
z7o2rs@r0Z~ROZ`7_>p0TJy|b~ecKc8e`tYh)0+)7%a}-X*^QC5GxBa0$ry4GI@8>w
z{@$s-(E684d8J<d4-s?6pzi}lsk4zC4KzxgPeFh}q)6JyZrPbTW~?*vc>Btz#z2?|
ztKkh}c>Cr+!kIHp0lYSyAO5^(1F;B-U$yht(L!g+OH(87Wk#2dj8&f(!+mWUJKHYw
zRN2;WtZ#Qk>EjfE(Is0%gALl;9+zmLvou_JMFqOqjD7s-LNE{`ZKljS%u}9&qo`-t
z;*sL0yhAWfCLv+1%{RkY60XWzfYo?6tO}t*eG#nd@>pGOFUMQ5o@h7xli;QV^AwD_
zz9{Dy=cVXRGkmp3oS?B1Zt#8XYtyskR-Hmag#`C<Sw3YSPRKAKZ`1aM$@EK{1X6J9
zn=_t>8{xf$Nr!uYwk>Lc?kmdCLNv^kAioDZ>TJJmsS+PE<FsJ;&n%#XkAGokB|+Zw
zV^wMAm)yG+x5VI;2Bo(bULgAQ7M!4RY1xAISZ7-letd5ei4tP>Fr;j^G~43TJ)^y&
z64^2whSFNu9&1?K)yt4t1=;{bOS9!~%DqD$Cr63ElZeY?*4-ED(x<gI5kDp3VK+`e
zsWsl!0iz(6917b4&OUGwYlxqrK{Ml$w-8Z5w5J8lwPF1Qm+47Dje&C{5#az9%6lvc
z`?eCnOBnBE0d2~s>i`ny*Ru-cwI%yDw&O<taA#14u_UvUHJQ@47M&?IXY>)b+0rP_
zG2<CF#-y>8u{NT)0U@g_*LTYPYE>jgSuO8jv7M}-X!Q4t()t2Vn1*pnAX}ba%!-$0
zIMU@LUM;3iHYOW$U@j;$=DRB(dCxItAw1Qngt@q!^0o=H0-OnGcbj6&3wVqf_*a1(
zzq#71#_trAD>vrD(j0{RS&U~2@_`Xj`xHYyfM({Ubid2lXx(BQabB`&l~?Ad)~ZFO
zOh!>uIff?k2-bhRecw{{Pcgk9?XzcWoY^G7OBB_TR5aE86d~*~YWNiyk6X44KzrQn
zeC_yxj2A`uN0SH*y9e<k?gU7}d<zqxnSu$IGI2F;FhTQIVw^Ya>*-A5BrWAr^4Ur7
zB^U=GO}V`&2QvbVhZKc<tPeTMUjZzw1ymGjTx(#NT<=0@Z#9dLEv$<=`<tpVzWs;5
zDQ@KSs-7d@qIuE<^Q`K!@~*gnsduu_LWcg0fc7D^?uPw{6g{0}p%Lrshgj%Z08n~&
zL(XWfC`yl43G>xFG@@d}?X0=lIjD^Asx^be#@5II0F1Xf;Zgg6N6$k)5%A4fUmH&0
znlDgX1S;BYr}mm<P`6N4AtzpRn<icKQ%&N2W-4r?2+8@AY?CAc{jRZ7JbaV`DoGqA
zjTfR7E>No9<M;W!6^8sIYL2DGkn^K!GIgpQBgjTVd4i1R2=aZA8+MP-oZ^(@#{PdH
z!BcD)1s#{s7j~|XRyjq3!9boOot~4m+b>3Wn?_OI^0AZ$Qt~LegXtQCrK!&69B8V^
zC^AxXfu7TK`(hL{t^+wIw7*I`@<GnQD6#?e5vobDk3J+hz0t(Mhgej#;eRC=`3aUC
zH~gILUSd?LEErbbEH!RZ9xT=b#25U3rk&jBrRFY)PQEWsK5oC)?*aKdj2Gp}hq*()
zX!svvBiYA={OEcP3A=)(DSB4j$`jO0cmp^b#^`yw+PN<>)j4F(UC%XkA*%Pn+HP?>
zU$p0sBBWp0wL33HNZUpc((18<^lmH+e5qEGQl2CI|ICd3Iz4>=W99s|>6VCkJ>~ro
zRX{$1Nm08hm=XDqK0Wk}V@8?-o@7~t<RUb7&h7l5c}+&)qad%4Ds@BKY_>6tf?V@q
z_4j=o%ggOdEFs!*fUE~+hrO6^ra%_M$!xei+LfSDbAMea&!dy$_pR%RfqXCPjtU_&
zdl=tKpk%XqkfYDil^n%wTmFFH0011$6lp18;9XfI=~EArVyS9##BI8<)Jek3hT5Yc
zv=+Bva1?*N+nE%bXFN8vCH8P=d5)6F@rLE@>xdU8Cl{<(%H5xIJ2l-qSCau@VrYGm
z(*ZCx)OxfgqYPlwqk!>!C<tYVe#Y?r`Dpi=m6VE6oX6!j5sEi|=iB^Sx(Vo*rUpE(
zEa@pa;mo^pc963W?M1N?XXWX`^)oe|*{tIs@Dlh4VDxZFKLL*JTrx-yA}Ay%A}DrV
zKD}Uzn&C$oryzpnD%Qs-=VdyVGzq35D{IV$(r4nTjcT{ABP!eAe^5gGi3ew8&apE^
zYx2~GMjqvGkR~DT7YmqeJ0a7Lz6PXN28Z;&o%_#Bn@Twa)zrIHvvRtHYSBFB-|m`k
z$~wo_v%si!B6m;DgweFVz-;IPI4f^VS|OU=qYbh4Sh81EAv0TWP#*)!DTo2}40h&T
z{hA&4kvka%7aI=mQ1(kW0N@A?<=G-RjdhJI?gk{`w;}t=aRGqJrt|HdL%ZiP^YT%r
z%Ol&+SyrYrgJZT>V^YgB>HtnHX0_M8!ZK=IJmR$WF3Nn6@x`nYqU$dlV7St>Cp~?9
zDc+w};=MoGiE~2Ur1=Z$EMU)E1Pck2(7l);CB-ack&(@iT$g=>!`;sq^aEp}E81uW
zvp)7xcFE5;;lAp~P0S5K@oR+I5rg8ehzU7&1tq;pwTjiU-N4sf>^!n;2deImCwJg-
zmQeS%<VRU{m2-dJq|Nd^hSvbl(HoBcu|%w)9S$UXA)JHY8;olhN`;XhV^PICmNRrO
z!KF^gwNo>GHcV+WO6AcSjrOpKB5RzRP2EDxT@+)vlKFY{U2;a{%>>juv`4E?pV3$a
zT}tUl3}3r|N09R6$c6*M+Q=#GjzKw8dowHQea89zwKHbLk(l-KIvE8yxtA~5M{qg8
z0KpNbCbhhlgr5!R>v-z8LMu>;zwdHdQ}eW&z9FAuodwPVsafWO?x7b_9%K4%m_&0*
zV+eqMB=RCWTk#~GEFbz+y4kyKVKSJ4N&Bpk^bah;vlzBVcbmI)3EcVC_QTu_GcN9d
zrd*Mw=L&ir!;u@x!eCr*_4%O>p1t411;6<>WdJ{Q(~pn;EH$yUD!C9<Yap_@1kM*)
zLpb_3;ncB1X#`7CrIV-|gv3rL9{6{onbKa2@Ir)He3#frQ>WOJUXiM7odR3!WvLZn
zx7g#ne)GiIsY!q9`aI?t$(3R~jOYtf1-LyObj(}k8q1xHx79j3ZxPeghdxqushRom
zp&F<x;_`klMqCJ9?tJ@}nY9PR70vLNKv&)AMJTD-mkP8MoGX-*6|TBl+XEWcVsW{1
z{?@~xD|ulQevfsjg-+M4l}^QNlL`isbI_-Q&dS@CWybXtUVgb=jl=&AdRQwC4R}ON
zvOcvkwU8<m91sn0c(|mvsv5Ur`6-g*CYCu1@fRM5$Jhtd-kdKt2XC>&aoTP7kIBdl
zUqAK*OwQN}qw>yqX=?{MJ!!^(DNf*%1jUIH&a4m2tv{E~sS&f$gqmvH`4?BGt9#y`
zx+I?$YIB_{mX)b$Q3LaLb?>Ysj#82Z;u`0+56lj|qa9Gr6`i@udFBHPGI!y96|I{^
zWVqofHfL1ql8chk8bv6qwUNC_ZpG2~<g+}6Sesr40x1CnpKOE9ib*xo>r-n<hM;K=
z+V&^JYILkT7mczzI@3A*!BfUb=dBN}*xZN4;Pt7E$!(zJja9DHRy=TQIctI@Bgsvv
zVxqx!XByO8Y)oxMmpArY1&5X>4tI7tC;Dr>kDn_!b*FP(|D>urVaTrmhWgY@nE8=<
z#MRTDT<ARBKfUCw)b`XCoRO5uV&uH#{Iq|P@s{(4{@To2_(DL*sr4AB%GT|okzbNg
zF0K<>cko*>qD92GTF}-V_>RoPTy7;|{HDxs?!?xeq64LM-AcqdkvoiO*XXq>;tGF7
zt4hrbFX3ErC6e7(+`7=j*2^#uihHOT>B?0_syIb^JesvzKx`S<Bep4P|HN_e%8OyZ
zmuk!#T=-gzc|)x$nP*Qt*5Q-hk=ik^7tM$`OK+cEyI1UOHgk;tQPSI4y{ckhA4+uI
zKFP41TW+77xk$O>f&3U;DY_9mog3SO1Ba4_ab$3Y?2qeg;0nBNj;Lu?$8(sAtHjkD
zh-UXNj%3-u5gt9jHMv~&B;BHCzyqgbuu2eUS?kd)q8ADI)dep%4?LEAqEEfji%Pcu
z#WA9Q$8y<#rI3NHQ#GlzsdcIP<Z<V~hbJ1>I!zy5t|isniRv%F%!g&H+YuhGTT-mL
z%;UWJ;p1i<MKpE!>2Em^xWa?og6cf3GQC$ZSAK?@&m*YRj7R-Jx_X?6PY_)3S-@P@
zf7Iu_nC0h~IZ7&)PqN&56pQC%hEN)9R5Bxl*4!IjA3kxw`NHS<hA$9&aiqY{>gCgn
zed9mpow)t$|57ji-wD*kz^Wx`;txBE!fSW_txCmBBcnZRB*VzvNn8j1R$PN4xPaqJ
zXw=&_<^o|gLeDt=`jG`0Dh=6^oDc~+e5A1*w&$N9)hf}lOjn#{_!xxiO7`o$=m_%5
zC<!h9E3BYY8GLR+YmvD*@odt8!jnfgVal|u@7#UYygJDf1k}BqHPqpOYVB_s{iSYx
zUJ=+9Qh7>2DR<W2Jws~*MPU?<{pK%!`3q22J2&0EWC_&(`Bj4Vqyk{PCpq`8p0;#E
zkI}b~|H_p(v~88w*?Z5<ZEB)MISpQHO|F7LEVr!W^Jo`ZI!eosUvqwV&y^XaRh-a@
z{)lWH#haAh&3GwW=^RR*P@08ODz0H+>a_a$T?<1?I1XCdVDuY|ev{w@f)@#t0+OKy
zB41*>k0ELv@?}O}AsDA{e3EHYDfG$6w;8&a{Z?rrR3PNnSwm?;4<ZVk2DVh3XlRe=
zWrpkc1zOCjtno7z#tm{Dr6tr6hWrkbzeVtEg3l27YYb6wk$=Y!9R%eB0;Q?E;e73*
z6Ei0n*E$Hxe3#&R1bJOV{yn19hWObDt$dK<Y7JyLi)!VAp&t_bh=5v%{4qnfF{JDo
zy$sQ$frCNl%<z4WtvC%UMBUEu7LIpSL~a2_=Y@|=o2&fBmBpL72`>^T4`ii7?B{EF
z>LUNqdGn(aGH)_YjYDhnKSvbnUF#mXqXw;lSVh!EsB0#mN!ymcbJ(eWfGCa-sMEqr
ziL8y_ADDKMAa8WncCO2rHb@Zt1)<VUcCpIqjQ%shzYzQ@0hhJ>B|}QH`8YyW0X{K8
z^+>+O^tTCqMQ{{gSo2gq@RsxE`=-NnKYahZ7ORrq1Bv6py!J(F=%x9;*|u<n0M`sh
zZW@0dZ?tffuJ09S1?|^FpI`6V*835p?Iat01OXS^w#spEG}Uze?*1PbyPSv4&owqV
zubqFOM6zIG0m|^N*co`RHUuZlTOhE&VX?sZ%BL^O_<LSY7BB@q9N&brO23>yI@I%p
zB>mu-`YEKsQyyAaK)^wisgPxZ5NWU{rr_BFpVTLj8%z}`!Ndd4jv`(GLyiY-ltD<p
zp`>R3u1KjMQ|U)T$oM5lKw>JSqdFVmf)omfkdh-qA_yrqgjc`;MFf@0MWOTZCu&Q4
zqOi@4SL{MwNy<0$AD{3UUVceMz{6mM6Mnp6&S<XpZj#8=Ug!8jqC^!nB$hb&-cjeZ
zhiWG7KtvD!0t}(LMIG@ve}1Skqu4zZ8{DQW(+1G>EWz4YPn+GF&<01f))jwKJWuhr
zGc35DKrys^46WdccoI6YE{DDScQUn~g?x~B^!1FNA7-M$@qaP&CbN}EKrJ`5vbY#{
z!&P0XP@O=x6bjWVc8IsfJKztXya5siz0OY`UTOq}ZvJ${oao}q3<8tDP2h1}fAms)
zjWf5CyxW+*&pGyZ%dkW@Jfv)-7kI|BLlJj=bl*Q}2SUeOBZgo|2psd+AF+eVsm#)&
zUp{X1WdU2NcAd<)(SMjj^*2pH>3B9cG6mxX{}^$P@H)kl*3Dxx^OGU{Y*~&X$4vyj
z!=Wc4cM{U={PBr(M$lRP+3pN?M#vsd5m^OaqiB!ZVmxQ~FXZ0O5hSm7;b^!5=Qfh5
zd7tzOWQ6c2d&?h@lk$YLiVZEk%9Zvdp=Lt~FsZZPj*Y7o-ciVcS0H^BUD|80*@Mgg
z4=#q|P=bO<uW%3eQo)>)T`*N3jP($B(;+;Cc#7cRUJ1?u|8@0^H-EZVc+w@pNSF3}
zD$jAwr_103SJdW{=TDtWc{7PW!RgRw_2Xqlsj_BRBXBs3dm<{$xJVjY;V&rX9-BrF
z!*sB7ax$Q>1oF8iF*6TrtOTY5LuO{B%F?A!=a6vLOqZNG2i3O-cT7sj@6cx-`V0?s
zHDv1RWJoJO8u%eqE#$$Iow@75=uodj@ySAjiUtbl?~dC>J;3S3%YHV})!|M?(R1{2
zIEonhH!6mn>6q$MIbQC^%ge_!$QUR};hVVW2{21cAipJSGT|Q^fPbvBXY<s6ycNWt
z#QHyMS(xdYzc1Dq>*($XbB_wczeD!oUI#48$~!zfPud;G5jbI~1EW%ICu+@6(!=*%
zlF`OCS&FcX5;PDfq3S;vDkC^aK)H)P>~NA+7vN%@brsoy#*F&H&@zNnVEJ%h#c*Jy
z{7=^O>2H3?|01scL+}TJMgq$6_~sEm5SG9N|9E@1+JA9BB_6@;Z_O_MZ<hQc!8-s}
zNMHGG&+t1(utMX4+qgBBora?V5owOcMeaLDk#5+hVb+h}5Xu(4pWyo20%R_A_b+Ia
zoB9r*2RB4KaM8do95#5N%41mOWUNxggK*`5>opud%1!ucZ^CV>peZYz4KGxMPGDwo
zN}>~aVgBq&#QUZ-;x0~SjNfCTo4WoqR6>31Y;kKM-s$weP?ec{(LyTUDmZv#Yh>y2
zWmXYP+7Wyy4VLW4vQ=ZFD=hDxgG-mL7@K<OaOLGzXzRgU5o`0d$kJtNtl-YQ`yvM+
zBQ4e8fQl?zW=(glc)7wk`LdT!CK8G0+d4Q#fLA;(z$TiXo_gc2uXy@=ec!f=`TuK=
z34W^~vNVE=KFV<o7pC~}DRpm(P^*S0le<LqE~qm3a*<{41y|p$i?*m8nnC^o<a=wC
z#&mq_)h{%BEUNd^sJ90|VOXwu8Q_a@kdCv)8Rxg3wKH!aYve;2_@>EknD8rt4gw`8
zv@xVM-~0}o`_P5{9-r3-I>rM!t}x~qavC<TMm*||`wIn2wTn=Crti9dK^`|SO(#9n
zJI#lybYzb%T;ov}0dRQ<n+90iElBe_MNd_E!O4WrDSK*y;deGXH61YWRP9LWH0D!~
zR5j_&qkR9#lV4Ui{xer+vh&-gyr?}AU-tRMnZT$cU&Wu4-lv$5(!0o*aQNHVs}wSv
ztbYgzQfU-_;(}A64P}2T1d(5RcH&G?4ITylaf2uS18E@+8h3q{r>+9yMkg-b!dMBn
ztO_`3EjQ*EtMO|Z^Wes1B8F;DexcJ?qBP`@X5PfMD3MXMm2489smQm`a_;|P%?Z^^
zC37eh<xSRM^YnLisHhBo%a6-tTzj)MdKu(g-oQk7wid^U77Q+PeNM&Gvy7QTmp<KU
zWYoaTA@;G&8q=aOT+D07y&T*K&fVj+g8J$pH8`sXoV*i0na%om=_0L`GeBMX((qWa
zEYB7<<Lqq&**c%I@|jt_S)@p>v-v9(j(p~Om5*ZlM|=bfWxnZ~tM2B<zg(FaCx0u^
zd!$cF2)~V8yqutcAg>(P1*`x`K&f`4{0o&%tV$`OP^7brHPv)n&<Q#LviBO#KA%(b
zm9R0}S^t&a*DPYKYIFd;mZzZ8L%RcQHi{;>*tzW4+BI3DOI=T~yj@^<?c=9%V)6_W
zWOY@GU#yAaa26ABjn8@D*=l3S(9_R4CO+1%_xY8^)6VtJPsn&U(iHG>X1Bh=tvRg|
z`*~O$LaXEG=HFaFw*0xTHcumSp8C6bC3M$K7-gtY3+fWKe?$mS-C4qyyn;UI(NqJz
zmeF-IrfKEnhpabl>^Y#g*trLV_>I2YEqh&_<BT!DUG!!8XSrO=*K8u564nT@26h5E
zG>;#q2{|FVn4bv3&8`l5a>snmz}FT|;q}Zu-~MBV7cE^qyu0bmZ`n%cl^145zKd79
zJqZAg27Zb{op=L1-#e518ltRWJ0=s5DavQi80b05*MeX;%yjv0Z0VS7_Yx=k;zF};
zp>zF9wayDKmO9g4oPaO(;pVM#?~C^~$eW1nz4`Qi-v=IqcEC2S6xKd`1kmmZFRk;9
zI!L~m&7kQ<{vU!rXh(_AOLOc0mRJ2?<W;%1EOY+&;w66*jvrotV>xkb`CHoYqYH4X
zfW=n%exOU4zXjEgM^M4vL@n--v1-Ox^}&T6lP>aCpnuQlCUwXbW>4shhn3XDPyU5*
zk*_<`p3EH={CmCSKiC#i-)ift+81r3K%E!SRAGFdb2Uk^cJY$o&jj#$6u3X3#XNdz
zEbi)CyLgG{#mC!XjS<LPv0eN2MH*#1kpOexU5hbdO=gsg#0%zqZU{?#6Zh&%j$HA{
zc)Trq?b^jF(|uFMRUTCrQq}j^P3pFmkJWLXQfUo7ZH4^_bKfNRF~QFXKH$Q|GIhT*
z=QQ+aee<;ykMDf5@6LMXmoLwPk*x5QI`fzMiu>+YrVbkoXconyIZ7moy^iIlinKog
z8D#hvNPlcb&6rgPbtP%{Xlv|nEa!H?ZQ|`hjhq}qf9tafA*?U3;1am)sH?9STN{+)
z{MRdU12H<Tpfz)y>95Y*_6XagtPeZt@O>GE4yl8s9YDd@@Q*3OD7vH0JcLIsC38^O
zugZnaZLdygksn6ktS{$@Ps;=&Er||%>J2%tH`>vLxa^eVn{Z=;Rq0{A+fAVMh+=(J
zAdVaL9r!qngz^CwW8O=_oWdOn)C_+sXUexG&7wtp%mn~?*F&J8ESL3AWsvW8j(uzK
z+5qbW31}#g94lEwP)txlP)a~8&hmBPvJgH<vYG959{JX+r8G$4W6s#IyTF8#Y#P7#
zr|i;SN8Tu3qnA{rKYeR@=I@E{7epAK8JJu(p-txh#Ke*YeDp-6;*tO_4a{cwSp*9S
zmJlo_SV8a{!5aksP4FjzKNCPAb;&~n7D0r-2FR$3`|5J5y1%I$=%qS5c$j4#BT&vO
zw5iGfXEI1Y+ls7Y=$pjf!Dv0fCV~V(1#%&g>CZOf8@i1xc4J#kUDLwE27IDkE6Q*w
zf+7hGM*B>uWz;~#C45CGK^Z|gf!c4WyW(TIvCiR-Ve5+-o>Seiq>|Z{Tq<X>0tF^H
zn?UW&-14l#&N!6W476E%<20fDmy`wV=X~ud7GA{=<x{J)D~>PD=`VNduT4MC<lM)+
z^|wORHzwpiFujoK+ORaJKBNgfT-D|`AODQOvC5TgJmEb2ouUj4@WY>*!KZav6LLCR
zSO8bDNUxo=ilg{KUcNp)q)JK|kqWJ<v2i74;qE$Q&*75=5yN+BbUn&p(L`?tzD@OW
zzPBv8opd)^NVO^zDUM2t@N=073*REA(hZ}24w(+)@Z_}R^t?8Go>e;XIUA(rU#h6M
z31^uJ`QLo0^X6-l@%bX6)BVyCrn1krVg&Vd6J<-xi~jlqLn$DX|1AhLT13(<QoXi{
zb1lMI15sLqaBhp8&b3+o59(>(yLN+jbKXEaYvewAt>k1XhhQVLD{#Ys8Cib^AX|tl
z>rh|$C3Aj>376{O4A&lo0vjXJAF<J3`7i--*UEh%H^J03B&8OFtU_(H=}}h<<>5=1
z(ZeF^h~jaEo*?)%!KVmrBT$05+K}fnIsz6Tn)R>8ZTaqYn|#xi_`S32^+oXiioCuu
zBm6b~nf@w&t$%~R*uTSH>7U@A5zM$|PtKt8#lacDdH#anA{bs^iqZFPO8IMg{P+%z
zx?AH>N*q%?xRay*QN6=IcNls4Lm=gd6HB?_bK~a26n}-kJXki%9|~3naYi#5IldQS
z`jk%S!G|J?4d<t?ugSRiEs<ig5)c+3*ovy%rn=bUgQmy}K0c$=x6si4mI9wAUWmH<
z_DBV+M}9PpVWtqb$5^VYO~ptlHY$P713V@y!de$CRQ9RKrXTjHLN>e<U@{^!S21bU
KM5pPEdH)Xq?w?Tr
diff --git a/__pycache__/index_mice.cpython-37.pyc b/__pycache__/index_mice.cpython-37.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..dc0d0f83258d0e8e4559bfb0f57d89d9f097425c
GIT binary patch
literal 1756
zcma)6&2Az`5bo|77>qH-KkEhmXXAsdwGGDWtXHyD`B^EOXzhy=Ss6KL4RkGxfnjzU
zY(Tkj%O!q*D1zLQLmncBJVVslTck*SgD9s|58&P9U@fCAs=leN`s(Y}d^$B1C;0oK
z{G<3|f{?#;F#1{Oe2h!}zz7jU&?>P9|J0^o%xor@FoY>0A}TBqD_OST65GTTY5fDF
z7ZnqI`iy-0mu-R7i-~xj%9B@5f%PWDWS`n`i1#MJLN5V{-jqn7mkhlr^rk~EiC!x7
zrqP=Ty%c)0FfC>(EcgXXd(KN>G>g&i5V7YV33C-I__IVxn(^FE1|W7sU>;_>$Q|?k
z)kubj`3&jZ1FT2M5V;p3Ga5mr3qy1vGu+)}BjgC*pOMyYVD|3!2~PR5orVmUkjA=L
z1v}lZuxl#68bi|>4K_-|Uc|N3`VyAk>>7H=tIM$F_Bpd>aQ=BnVZBTxf;(HVh@3-)
zury+SO<%Rn2rP~&;j~MmT1N!mZk(YS_6n|5Tx+e1Ti)Lx?{ygC-esuByVGdiEJ(u&
z_DpMiyg00md`R~`5Q}|^{9eP_2*tjyk<~h^cD2ImU7{_lyc4a35$!!dz8{JuvD{^S
z-CZLQL{>(Kpr>SO^+w;&lg`7>PM^QX9X#PL-A1GB1NQ?cxQ^tPPQY`!xqZHm>(K$<
z+x0tnw@~F*m<S%<`{=GM2mJ6b#~<OI95&$gC(n2;f9H-vez5--JNoZWNhfc=2md^O
z@;E;zee%uMJGnb~<oIEJFHGJXOCFcyP3H1$R<u2*<`((qFOJ>k$H%-F)TO+1sxsSQ
zTDi^5^Q-T7wjOTXs)hhYw=!>@3k=_~Hxv)(#?Z#pCZbIgjk1Pkq)bg|%TN&sAf$@n
z8FH^#lWL+^c3e;U0e+B3pubPSq%M%tW*v=tjrCOw`SMe@#yj&x$E$m8(W!FRds&ye
z7qadtdezB{z2LfO_;zbUvAR@dn8^ZAW&>oyo*rEx*(~T%b|O@67%8jj%7$<QO}7fD
z8mJg}LgKwDE2~(q`3S3M03~SoDk{MVie(+kX0@RrrJ&yQm4($Vpx7A9tS_!IwFt@@
zY7kDV;?=qyNVu43xL!j=9pA@h%BVRlWmcdq6_v_B9(KZmD}6Yjg7$=t2C^C**H=ct
zE!h?-gyWa&gha`P6r_syfvz2w4JT*}<}gLG<}0&~&3+}p4(c76|MJ;NyO4ECzR1c#
z&~z_6nRUyx8ZSWEx#U;@CCs{>fELH=1iOA)#g1xq(X7H#eIeupni$2$h|&ZzsmW3_
zP8qJZ9u}if%rSALXoAktdG^3a(rrXcx<Qke8{&qJO`0?o=mJ)YV{U4=7x%MS`ak`=
B?#=)J
literal 0
HcmV?d00001
diff --git a/__pycache__/mcmc_chain_analysis.cpython-37.pyc b/__pycache__/mcmc_chain_analysis.cpython-37.pyc
index dfcc898f4df0897c89643b35287844e38e54193b..852f09637cc815a928e81c03093be6c9e08945c3 100644
GIT binary patch
delta 1170
zcmZuw&2G~`5Z>81juXd8o2D)80n#7Z2El<AIG_sgTOlDhz@Z`%vO?BP+$MIryFqPO
zD~EDJNUeHD$(0)y9)L&S=o4370dZi~p@m9Vp4pjiX1<+Suisz$a<g#XadHHowL3lb
z_~LQl0t}9qSL_9rX68GnW!U0z?d8m6c9C5|cY-akIW~`0nb8-dy8QX;uX6u~g#-;#
z)8JD^u~8nNh@U`+GHjyAHeoPBrb%Q<7*J(0I!%TUWlkfJg?LI41&XRv)4*6(8h`ym
zh;lN`QUg=vhHwZ#R8D0Wtdo9|L^etu#W{lgDVb5`Av`2U@XCq`%!rB%qOtu?h+7(e
z^)VV(<I34F`NzE@l@`CGAz|CX%^Wq}glTM*r9Y$+4E}@Bb69_rm^ugI1iIs3R%R2L
z>{*m*Y8}J=i#ymPUsU7Km`wab8Ws1BHnEX`t;U8pCC+T(%tcj5PK@fw5z8RvjLxBt
z!OSSUPLz#>MO~=!XgP!Y5RVzi-aeA1akKA#XFWgVKOWBakB5Veykw*jZf@b`rL%_m
zx+54k%EoUtq4H{C2HVhjOW%>~Gme*l@HDDkdDEAzuq!Lv>x6IWZE?HuHu|xOp(=2F
zgL<=1{lM>ZI~8AGsIrnIdJOGq&wj3}Jf@3!S56cA2=R)C>Oi*qp!uK5YJtz`9x|xn
zbkF0Vh|@yW8$0~2#-D2SBnk3*(DY)Hd2LU6{GP^A9Y5gf+A<p9?m=wUc@v>(InFd%
zb-y2o>mw@=tEa}=sc|FDa*uTzo-4!H6n<y7?ZwuJ-VI^{i9~D#u4W1ZU(}sk)#?ta
zT2_SZuJpq|)Jz$=SQdkijBcwgT`Y}LBJ8)h{)KDy2yKUrdOJ2=dqVQOCQ6ipR(xs8
z5188wL-uDow;PDvwlCc#_dIs+9+Sj<R0PJL068d8y;{@&1EmDGY9A=BH9CtK=-+}I
zEhC(#y1<0HP^LM`TS!x-$+zFA0>1pw2zidnMs@P!+oXWj$-Sv1m_KPvEzn%~H=9E&
A8UO$Q
delta 97
zcmdmH+O5my#LLUY00dg^pC?}vpU5Y{IBlZ32Dd^AYYJNnV-$Oe;>4&(ex($-6nP+5
vIE5*gK~s6-nLi>TnoLEEK#gvij6hW63S?@EO#UWWJULHFiBVzlL@8wezC;*T
diff --git a/behavioral_state_data.py b/behavioral_state_data.py
index be9aef83..11a8cc41 100644
--- a/behavioral_state_data.py
+++ b/behavioral_state_data.py
@@ -73,7 +73,7 @@ to_introduce = [2, 3, 4, 5]
# "ibl_witten_06", "ibl_witten_07", "ibl_witten_12", "ibl_witten_13", "ibl_witten_14", "ibl_witten_15",
# "ibl_witten_16", "KS003", "KS005", "KS019", "NYU-01", "NYU-02", "NYU-04", "NYU-06", "ZM_1367", "ZM_1369",
# "ZM_1371", "ZM_1372", "ZM_1743", "ZM_1745", "ZM_1746"] # zoe's subjects
-subjects = ['ibl_witten_14']
+subjects = ['ZFM-05236']
data_folder = 'session_data'
# why does CSHL058 not work?
@@ -149,7 +149,8 @@ for subject in subjects:
contrast_set = {0, 1, 9, 10}
rel_count = -1
- for i, (eid, extra_eids) in enumerate(zip(fixed_eids, additional_eids)):
+ quit()
+ for i, (eid, extra_eids, date) in enumerate(zip(fixed_eids, additional_eids, fixed_dates)):
try:
trials = one.load_object(eid, 'trials')
@@ -182,6 +183,8 @@ for subject in subjects:
df = pd.concat([df, df2], ignore_index=1)
print('new size: {}'.format(len(df)))
+ pickle.dump(df, open("./sofiya_data/{}_df_{}_{}.p".format(subject, rel_count, date), "wb"))
+
current_contrasts = set(df['signed_contrast'])
diff = current_contrasts.difference(contrast_set)
for c in to_introduce:
diff --git a/behavioral_state_data_easier.py b/behavioral_state_data_easier.py
index e526d7cc..18a9de78 100644
--- a/behavioral_state_data_easier.py
+++ b/behavioral_state_data_easier.py
@@ -34,14 +34,12 @@ misses = []
to_introduce = [2, 3, 4, 5]
amiss = ['UCLA034', 'UCLA036', 'UCLA037', 'PL015', 'PL016', 'PL017', 'PL024', 'NR_0017', 'NR_0019', 'NR_0020', 'NR_0021', 'NR_0027']
-subjects = ['NYU-21']
+subjects = ['ZFM-04019', 'ZFM-05236']
fit_type = ['prebias', 'bias', 'all', 'prebias_plus', 'zoe_style'][0]
if fit_type == 'bias':
loading_info = json.load(open("canonical_infos_bias.json", 'r'))
- r_hats = json.load(open("canonical_info_r_hats_bias.json", 'r'))
elif fit_type == 'prebias':
loading_info = json.load(open("canonical_infos.json", 'r'))
- r_hats = json.load(open("canonical_info_r_hats.json", 'r'))
already_fit = list(loading_info.keys())
remaining_subs = [s for s in subjects if s not in amiss and s not in already_fit]
@@ -70,13 +68,14 @@ for subject in subjects:
print('_____________________')
print(subject)
- if subject in already_fit or subject in amiss:
- continue
+ # if subject in already_fit or subject in amiss:
+ # continue
trials = one.load_aggregate('subjects', subject, '_ibl_subjectTrials.table')
# Load training status and join to trials table
training = one.load_aggregate('subjects', subject, '_ibl_subjectTraining.table')
+ quit()
trials = (trials
.set_index('session')
.join(training.set_index('session'))
diff --git a/canonical_infos.json b/canonical_infos.json
index 12a3fa3a..99ae626e 100644
--- a/canonical_infos.json
+++ b/canonical_infos.json
@@ -1 +1 @@
-{"SWC_022": {"seeds": ["412", "401", "403", "413", "406", "407", "415", "409", "405", "400", "408", "404", "410", "411", "414", "402"], "fit_nums": ["347", "54", "122", "132", "520", "386", "312", "59", "999", "849", "372", "300", "485", "593", "358", "550"], "chain_num": 19}, "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}, "ZFM-05236": {"seeds": ["404", "409", "401", "412", "406", "411", "410", "402", "408", "405", "415", "403", "414", "400", "413", "407"], "fit_nums": ["106", "111", "333", "253", "395", "76", "186", "192", "221", "957", "989", "612", "632", "304", "50", "493"], "chain_num": 14, "ignore": [3, 0, 8, 2, 6, 15, 14, 10]}, "ibl_witten_15": {"seeds": ["408", "412", "400", "411", "410", "407", "403", "406", "413", "405", "404", "402", "401", "415", "409", "414"], "fit_nums": ["40", "241", "435", "863", "941", "530", "382", "750", "532", "731", "146", "500", "967", "334", "375", "670"], "chain_num": 19}, "ibl_witten_13": {"seeds": ["401", "414", "409", "413", "415", "411", "410", "408", "402", "405", "406", "407", "412", "403", "400", "404"], "fit_nums": ["702", "831", "47", "740", "251", "929", "579", "351", "515", "261", "222", "852", "754", "892", "473", "29"], "chain_num": 19}, "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}, "ibl_witten_19": {"seeds": ["412", "415", "413", "408", "409", "404", "403", "401", "405", "411", "410", "406", "402", "414", "407", "400"], "fit_nums": ["234", "41", "503", "972", "935", "808", "912", "32", "331", "755", "117", "833", "822", "704", "901", "207"], "chain_num": 19}, "CSH_ZAD_017": {"seeds": ["408", "404", "413", "406", "414", "411", "400", "401", "415", "407", "402", "412", "403", "409", "405", "410"], "fit_nums": ["928", "568", "623", "841", "92", "251", "829", "922", "964", "257", "150", "970", "375", "113", "423", "564"], "chain_num": 19}, "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_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}, "SWC_021": {"seeds": ["404", "413", "406", "412", "403", "401", "410", "409", "400", "414", "415", "402", "405", "408", "411", "407"], "fit_nums": ["840", "978", "224", "38", "335", "500", "83", "509", "441", "9", "135", "890", "358", "460", "844", "30"], "chain_num": 19}, "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}, "ZFM-05245": {"seeds": ["400", "413", "414", "409", "411", "403", "405", "412", "406", "410", "407", "415", "401", "404", "402", "408"], "fit_nums": ["512", "765", "704", "17", "539", "449", "584", "987", "138", "932", "869", "313", "253", "540", "37", "634"], "chain_num": 11}, "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}, "ZFM-04019": {"seeds": ["413", "404", "408", "403", "415", "406", "414", "410", "402", "405", "411", "400", "401", "412", "409", "407"], "fit_nums": ["493", "302", "590", "232", "121", "938", "270", "999", "95", "175", "576", "795", "728", "244", "32", "177"], "chain_num": 14, "ignore": [11, 12, 6, 3, 5, 2, 15, 1]}, "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}, "KS003": {"seeds": ["405", "401", "414", "415", "410", "404", "409", "413", "412", "408", "411", "407", "402", "406", "403", "400"], "fit_nums": ["858", "464", "710", "285", "665", "857", "990", "438", "233", "177", "43", "509", "780", "254", "523", "695"], "chain_num": 19}, "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}, "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
+{"NYU-45": {"seeds": ["513", "503", "500", "506", "502", "501", "512", "509", "507", "515", "510", "505", "504", "514", "511", "508"], "fit_nums": ["768", "301", "96", "731", "879", "989", "915", "512", "295", "48", "157", "631", "666", "334", "682", "714"], "chain_num": 14}, "UCLA035": {"seeds": ["512", "509", "500", "510", "501", "503", "506", "513", "505", "504", "507", "502", "514", "508", "511", "515"], "fit_nums": ["715", "834", "996", "656", "242", "883", "870", "959", "483", "94", "864", "588", "390", "173", "967", "871"], "chain_num": 14}, "NYU-30": {"seeds": ["505", "508", "515", "507", "504", "513", "512", "503", "509", "500", "510", "514", "501", "502", "511", "506"], "fit_nums": ["885", "637", "318", "98", "209", "171", "472", "823", "956", "89", "762", "260", "76", "319", "139", "785"], "chain_num": 14}, "CSHL047": {"seeds": ["509", "510", "503", "502", "501", "508", "514", "505", "507", "511", "515", "504", "500", "506", "513", "512"], "fit_nums": ["60", "589", "537", "3", "178", "99", "877", "381", "462", "527", "6", "683", "771", "950", "294", "252"], "chain_num": 14}, "NYU-39": {"seeds": ["515", "502", "513", "508", "514", "503", "510", "506", "509", "504", "511", "500", "512", "501", "507", "505"], "fit_nums": ["722", "12", "207", "378", "698", "928", "15", "180", "650", "334", "388", "528", "608", "593", "988", "479"], "chain_num": 14}, "NYU-37": {"seeds": ["508", "509", "506", "512", "507", "503", "515", "504", "500", "505", "513", "501", "510", "511", "502", "514"], "fit_nums": ["94", "97", "793", "876", "483", "878", "886", "222", "66", "59", "601", "994", "526", "694", "304", "615"], "chain_num": 14}, "KS045": {"seeds": ["501", "514", "500", "502", "503", "505", "510", "515", "508", "507", "509", "512", "511", "506", "513", "504"], "fit_nums": ["731", "667", "181", "609", "489", "555", "995", "19", "738", "1", "267", "653", "750", "332", "218", "170"], "chain_num": 14}, "UCLA006": {"seeds": ["509", "505", "513", "515", "511", "500", "503", "507", "508", "514", "504", "510", "502", "501", "506", "512"], "fit_nums": ["807", "849", "196", "850", "293", "874", "216", "542", "400", "632", "781", "219", "331", "730", "740", "32"], "chain_num": 14}, "UCLA033": {"seeds": ["507", "505", "501", "512", "502", "510", "513", "514", "506", "515", "509", "504", "508", "500", "503", "511"], "fit_nums": ["366", "18", "281", "43", "423", "877", "673", "146", "921", "21", "353", "267", "674", "113", "905", "252"], "chain_num": 14}, "NYU-40": {"seeds": ["513", "500", "503", "507", "515", "504", "510", "508", "505", "514", "502", "512", "501", "509", "511", "506"], "fit_nums": ["656", "826", "657", "634", "861", "347", "334", "227", "747", "834", "460", "191", "489", "458", "24", "346"], "chain_num": 14}, "NYU-46": {"seeds": ["507", "509", "512", "508", "503", "500", "515", "501", "511", "506", "502", "505", "510", "513", "514", "504"], "fit_nums": ["503", "523", "5", "819", "190", "917", "707", "609", "145", "416", "376", "603", "655", "271", "223", "149"], "chain_num": 14}, "KS044": {"seeds": ["513", "503", "507", "504", "502", "515", "501", "511", "510", "500", "512", "508", "509", "506", "514", "505"], "fit_nums": ["367", "656", "73", "877", "115", "627", "610", "772", "558", "581", "398", "267", "353", "779", "393", "473"], "chain_num": 14}, "NYU-48": {"seeds": ["502", "507", "503", "515", "513", "505", "512", "501", "510", "504", "508", "511", "506", "514", "500", "509"], "fit_nums": ["854", "52", "218", "963", "249", "901", "248", "322", "566", "768", "256", "101", "303", "485", "577", "141"], "chain_num": 14}, "UCLA012": {"seeds": ["508", "506", "504", "513", "512", "502", "507", "505", "510", "503", "515", "509", "511", "501", "500", "514"], "fit_nums": ["492", "519", "577", "417", "717", "60", "130", "186", "725", "83", "841", "65", "441", "534", "856", "735"], "chain_num": 14}, "KS084": {"seeds": ["515", "507", "512", "503", "506", "508", "510", "502", "509", "505", "500", "511", "504", "501", "513", "514"], "fit_nums": ["816", "374", "140", "955", "399", "417", "733", "149", "300", "642", "644", "248", "324", "830", "889", "286"], "chain_num": 14}, "CSHL052": {"seeds": ["502", "508", "509", "514", "507", "515", "506", "500", "501", "505", "504", "511", "513", "512", "503", "510"], "fit_nums": ["572", "630", "784", "813", "501", "738", "517", "461", "690", "203", "40", "202", "412", "755", "837", "917"], "chain_num": 14}, "NYU-11": {"seeds": ["502", "510", "501", "513", "512", "500", "503", "515", "506", "514", "511", "505", "509", "504", "508", "507"], "fit_nums": ["650", "771", "390", "185", "523", "901", "387", "597", "57", "624", "12", "833", "433", "58", "276", "248"], "chain_num": 14}, "KS051": {"seeds": ["502", "505", "508", "515", "512", "511", "501", "509", "500", "506", "513", "503", "504", "507", "514", "510"], "fit_nums": ["620", "548", "765", "352", "402", "699", "370", "445", "159", "746", "449", "342", "642", "204", "726", "605"], "chain_num": 14}, "NYU-27": {"seeds": ["507", "513", "501", "505", "511", "504", "500", "514", "502", "509", "503", "515", "512", "510", "508", "506"], "fit_nums": ["485", "520", "641", "480", "454", "913", "526", "705", "138", "151", "962", "24", "21", "743", "119", "699"], "chain_num": 14}, "UCLA011": {"seeds": ["513", "503", "508", "501", "510", "505", "506", "511", "507", "514", "515", "500", "509", "502", "512", "504"], "fit_nums": ["957", "295", "743", "795", "643", "629", "142", "174", "164", "21", "835", "338", "368", "341", "209", "68"], "chain_num": 14}, "NYU-47": {"seeds": ["515", "504", "510", "503", "505", "506", "502", "501", "507", "500", "514", "513", "508", "509", "511", "512"], "fit_nums": ["169", "941", "329", "51", "788", "654", "224", "434", "385", "46", "712", "84", "930", "571", "273", "312"], "chain_num": 14}, "CSHL045": {"seeds": ["514", "502", "510", "507", "512", "501", "511", "506", "508", "509", "503", "513", "500", "504", "515", "505"], "fit_nums": ["862", "97", "888", "470", "620", "765", "874", "421", "104", "909", "924", "874", "158", "992", "25", "40"], "chain_num": 14}, "UCLA017": {"seeds": ["510", "502", "500", "515", "503", "507", "514", "501", "511", "505", "513", "504", "508", "512", "506", "509"], "fit_nums": ["875", "684", "841", "510", "209", "207", "806", "700", "989", "899", "812", "971", "526", "887", "160", "249"], "chain_num": 14}, "CSHL055": {"seeds": ["509", "503", "505", "512", "504", "508", "510", "511", "507", "501", "500", "514", "513", "515", "506", "502"], "fit_nums": ["957", "21", "710", "174", "689", "796", "449", "183", "193", "209", "437", "827", "990", "705", "540", "835"], "chain_num": 14}, "UCLA005": {"seeds": ["507", "503", "512", "515", "505", "500", "504", "513", "509", "506", "501", "511", "514", "502", "510", "508"], "fit_nums": ["636", "845", "712", "60", "733", "789", "990", "230", "335", "337", "307", "404", "297", "608", "428", "108"], "chain_num": 14}, "CSHL060": {"seeds": ["510", "515", "513", "503", "514", "501", "504", "502", "511", "500", "508", "509", "505", "507", "506", "512"], "fit_nums": ["626", "953", "497", "886", "585", "293", "580", "867", "113", "734", "88", "55", "949", "443", "210", "555"], "chain_num": 14}, "UCLA015": {"seeds": ["503", "501", "502", "500"], "fit_nums": ["877", "773", "109", "747"], "chain_num": 14}, "KS055": {"seeds": ["510", "502", "511", "501", "509", "506", "500", "515", "503", "512", "504", "505", "513", "508", "514", "507"], "fit_nums": ["526", "102", "189", "216", "673", "477", "293", "981", "960", "883", "899", "95", "31", "244", "385", "631"], "chain_num": 14}, "UCLA014": {"seeds": ["513", "509", "508", "510", "500", "504", "515", "511", "501", "514", "507", "502", "505", "512", "506", "503"], "fit_nums": ["628", "950", "418", "91", "25", "722", "792", "225", "287", "272", "23", "168", "821", "934", "194", "481"], "chain_num": 14}, "CSHL053": {"seeds": ["505", "513", "501", "508", "502", "515", "510", "507", "506", "509", "514", "511", "500", "503", "504", "512"], "fit_nums": ["56", "674", "979", "0", "221", "577", "25", "679", "612", "185", "464", "751", "648", "715", "344", "348"], "chain_num": 14}, "NYU-12": {"seeds": ["508", "506", "509", "515", "511", "510", "501", "504", "513", "503", "507", "512", "500", "514", "505", "502"], "fit_nums": ["853", "819", "692", "668", "730", "213", "846", "596", "644", "829", "976", "895", "974", "824", "179", "769"], "chain_num": 14}, "KS043": {"seeds": ["505", "504", "507", "500", "502", "503", "508", "511", "512", "509", "515", "513", "510", "514", "501", "506"], "fit_nums": ["997", "179", "26", "741", "476", "502", "597", "477", "511", "181", "233", "330", "299", "939", "542", "113"], "chain_num": 14}, "CSHL058": {"seeds": ["513", "505", "514", "506", "504", "500", "511", "503", "508", "509", "501", "510", "502", "515", "507", "512"], "fit_nums": ["752", "532", "949", "442", "400", "315", "106", "419", "903", "198", "553", "158", "674", "249", "723", "941"], "chain_num": 14}, "KS042": {"seeds": ["503", "502", "506", "511", "501", "505", "512", "509", "513", "508", "507", "500", "510", "504", "515", "514"], "fit_nums": ["241", "895", "503", "880", "283", "267", "944", "204", "921", "514", "392", "241", "28", "905", "334", "894"], "chain_num": 14}}
\ No newline at end of file
diff --git a/combine_canon_infos.py b/combine_canon_infos.py
new file mode 100644
index 00000000..a31fe355
--- /dev/null
+++ b/combine_canon_infos.py
@@ -0,0 +1,33 @@
+"""
+ Script for combining local canonical_infos.json and the one from the cluster
+"""
+import json
+
+dist_info = json.load(open("canonical_infos.json", 'r'))
+local_info = json.load(open("canonical_infos_local.json", 'r'))
+
+cluster_subs = ['KS045', 'KS043', 'KS051', 'DY_008', 'KS084', 'KS052', 'KS046', 'KS096', 'KS086', 'UCLA033', 'UCLA005', 'NYU-21', 'KS055', 'KS091']
+
+for key in cluster_subs:
+ print('ignore' in dist_info[key])
+
+quit()
+
+for key in cluster_subs:
+ if key not in local_info:
+ print("Adding all of {} to local info".format(key))
+ local_info[key] = dist_info[key]
+ continue
+ else:
+ for sub_key in dist_info[key]:
+ if sub_key not in local_info[key]:
+ print("Adding {} into local info for {}".format(key))
+ local_info[key][sub_key] = dist_info[key][sub_key]
+ else:
+ if local_info[key][sub_key] == dist_info[key][sub_key]:
+ continue
+ else:
+ assert len(dist_info[key][sub_key]) == 16
+ for x in dist_info[key][sub_key]:
+ assert x in local_info[key][sub_key]
+ local_info[key][sub_key] = dist_info[key][sub_key]
diff --git a/dyn_glm_chain_analysis.py b/dyn_glm_chain_analysis.py
index e20781d8..1b806222 100644
--- a/dyn_glm_chain_analysis.py
+++ b/dyn_glm_chain_analysis.py
@@ -19,7 +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
+from mcmc_chain_analysis import state_num_helper, ll_func, r_hat_array_comp, rank_inv_normal_transform
import pandas as pd
from analysis_pmf import pmf_type, type2color
@@ -670,10 +670,12 @@ def control_flow(test, indices, trials, func_init, first_for, second_for, end_fi
continue
first_for(test, results)
+ counter = -1
for i, m in enumerate([item for sublist in test.results for item in sublist.models]):
if i not in indices:
continue
- second_for(m, j, session_trials, trial_counter, results)
+ counter += 1
+ second_for(m, j, counter, session_trials, trial_counter, results)
end_first_for(results, indices, j, trial_counter=trial_counter, session_trials=session_trials)
trial_counter += len(only_for_length)
@@ -681,181 +683,6 @@ def control_flow(test, indices, trials, func_init, first_for, second_for, end_fi
return results
-def create_mode_indices(test, subject, fit_type):
- dim = 3
-
- try:
- xy, z = pickle.load(open("multi_chain_saves/xyz_{}_{}.p".format(subject, fit_type), 'rb'))
- except Exception:
- print('Doing PCA')
- ev, eig, projection_matrix, dimreduc = test.state_pca(subject, pca_type='dists', dim=dim)
- xy = np.vstack([dimreduc[i] for i in range(dim)])
- from scipy.stats import gaussian_kde
- z = gaussian_kde(xy)(xy)
- pickle.dump((xy, z), open("multi_chain_saves/xyz_{}_{}.p".format(subject, fit_type), 'wb'))
-
- print("Mode indices of " + subject)
-
- threshold_search(xy, z, test, 'first_', subject, fit_type)
-
- print("Find another mode?")
- if input() not in ['yes', 'y']:
- return
-
- threshold_search(xy, z, test, 'second_', subject, fit_type)
-
- print("Find another mode?")
- if input() not in ['yes', 'y']:
- return
-
- threshold_search(xy, z, test, 'third_', subject, fit_type)
- return
-
-def threshold_search(xy, z, test, mode_prefix, subject, fit_type):
- happy = False
- conds = [0, None, None, None, None]
- x_min, x_max, y_min, y_max = None, None, None, None
- while not happy:
- print()
- print("Pick level")
- prob_level = input()
- if prob_level == 'cond':
- print("x > ?")
- resp = input()
- if resp not in ['n', 'no']:
- x_min = float(resp)
-
- print("x < ?")
- resp = input()
- if resp not in ['n', 'no']:
- x_max = float(resp)
-
- print("y > ?")
- resp = input()
- if resp not in ['n', 'no']:
- y_min = float(resp)
-
- print("y < ?")
- resp = input()
- if resp not in ['n', 'no']:
- y_max = float(resp)
-
- print("Prob level")
- prob_level = float(input())
- conds = [prob_level, x_min, x_max, y_min, y_max]
- print("Condtions are {}".format(conds))
- else:
- prob_level = float(prob_level)
- conds[0] = prob_level
-
- print("Level is {}".format(prob_level))
-
- mode = conditions_fulfilled(z, xy, conds)
- print("# of samples: {}".format(mode.sum()))
- mode_indices = np.where(mode)[0]
- if mode.sum() > 0:
- print(xy[0][mode_indices].min(), xy[0][mode_indices].max(), xy[1][mode_indices].min(), xy[1][mode_indices].max())
- print("Happy?")
- happy = 'yes' == input()
-
- print("Subset by factor?")
- if input() == 'yes':
- print("Factor?")
- print(mode_indices.shape)
- factor = int(input())
- mode_indices = mode_indices[::factor]
- print(mode_indices.shape)
- loading_info[subject]['mode prob level'] = prob_level
-
- pickle.dump(mode_indices, open("multi_chain_saves/{}mode_indices_{}_{}.p".format(mode_prefix, subject, fit_type), 'wb'))
- consistencies = test.consistency_rsa(indices=mode_indices)
- pickle.dump(consistencies, open("multi_chain_saves/{}mode_consistencies_{}_{}.p".format(mode_prefix, subject, fit_type), 'wb'))
-
-
-def conditions_fulfilled(z, xy, conds):
- works = z > conds[0]
- if conds[1]:
- works = np.logical_and(works, xy[0] > conds[1])
- if conds[2]:
- works = np.logical_and(works, xy[0] < conds[2])
- if conds[3]:
- works = np.logical_and(works, xy[1] > conds[3])
- if conds[4]:
- works = np.logical_and(works, xy[1] < conds[4])
-
- return works
-
-
-def state_set_and_plot(test, mode_prefix, subject, fit_type):
- mode_indices = pickle.load(open("multi_chain_saves/{}mode_indices_{}_{}.p".format(mode_prefix, subject, fit_type), 'rb'))
- consistencies = pickle.load(open("multi_chain_saves/{}mode_consistencies_{}_{}.p".format(mode_prefix, subject, fit_type), 'rb'))
- session_bounds = list(np.cumsum([len(s) for s in test.results[0].models[-1].stateseqs]))
-
- import scipy.cluster.hierarchy as hc
- consistencies /= consistencies[0, 0]
- linkage = hc.linkage(consistencies[0, 0] - consistencies[np.triu_indices(consistencies.shape[0], k=1)], method='complete')
-
- # R = hc.dendrogram(linkage, truncate_mode='lastp', p=150, no_labels=True)
- plt.savefig("peter figures/{}tree_{}_{}".format(mode_prefix, subject, 'complete'))
- plt.close()
-
- session_bounds = list(np.cumsum([len(s) for s in test.results[0].models[-1].stateseqs]))
-
- fig, ax = plt.subplots(ncols=5, sharey=True, gridspec_kw={'width_ratios': [10, 1, 1, 1, 1]}, figsize=(13, 8))
- from matplotlib.pyplot import cm
- for j, criterion in enumerate([0.95, 0.8, 0.5, 0.2]):
- clustering_colors = np.zeros((consistencies.shape[0], 100, 4))
- a = hc.fcluster(linkage, criterion, criterion='distance')
- b, c = np.unique(a, return_counts=1)
- print(b.shape)
- print(np.sort(c))
-
- if criterion in [0.95]:
- state_sets = []
- for x, y in zip(b, c):
- state_sets.append(np.where(a == x)[0])
- print("dumping state set")
- pickle.dump(state_sets, open("multi_chain_saves/{}state_sets_{}_{}.p".format(mode_prefix, subject, fit_type), 'wb'))
- state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, save_append='_{}{}'.format(mode_prefix, criterion), show=True, separate_pmf=True, type_coloring=True)
-
- # I think this is about finding out where states start and how long they last
- # state_id, session_appears = np.where(states)
- # for s in np.unique(state_id):
- # state_appear_mode.append(session_appears[state_id == s][0] / (test.results[0].n_sessions - 1))
-
- # cmap = cm.rainbow(np.linspace(0, 1, 17))#len([x for x, y in zip(b, c) if y > 50])))
- # rank_to_color_place = dict(zip(range(17), [0, 16, 8, 4, 12, 2, 6, 10, 14, 1, 3, 5, 7, 9, 11, 13, 15])) # handcrafted to maximise color distance, I think
- i = -1
- b = [x for _, x in sorted(zip(c, b))][::-1]
- c = [x for x, _ in sorted(zip(c, b))][::-1]
- for x, y in zip(b, c):
- if y > 50:
- i += 1
- # clustering_colors[a == x] = cmap[rank_to_color_place[i]]
- clustering_colors[a == x] = cm.rainbow(np.mean(np.where(a == x)[0]) / a.shape[0])
-
- ax[j+1].imshow(clustering_colors, aspect='auto', origin='upper')
- for sb in session_bounds:
- ax[j+1].axhline(sb, color='k')
- ax[j+1].set_xticks([])
- ax[j+1].set_yticks([])
- ax[j+1].set_title("{}%".format(int(criterion * 100)), size=20)
-
- ax[0].imshow(consistencies, aspect='auto', origin='upper')
- for sb in session_bounds:
- ax[0].axhline(sb, color='k')
- ax[0].set_xticks([])
- ax[0].set_yticks(session_bounds[::-1])
- ax[0].set_yticklabels(session_bounds[::-1], size=18)
- ax[0].set_ylim(session_bounds[-1], 0)
- ax[0].set_ylabel("Trials", size=28)
- plt.yticks(rotation=45)
-
- plt.tight_layout()
- plt.savefig("peter figures/{}clustered_trials_{}_{}".format(mode_prefix, subject, 'criteria comp').replace('.', '_'))
- plt.close()
-
-
def state_pmfs(test, trials, indices):
def func_init(): return {'pmfs': [], 'session_js': [], 'pmf_weights': []}
@@ -863,7 +690,7 @@ def state_pmfs(test, trials, indices):
results['pmf'] = np.zeros(test.results[0].n_contrasts)
results['pmf_weight'] = np.zeros(4)
- def second_for(m, j, session_trials, trial_counter, results):
+ def second_for(m, j, counter, session_trials, trial_counter, results):
states, counts = np.unique(m.stateseqs[j][session_trials - trial_counter], return_counts=True)
for sub_state, c in zip(states, counts):
results['pmf'] += weights_to_pmf(m.obs_distns[sub_state].weights[j]) * c / session_trials.shape[0]
@@ -884,7 +711,7 @@ def state_weights(test, trials, indices):
def first_for(test, results):
results['weights'] = np.zeros(test.results[0].models[0].obs_distns[0].weights.shape[1])
- def second_for(m, j, session_trials, trial_counter, results):
+ def second_for(m, j, counter, session_trials, trial_counter, results):
states, counts = np.unique(m.stateseqs[j][session_trials - trial_counter], return_counts=True)
for sub_state, c in zip(states, counts):
results['weights'] += m.obs_distns[sub_state].weights[j] * c / session_trials.shape[0]
@@ -907,7 +734,7 @@ def lapse_sides(test, state_sets, indices):
def first_for(test, results):
results['pmf'] = np.zeros(test.results[0].n_contrasts)
- def second_for(m, j, session_trials, trial_counter, results):
+ def second_for(m, j, counter, session_trials, trial_counter, results):
states, counts = np.unique(m.stateseqs[j][session_trials - trial_counter], return_counts=True)
for sub_state, c in zip(states, counts):
results['pmf'] += weights_to_pmf(m.obs_distns[sub_state].weights[j]) * c / session_trials.shape[0]
@@ -1102,6 +929,8 @@ def state_development(test, state_sets, indices, save=True, save_append='', show
if test.results[0].name.startswith('GLM_Sim_'):
print("./glm sim mice/truth_{}.p".format(test.results[0].name))
truth = pickle.load(open("./glm sim mice/truth_{}.p".format(test.results[0].name), "rb"))
+ # truth['state_posterior'] = truth['state_posterior'][:, [0, 1, 3, 4, 5, 6, 7]] # 17, mode 2
+ # truth['weights'] = [w for i, w in enumerate(truth['weights']) if i != 2] # 17, mode 2
states_by_session = np.zeros((len(state_sets), test.results[0].n_sessions))
trial_counter = 0
@@ -1172,7 +1001,7 @@ def state_development(test, state_sets, indices, save=True, save_append='', show
counter += 1
pmfs_to_score.append(np.mean(pmfs))
# test.state_mapping = dict(zip(range(len(state_sets)), np.argsort(np.argsort(pmfs_to_score)))) # double argsort for ranks
- test.state_mapping = dict(zip(range(len(state_sets)), np.flip(np.argsort((states_by_session != 0).argmax(axis=1)))))
+ test.state_mapping = dict(zip(np.flip(np.argsort((states_by_session != 0).argmax(axis=1))), range(len(state_sets))))
for state, trials in enumerate(state_sets):
cmap = matplotlib.cm.get_cmap(cmaps[state]) if state < len(cmaps) else matplotlib.cm.get_cmap('Greys')
@@ -1218,7 +1047,7 @@ def state_development(test, state_sets, indices, save=True, save_append='', show
test.state_mapping[state] + states_by_session[state] - 0.5, color=state_color)
else:
- n_points = 150
+ n_points = 400
points = np.linspace(1, test.results[0].n_sessions, n_points)
interpolation = np.interp(points, np.arange(1, 1 + test.results[0].n_sessions), states_by_session[state])
@@ -1229,7 +1058,7 @@ def state_development(test, state_sets, indices, save=True, save_append='', show
ax1.annotate(len(state_sets) - test.state_mapping[state], (test.results[0].n_sessions + 0.1, test.state_mapping[state] - 0.15), fontsize=22, annotation_clip=False)
if test.results[0].name.startswith('GLM_Sim_'):
- ax1.plot(range(1, 1 + test.results[0].n_sessions), truth['state_map'][test.state_mapping[state]] + truth['state_posterior'][:, state] - 0.5, color='r')
+ ax1.plot(range(1, 1 + test.results[0].n_sessions), state + truth['state_posterior'][:, state] - 0.5, color='r')
alpha_level = 0.3
ax2.axvline(0.5, c='grey', alpha=alpha_level, zorder=4)
@@ -1259,7 +1088,7 @@ def state_development(test, state_sets, indices, save=True, save_append='', show
if test.results[0].name.startswith('GLM_Sim_'):
sim_pmf = weights_to_pmf(truth['weights'][state])
- ax2.plot(np.where(defined_points)[0] / (len(defined_points)-1), sim_pmf[defined_points] - 0.5 + truth['state_map'][test.state_mapping[state]], color='r')
+ ax2.plot(np.where(defined_points)[0] / (len(defined_points)-1), sim_pmf[defined_points] - 0.5 + state, color='r')
if not test.state_mapping[state] in dont_plot:
ax2.annotate("Type {}".format(1 + pmf_type(pmf[defined_points])), (1.05, test.state_mapping[state] - 0.37), rotation=90, size=13, color=type2color[pmf_type(pmf)], annotation_clip=False)
@@ -1354,15 +1183,71 @@ def state_development(test, state_sets, indices, save=True, save_append='', show
plt.tight_layout()
if save:
- # print("saving with {} dpi".format(dpi))
+ print("saving with {} dpi".format(dpi))
plt.savefig("dynamic_GLM_figures/meta_state_development_{}_{}{}.png".format(test.results[0].name, separate_pmf, save_append), dpi=dpi)
if show:
plt.show()
else:
plt.close()
+ if subject.startswith("GLM_Sim_"):
+ plt.figure(figsize=(16, 9))
+ for state, trials in enumerate(state_sets):
+ dur_params = dur_hists(test, trials, indices)
+
+ plt.subplot(4, 4, 1 + 2 * test.state_mapping[state])
+ plt.hist(dur_params[:, 0])
+ plt.axvline(truth['durs'][state][0], color='red')
+
+ plt.subplot(4, 4, 2 + 2 * test.state_mapping[state])
+ plt.hist(dur_params[:, 1])
+ plt.axvline(truth['durs'][state][1], color='red')
+
+ plt.tight_layout()
+ plt.savefig("dur hists")
+ plt.show()
+
+ from scipy.stats import nbinom
+ points = np.arange(900)
+ plt.figure(figsize=(16, 9))
+ for state, trials in enumerate(state_sets):
+ dur_params = dur_hists(test, trials, indices)
+
+ plt.subplot(2, 4, 1 + test.state_mapping[state])
+ plt.plot(nbinom.pmf(points, np.mean(dur_params[:, 0]), np.mean(dur_params[:, 1])))
+ plt.plot(nbinom.pmf(points, truth['durs'][state][0], truth['durs'][state][1]), color='red')
+ plt.xlabel("# of trials")
+ plt.ylabel("P")
+
+ plt.tight_layout()
+ plt.savefig("dur dists")
+ plt.show()
+
return states_by_session, all_pmfs, all_pmf_weights, durs, state_types, contrast_intro_types, smart_divide(introductions_by_stage, np.array(durs)), introductions_by_stage, states_per_type
+def dur_hists(test, trials, indices):
+ def func_init(): return {'dur_params': np.zeros((len(indices), 2))}
+
+ def first_for(test, results):
+ pass
+
+ def second_for(m, j, counter, session_trials, trial_counter, results):
+ states, counts = np.unique(m.stateseqs[j][session_trials - trial_counter], return_counts=True)
+ for sub_state, c in zip(states, counts):
+ try:
+ p = m.dur_distns[sub_state].p_save # this cannot be accessed after deleting data, but not every model's data is deleted
+ except:
+ p = m.dur_distns[sub_state].p
+ results['dur_params'][counter] += np.array([m.dur_distns[sub_state].r * c / session_trials.shape[0], p * c / session_trials.shape[0]])
+
+ def end_first_for(results, indices, j, **kwargs):
+ pass
+
+ results = control_flow(test, indices, trials, func_init, first_for, second_for, end_first_for)
+ results['dur_params'] = results['dur_params'] / len(indices)
+ return results['dur_params']
+
+
def compare_pmfs(test, states2compare, states_by_session, all_pmfs, title=""):
"""
Take a set of states, and plot out their PMFs on all sessions on which they occur.
@@ -1645,14 +1530,19 @@ if __name__ == "__main__":
fit_type = ['prebias', 'bias', 'all', 'prebias_plus', 'zoe_style'][0]
if fit_type == 'bias':
loading_info = json.load(open("canonical_infos_bias.json", 'r'))
- r_hats = json.load(open("canonical_info_r_hats_bias.json", 'r'))
elif fit_type == 'prebias':
loading_info = json.load(open("canonical_infos.json", 'r'))
- r_hats = json.load(open("canonical_info_r_hats.json", 'r'))
no_good_pcas = ['NYU-06', 'SWC_023']
subjects = list(loading_info.keys())
- # subjects = ['SWC_021', 'ibl_witten_15', 'ibl_witten_13', 'KS003', 'ibl_witten_19', 'SWC_022', 'CSH_ZAD_017']
- subjects = ['KS021']
+
+ new_subs = ['KS044', 'NYU-11', 'DY_016', 'SWC_061', 'ZFM-05245', 'CSH_ZAD_029', 'SWC_021', 'CSHL058', 'DY_014', 'DY_009', 'KS094', 'DY_018', 'KS043', 'UCLA014', 'SWC_038', 'SWC_022', 'UCLA012', 'UCLA011', 'CSHL055', 'ZFM-04019', 'NYU-45', 'ZFM-02370', 'ZFM-02373', 'ZFM-02369', 'NYU-40', 'CSHL060', 'NYU-30', 'CSH_ZAD_019', 'UCLA017', 'KS052', 'ibl_witten_25', 'ZFM-02368', 'CSHL045', 'UCLA005', 'SWC_058', 'CSH_ZAD_024', 'SWC_042', 'DY_008', 'ibl_witten_13', 'SWC_043', 'KS046', 'DY_010', 'CSHL053', 'ZM_1898', 'UCLA033', 'NYU-47', 'DY_011', 'CSHL047', 'SWC_054', 'ibl_witten_19', 'ibl_witten_27', 'KS091', 'KS055', 'CSH_ZAD_017', 'UCLA035', 'SWC_060', 'DY_020', 'ZFM-01577', 'ZM_2240', 'ibl_witten_29', 'KS096', 'SWC_066', 'DY_013', 'ZFM-01592', 'GLM_Sim_17', 'NYU-48', 'UCLA006', 'NYU-39', 'KS051', 'NYU-27', 'NYU-46', 'ZFM-01936', 'ZFM-02372', 'ZFM-01935', 'ibl_witten_26', 'ZFM-05236', 'ZM_2241', 'NYU-37', 'KS086', 'KS084', 'ZFM-01576', 'KS042']
+
+ miss_count = 0
+ for s in new_subs:
+ if s not in subjects:
+ print(s)
+ miss_count += 1
+ quit()
print(subjects)
fit_variance = [0.03, 0.002, 0.0005, 'uniform', 0, 0.008][0]
@@ -1660,8 +1550,6 @@ if __name__ == "__main__":
# fig, ax = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
- thinning = 25
- summary_info = {"thinning": thinning, "contains": [], "seeds": [], "fit_nums": []}
pop_state_starts = np.zeros(20)
state_appear_dist = np.zeros(10)
state_appear_mode = []
@@ -1703,227 +1591,173 @@ if __name__ == "__main__":
for subject in subjects:
- if subject.startswith('GLM_Sim_') or subject == 'ibl_witten_18':
+ # NYU-11 is quite weird, super errrativ behaviour, has all contrasts introduced at once, no good session at end
+ if subject.startswith('GLM_Sim_'):
continue
print()
print(subject)
- try:
- test = pickle.load(open("multi_chain_saves/canonical_result_{}_{}.p".format(subject, fit_type), 'rb'))
+ test = pickle.load(open("multi_chain_saves/canonical_result_{}_{}.p".format(subject, fit_type), 'rb'))
- # mode_specifier = ''
+ mode_specifier = 'first'
+ try:
+ mode_indices = pickle.load(open("multi_chain_saves/{}_mode_indices_{}_{}.p".format(mode_specifier, subject, fit_type), 'rb'))
+ state_sets = pickle.load(open("multi_chain_saves/{}_state_sets_{}_{}.p".format(mode_specifier, subject, fit_type), 'rb'))
+ except:
try:
- mode_indices = pickle.load(open("multi_chain_saves/first_mode_indices_{}_{}.p".format(subject, fit_type), 'rb'))
- state_sets = pickle.load(open("multi_chain_saves/first_state_sets_{}_{}.p".format(subject, fit_type), 'rb'))
+ 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'))
except:
- try:
- 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'))
- except:
- continue
+ continue
- # lapse differential
- # lapse_sides(test, [s for s in state_sets if len(s) > 40], mode_indices)
-
- # training overview
- # _ = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, save_append='step 0', show=1, separate_pmf=1, type_coloring=False, dont_plot=list(range(8)), plot_until=-1)
- # _ = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, save_append='step 1', show=1, separate_pmf=1, type_coloring=False, dont_plot=list(range(7)), plot_until=2)
- # _ = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, save_append='step 2', show=1, separate_pmf=1, type_coloring=False, dont_plot=list(range(6)), plot_until=7)
- # _ = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, save_append='step 3', show=1, separate_pmf=1, type_coloring=False, dont_plot=list(range(4)), plot_until=13)
- states, pmfs, pmf_weights, durs, state_types, contrast_intro_type, intros_by_type, undiv_intros, states_per_type = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, show=0, separate_pmf=1, type_coloring=True)
- # state_nums_5.append((states > 0.05).sum(0))
- # state_nums_10.append((states > 0.1).sum(0))
- # continue
- a = compare_params(test, 26, [s for s in state_sets if len(s) > 40], mode_indices, [3, 5])
- compare_pmfs(test, [3, 5], states, pmfs, title="{} convergence pmf".format(subject))
- quit()
- new = type_2_appearance(states, pmfs)
-
- if new == 2:
- print('____________________________')
- print(subject)
- print('____________________________')
- if new == 1:
- new_counter += 1
- if new == 0:
- transform_counter += 1
- print(new_counter, transform_counter)
-
- consistencies = pickle.load(open("multi_chain_saves/consistencies_{}_{}.p".format(subject, fit_type), 'rb'))
- consistencies /= consistencies[0, 0]
- temp = contrasts_plot(test, [s for s in state_sets if len(s) > 40], dpi=300, subject=subject, save=True, show=True, consistencies=consistencies, CMF=False)
- continue
+ # lapse differential
+ # lapse_sides(test, [s for s in state_sets if len(s) > 40], mode_indices)
- state_dict = write_results(test, [s for s in state_sets if len(s) > 40], mode_indices)
- pickle.dump(state_dict, open("state_dict_{}".format(subject), 'wb'))
- quit()
- # abs_state_durs.append(durs)
- # continue
- # all_pmf_weights += pmf_weights
- # all_state_types.append(state_types)
- #
- # state_types_interpolation[0] += np.interp(np.linspace(1, state_types.shape[1], 150), np.arange(1, 1 + state_types.shape[1]), state_types[0])
- # state_types_interpolation[1] += np.interp(np.linspace(1, state_types.shape[1], 150), np.arange(1, 1 + state_types.shape[1]), state_types[1])
- # state_types_interpolation[2] += np.interp(np.linspace(1, state_types.shape[1], 150), np.arange(1, 1 + state_types.shape[1]), state_types[2])
- #
- # all_first_pmfs_typeless[subject] = []
- # for pmf in pmfs:
- # all_first_pmfs_typeless[subject].append((pmf[0], pmf[1][0]))
- # all_pmfs.append(pmf)
- #
- # first_pmfs, changing_pmfs = get_first_pmfs(states, pmfs)
- # for pmf in changing_pmfs:
- # if type(pmf[0]) == int:
- # continue
- # all_changing_pmf_names.append(subject)
- # all_changing_pmfs.append(pmf)
- #
- # all_first_pmfs[subject] = first_pmfs
- # continue
- #
- # consistencies = pickle.load(open("multi_chain_saves/consistencies_{}_{}.p".format(subject, fit_type), 'rb'))
- # consistencies /= consistencies[0, 0]
- # contrasts_plot(test, [s for s in state_sets if len(s) > 40], dpi=300, subject=subject, save=True, show=False, consistencies=consistencies, CMF=False)
-
-
- # b_flips = bias_flips(states, pmfs, durs)
- # all_bias_flips.append(b_flips)
-
- # regression, diffs = pmf_regressions(states, pmfs, durs)
- # regression_diffs += diffs
- # regressions.append(regression)
- # continue
- # compare_pmfs(test, [s for s in state_sets if len(s) > 40], mode_indices, [4, 5], states, pmfs, title="{} convergence pmf".format(subject))
- # compare_weights(test, [s for s in state_sets if len(s) > 40], mode_indices, [4, 5], states, title="{} convergence weights".format(subject))
- # quit()
- # all_intros.append(undiv_intros)
- # all_intros_div.append(intros_by_type)
- # if states_per_type != []:
- # all_states_per_type.append(states_per_type)
- #
- # intros_by_type_sum += intros_by_type
- # for pmf in pmfs:
- # all_pmfs.append(pmf)
- # for p in pmf[1]:
- # all_pmf_diffs.append(p[-1] - p[0])
- # all_pmf_asymms.append(np.abs(p[0] + p[-1] - 1))
- # contrast_intro_types.append(contrast_intro_type)
-
- # num_states.append(np.mean(test.state_num_dist()))
- # num_sessions.append(test.results[0].n_sessions)
-
- # a, b = np.where(states)
- # for i in range(states.shape[0]):
- # state_appear.append(b[a == i][0] / (test.results[0].n_sessions - 1))
- # state_dur.append(b[a == i].shape[0])
-
- # state_development_single_sample(test, [mode_indices[0]], show=True, separate_pmf=True, save=False)
-
- # session overview
- # consistencies = pickle.load(open("multi_chain_saves/consistencies_{}_{}.p".format(subject, fit_type), 'rb'))
- # consistencies /= consistencies[0, 0]
- # c_n_a, rt_data = contrasts_plot(test, [s for s in state_sets if len(s) > 40], dpi=300, subject=subject, save=True, show=False, consistencies=consistencies, CMF=True)
- # compare_performance(cnas, (1, 1), title="{} contrast {} performance test".format(subject, (1, 1)))
- # compare_performance(cnas, (0, 0.848))
-
- # duration of different state types (and also percentage of type activities)
- # continue
- # simplex_durs = np.array(durs).reshape(1, 3)
- # print(simplex_durs / np.sum(simplex_durs))
- # from simplex_plot import projectSimplex
- # print(projectSimplex(simplex_durs / simplex_durs.sum(1)[:, None]))
-
- # compute state type proportions and split the pmfs accordingly
- # ret, trans, pmf_types = state_cluster_interpolation(states, pmfs)
- # state_trans += trans # matrix of how often the different types transition into one another
- # for i, r in enumerate(ret):
- # state_trajs[i] += np.interp(np.linspace(1, test.results[0].n_sessions, n_points), np.arange(1, 1 + test.results[0].n_sessions), r) # for interpolated population state trajectories
- # plot_pmf_types(pmf_types, subject=subject, fit_type=fit_type)
- # continue
- # plt.plot(ret.T, label=[0, 1, 2])
- # plt.legend()
- # plt.show()
-
- # Analyse single sample
- # xy, z = pickle.load(open("multi_chain_saves/xyz_{}.p".format(subject), 'rb'))
- #
- # single_sample = [np.argmax(z)]
- # for position, index in enumerate(np.where(np.logical_and(z > 2.7e-7, xy[0] < -500))[0]):
- # print(position, index, index // test.n, index % test.n)
- # states, pmfs = state_development_single_sample(test, [index], save_append='_{}_{}'.format('single_sample', position), show=True, separate_pmf=True)
- # if position == 9:
- # quit()
- # quit()
-
- # all_state_starts = test.state_appearance_posterior(subject)
- # pop_state_starts += all_state_starts
- #
- # a, b, c = test.state_start_and_dur()
- # state_appear += a
- # state_dur += b
- # state_appear_dist += c
-
-
- # all_state_starts = test.state_appearance_posterior(subject)
- # plt.plot(all_state_starts)
- # plt.savefig("temp")
- # plt.close()
-
- print('Computing sub result')
- create_mode_indices(test, subject, fit_type)
- state_set_and_plot(test, 'first_', subject, fit_type)
- print("second mode?")
- if input() in ['y', 'yes']:
- state_set_and_plot(test, 'second_', subject, fit_type)
-
- except FileNotFoundError as e:
- print(e)
- print('no canoncial result')
- continue
- print(r_hats[subject])
- if r_hats[subject] >= 1.05:
- print("Skipping")
- continue
- else:
- print("Making canonical result")
- results = []
- for j, (seed, fit_num) in enumerate(zip(loading_info[subject]['seeds'], loading_info[subject]['fit_nums'])):
- if j in loading_info[subject]['ignore']:
- continue
- summary_info["contains"].append(j)
- summary_info["seeds"].append(seed)
- summary_info["fit_nums"].append(fit_num)
- info_dict = pickle.load(open("./session_data/{}_info_dict.p".format(subject), "rb"))
-
- samples = []
- mini_counter = 1 # start from one here, first 4000 as burn-in
- 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
- save_id = "{}_fittype_{}_var_{}_{}_{}.p".format(subject, fit_type, fit_variance, seed, fit_num).replace('.', '_')
-
- print("Loaded")
-
- result = MCMC_result(samples[::50], infos=info_dict,
- data=samples[0].datas, sessions=fit_type, fit_variance=fit_variance,
- dur=dur, save_id=save_id)
- results.append(result)
- test = MCMC_result_list(results, summary_info)
- pickle.dump(test, open("multi_chain_saves/canonical_result_{}_{}.p".format(subject, fit_type), 'wb'))
-
- test.r_hat_and_ess(state_num_helper(0.05), False)
- test.r_hat_and_ess(state_num_helper(0.02), False)
- test.r_hat_and_ess(state_size_helper(), False)
- test.r_hat_and_ess(state_size_helper(1), False)
- test.r_hat_and_ess(gamma_func, True)
- test.r_hat_and_ess(alpha_func, True)
+ # training overview
+ # _ = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, save_append='step 0', show=1, separate_pmf=1, type_coloring=False, dont_plot=list(range(8)), plot_until=-1)
+ # _ = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, save_append='step 1', show=1, separate_pmf=1, type_coloring=False, dont_plot=list(range(7)), plot_until=2)
+ # _ = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, save_append='step 2', show=1, separate_pmf=1, type_coloring=False, dont_plot=list(range(6)), plot_until=7)
+ # _ = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, save_append='step 3', show=1, separate_pmf=1, type_coloring=False, dont_plot=list(range(4)), plot_until=13)
+ states, pmfs, pmf_weights, durs, state_types, contrast_intro_type, intros_by_type, undiv_intros, states_per_type = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, show=0, separate_pmf=1, type_coloring=True)
+ quit()
+ # state_nums_5.append((states > 0.05).sum(0))
+ # state_nums_10.append((states > 0.1).sum(0))
+ # continue
+ # a = compare_params(test, 26, [s for s in state_sets if len(s) > 40], mode_indices, [3, 5])
+ # compare_pmfs(test, [3, 2, 4], states, pmfs, title="{} convergence pmf".format(subject))
+ consistencies = pickle.load(open("multi_chain_saves/first_mode_consistencies_{}_{}.p".format(subject, fit_type), 'rb'))
+ consistencies /= consistencies[0, 0]
+ temp = contrasts_plot(test, [s for s in state_sets if len(s) > 40], dpi=300, subject=subject, save=True, show=True, consistencies=consistencies, CMF=False)
+ continue
+ quit()
+ new = type_2_appearance(states, pmfs)
+
+ if new == 2:
+ print('____________________________')
+ print(subject)
+ print('____________________________')
+ if new == 1:
+ new_counter += 1
+ if new == 0:
+ transform_counter += 1
+ print(new_counter, transform_counter)
+
+
+ state_dict = write_results(test, [s for s in state_sets if len(s) > 40], mode_indices)
+ pickle.dump(state_dict, open("state_dict_{}".format(subject), 'wb'))
+ quit()
+ # abs_state_durs.append(durs)
+ # continue
+ # all_pmf_weights += pmf_weights
+ # all_state_types.append(state_types)
+ #
+ # state_types_interpolation[0] += np.interp(np.linspace(1, state_types.shape[1], 150), np.arange(1, 1 + state_types.shape[1]), state_types[0])
+ # state_types_interpolation[1] += np.interp(np.linspace(1, state_types.shape[1], 150), np.arange(1, 1 + state_types.shape[1]), state_types[1])
+ # state_types_interpolation[2] += np.interp(np.linspace(1, state_types.shape[1], 150), np.arange(1, 1 + state_types.shape[1]), state_types[2])
+ #
+ # all_first_pmfs_typeless[subject] = []
+ # for pmf in pmfs:
+ # all_first_pmfs_typeless[subject].append((pmf[0], pmf[1][0]))
+ # all_pmfs.append(pmf)
+ #
+ # first_pmfs, changing_pmfs = get_first_pmfs(states, pmfs)
+ # for pmf in changing_pmfs:
+ # if type(pmf[0]) == int:
+ # continue
+ # all_changing_pmf_names.append(subject)
+ # all_changing_pmfs.append(pmf)
+ #
+ # all_first_pmfs[subject] = first_pmfs
+ # continue
+ #
+ # consistencies = pickle.load(open("multi_chain_saves/consistencies_{}_{}.p".format(subject, fit_type), 'rb'))
+ # consistencies /= consistencies[0, 0]
+ # contrasts_plot(test, [s for s in state_sets if len(s) > 40], dpi=300, subject=subject, save=True, show=False, consistencies=consistencies, CMF=False)
+
+
+ # b_flips = bias_flips(states, pmfs, durs)
+ # all_bias_flips.append(b_flips)
+
+ # regression, diffs = pmf_regressions(states, pmfs, durs)
+ # regression_diffs += diffs
+ # regressions.append(regression)
+ # continue
+ # compare_pmfs(test, [s for s in state_sets if len(s) > 40], mode_indices, [4, 5], states, pmfs, title="{} convergence pmf".format(subject))
+ # compare_weights(test, [s for s in state_sets if len(s) > 40], mode_indices, [4, 5], states, title="{} convergence weights".format(subject))
+ # quit()
+ # all_intros.append(undiv_intros)
+ # all_intros_div.append(intros_by_type)
+ # if states_per_type != []:
+ # all_states_per_type.append(states_per_type)
+ #
+ # intros_by_type_sum += intros_by_type
+ # for pmf in pmfs:
+ # all_pmfs.append(pmf)
+ # for p in pmf[1]:
+ # all_pmf_diffs.append(p[-1] - p[0])
+ # all_pmf_asymms.append(np.abs(p[0] + p[-1] - 1))
+ # contrast_intro_types.append(contrast_intro_type)
+
+ # num_states.append(np.mean(test.state_num_dist()))
+ # num_sessions.append(test.results[0].n_sessions)
+
+ # a, b = np.where(states)
+ # for i in range(states.shape[0]):
+ # state_appear.append(b[a == i][0] / (test.results[0].n_sessions - 1))
+ # state_dur.append(b[a == i].shape[0])
+
+ # state_development_single_sample(test, [mode_indices[0]], show=True, separate_pmf=True, save=False)
+
+ # session overview
+ # consistencies = pickle.load(open("multi_chain_saves/consistencies_{}_{}.p".format(subject, fit_type), 'rb'))
+ # consistencies /= consistencies[0, 0]
+ # c_n_a, rt_data = contrasts_plot(test, [s for s in state_sets if len(s) > 40], dpi=300, subject=subject, save=True, show=False, consistencies=consistencies, CMF=True)
+ # compare_performance(cnas, (1, 1), title="{} contrast {} performance test".format(subject, (1, 1)))
+ # compare_performance(cnas, (0, 0.848))
+
+ # duration of different state types (and also percentage of type activities)
+ # continue
+ # simplex_durs = np.array(durs).reshape(1, 3)
+ # print(simplex_durs / np.sum(simplex_durs))
+ # from simplex_plot import projectSimplex
+ # print(projectSimplex(simplex_durs / simplex_durs.sum(1)[:, None]))
+
+ # compute state type proportions and split the pmfs accordingly
+ # ret, trans, pmf_types = state_cluster_interpolation(states, pmfs)
+ # state_trans += trans # matrix of how often the different types transition into one another
+ # for i, r in enumerate(ret):
+ # state_trajs[i] += np.interp(np.linspace(1, test.results[0].n_sessions, n_points), np.arange(1, 1 + test.results[0].n_sessions), r) # for interpolated population state trajectories
+ # plot_pmf_types(pmf_types, subject=subject, fit_type=fit_type)
+ # continue
+ # plt.plot(ret.T, label=[0, 1, 2])
+ # plt.legend()
+ # plt.show()
+
+ # Analyse single sample
+ # xy, z = pickle.load(open("multi_chain_saves/xyz_{}.p".format(subject), 'rb'))
+ #
+ # single_sample = [np.argmax(z)]
+ # for position, index in enumerate(np.where(np.logical_and(z > 2.7e-7, xy[0] < -500))[0]):
+ # print(position, index, index // test.n, index % test.n)
+ # states, pmfs = state_development_single_sample(test, [index], save_append='_{}_{}'.format('single_sample', position), show=True, separate_pmf=True)
+ # if position == 9:
+ # quit()
+ # quit()
+
+ # all_state_starts = test.state_appearance_posterior(subject)
+ # pop_state_starts += all_state_starts
+ #
+ # a, b, c = test.state_start_and_dur()
+ # state_appear += a
+ # state_dur += b
+ # state_appear_dist += c
+
+
+ # all_state_starts = test.state_appearance_posterior(subject)
+ # plt.plot(all_state_starts)
+ # plt.savefig("temp")
+ # plt.close()
# pickle.dump(all_first_pmfs, open("all_first_pmfs.p", 'wb'))
diff --git a/dynamic_GLMiHMM_consistency.py b/dynamic_GLMiHMM_consistency.py
index 8bf1924c..a2aa0174 100644
--- a/dynamic_GLMiHMM_consistency.py
+++ b/dynamic_GLMiHMM_consistency.py
@@ -87,6 +87,7 @@ for subject in subjects:
mega_data[:, 4] = 1
mega_data[:, 5] = data[~bad_trials, 1] - 1
mega_data[:, 5] = (mega_data[:, 5] + 1) / 2
+ print(mega_data.sum(0))
posteriormodel.add_data(mega_data)
import pyhsmm.util.profiling as prof
diff --git a/dynamic_GLMiHMM_fit.py b/dynamic_GLMiHMM_fit.py
index bf1d046d..a16e7376 100644
--- a/dynamic_GLMiHMM_fit.py
+++ b/dynamic_GLMiHMM_fit.py
@@ -19,15 +19,6 @@ import json
import sys
-def crp_expec(n, theta):
- """
- Return expected number of tables after n customers, given concentration theta.
-
- From Wikipedia
- """
- return theta * (digamma(theta + n) - digamma(theta))
-
-
def eleven2nine(x):
"""Map from 11 possible contrasts to 9, for the non-training phases.
@@ -87,13 +78,7 @@ subjects = ['ibl_witten_15', 'ibl_witten_17', 'ibl_witten_18', 'ibl_witten_19',
# test subjects:
subjects = ['KS014']
-# subjects = ['KS021', 'KS016', 'ibl_witten_16', 'SWC_022', 'KS003', 'CSHL054', 'ZM_3003', 'KS015', 'ibl_witten_13', 'CSHL059', 'CSH_ZAD_022', 'CSHL_007', 'CSHL062', 'NYU-06', 'KS014', 'ibl_witten_14', 'SWC_023']
-
-# subjects = [['GLM_Sim_15', 'GLM_Sim_14', 'GLM_Sim_13', 'GLM_Sim_11', 'GLM_Sim_10', 'GLM_Sim_09', 'GLM_Sim_12'][2]]
-# (0.03, 0.3, 5, 'contR', 'contL', 'prevA', 'bias', 1, 0.1):
cv_nums = [15]
-# conda activate hdp_pg_env
-# python dynamic_GLMiHMM_fit.py
cv_nums = [200 + int(sys.argv[1]) % 16]
subjects = [subjects[int(sys.argv[1]) // 16]]
@@ -267,12 +252,6 @@ for loop_count_i, (s, cv_num) in enumerate(product(subjects, cv_nums)):
print("meh, skipped session")
continue
- # if j == 15:
- # import matplotlib.pyplot as plt
- # for i in [0, 2, 3,4,5,6,7,8,10]:
- # plt.plot(i, data[data[:, 0] == i, 1].mean(), 'ko')
- # plt.show()
-
if params['obs_dur'] == 'glm':
for i in range(data.shape[0]):
data[i, 0] = num_to_contrast[data[i, 0]]
@@ -291,10 +270,7 @@ for loop_count_i, (s, cv_num) in enumerate(product(subjects, cv_nums)):
elif reg == 'cont':
mega_data[:, i] = data[mask, 0]
elif reg == 'prevA':
- # prev_ans = data[:, 1].copy()
new_prev_ans = data[:, 1].copy()
- # prev_ans[1:] = prev_ans[:-1]
- # prev_ans -= 1
new_prev_ans -= 1
new_prev_ans = np.convolve(np.append(0, new_prev_ans), params['exp_filter'])[:-(params['exp_filter'].shape[0])]
mega_data[:, i] = new_prev_ans[mask]
@@ -334,9 +310,6 @@ for loop_count_i, (s, cv_num) in enumerate(product(subjects, cv_nums)):
posteriormodel.add_data(mega_data)
- # for d in posteriormodel.datas:
- # print(d.shape)
-
# if not os.path.isfile('./{}/data_save_{}.p'.format(data_folder, params['subject'])):
pickle.dump(data_save, open('./{}/data_save_{}.p'.format(data_folder, params['subject']), 'wb'))
quit()
diff --git a/index_mice.py b/index_mice.py
index 125425e2..96d35de5 100644
--- a/index_mice.py
+++ b/index_mice.py
@@ -15,6 +15,8 @@ for filename in os.listdir("./dynamic_GLMiHMM_crossvals/"):
regexp = re.compile(r'((\w|-)+)_fittype_(\w+)_var_0.03_(\d+)_(\d+)_(\d+)')
result = regexp.search(filename)
subject = result.group(1)
+ if subject == 'ibl_witten_26':
+ print('here')
fit_type = result.group(3)
seed = result.group(4)
fit_num = result.group(5)
@@ -49,8 +51,8 @@ big = []
non_big = []
sim_subjects = []
for s in prebias_subinfo.keys():
- assert len(prebias_subinfo[s]["seeds"]) in [16, 32]
- assert len(prebias_subinfo[s]["fit_nums"]) in [16, 32]
+ # assert len(prebias_subinfo[s]["seeds"]) in [16, 32], s + " " + str(len(prebias_subinfo[s]["seeds"]))
+ # assert len(prebias_subinfo[s]["fit_nums"]) in [16, 32]
print(s, len(prebias_subinfo[s]["fit_nums"]))
if len(prebias_subinfo[s]["fit_nums"]) == 32:
big.append(s)
@@ -58,8 +60,6 @@ for s in prebias_subinfo.keys():
non_big.append(s)
if s.startswith("GLM_Sim_"):
sim_subjects.append(s)
- else:
- non_big.append(s)
print(non_big)
print()
diff --git a/multi_chain_saves/.gitignore b/multi_chain_saves/.gitignore
new file mode 100644
index 00000000..b84a8f82
--- /dev/null
+++ b/multi_chain_saves/.gitignore
@@ -0,0 +1,23 @@
+
+old_ana_code/
+figures/
+dynamic_figures/
+iHMM_fits/
+dynamic_iHMM_fits/
+overview_figures/
+beliefs/
+session_data/
+WAIC/
+*.png
+*.p
+*.npz
+*.pdf
+*.csv
+*.zip
+*.json
+peter_fiugres/
+consistency_data/
+dynamic_GLM_figures/
+dynamic_GLMiHMM_fits2/
+glm sim mice/
+dynamic_GLMiHMM_crossvals/
diff --git a/multi_chain_saves/work.py b/multi_chain_saves/work.py
deleted file mode 100644
index 25f82f9d..00000000
--- a/multi_chain_saves/work.py
+++ /dev/null
@@ -1,15 +0,0 @@
-import os
-
-i = 0
-for filename in os.listdir("./"):
- if not filename.endswith('.p'):
- continue
- if 'bias' in filename:
- continue
-
- if not filename.endswith('bias.p'):
- i += 1
- print("Rename {} into {}".format(filename, filename[:-2] + '_prebias.p'))
- os.rename(filename, filename[:-2] + '_prebias.p')
-
-print(i)
diff --git a/raw_fit_processing.py b/raw_fit_processing_part1.py
similarity index 60%
rename from raw_fit_processing.py
rename to raw_fit_processing_part1.py
index 7a6c08b1..82807f9c 100644
--- a/raw_fit_processing.py
+++ b/raw_fit_processing_part1.py
@@ -5,7 +5,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
-
+ This cannot be run in parallel (because the loading_info dict gets changed and dumped)
"""
import numpy as np
import pyhsmm
@@ -13,8 +13,9 @@ import pickle
import json
from dyn_glm_chain_analysis import MCMC_result
import time
-from mcmc_chain_analysis import state_size_helper, state_num_helper, find_good_chains_unsplit_greedy
+from mcmc_chain_analysis import state_size_helper, state_num_helper, find_good_chains_unsplit_greedy, gamma_func, alpha_func
import index_mice # executes function for creating dict of available fits
+from dyn_glm_chain_analysis import MCMC_result_list
fit_type = ['prebias', 'bias', 'all', 'prebias_plus', 'zoe_style'][0]
@@ -26,9 +27,13 @@ elif fit_type == 'prebias':
loading_info = json.load(open("canonical_infos.json", 'r'))
r_hats = json.load(open("canonical_info_r_hats.json", 'r'))
-subjects = ['ZFM-04019', 'ZFM-05236'] # list(loading_info.keys())
+# done: 'NYU-45', 'UCLA035', 'NYU-30', 'CSHL047', 'NYU-39', 'NYU-37',
+# can't: 'UCLA006'
+subjects = ['NYU-40', 'NYU-46', 'KS044', 'NYU-48']
+# 'UCLA012', 'CSHL052', 'NYU-11', 'UCLA011', 'NYU-47', 'CSHL045', 'UCLA017', 'CSHL055', 'UCLA005', 'CSHL060', 'UCLA015', 'UCLA014', 'CSHL053', 'NYU-12', 'CSHL058', 'KS042']
+
-thinning = 25
+thinning = 50
fit_variance = 0.03
func1 = state_num_helper(0.2)
@@ -37,13 +42,15 @@ func3 = state_size_helper()
func4 = state_size_helper(1)
dur = 'yes'
-m = len(loading_info[subjects]["fit_nums"])
-assert m == 16
for subject in subjects:
+ results = []
+ summary_info = {"thinning": thinning, "contains": [], "seeds": [], "fit_nums": []}
+ m = len(loading_info[subject]["fit_nums"])
+ assert m == 16
print(subject)
n_runs = -1
counter = -1
- n = (loading_info[subject]['chain_num'] + 1) * 4000 // thinning
+ n = (loading_info[subject]['chain_num']) * 4000 // thinning
chains1 = np.zeros((m, n))
chains2 = np.zeros((m, n))
chains3 = np.zeros((m, n))
@@ -53,11 +60,27 @@ for subject in subjects:
print(seed)
info_dict = pickle.load(open("./session_data/{}_info_dict.p".format(subject), "rb"))
samples = []
- for mini_counter in range(m):
- 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 # start at 1, discard first 4000 as burnin
+ 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('.', '_')
print("loaded seed {}".format(seed))
@@ -65,7 +88,8 @@ for subject in subjects:
result = MCMC_result(samples[::thinning],
infos=info_dict, data=samples[0].datas,
sessions=fit_type, fit_variance=fit_variance,
- dur=dur, save_id=save_id))
+ dur=dur, save_id=save_id)
+ results.append(result)
print("Making result {} took {:.4}".format(counter, time.time() - lala))
res = func1(result)
@@ -82,8 +106,6 @@ for subject in subjects:
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'))
- r_hats = {}
-
# R^hat tests
# test = MCMC_result_list([fake_result(100) for i in range(8)])
# test.r_hat_and_ess(return_ascending, False)
@@ -94,10 +116,6 @@ for subject in subjects:
continue
print()
print("Checking R^hat, finding best subset of chains")
- chains1 = chains1[:, 160:]
- chains2 = chains2[:, 160:]
- chains3 = chains3[:, 160:]
- chains4 = chains4[:, 160:]
# mins = find_good_chains(chains[:, :-1].reshape(32, chains.shape[1] // 2))
sol, final_r_hat = find_good_chains_unsplit_greedy(chains1, chains2, chains3, chains4, reduce_to=chains1.shape[0] // 2)
@@ -110,4 +128,27 @@ for subject in subjects:
json.dump(r_hats, open("canonical_info_r_hats_bias.json", 'w'))
elif fit_type == 'prebias':
json.dump(loading_info, open("canonical_infos.json", 'w'))
- json.dump(r_hats, open("canonical_info_r_hats.json", 'w'))
\ No newline at end of file
+ json.dump(r_hats, open("canonical_info_r_hats.json", 'w'))
+
+ if r_hats[subject] >= 1.05:
+ print("Skipping canonical result")
+ continue
+ else:
+ print("Making canonical result")
+
+ # subset data
+ summary_info['contains'] = [i for i in range(m) if i not in sol]
+ summary_info['seeds'] = [loading_info[subject]['seeds'][i] for i in summary_info['contains']]
+ summary_info['fit_nums'] = [loading_info[subject]['fit_nums'] for i in summary_info['contains']]
+ results = [results[i] for i in summary_info['contains']]
+
+ test = MCMC_result_list(results, summary_info)
+ pickle.dump(test, open("multi_chain_saves/canonical_result_{}_{}.p".format(subject, fit_type), 'wb'))
+
+ test.r_hat_and_ess(state_num_helper(0.2), False)
+ test.r_hat_and_ess(state_num_helper(0.1), False)
+ test.r_hat_and_ess(state_num_helper(0.05), False)
+ test.r_hat_and_ess(state_size_helper(), False)
+ test.r_hat_and_ess(state_size_helper(1), False)
+ test.r_hat_and_ess(gamma_func, True)
+ test.r_hat_and_ess(alpha_func, True)
diff --git a/raw_fit_processing_part2.py b/raw_fit_processing_part2.py
new file mode 100644
index 00000000..4a0835c4
--- /dev/null
+++ b/raw_fit_processing_part2.py
@@ -0,0 +1,238 @@
+import json
+import numpy as np
+from dyn_glm_chain_analysis import MCMC_result_list
+import pickle
+import matplotlib.pyplot as plt
+import os
+
+
+def create_mode_indices(test, subject, fit_type):
+ dim = 3
+
+ try:
+ xy, z = pickle.load(open("multi_chain_saves/xyz_{}_{}.p".format(subject, fit_type), 'rb'))
+ except Exception:
+ return
+ print('Doing PCA')
+ ev, eig, projection_matrix, dimreduc = test.state_pca(subject, pca_type='dists', dim=dim)
+ xy = np.vstack([dimreduc[i] for i in range(dim)])
+ from scipy.stats import gaussian_kde
+ z = gaussian_kde(xy)(xy)
+ pickle.dump((xy, z), open("multi_chain_saves/xyz_{}_{}.p".format(subject, fit_type), 'wb'))
+
+ print("Mode indices of " + subject)
+
+ threshold_search(xy, z, test, 'first_', subject, fit_type)
+
+ print("Find another mode?")
+ if input() not in ['yes', 'y']:
+ return
+
+ threshold_search(xy, z, test, 'second_', subject, fit_type)
+
+ print("Find another mode?")
+ if input() not in ['yes', 'y']:
+ return
+
+ threshold_search(xy, z, test, 'third_', subject, fit_type)
+ return
+
+
+def threshold_search(xy, z, test, mode_prefix, subject, fit_type):
+ happy = False
+ conds = [0, None, None, None, None, None, None]
+ x_min, x_max, y_min, y_max, z_min, z_max = None, None, None, None, None, None
+ while not happy:
+ print()
+ print("Pick level")
+ prob_level = input()
+ if prob_level == 'cond':
+ print("x > ?")
+ resp = input()
+ if resp not in ['n', 'no']:
+ x_min = float(resp)
+ else:
+ x_min = None
+
+ print("x < ?")
+ resp = input()
+ if resp not in ['n', 'no']:
+ x_max = float(resp)
+ else:
+ x_max = None
+
+ print("y > ?")
+ resp = input()
+ if resp not in ['n', 'no']:
+ y_min = float(resp)
+ else:
+ y_min = None
+
+ print("y < ?")
+ resp = input()
+ if resp not in ['n', 'no']:
+ y_max = float(resp)
+ else:
+ y_max = None
+
+ print("z > ?")
+ resp = input()
+ if resp not in ['n', 'no']:
+ z_min = float(resp)
+ else:
+ z_min = None
+
+ print("z < ?")
+ resp = input()
+ if resp not in ['n', 'no']:
+ z_max = float(resp)
+ else:
+ z_max = None
+
+ print("Prob level")
+ prob_level = float(input())
+ conds = [prob_level, x_min, x_max, y_min, y_max, z_min, z_max]
+ print("Condtions are {}".format(conds))
+ else:
+ try:
+ prob_level = float(prob_level)
+ except:
+ print('mistake')
+ prob_level = float(input)
+ conds[0] = prob_level
+
+ print("Level is {}".format(prob_level))
+
+ mode = conditions_fulfilled(z, xy, conds)
+ print("# of samples: {}".format(mode.sum()))
+ mode_indices = np.where(mode)[0]
+ if mode.sum() > 0:
+ print(xy[0][mode_indices].min(), xy[0][mode_indices].max(), xy[1][mode_indices].min(), xy[1][mode_indices].max(), xy[2][mode_indices].min(), xy[2][mode_indices].max())
+ print("Happy?")
+ happy = 'yes' == input()
+
+ print("Subset by factor?")
+ if input() == 'yes':
+ print("Factor?")
+ print(mode_indices.shape)
+ factor = int(input())
+ mode_indices = mode_indices[::factor]
+ print(mode_indices.shape)
+ if subject not in loading_info:
+ loading_info[subject] = {}
+ loading_info[subject]['mode prob level'] = prob_level
+
+ pickle.dump(mode_indices, open("multi_chain_saves/{}mode_indices_{}_{}.p".format(mode_prefix, subject, fit_type), 'wb'))
+ # consistencies = test.consistency_rsa(indices=mode_indices) # do this on the cluster from now on
+ # pickle.dump(consistencies, open("multi_chain_saves/{}mode_consistencies_{}_{}.p".format(mode_prefix, subject, fit_type), 'wb', protocol=4))
+
+
+def conditions_fulfilled(z, xy, conds):
+ works = z > conds[0]
+ if conds[1]:
+ works = np.logical_and(works, xy[0] > conds[1])
+ if conds[2]:
+ works = np.logical_and(works, xy[0] < conds[2])
+ if conds[3]:
+ works = np.logical_and(works, xy[1] > conds[3])
+ if conds[4]:
+ works = np.logical_and(works, xy[1] < conds[4])
+ if conds[5]:
+ works = np.logical_and(works, xy[2] > conds[5])
+ if conds[6]:
+ works = np.logical_and(works, xy[2] < conds[6])
+
+ return works
+
+
+def state_set_and_plot(test, mode_prefix, subject, fit_type):
+ mode_indices = pickle.load(open("multi_chain_saves/{}mode_indices_{}_{}.p".format(mode_prefix, subject, fit_type), 'rb'))
+ consistencies = pickle.load(open("multi_chain_saves/{}mode_consistencies_{}_{}.p".format(mode_prefix, subject, fit_type), 'rb'))
+ session_bounds = list(np.cumsum([len(s) for s in test.results[0].models[-1].stateseqs]))
+
+ import scipy.cluster.hierarchy as hc
+ consistencies /= consistencies[0, 0]
+ linkage = hc.linkage(consistencies[0, 0] - consistencies[np.triu_indices(consistencies.shape[0], k=1)], method='complete')
+
+ # R = hc.dendrogram(linkage, truncate_mode='lastp', p=150, no_labels=True)
+ # plt.savefig("peter figures/{}tree_{}_{}".format(mode_prefix, subject, 'complete'))
+ # plt.close()
+
+ session_bounds = list(np.cumsum([len(s) for s in test.results[0].models[-1].stateseqs]))
+
+ plot_criterion = 0.95
+ a = hc.fcluster(linkage, plot_criterion, criterion='distance')
+ b, c = np.unique(a, return_counts=1)
+ state_sets = []
+ for x, y in zip(b, c):
+ state_sets.append(np.where(a == x)[0])
+ print("dumping state set")
+ pickle.dump(state_sets, open("multi_chain_saves/{}state_sets_{}_{}.p".format(mode_prefix, subject, fit_type), 'wb'))
+ state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, save_append='_{}{}'.format(mode_prefix, plot_criterion), show=True, separate_pmf=True, type_coloring=True)
+
+ fig, ax = plt.subplots(ncols=5, sharey=True, gridspec_kw={'width_ratios': [10, 1, 1, 1, 1]}, figsize=(13, 8))
+ from matplotlib.pyplot import cm
+ for j, criterion in enumerate([0.95, 0.8, 0.5, 0.2]):
+ clustering_colors = np.zeros((consistencies.shape[0], 100, 4))
+ a = hc.fcluster(linkage, criterion, criterion='distance')
+ b, c = np.unique(a, return_counts=1)
+ print(b.shape)
+ print(np.sort(c))
+
+ cmap = cm.rainbow(np.linspace(0, 1, 17))
+ rank_to_color_place = dict(zip(range(17), [0, 16, 8, 4, 12, 2, 6, 10, 14, 1, 3, 5, 7, 9, 11, 13, 15])) # handcrafted to maximise color distance, I think
+ i = -1
+ b = [x for _, x in sorted(zip(c, b))][::-1]
+ c = [x for x, _ in sorted(zip(c, b))][::-1]
+ plot_above = 50
+ while len([y for y in c if y > plot_above]) > 17:
+ plot_above += 1
+ for x, y in zip(b, c):
+ if y > plot_above:
+ i += 1
+ clustering_colors[a == x] = cmap[rank_to_color_place[i]]
+ # clustering_colors[a == x] = cm.rainbow(np.mean(np.where(a == x)[0]) / a.shape[0])
+
+ ax[j+1].imshow(clustering_colors, aspect='auto', origin='upper')
+ for sb in session_bounds:
+ ax[j+1].axhline(sb, color='k')
+ ax[j+1].set_xticks([])
+ ax[j+1].set_yticks([])
+ ax[j+1].set_title("{}%".format(int(criterion * 100)), size=20)
+
+ ax[0].imshow(consistencies, aspect='auto', origin='upper')
+ for sb in session_bounds:
+ ax[0].axhline(sb, color='k')
+ ax[0].set_xticks([])
+ ax[0].set_yticks(session_bounds[::-1])
+ ax[0].set_yticklabels(session_bounds[::-1], size=18)
+ ax[0].set_ylim(session_bounds[-1], 0)
+ ax[0].set_ylabel("Trials", size=28)
+ plt.yticks(rotation=45)
+
+ plt.tight_layout()
+ plt.savefig("peter figures/{}clustered_trials_{}_{}".format(mode_prefix, subject, 'criteria comp').replace('.', '_'))
+ plt.close()
+
+
+fit_type = ['prebias', 'bias', 'all', 'prebias_plus', 'zoe_style'][0]
+if fit_type == 'bias':
+ loading_info = json.load(open("canonical_infos_bias.json", 'r'))
+elif fit_type == 'prebias':
+ loading_info = json.load(open("canonical_infos.json", 'r'))
+subjects = list(loading_info.keys())
+# error: KS043, KS045, 'NYU-12', ibl_witten_15, NYU-21, CSHL052, KS003
+# done: NYU-46, NYU-39, ibl_witten_19, NYU-48
+subjects = ['DY_013', 'ZFM-01592', 'NYU-39', 'NYU-27', 'NYU-46', 'ZFM-01936', 'ZFM-02372', 'ZFM-01935', 'ibl_witten_26', 'ZM_2241', 'KS084', 'ZFM-01576']
+
+for subject in subjects:
+ test = pickle.load(open("multi_chain_saves/canonical_result_{}_{}.p".format(subject, fit_type), 'rb'))
+ if os.path.isfile("multi_chain_saves/{}mode_indices_{}_{}.p".format('first_', subject, fit_type)):
+ print("It has been done")
+ continue
+ print('Computing sub result')
+ create_mode_indices(test, subject, fit_type)
+ # state_set_and_plot(test, 'first_', subject, fit_type)
+ # print("second mode?")
+ # if input() in ['y', 'yes']:
+ # state_set_and_plot(test, 'second_', subject, fit_type)
diff --git a/remove_duplicates.py b/remove_duplicates.py
new file mode 100644
index 00000000..97f6683b
--- /dev/null
+++ b/remove_duplicates.py
@@ -0,0 +1,86 @@
+import os
+import re
+
+memory = {}
+name_saves = {}
+seed_saves = {}
+
+do_it = True
+
+for filename in os.listdir("./dynamic_GLMiHMM_crossvals/"):
+ if not filename.endswith('.p'):
+ continue
+ regexp = re.compile(r'((\w|-)+)_fittype_(\w+)_var_0.03_(\d+)_(\d+)_(\d+)')
+ result = regexp.search(filename)
+ subject = result.group(1)
+ fit_type = result.group(3)
+ seed = result.group(4)
+ fit_num = result.group(5)
+ chain_num = result.group(6)
+
+ if fit_type == 'prebias':
+
+ if subject not in seed_saves:
+ seed_saves[subject] = [seed]
+ else:
+ if seed not in seed_saves[subject]:
+ seed_saves[subject].append(seed)
+
+ if (subject, seed) not in name_saves:
+ name_saves[(subject, seed)] = []
+
+ if fit_num not in name_saves[(subject, seed)]:
+ name_saves[(subject, seed)].append(fit_num)
+
+ if (subject, seed, fit_num) not in memory:
+ memory[(subject, seed, fit_num)] = {"chain_num": int(chain_num), "counter": 1}
+ else: # if this is the first file of that chain, save some info
+ memory[(subject, seed, fit_num)]["chain_num"] = max(memory[(subject, seed, fit_num)]["chain_num"], int(chain_num))
+ memory[(subject, seed, fit_num)]["counter"] += 1
+
+total_move = 0
+nyu_11_move = 0
+dicts_removed = 0
+moved = []
+completed = []
+incompleted = []
+for key in name_saves:
+ subject = key[0]
+ seed = key[1]
+ complete = False
+ save_fit_num = -1
+ for fit_num in name_saves[key]:
+ if memory[(subject, seed, fit_num)]['chain_num'] == 14 and memory[(subject, seed, fit_num)]['counter'] == 15:
+ save_fit_num = fit_num
+ complete = True
+ if len(seed_saves[subject]) == 16:
+ if subject not in completed:
+ completed.append(subject)
+ else:
+ if subject not in incompleted:
+ incompleted.append(subject)
+ if complete and len(name_saves[key]) > 1:
+ assert save_fit_num != -1
+ for fit_num in name_saves[key]:
+ if fit_num != save_fit_num:
+ for i in range(15):
+ if do_it:
+ if os.path.exists("./dynamic_GLMiHMM_crossvals/{}_fittype_prebias_var_0.03_{}_{}_{}.p".format(subject, seed, fit_num, i)):
+ os.rename("./dynamic_GLMiHMM_crossvals/{}_fittype_prebias_var_0.03_{}_{}_{}.p".format(subject, seed, fit_num, i),
+ "./del_test/{}_fittype_prebias_var_0.03_{}_{}_{}.p".format(subject, seed, fit_num, i))
+ if subject not in moved:
+ moved.append(subject)
+ total_move += 1
+ else:
+ if os.path.exists("./dynamic_GLMiHMM_crossvals/{}_fittype_prebias_var_0.03_{}_{}_{}.p".format(subject, seed, fit_num, i)):
+ print("I would move ")
+ print("./dynamic_GLMiHMM_crossvals/{}_fittype_prebias_var_0.03_{}_{}_{}.p".format(subject, seed, fit_num, i))
+ print(" to ")
+ print("./del_test/{}_fittype_prebias_var_0.03_{}_{}_{}.p".format(subject, seed, fit_num, i))
+ total_move += 1
+ nyu_11_move += subject == "NYU-11"
+
+print(moved)
+print(completed)
+print(incompleted)
+print("Would move {} in total, and {} of NYU-11".format(total_move, nyu_11_move))
diff --git a/test_codes/bias_shift_test.py b/test_codes/bias_shift_test.py
new file mode 100644
index 00000000..597c7071
--- /dev/null
+++ b/test_codes/bias_shift_test.py
@@ -0,0 +1,38 @@
+"""
+ Studying how much easier it is for small changes in the regressors to affect the pmf around 0 versus further from zero
+
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+weights = list(np.zeros((17, 3)))
+
+for i, weight in enumerate(weights):
+ weight[-1] = i * 0.2 - 1.6
+
+contrasts_L = np.array([1., 0.987, 0.848, 0.555, 0.302, 0, 0, 0, 0, 0, 0])
+contrasts_R = np.array([1., 0.987, 0.848, 0.555, 0.302, 0, 0, 0, 0, 0, 0])[::-1]
+bias = np.ones(11)
+
+predictors = np.vstack((contrasts_L, contrasts_R, bias)).T
+
+for weight in weights:
+ plt.plot(1 / (1 + np.exp(- np.sum(weight * predictors, axis=1))))
+
+plt.ylim(0, 1)
+plt.show()
+
+
+
+weights = list(np.zeros((17, 3)))
+
+for i, weight in enumerate(weights):
+ weight[-1] = i * 0.2 - 1.6
+ weight[0] = 2
+
+for weight in weights:
+ plt.plot(1 / (1 + np.exp(- np.sum(weight * predictors, axis=1))))
+
+plt.ylim(0, 1)
+plt.show()
\ No newline at end of file
--
GitLab