From 1330eee0e9cd721d547de0f07c080f9f692e5783 Mon Sep 17 00:00:00 2001
From: SebastianBruijns <>
Date: Wed, 12 Jul 2023 16:32:39 +0200
Subject: [PATCH] code update
---
.gitignore | 2 +-
__pycache__/analysis_pmf.cpython-37.pyc | Bin 9525 -> 11418 bytes
.../dyn_glm_chain_analysis.cpython-37.pyc | Bin 53646 -> 53918 bytes
__pycache__/simplex_plot.cpython-37.pyc | Bin 3111 -> 3129 bytes
analysis_pmf.py | 191 +++++++-----
...ght_analysis.py => analysis_pmf_weights.py | 109 +++++++
analysis_regression.py | 12 +-
analysis_state_intros.py | 8 +-
behavioral_state_data_easier.py | 77 ++++-
canonical_infos.json | 2 +-
dyn_glm_chain_analysis.py | 276 +++++++++++-------
index_mice.py | 6 +-
simplex_animation.py | 17 +-
simplex_plot.py | 8 +-
14 files changed, 489 insertions(+), 219 deletions(-)
rename pmf_weight_analysis.py => analysis_pmf_weights.py (61%)
diff --git a/.gitignore b/.gitignore
index b84a8f82..e0d406ba 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,4 @@
-
+summary_figures/
old_ana_code/
figures/
dynamic_figures/
diff --git a/__pycache__/analysis_pmf.cpython-37.pyc b/__pycache__/analysis_pmf.cpython-37.pyc
index f9fd33bf94c4befd1c5937c83c151a084e749d91..899a197bee3c2ffcdbb2f048e7d0088d934a594e 100644
GIT binary patch
literal 11418
zcmd5?X>1%vcJ5<lI2;~BQWOtSmPfV@(koFUMa#PE(czV?iz}Nt+mgpkR&m(moYhTn
zIF>gKC?gJFJM&``IE(xk0)~N204E6iBLST3A^{djkQhc1Ae%-K*x3ZgU=74SasZ1U
z$@i*zhBG{5C6)~}(cM+`>Ui&c@4b3m!&7Z-5d~NAw?1Hh@Uo)(lq%J$4TTf9?aY9p
zFomfp#qi#0S~GOjD+i3Aq=k$orZJrbSnxM=Bg}*4W)@;iEX<l&1U*_<D{EuzEShXt
zQH_XC4}qd2188Ys9jx=NVzf$17u&+R*;dv=^y;4+*>-O_Q}GviL+@9NcJ>Rbm-XFM
zjVSA9+c0+r8(`b<?&MwUk$J5&Am+GI-om@t4)*A>1}?U8{FQrFlmUgsx)oz9+iCQ2
zeXb9y_JfZodu&DLTba_Y@Gw2u<8MO~=&_CM;=OG5+nTHmusy4_@^*}jtq`TD{K$$z
zJ<2=$_oHm@iYi-TUdy_*qH--fAHZ6JU$WMaTno<!u+|d|Yk7M9mapfXzTWrw^~e1B
ze!u=WG`b7g9%cvF$lLl&&Dd@1DSX6tveAC!S>@L6pk=Qy$oKFLKDey*D^hDi9J-X)
z4D!7sNtMdRBuc^Z6TFrBv!NA@9h?u9bkSNj$G3oTeUEWI$U@M?lYPpJ#ts2pV4WGc
zj}cGx){Q7>`FE}<h2Q1-_?~4I+`kX*_Z!1}c<ul@%#I|rex*FZj!JzEmq*FYq_ym4
z%DfI9zsHBmV@buQAM{F@o9vh|jycA_^^@%QipnE=oITB+c{>1m8lF3ddC249@Hz6)
zS{c2hB<OR9onR+@s~w(uik+&~N_zgnwMIKVU5{r{F%B2r;yWSDUE>HiI0~Jr#<7CI
zck&|=qX{3QVEH(gc!2&if0{izA1DO`<Hwj1BoB%{rRd6WC9@5*Z_g{Bje!=er<DTv
zF>rIZ{0vr#^JDy2`2_jsib_0zx=r$Wk{|YYJjKrVrETmvzcj#}m!*w-@T}iQ`XD=J
zoMz{ZXJKs@U`5Y@NAma88hlYxr2To3VN<)~r`bh`hAv}(hpJxwHqe-m2<zoDKCXKC
zxfP{`@+yyLg_V2<7V<p1<a6n@z0lCM8c$waBMZqf*&s)?_HfH_xj~L$iN}=&IjXgX
zTaK#@azv#Z*Ba!g)*fy-t~bbWNy>4fL5^ze;g-W_kmH<`qv{c!pI2-DXF0y9s4;g2
z5hR+_*h|%ngQCaW7+P<#iKLFqi|D$fQ+~1<Py9GB_b9u?US<g|UnzTW?}6GT`-VLB
zf#~9ub-MV5?E8UmZ?5Bh3VZKa;~ZsC{477$kVR1*va>D^Il}kzvuvnaK}OsnKEI?O
zzuC`*yA^(Jn^G<HAfG{&b%-Bf)1J*0SZ$}jUdzV3aWW_Ctsh@&@#mp@dedFg%hpP*
z`R&(~QcM1?kvn}<%bh6m+i#pNyoZ+aQi^xca^ZfmDM%UZql~K6cKHRFTeXTWHq5O~
zF!t4QD{%jPaNnALb4?b-8=}v>aQdejPCt<{ui~VX+Qhr1w)~yU4vF0*U`K5qR*~#e
zjQtr6NUkE=CVu?TAlsJz59kK`KX7E*?GJrq+r^)KO~|(A9dO^CpIU?L%vVI#Loq5u
zDVo3PBYR#V8x`+Vac)4CMrfs|C_gmFM)RLQH{kz)BimuW<0IQ4-us%6?Z{7n`;Pqa
zHOMBvBC;Neo#H~NGe6=ZJ1ddx6tfap^k_hqMrfr@QG95S?acosbOZh$II>;#tdDG$
z`0m$)Y*&5|+;`=B)*x%uBkP|*%=u>!7ybIYU!TBk`W?BGdgnE#Sc*oM)WSF-^(B6h
zPb_P+BbC%!F^br@2aOkS*7G8=vUAu;qrJ+CJS{4p;xEwH!U>!(Aj`y8^tdsJ`MrH{
zkiIwllZ;99Kn-~|on}nZS%-HvQ`(Y$MM^`bOWv1dHOrGLu{3tNE#fCwzy8b8t5wes
z7}gxXSAV*iN8%jo#A<DmtTKb%o*a~QZhFqZs{K5EvaskaV&5eNobOwML#`eNoGn~t
zig5+^RovI0;Sr9L2IzH<afIb>28`>*4Z|>A@-h(PX5nLg6FU4cPJ{-bbvhZ683*=w
zt+Z8a<Co>h0nQ;vS4^-_GC&%ZtIG3<#&*^|A1rl?pR<>Z8^~1FjRRe!c|pJKuN1{f
z-H1#Z_YR`>^XUE3CUlj?f^L+$^Y1}7SLjs4xWcdU>lk^L8>OwexA=|H*6dsSB~Y&N
zYedPuPZT_f^1h^8!wh#Z!!2f)U*=a<)Mc{PTg*WXU!q=(Uz2aWe3DncvJ;Xs)O;l9
zjG}s8@d}+x(rFVrX_&rOEcjkA#co66o8+2oF{819G0mr00s83?m)WRyhG(h9&O(_`
zz4b26kK~yVaC}qZ_<9|VanJreKUvlc#(Pi}@CSJ|hZYXbXXL#X=S+M)5LH$*l-RsN
zC1Y>hdENcwojv&!Vnt8hfEM?V@70~i?ICY%sN<}(p^pBWoXE+so6<1PrZk#!Q!lHX
z=6UnqPd{4ge)`e6>rW`XbvA0Jc{ge&9VTH<-kF!D!}>buu)a>~;$)JWepHeu(CK8Y
zrT+YUiBBW0dNEmAs&8$5;v+u0OqQz{7AwL^=Zty!B5F8SPT^bPzoZ@hn{lM@dF`9#
z(}p+_MJ(+Xky2m&3s_pZW|KG_?uDQ9QM7&)F`-`^=9ctbiwzkWoE2xizAWn5E51Jh
z%c_y#b38-Q3BEkX=TS!dKzW`&2zx)zGtl4$w$ASO-&T+VHr|*2N7(rJH8$?mSMiYb
z4RH0^zub<zejA@sd$fA2&#&<?uU^j$d7!_4fp^q4<c-=lTa#$@)^whc1+JdEd-eJ?
ziN3n##&#}8ZL;Pq_8Q_ZA#JVbe(YU6vHfnn^)hu*8G5}!zrszp!KC0k;|8Wu3*6NO
zXG6MBL8-`xMy3D3m89N^^63x#!Y9?jkA9fEe&oYHx&4Q~T`R=>K0m4!KB*Qyjf)^=
znWnkvqBs3KcpA4I!^2U^DsIiu@8IV^@8>?+NXh+Hi0$Bn<W5dg`EM{LnlNlA;A)wi
z8%T3A<Lc&=?RMq3n8}K?nVIH^T-M4s_K@o8x%A8r6w!_$KOa9Zzch8gOy<}DdrB-=
za~b=9HJeT+rueLRJ7GJf!|elR#!M~QmYu-J;oOqjf_K8nCTRLrOk<OV6isyC>5KhE
zJbB>ncZGiOSN}Dh{N<Cc{PQP|{VbmR<9{E#`73|)Sv>hmJ%9W0e@^~Yd|D2~C>=K`
zMe5)X3M6?yG|^2Z9Emt;H{`#gqJ|$Gsr;G-&1jBR*6-ja#a{eO0nG&JG*eYd%XBv5
z2-9}##EWiS<$`F&LaxpWxgp(cTI5!8*0J5@Ma!8@Oj)K4P1%;zllO>H_jbRZls`SQ
z1#NZ!4@KRjie9{HxRPviOPK~ya3IUF`kHzNKc3DjqMK=xkcjCh$5>!eVL{v^MQ9><
zP+mWcOQ4e<zxkK{@M%lMGaM1ZTtlkV_K+s%_hSP#ySSFqutaGo$B8SOa8lGT{sMGs
zq7U^NuVjNY{lAWir+-iPeb700>?3|j2Eu~{ydHhZvc`h*0sG9F5h08y_<cba&guy2
zFC0f@O=ay4Y?>UE=9W4L`#NZzkobA>fB)C%`cqH-<Ok1vR_y-ZZ{kU>`LlR&D6x87
z7TJZV6t}Z!9=mk$%+c7ZudNajp7L~8Pn!idm=o{`Hza+;fl5i2Yt`6VRbwKXVmr~n
zCNo#mZYXhB-LAGnp-|^`w7KC#B5hikM4~t#T_@o{{dNM2aFNTVOvlP*hI6harYa2m
zD_!*XH9X*!!Rc%&D~bcVW7(M)jf<JnB5T{RTsF7h4MnejZ54U3|7>kW&SEyVZA=_<
ztTZn^GJL>ZNT*G)l$f!S3j*#dC!4iwCpLO-@|e`lz55(2_M52`<`%X?@xYccLf`gq
zt{5Q_HKN=>Wpbm$Ov=ikr9w_Vq*X|L6>>d@(*zz#A|rTya%6NIE$Z};mVD!9FZ|w)
zHxMyCm+wkEIU}!~Cw$C*5U-HuC133(WB6X&4GQvJS9h|xGZhMLl@@yJr)4zpI9_g8
zWD(AwC_hB~OT5^2X;6+^GzGJVD%;S_Xg*yH7e-{>22)a?nyn@O;?L>&yLhqnqE&$Q
zNkncJn;$3K($+LDj-;*3;O<z8&o~|x5XnHgFSa`-*K2GZD(+nmT25qhSwZY3vd%0Q
zQk4~o^A+-cqv`bI^A|=AjY2w#Ev}YYtk5pv258+1g(NpLm35qKy4Yr0nPiH`c)?7=
z20{}zuIwK<)Yv>`v{_TB1cIN#Gl|jhM(RN$OzR~^Mh`a<#v6w{)kuA^g8ykzDim3Y
z?X-FhTu5^qy$q~d5SBA`&J;Z3tf7yssUL)#M=*^TqX|3dK@qJuuwHX!#!8zh8-m`0
zICd}M!Akf`<{mFbr)SMf5_+xLdv|?9B4egGc|q4@^5iOs0+L<n!Q^5$i9mK)e004C
zwf^KwLxD56gk-5Zvyh(RqH1a}%SL!B9)jg9&RWy6sKl^uBXq=ag5S2X3wCM=u%lwz
zzOBYAL<(YR%yZFIPae)e$%K{IMSL~Gfn`s-*do2tuXL^F7vpFLsSnn_A9vNM;?a61
zAYs7azD6IYwPyGN>>(BN8;X0@ufwdF8BW$>AusWzK-eaphjMNxXHCzixZ9aFom?vG
zq^v1qD^vo>d+WKB;|8I`+>*EZH`&~r$_gY*cE)xiq?j=oh~1zmgt_FB77~;}*n$jB
z?z)2Fwjd*vJHFgXYa0L|pGd$BvnF!`*&Oo9{!KRnd|;FB&e(21@)=H9<Pg)`4K1b2
zDV}me1@9eLN?B>QiE%q;Ww_hwP`pc|%%$vt<2Km<9ibGPNGfY{x0je+1-2HV-7pbI
z!EQU1R>QL!wil*|1>02<$Q}#m!7ycpGe2{bo1+_CL@eS6x7mdZ!*rPp1y$}gb10Py
zDyo*Nnr=w)X}c{pcM`WP8^MKR5mclg5h@iN00yngR#`Th>6M@@!3~*+4oU8lJuP^{
zWZCGP;fsiiUTDP-hSYtYFS)@Y$G&N_d2v9>$Zg@pMw!-y%p6c*$N?<EeIe}zXm@mj
zj^(7dI0#6LW)Ct5<pv6*`OuQbte%`U-8SMFQ1L~JIkRpklMvZOo8GV;+trX;sB><o
zG@be|;r6U;@&W}?5V=>odeK4#PbLIRJwRMFXN0^>#%N2!qd68?p#9&7!Y+}PBVapH
zOjnyt8zFMJ#Hi6EcN!GDhOuhHn5z~9?R0`-v)g1AZj&U)EM~KdZba(V13;R{sudX#
zva*C!lF<YsN<e>p0Rp$AvPo;&OeIV*3<)pkEs$TZV4zi&oBcu}Wz>qo)giv14ciEj
ze3(r~&QzC)(T*-;{vAG^&a#CRe})2%eH=MOSZ$@eR&P~fT9+0M#MCFWNFW>uw@0)N
zwL=T2k0D8G_j;)jHLQJE9iuR=TkTSZ)gE<BJA_g<#z*kiiuwAG>$jsNOtZBHWP2A{
zBIp;vh#j(pdQ-lQ=ea<DxER!0<t*)J?XHg3H>Lc4#%w-YP_Bby)#pR+Qe%=oa3ggz
zcrO}pUOz%72}ov5-GjG*+tX1W<n?qk;>oU6dxp@vUDF<>F*@N)xW>SJ1e}e7_b~B^
zc6!Eu8R<u&2;F*ro-Gl&?f3xLj#j#<OeJdHrp9y{9}Y%fOR!HlwgVK;YFahw6~Vlu
zx0rgAEERt=+jezO{`JAm`ebVd^^$xMA0f;|^N_XS;{)~{Rzq4JR%z9|*@oq)KFqQm
P<7nP?%(SLV>#P3@FI;SK
literal 9525
zcmd5>Z){vwm4AQ6<MG&@aokBfZW6~0b+FUcaU3^kLz}UikcOl+PUAM|qviG4_j>&7
znKv`{J;!5{u0UylU8!O|0NSniFh+m`9{?eQKtfB!e(6dpRCSSd_o=&T1-oJd0)&KE
ztq|}#_r3XNVu-d$l^MPF?m73|d(ZiubMHOxzQ=ld;|hL1e)9tR@-vF^6C&!rUKFNr
zxvvf=3R9SxQ;guL<~2iCgL1^^kaba`lW9z65!Ueq-H7pyY8Q*LP8MTbEDnio*28*P
zf+e%vTdENc>!YYqvJr5)Ss&Z=l4A78ntpZ{OR?Q-57pQIvSsgIu4gL#(U<i>#YixR
z4Y2e}s*z-a>~8ek$M&*)c<$o;Z2ziO84#c6O7$+DVngh}bq!kV=J;3lY$<ycHk?w7
z-RvG?Ae1>6%A`a3eIY#<()Wh+-CSSU%aUw_9b)f)MSo5+_8I%j-{)yII;fmcp8qB|
zL&gEVpZD<t*VRE~Wf(Ty!-uw1rVJ|l03X6VR<C7auK;UE9OOMD&kk>C?7r1#MHR2N
z^zpuJ68G~1EIO#L2hz%t#>UwZ>cyG3i53$BEiEc)>HKA-{4IVj-+x_&_FsYaBgP?q
zXyyIvD4Wb`gGzOj9otf&<Du#pe3jhNP#$5*s$S8=9Jq(GN?3niP;#DQQ^x)1<1n;-
zfIYaS@;JYr9cL$AiBLZ)_n{w(mWL286km-pq@*Q~8D|f(lb~p{Ln}wv2kKNsE5$A=
zT(3@SDX7z{N7*BCe7!nJ9QLSjjGZ#3;E@cFo&rjQdy`5IDw-mFr2l2#@MCORYOGfu
zY*?;WkE1*zQP8U=LfiG~Lt9D%iO%Z7p>E&;$;W`=Np?EaGT@zQ<27yh;5J*xj<aoc
z)Tuk$j>p^V7?J)s*JekZy0h*0P@5e|X~+3CJL=S(ZO4T+J5EbGo@le9PTkpdJQ>>Y
zfpD(h!uj!Dtr{zjU~QWX*EW6SQTE{$>IAJcS&hxr>uF_KISlSaHj~w{>r%Zg>uKj`
zt~Jdy<f%4{ZfVT2_h=^Hc8o5yVe}r&hudz|e4ABIV>dr#WGG|dr+B6<W1&4i-LmIT
z^AUcE9Zo68oZ=!iry!RYVIwJpXYQu`UY7PCcR)rn&ZpTkL41}OPxEP+`CMsaK7q8%
ze+F9HH#lKGH<1}t;-xc}l}cCX_sCPe-^f!a%Nj9e%3lI!Mq2Ur;5>FSnG38;rYU>r
zb~XN(%wxJm`i6PTG`sxn@)&6ULuel_eREsJa^<G|^wWl){$1MK^pjHQ7Oz#hOJ8z^
zr0-6{cf@<E?#Z-bj5r!R+4O9;_~jkr*`Cr6uz~)!+_OFITcKxr#BXksXM0Ltfc8D5
zzuV^7XWyk~126W9S1Y}xzX?5iN_w_ee5dZ`HqTNEt<o!AzhgYxTl!yM1O0EgXA|ys
zLeD0|zuYFzCQ4VKeWJv+dDeKBo(;U16ka7+x)OT!Bhs@;QIVd7M4M-+g;q(5FWxbp
zO_tsSHqif;d$!N5gr4maU%5@5?JLbe`@YhtZJvFu)wALI{`2Ac{^^kZXh@&IZu&vF
zlhQk}UOmf9YF$ywQ&_1V=coCZ>l*Dy74`XHTCs5*FwQ~pL&)GB#ZH<WP_`&@LtB21
z+LovBHjk_u@846#dGsIbiwEevH9R#q4+&7n*XiWoJe@1B7wDb8va9rIX-)GDwVn~m
zk(h;@u3yk-t#7CSRexV+81qrQVO#Z_Q|2iRYKIJT0n&jTl%ekUPS5Icfj*0f=v`ud
zjx?X&<_Fg52RcJwCyfgs{yfD03AO}8cgSkkvZFEEc#=QKY>d+{{)dg?{TU}SwrZry
z%RCK?UofVFoW+<e&+%Ds&ehLr$O|oJHmx3z9%}P~<Va7*Gn`6)=?{pm{?adC)qe6;
z%X!Iu>bZ?Ryp4wZ4(BDZ?T&TmXU96~b4R~4&RT;0Z$^&Bx*0hd_s`Fdf^lxdr+#n5
zCk#K_IBf|qY^4dqR+`4exe>n*#<BDQogFneEvH7j?M*U_nBAhE&t0ceB*nN$d2`!I
z-&4pqe=1|>En~X;<3{HDIwF$J{A7O1l9+M3#g)qL(kYyOT?%3bcbwM^;N2}vb9!Y}
zh0G`UMV35)k=S8l-bwSP16h^~q>-gYuB*m8f2Pp~#iwRm=6k?d#cCjX1?LJzdRAUT
zIAOf98c8Z!8cLT}6)G7+e1^{hof~A_EEk+XlD>%F-K7b{@65LN4d^D;lD>h@0{&l*
zkwEUo=Yv6As~=~!&7Od6)rMwQxSBzCP&d?#_&beBb*+&y!zj9%C!PV_IwsXO$GmY2
z%T_y#Y4D5JgFU4A=ALg`0o6;&SfNHI_OpIRR&eh4k%inE_jSS9m@ZV*)Z|U0(jW25
zlC+|{{e!UZM!oRb53*Mt{`Tuv|M44*LMD`XtzLMeUU)krI?&4^^-Vu|(x-&exZGjf
zJf*7Q(meedJ{<&~2C@+f@GW(V)b5y-+R1Kee}H156V1jVzUCDDNS<4cuUiYQ-(Tcn
zsUY%}v&hY2!FD`%O!f6*e(9@<NTA98O&wX?SU6&3i|mNIAlB>^$30>%=kw+QU$(BA
zu4j4NJz_alZo{=*6D`M!8-6#QrdKej`yO=T5<-e5`fv};{&gmM<Uj9?e&Tb#%4Gj`
z=7s-$<KWLT*?;=&=yRX>;X9e^r}q5EkA8RQUo(rcAzJCU2o#CIaTLh%K_HQ$5?+fv
zH3j>xsi5$Ig32#yKt}Vls(uZh9|rK*0mwvJkg0>xy;yKOVY!}bKIOL%PlzrI<m<d#
z9Mk>Ib#7;uJ=gDAx4mU^!M0o=<=PUb;FcuuL?|ci&oA8to*ThUQSVg+zV<ZzO19A@
zZ5l<vgDuz97u9R{(6XwC6w@xjBBrA}%p#W*)`5$xh|Xq@$=^>h($U!;efrJ+db_(8
z1dfQJuQ64kJ*Elz)?^@N*SBjLlU~{=a?;8rKPhSq{|H^0NP}*R!3O>xM@@i#fPEU6
zL&r4fOE!={SR{~0W65M4s}c9?wiZ#eD2KAB7hbak$qUa@S!Y$dhM1;_N_|TVVj%`y
zv(kTo{jdH5{r)+yKl}B^-`Pxk{ogX#fc;KpbIff1=0#y`A;;Z9o)6DGdG_S+$6jn!
zngMvaujj3@-%%8p6Mj_A5f6|OmfK<Msl#X%a_k-mxD@7!ngWR9>OM6AKmq4I@cfu*
z<}KSX&CR`X>X;tjcTEh!MX`{xJiFkG7ky1E)EMMz{kUxoTDhFLWDD1$`Q^$s(BB;|
zZpNvKm?-yAnZnXs%GpJ5YBXuav>K^Xqfo}&?p(67uD!|qjtiG2j!uE2E(RuR*_VHQ
z@w-DWV;cTKK5LomlKkB>9Rm5)OpU@BjUv=X_G2kR`m-6oL(tUlb+1r7TcdTO)=jqu
zNLUjT+I~zFuvWmlaQ^=+v)Ma0D%-7Fg1KY0yCG(DovzPCBfjcFLup-&maiIxwb6O{
z(zz!ljvqz8hp1mIw_c;g(U0W#l2@ad=0_I_o>$0k_PVx{&GBJgw(`SAH>0!9KK;PN
z@pks4(Q7Z{Ow4JIJLb`;_S$1cjK(u3jy}|0G1cDeM0@RnwPDgDRyU@$7i5j*&t`&@
zDMBgoJd|puC%4hZVCG3Q6Gy3oU38;)y4htpcHYXl=%7)YY!;8zVr^1Wo5{sx%gF+c
zdbFfk8Kz_9ImJ-_JZV+0qUk{!LdRkuR}i&5=u;dd9fqj)-i8<jU=V)IC><|Cb+TB{
zl$eKZxtz4w-Rzd37Z71&eJ>+OV_#j^Ofh?Di3{#{HujbPrb$da>TB#Heza&WuI9ME
zD{pzlT*1rP3)pd}gsR}F7jvH9;jHD08^H?R>2hnKAh7AVj_bz>`AKP#-(d-1ZTO^|
zNxOh62wu5z3YvR@AeXDQTn-zPnD8=9#9YB*exy)DZZNpxgcgJ7#MhQwe}KkmV$_()
zevB#zMn6HN<}CJO?%D!b=lZINO}h*U#+Ph|vv9+fYmy(42z0D31365`!kWX-spHV3
z<9-(h8eC9OwSCp{qizvhxPG_GJ@cyVVhV5!1!5y-=lwX9${vsgT#_rxMi)Iyxqd_~
zfq>#+b<SGgxv*JK?(t|zFmu*MVa@Y9UF#~e&LWg@1(%CSAYya{h(R?!B8QEXY4GSq
zFtDCowERxSrK)}}sgZ-y>o)V2u>dU>)*%`U7A=QQ#}tKi*N-~zW)A!FD&O!mEad8n
zzf1bKbr$%0wzGnHfc4(at`Ux&?Ex2k(?-@rE{7d~H+;3|$K{v-8q#Y`#~N|6*9@Ag
z<p5b}T%@#S1CB*F2ur@NV<bUqr4e0oxd}~m+Hmx|?HHYi2eVl~U3adKwHK|NX*o>L
za@n#K(W)zC0up;9UtnuFev0x&_iM;7Vrma<?0Sznto3WL$gp}?i$`MdSR$_VseM{x
zAMSXj)P8kb-J?!w$5Be5RUH2wEuyB8yd=PhfsSP(k%)Rw9o2fIL;~DYy-we;=Ks@X
z=iWNxIM`Y1dL*D<zdB5Il8(?BHGi#MTClC9PS*6JH%&c+r-3WLBHh4ZA{?VG8HIF0
z)9$1;I{B0QHw^9L&}<5RjFFz;(|s8JgXh2lTEIDi2q1fi?$N-*bC#k4|9$d53Oq*<
u3G@&6KqhDp*(mx*;Yr9v5MePjs--o;ZWu8|@@cf*hZa;$fNC!Xe)|*Ak0|#5
diff --git a/__pycache__/dyn_glm_chain_analysis.cpython-37.pyc b/__pycache__/dyn_glm_chain_analysis.cpython-37.pyc
index 797f6ea4ca71054cc38cb50347dc59cdea7fef0b..06557f92116c59895d3c014e28181f1e3c4ae7f7 100644
GIT binary patch
delta 17789
zcmb_@31F1fwf}sR$uL<75SFkRAS_`C`woJvDiGNe8NkVq`y~uXX2P8bS&SVOT&g00
z>xP0R;=XI@Qme0OrFL;^Uv2xoXJ7w)wYFB<KD$3_?f>^X-z1Yn{8#_)1-|p$?=I(_
zd(OG%oV(2RCo^upJ0pKZUS75be>?7fM15R!B!8IqsjSB*`K+Y!D64j2aH8#NpQQ5f
z&!armz78u`s`{&fgP!1I+qZX$8lVQML5IA-srZ_v3f164mgP}Hf@St}dxkwz4Lx9~
zVX8<CS0mI&HA;<EXQ*N|Mvbk`Ip95D1!n>8I8~y?s|jkNnxsnAWHm)iRmo|`bAq$c
z`5b%p-nr;@UW$+|RGw~B<EnE|J|E=^)X)wKq@K1;l{dIh>P(S3E<s&b##Lvbe3A4q
zR+XvgNi{>w<e0nuNhgcb6;4;KEd0yABBRI?T%zt)v()T^R&eQBkD8<A9`vZW2fV>E
z)jTyH+*^j+1!^I37b17Lx<@Tii&6V5<S$W64<<coDGFAoGu1LwU#XU>v+!P}?p3SR
z8dO}Z)~d6Svqqhx&c%DJTBp|IeYV=5%JDu&ou@Y9eXiQ1Hsifc-KVywt!TL(Vp@e~
zrCOm1t3B%cs(`vcZ9~Zhb)njhy5;KrqzbAX$U09wpmwTVNNrRPs*6+zsZFXvg;8y@
z+N~<_-lCMU@!qPcR5jk`t34`$_XTRNs=<4kdPvo(I<&k{MO6$r+ttHrpVG()sz+2@
zC6L--?^F#n-nK0Lgxar?jjBmCtCs3a<pDc)*}GID3N!S5c2WC9c1T^US`T<3_loFn
zl<k+Ymr`Yybd^<mf?<1N`)>6m)U0gUX<O||B|EHDo?YcZ+qT=Q)TJF3u)d9SwYrSy
z_{VJ-``%94mF;_M8_gniL|HX~HlHqJfgNeztHvR}#tygF+BF@}cjj;N*mZW)4j-@%
za0Xni4yc3b3TV%j)qWX5OkEYMhyLtSSEIF-6x1~s2ayRp4S4qBX~ff1o%95orOJ3z
zX;EJeUhIC^7gurAt@3pZQkSU0;Fr)O8Sw-!#d8_bm*Y8r=OCUd@LY-KDm+)?`7)ks
z@LY@MIy{H)T#x54o;G#uPH*r?@Tl!eg3jxJ-XQx1JU8OGNqK@dxA~e{?eHp3n{Sr}
z4{E@|LttBGe7IesuCMX8W$Hz3nfr#@bq72ac~V1u*jXa9!|k`&HTJF!uWA$F9aBew
zw>B?Gl^#Xutx}$$l5U0DQbndmky~<nn<rXtp$p5?JSSD_hCb-KTHn|wpKZ(BwsxJm
zNeHY{H?vlUh|jJ|y4h`+jl;KlqTUNVF0?n;$J8zA*g=nde3j>bW#5KLL<?4^Zavg*
zXIAh;@b=&x!S-NB@Xp}L;9bGH?KWu19T3+&_6gdDE4<26<t_4bSVgccMV>xo`KdC@
z|NJ6U6!juIlR(qXFZ*USZaxH&8uM-Kx2WTwyNjjwsuT9T5QATT-#&~s_(mu9s@sb^
z!_cYaiB3-C-_bjtqC&k?ZnDU8AOo7su8s0!l2Tv$eR%Fy<Lr*T52$gV|H~a2$XS5#
zcjpcDba?F!j;teNz9+-8#Dn(ZzHs{_(Dy)4@`VwU^k~1;K4ISuYx<!5peli-=%+t`
z+}mBbrxdj7`$091HIn)h`-m|52pV~i70?CtQOP+9RvodA9`HgKhg64X!yVKHh~<!d
z1DM`HMpt|6%j}vZo~C*B?JziH!H4)=<U#pEBKe1ehv1YeX8RyZMBssa$yG2Y+g!Ra
zCm#+z(zL*S1ZZXe&7&>?sNZ8k*j;Ubeb?AmqPdtO`-&cjt^uBh?Mq3B{kTlYwZv+-
zAIJdSCxV^!wRR^6zt+A6t)37KyQ;l5Y2PKutC4w9Qg>H-YBFG*huB|6-rbVdVxL6Q
zmZse@!XQZN9z)aFAQD0Z$xDJ?u`frPuTZ4{F=AI0?E}d5wLj_Rx#dGKn^3yeVh3$m
z`<84YEWs$6QTa~Eg;L@>tNkgv4W8gmbrQye!l+YsVGiHy&S7P3@svK=SrI6-cE0xy
z{8d^EbnwSFgPxUm;>7?7PrHT3o5*NI)T6vAL;0$GMIPuf$N-wma?M8X3HSuKbRj-V
zeU3L;ZzkpzXIoB2eO=Z49zD|3=d2o;gO5-m7OJ;()hN`AdjULNJsQv?vTpL~!6rBR
zMsrX0CH}$4(}kuq=Zs7}l;!Jlip%s!<_chYC%@5-RPRZ|1!-N?Qsl*b03NHv(z82n
z&K+uvn~P69kDyQEKqwT6MiQaWD)YU($C3+B#mPS_9=CNO5{s_ZI;NeBNHn2|!SUC{
zlwBLw^ZAii8Ec3pLg9EkQXP$#X6VIIr7B#RU?XkE`*!A!H|&nviE~jeo(tgd4zvoa
zT#K+}=!TUWj`shmHP}4ce^DMs>twdrIu_S+ja6`NU?Fl!`{~7IOTl<+g=sDrTiAOn
zdL3$PF^?8}CCJvA-RMOGGOA1+j#k?`pBWrX7s1)_L^xrGq78L%$BOhNIlsF!q$2S|
zMd$Ma=KB^>E~Q?*tqXan*YTFk(A$}92_Mq(`6zoyEm~H_>gvyIFw=%yI*DYX^$xtq
zTUBZ~zS>A#id~ec-U+bZEE-yDmYJfd(!kQf0k@sr#aeyF&qj@Cw5miJq<Zky_8pX9
zZ(_yqJu2$ZXP|60Im4#rP^m@Hp`I&-g|kZdJp{l!G;G3>(=j~D22X7b<_f=EZY@I#
z#~-e*x1-9*sEb7P!p{7nftIz{j2S+AfX!lSuj7aIG}I>ac4H5pXYDe_ho9X=>qYFC
zs-{Kz^K@uW7+M%rAv+$g>dYGPsyBHx(d1OcYL%@*kb8~Z!v>t9S_l^bg6+E_&|FKO
zje_#hJk9QPgh1+xd_{veJ|}Ar<QS{g;W}N<I_sI)4?1068%sDomg;h;=!@@(HI@dT
zth+{!ol3QE{4!ghVj+DE8~2u-;}_iVlbsKae$JXFVzD%t1vM2#VL#L(#}`T1b$yzP
zL{FI$#cMEq-YG86Ny}rdSvqEzzk*#>njK?CTH8$1m~jX@IRMi!rew%QBr%NeemfLV
z`g(ko=792lA5*ga7>lwJkvcoXa#3)ZQjTBd1%?xMGOM7&pPK9aD`U|(=(VGj%^@8R
zKWo;F-C~t=-Z9pfIs8HPZUEkUfGLp+MWRd0^An4_Q$6N;_YX4fPaHk7n@y<+lbT~r
zKt<}btr{xR6HS=#u~{@}!-%oOH3a}ut6x>5Iv#1Uk+drHcys5Zr~Q*zKEZ@a>le`!
z=t=^2u<m>Z!TJ_tI{l!Cd&2d$lPND9W=-MlttM;oY)hM^lZPg$C{9kiGE(1MMg(#F
z80(3G_|boZw!z`bZ~}8dU&qS+rrPjsyVmhHyKi3;EP!Z_ZfDc;2t-I7Oi@PqFu^7Q
z8Z=!%K%3y0f1NznFH>T($)7T2STjq`VTl!X0`-;Q(0;ozuD@#5PARb-GqqDDo-3NZ
zkM$@9;i{7XN*#YRq+^Y7#~Z10`c>`@>rl8!c#LWCD67h>%W~NS)K2G5rtG);6InFQ
zY@fDZ>Nq4(`Aoi2*6Kc%c8gL_h_qIi=clb6DU-3E%T+PjA`{_hp*d@2VZKZR(Yw+t
zXy)>=h4WSrV-tYmjfe55KIol6E`@;FwGOI5&o_^jmFHdU5+-)wL6bdw@|q<?MvpHp
zL<gC=ir@wUQPH=VdWRrA3`vQ(o`Zswre*pFYfk4))3;kgMVUpdS*e>1#y_*Tf4bRD
zl<hQgXU<PjjQTbL>Xs9LT@rhtZ)Kh)s9=G=AsX4&0HcZz@Iq|UuUXJP2G6SZ!`#WZ
z<NBw}O?j^>4WA5mo)UB1-Unrcwd|^U-x;Xv3fZfPPQOneL;DF*YfSd6Y1U6V7tC5<
z_1}yP$0~F6rjRqW#H=%Ur*l%DCH7kL_UyUVI5Tk0&I?5N8LseuL)5u29EngJTtp;V
z-Q^~&zVa3IiG#xuKKVVRLq$m(Uo;G#bJ)B*r+Dd|M3)(fR>k6uPlcfg!shwNc6lvB
z73RuZ&xXJ1V##qcX72h8%Tbjo#}+ik`fmgu12_RP5Il|RkC-QYWkARDC;Su{rv<$P
zc@>?{&h=W;yGv4rUEj$1H=U-@PUN+i_vbB1QYYMSMrc2dPbbffyF#v85$$j}kQvrG
z-0Ti;8yj+frJ3#^wljlT<78FU!fn!{=+lX5IqEc$z1i^a_sm~nP3`>o{6W^(_JUTc
z)w6fN2w5+-dRq||^Q+2>P&k{RaF;?X;7DNy$n)+AfOKEpaF~~;1=+I5Sc@l^Cl?M*
zO31N_O@2ijmVOVf>RXTT@BU?4G-qHBvmru6Bn2nwuS7W4Q0EFw=G6zp(|a~)>YF}6
zAd10BB{o~1WGYE!m*%>RyMQ_812_de4U2e1nT5Y)h2IefkN=CQ4++R!YN4`SJL-;K
z<`x6iB!@l;36D3|nqXxk9QIj501M5#iw1Z8ebM7qk|T9j6)q3x%4$jr+|h{=3;i<A
z0PywjQZj5RjOY)nyF&P&o!o?uG$8U<kxGQ=G6jlQ=~)6lP^8A~<fEdQwseHG$egou
z^=KIz&2C>QK8d`WIv-j(-b#vL6_O{gBBhAXN~TW*GlWdxby}~5ewqH+Fk}*HV`j}@
z^`{6#a6%V?DrP)PK*l3HvFqx4geT*Av%W8(&!ElMOnlj@U4KBT96{PIi0WN}UlJ@K
zSWZCkXbw!XQ4d7aY}Y#e9;px|h-u;pSlL4Zfxytfc>Vj%!sYWUzsm}<@vLQICCP>q
zUSCacyj3%QU>-kfV_BMcZe3kOc}j=fy`uK4x0pS`OkXk9nrb$$SU#(}N>{*WH9J-n
zCf;yC(2YeDZzt$@+q}GDw)Kel>x$_Mp2l~Ni6r5WVqO7(Sdkp096wfnHuPsVJ`g0@
zCbY72eV-^I%MHj8+B+HXSZ#U$y=WmBP5~pwP<1S(LX~^M2*qDDKUg{2nqUH}V&i4@
zd0lBQAQ{~wbXmUydA~BpSCuS0odDC7B_*St2qydzkGPD0r6-#IUNvRqB;<fJ(ZyGK
z_tV|fDbRMJK}WIJt;E6wj_HFgJXpDxn(Eby^8}ACg~vR)dT0S7u%1+#AFduf^g}lE
zqr~IoT-Q_ay+*&(nX@J%GcX0&uyc=_A?vp0O=fe>Bjh|{8rPMAF&*nB`DrY$De%U+
zS^b|t((z*DxCB`AQd6)#GFK*D4{b{ez@<evpD#`pyAX&$BX-q1vus0&`TqI^pR={w
ziW?5Zvk<NA!cLMk3UoM`b>Y1+4Ly9sJhNevRbqa&;YSPSvRG`DsME*%WPGe&A`tFR
zW@BI6ww0d3@l>#6q4}_UN>Z$c7>;w9NrUBNCK~E%ZGAg4o+hBiI(e>nPQ?$dB@eUo
z1WU8RwY5TXT>pug=Mj91CElvI{x-i(@hhVmk!}CDjp3=VZZ@}_x7>Q!ynWvA!5yr#
zk$@{DHyr<q$=O(9Ej2SXjvPEpn#b)rCxBQ#rV&Bw924GHT(Ft-w{*AMYL09iKZ2U=
z<iW^>Nf-80V4QuG(QsT3G~e5}c9aPAEHWlD8jHkj$5$7Q*XT9m&ssBg)1))fnx2-d
z4#%+{#AWGUP3xvBtm`_xo2Ob`+eve_jF~7x5UvU9g*pDPuC7D0&mp)=B4t{)Y$~~z
z-{QuwUIt~od(;NBTx34nGS@Pl6Sw}uDiXhqOjmnT4*S?LEbf+Fuq7$;B-=GAy$O&l
zFnx$aOIx6N;t2zk=IiTObeNzG05PoW+;Af;u8%NhU?l=D8{y+_JDzYuIsGNHg|0N$
z+gPfq$l5AYzJWNpho#?OAIq56duUF7?$v3#8jC=Q^_snyybSg|8Ro)mv#nb@uiNIg
z;Hsw;FOsDyrYcqSbH#QKm`^U8wnZkmSGSV<{8UUC5`BonH&p<cyL*zJ0V1%ANUeBH
z^02hOX0WF3V)e_+m$uK{M1*e3CS+oL5Up(vRmF6u){a&u_UPMKPTsijio`H|MlX@L
zj@5>$YGaKjj4wFF>M%2cH(7_xE5T9Y(zVkms-Qj&)K7y`@K61A)U6-EJX$3uo9uLV
zD@#+BQl_)mhJT<2?68bOGWuXW0jys&FYNeTQV`+`;{^nlgd-Bn(L^{IEs=UXo8?1M
z?nV*9rEGV>M25IS7bc`RB^q&*D6+8IlF}dA*j{=O!1wJy{JfdB>xRvujzZB%c1DWy
zod6h0m@O$I3|2-Y;t_H9YZ@`xJ~fr@K*MjCkr&NM5{0YkL<_a*;49#==Q0)wGO}8e
zGeVJOd$+6J%yUwkoZx@_De7=8sPb3=Z@?Gu2QmvSeLq@OW-u3j0gpK|m}7o?QF4;d
zkw%tIbxZm+SJKH)Vh43u*o}Q!kx5#*qP~5^@lS;P1_{LjLR*T>W#OSFuVbKj&(1Rk
zEA}Pt=1^y_LnOE{+*=9sJuL6j;pdjWzL%{(C3)mg@ZSefCJ}C0{QzqT4gXZ)`aag}
z!#;}ZM`q^k*{-;W?td)ll>bYe3A{=kb_+WrkQwkn*bkzkZeg>u^XI#NZEc!LMtYs>
z+Hidw8-5XG%ZlgIQ}q_IonryDGGAp}-5T)@I@Kf4hk3ysW{!tOo2~Z5#hfFbM)i4S
z{k!HV`$<>e-BbTlv)&c>|7h0#Kc(*Vyk|bBy7nTWKolJVbXPa}AtYT^dDAP&Qh%RG
zFE-w_yMp92%J5tBR`r6UprbGZr!1f(4%Dlemy6D-&zHZqmch}iVdf%IuoS@Y!C5s!
zF3f4c5lu8vEb=VA_>%I{0)*stqVL=kanimx8DXV#ay>C_AP6wCU$|*cZ6s>Dv6O^9
zjI{J$2sQy=tteX1gxa5Nmd<B&8Jd_c8bs{d?e3*VR3RrXC#`_-pPFHjQtOG%6_I;$
zW~Ee)w&Zhp*o+J_qG9a<nhKcXu_T}?!e*M0h#t?>M8Sv<t1?!LC?U!nQ2mm*wPAj;
zjEyaQ7G93GF;@NBtNfd>90-{wOizV4@40&t`bIYKy1Uv)QtM0Di*zO7Is2oe!p*2h
zNjd(On8tjoMgjibT$V!lZL@Lz!hz=zp`|zBrKc0jFh};!v?iO^_D@(LqJNmVL=C=(
z*`hA!uMtcmC`FNDP3Ea2j_+hJ{<^TPL5RA@N^Q0eFxwmR=67o)DBFQnn!1WLYre9v
zTEqTgeGJEeeb1`4n-h&^jC+Jty}}Y#BdMz(Tqru9bQGGOHcqfEF#e{J0Of(g)h2ZF
zaI4&GXd1bG73w(swz$!jK8LNnI2*I85}L~k;9iTD9M_;_cY2NR7pFO~#2(Zk893W`
zvZ?0WtweENk1#6jC1gwCHdMvlV5BnP+K+8)Dq)qlu3rl#f|tzR=4naU7b;?#9PSmS
z)<5pX{NI?3xw4;g#R-I?(HO&1%~g#K63C)v2U9x<b`gY$ypPqQP2^y(c=n=>69C3|
zeYn!jcKqPIt?TtGrnqI;iW97ujrbRf(nJJnp<&EOg~d?inAVT66Pc-XOhpM|1h^lA
z6_aUc8Ee_*L`z9|H8Qv@oCxiY#3Q?FvD6?=oS$H$5XU&WiWRUXN`#sdk;)o5w}^|0
z(7Rcb&7vk2#j$79|Kd3tIyg89<;7}rFPKWFx&}i+m@TZ-3<U0iiV|~j<92<RCx;<K
zEU5t1v|U_0o+|(ZXA<{g8{3WHj}imdC$6IX*t~V|*b)i*16=Mg=+Q}ba4Es%1oYR~
zB_7c_YZP_V$qCoi?+GI|-f!2+v4W_*KF92Bos^d{JEz!Yl)0m|#A{`Ao@(`3{=YG6
zl6m8jBCDnIy-RM&%sZOMZN)Vq+=1%1>tN?~2O50VC&qjA+ky7d*2%n}VBWcUqLpvX
zyk=Bq>6d+$b)_l0W|Q>;bLllFajIQ%?L=$8S%2*m>toY&?fje{xHpm7@)}#sv&d~V
z?_L`p^`EVa)CgR$$*A$Q<!gRb;A@2YW%gfpw^d<A9@>?8LtFoSE#^yyiUzc_=Bp7y
zJ*~@H$1U+l?tO<Q_J6%~b?XY;2kNi?h`iU$e;g{ZUN?U|G^*e0xLMJE-`LKP*Z<wJ
zE;7U0W?NBH-Zsj*(A2h_vTiZ6j{L!DGZT(#ciEHLKam}}e{$w?ule@To!$#*b>!e6
znPRU5%KxP8<A_kZHvB(_LD!v>&sPXuf_kKkt5}^UnfVk!@@c?S(I1yp&z8TRFVoMk
zu!r^6UuC`jHS@2Z1v(h3MO)YpEp>{9wC!J}pJUnc1pkt*@PC%+uk-60f10%O@Jqi~
zra#j?U{$*Ii~7uW(=6T!8!Nt(`%l2N?u+8Gd}(o6`UMdEo((ww?IQED8)sOX%)px_
zs1+RF=W$^ACDep7{RSI|^L?2qI=~d~e9PTL#f8gz`Jox+_)T-itnR9<KVj;$&S|_E
zdh?)xbl32t;#$RtE;UPT9?{EV<&^_R<OAIIK~JmAl{ZhHPVc2(A#gpH{ua}p>A1ks
zB6DKi8DH>7=ti^dma|uhC{nX!J53)mTp)4yPJG^&X21S6Na=QsOk4Up=0~?|_S4?#
zHZ%R$@+3pD9=<{hw_;u&uOU8E{G)hEagr6RCyP9~RQE*TBBGXHlCw|0%FovbUMF~i
z;3Ny^Q1zQu67VpS^rZTGEdD-09|y{ns-}x{H~nPUFN?iCUNp=qbcp)9Z1Fr&aL-tD
z$8^HoWD;jz!9sD$Thm_oUzzCz;I$37(YTyu>mM>3CfTE9`~4Ln{}Dgwi1m+|8pzZj
z0`a_#dGXf4$=&?^cNW_OZxj54;6DiHc=b;K=-qkoN<6Kc!s!2G*=~ZL6GRB!A@~JB
z9`U}*R0mVNovhxBJSPiBk{5+|xg;(Z5A-3{e32m23owGImkCY*!0{p~OyK-wKUROG
z=k$&1PvfigFHPrdGn3boidF(@qvKEQGKt5%jM+bCA0IIF8-m{w{Epzi2xKEl!o&~p
z>GV(e?bK$G>%N^ldU6I`@yZx3R0N`_eKL4+(H$Z(>_#@TnwWZ@dAAYV%TI|EWNm*r
z^JK^7_e`B4kp239W9knCe<b)50WF38Gl8({E~IdsBT^^P$zPbC^H;+ECTKIiIeu{c
zpV6HY;PQ)mv2Jwo5!+=$ARz`f*tHi4^q#oR05`Z%j6`DfxHg`++dO~b6NC)U-Ttn(
z*lT762brhuxMA@b6kp0v^_h`k94vP!Cg{Z^97p$W!6aZ??#lL+UDv4^NwduS-h{}m
z(vQq<+dmw0H>*f@NZ^$BQ(G-Z?;Np{u3xVpc$$@RdI&*;D?$(jk<R*vLqtE0BkDuN
z?ZD`m#yhu(m4xw^)^|~uBH;VjW1lwMR{6R4`<>I0BpOz)gOzBjzTgHLjHL7=?8^bJ
zAAo)Fd?4~zx!yoVAS;kjXkkNzWtExC$6qcU>=ELAbLaghH+wCgdGr1`lLsS52n9Vm
zn3jOG59K>i@Fx?z$CgVX=$_Vz4b1EJj7W;8q=gW~%6}n2=-lJ)af=G=grC#UxV&_L
zyPjFfK8eop)!5B(SB58;silNKZ1r_a$+?Ia=W>2ll7Y!8!fqqewdyZL_wLR-r>@&F
zz`P88(i*sNw`c)nhN+WU;BhP%onVl`aa_CrQilm5v48x14sa5(J^jot@12Ni2>tJS
z$QoS24#pEqAecz-1^HZL${wih3O>8|jDGYaAnQd{Z<ddMoKq|>$K4RBmF{M)o=i$A
zn8!^vm!lp0<l?7`Wm8y5c1yO1ehCAF<-)kM(Y%v54iku}Iz7K|mY3p}Sd~0e{oona
zw9cy@40-=I_%zdOe`L@92cKrA`Q&be$gvQoK}tU(AMvs@Wh;9W<)s5p%S^Ft$v&Ez
zX6|RIYaeYb=<CbUdYk6FsO%B)J1qn8D2r#1%x)FdBEi#(0~a4%aiq3jCzy{OU2K(^
zd5@{h9E+1Jx0-SF2D>i|kzmx`AogYKS0*l<M_;3s3!Nz^OKDm>xF=F*4Yzcz`QWk2
z9m1kKHkY|vgbKLuGlUrf2wI3wW}?_k(dTZZc5-;JIkd;F#YV#U?7=TR6k3xkY&1Ra
z_y^0yc2Rszzl}PEt${k%@nrTm#k+?c4P-}*8)0b0!i(-6F+Y4_+)Uvp{il-=wHx~k
z&l$i2y3ni;Zw5n+qm_2hOzteT=5}uFe9ubCu30x3qEFR8g+Hu@O7QC150<c$9;nBD
z4ynS{X0%!2v|gkQ;Y<`uER@OJ{M0e5%dj>MYor8p-5k?dXjNz4cyhXx)%lMn&$HHz
zXNPGZgrebRQ6<V;!8XFMFl*+d7-o(=ea2uRtD9jdOKM(wdbr%Ob9pDWa)cW#RMs}&
zrme<>9_|jsuwjFnX;+&;&y2U0bS`{myERKl5(<S|nW^&;($JflG$Aj4b=K0pG_d|6
z^6UoT3qQNr;^bFj1E?yZp-qg`j+@!fj`uG|?aEAZ)^mkMKl{||G#%uR9tV&X<ip5I
zZ*5z$k=QcPbd(mF*w=<6E7&6Kqr^Rl>dOekQFXA2xG2RGjjp>zo#kYKrjQ&1yOD9A
z3O{?`caJl9nIqI$glTqF*H)(leIq?!c+bJ|be-dLaXw;(KR3o&WtKemv#ByhVRPzK
zY8c5_O|-uFKE(G~_;xZXRh2F>=R7}pz63sRF*CD1v_D+icXA9)G$QB?*WtQD2scAx
z*m=bH_Senr&llyDgTMFVHzVv~aOca<U*#SAxTFctW$K8c#8kXcjBTdNUKo`0QGoQI
zFvMKqw&bo%#i=+XD3%4QEM2QV*!0ZkE8c(n>(ysm@jj9tE)zptI;h8@ZX0_LA1ISf
z#392aUGDl#Bq74ui_UaSl2B~!U)i{i+{902>q}6{$?A%q<8FP52Q0#Q|4jEaOAG7H
zA|8&!$zcl|+-h9W>C44IDZ`v6KfQ_RI$D|Um^)saW}V;p+KYd-7K{0p)r#0HYBqx}
zVwZ^^s@|>d_+bY(j*0cmcg!c>C`(FvX-aW`yN|Fb5#qAmT&88MoaR{t^JGCQjFSa(
zdJ!$QOe}f`D0*CU38Pupk-vlycM(rZ1#2x}l5uv|nXhba3;8m#J|V?GwmCp>kf8Sn
zoxDbk>w>re!5?G9e`az{&BPAi(o;ioZe;D}2%a~gQ<JfEdCjTen@cLg(O5K6iOVc*
zo6w|5I~!Z4l};%Q)oYu#N6PA33d^R~0bj(u^d(pK;?tc0pK9CUEvq*#oGPpS1IM(A
z;6?NZk#w6{q>+<415V|$j{cYxeK<?lwlXwh-rSa~;F=B7X3U*6rz<&U#>|CI9`?CI
zjS+5dhh{D?#owHUvx)L=KIC;u%%*P@nu>2_c0I+RmQh<`b#`J8RI7Mj1Nc(hsAEW@
z{T@sl^U*h#%RTNcl@Yhno1t=RIST*U6|?UI@$mm!I0i9Rza;T_m^&?l%BEj@3I4;V
zSH-X*uS}HHb7Nd?qJ(&z*~zi%_cY5bZWxt+LCb!kV9v-c;s=L9p-SY-Uk0EL&q{M_
z-nhlxgnh1)v=eEXTPB#xUKy6;4{mC)w516IVcy2wErsVpsCC%2FBTVC9EpQ1m<6fH
zFd2Q9`q$vH8xn5Fx2806tE9F{<{DJyqNUHvyLAnCdL9jmTUuWFbcO8MI9^=s?q^5$
zM<5BD+cnhHH_OV#82~zHp=%?%%j%ou2ReTDEyl~i$DVjyoosD!K@#2_(@{+Wg(Ffd
zk6tvt`pyRHRWt9qW36Lm`*+Vx?qgYQI2x{P#_w7<?DS@j<JGpiVKba#8x7!OSKHAl
z9ji+;B{ai&TpWgbtcPVv#Ompv<L<_e>`qDaBzECVcEcaf$i@M8J1{CIYrpCLYJT!h
z#2!Ws{R&easmoa<ciKYk*-IIVrV-30;B_Va6N29pOknGY1f>L139chB1h*1&5QtTh
zF#8F9{*546Oc`*otEY-I2+9Z;o$46`a|swA=_Leo|9TaHIE3?=x{%-(BxVIuF`}*q
zaB^#7VP38jmB$VHhTVJZ$^;}#8;i53B!3rHZ8z0B{dmh;MYIDij;$u%`L82}0B&c-
zIU(@_CB%`Hdvr7UaWboQtf5}~z7DhFN&<0PA2H=42{gA(L46FDiQV%h{&Xa+f5vQg
zSH{%5mRIt3e&@<h1>Cuq`mIIkREj?V7;cWgHf#_#k-84!Ak4jS++}>_HQV~oOnd$C
z@_`t-GnD&W?m7W`Sz(SgrX#qhmM20@P4>7m(2c5k&N=hK>yh<N0qrQKT-W76_cuda
zOgs7fv1TY<Z&yaDA~@)6<VeL`q(k$Q9Js`|L(HBx<|l_yAcxpGtK00jo1L=l`7}|?
zVI7Hy&SYvC!C@AhNzV9a8lCL)r8#|!eTbhfk_rm5w1=Y+Lxk80!dsFDk{fF6Gtt5m
zzs1aabJl!ZQX1$D^w0O^dVN@tjKs!EF@D%u3^+L>+n=BGBt4mdqWrq-5&4<8-_kF;
zDlj5YVs3eJ(O@5zNd?w0T&T~+Kl$?Jntyw9iRu5n0SJTgtsy?2*N-lSW&12I7F44x
zKi<R84RWn~JO$aayeqsuT-HS`U>W8eVin?hxrKi@+02)Je9Pk>Xk~ku9&Ht%hiv@1
i)$KXp^=B5aRwk**9-69_jc$fmQ>}vY2Kq3X-2V%=o!k5X
delta 17278
zcmcJ031F1fwf}rGlVP$F7J-B<B!CPM5`tj~5yGwlVHI!$hRNJ-k|DE%`%NGT!9jsq
zQ8q7bhzR0_D-e9Ht*>>fd$q6iz3=&4{{6Row6<!iuhrK2Kfm)$X2Pid>ifTcV9wm{
zF6W+m&bjBDd%inI?sDus?kJd@pPyr+zYT?Vii_7CC>UhFz<p|rQ%mU1jpvI3(OVSm
zv3tkRlSlLsefMaZP4x4Q)yL`M^>WdFhbGPw1H?cvNEC@;F<1-{L&Y#LyfJr&eTU|q
zAZ(&Uj1Z-wOpFwx#AtE87$ZDl?BQH*h3-6BsaI^AcuezFQF)RWm%?NdiEdOSjk(m}
z0_spL`XAHI>C~f7(C(egea3Jfiy&K)ppwRHYM;WB3>V`?xtJg-h_uu{o?~jy_QY2^
z2#x-Y+vO;-d8dhoMWvXyQ}fpFpm9W&^)H1|*cj<_Q6(nrw0Ub6*u({*dZ$fP@34C>
z6qCgis+~c#Q^ho@okq1Y#Ys^krW4k>1qqv|6&F&&h14)h%n&oF`)pAsX3_T?aZ1b+
z^Qq%pu|QlzHS@&9;u89vFBXbL^u0hV7WMRfkys*@()Yz;nOIKWmx$A1g;+@=7m_&V
z(r1pCE&4Rt#HB&ESS40d%ObHxT$Z5Di$$04ip#00UYrqY#X2f25s!%V!bhd0VuSEg
zuVtb^1n7IY5JIQ#6(T4a>3gMU5+VA&RBRN@^u0>_lL(6lja)6FB1Sc9#G_)9kW_P-
zcuW`~P9?8?xoB;+?@zd8rPwUCh&ItKIvTTtjo7|cUn{myqeDKV7ad)vuNPZI=MFpR
zhc8-8ZJW96C&{*L-EHCuk!ZAgH|T!dceLS{jTAE=wl@;gg&x=emA+D3C3c9Nq!GIs
zU6?E3)y4O`L6UBxxSA$t;*y6-p+=i`BYm3b6Q)muK2iF_xJN1VXc6D{ZnB>AO+k<+
z`RN)e8KRFjPQ6=0Ht%NoY@u=+ecI{ML7%Pk>7-9$8-XiOpwD*tTuGm+=(9uYUTgR6
z^zLGm*AmL_QSH_A`96KF*{`*kdSI??zqanD-rc&`Vbg51%H3F{p3OYF?&!7pZhh@B
zyV%3LyH4!&?rEQzY~4q#d-UB2>flIr*qdzHZ#DI3*|*vjU9Qm-nyvkUWUm7m=(|T>
zpHYvoh1S^p&Mi+|W{cX_*es|o(i7JSB@XVi>HFrAy!5>!nG(|08}_)?W_$O04|uQl
zD(^w>4c;5QH+gT?Z{c-X>$@>yqC-0U6T6FSLv7we`ds}waZ?(zxLMy#c>7R~+2WQx
zZtt!8otvJ}Z%Hu175aYt0Byj-`e8AGl*=W@=t3M~W_B~R`}wda;T|>mPEj(&=G{qS
zY+SQQ-^De%NDw>qT|4X~3NMjyYnlvw?+%CeHsatB{WenT!$%V&BeG<-k%Z8l*9Xd0
zcT4d|QaG(d(Cw_lHgPnZAGzIoJ4^g_y%h#)mtEL`_99#AsRU0j&i~h5BNE@xjfba@
zjmYeE^r$E~rls0>>YkHjw10(|RC~gry~hNk0ms6}B6oQ2Aj<C`?VqH#5{K>}=^bZA
zmW$FI8nlJ%8vSc@h&%UWQQux0kL#@#VsV$aTimn9LAAL@?>vUR(QD&f;$CqAj*ZH9
zcb8|=!#(0YaX&Rd0KGQeD;`Kc9X*j8e4ltw+)N}q#7KJ`y&v^LTVeA&-u1t-`)lg$
z+#Bs@H5JdhA?xgy^ry1(x4DYQ<H|C1UH)*VX^*z3lldji6E-<cb>_}Z6l~E$jZJaG
z7mxW`BEj+0*KiWBX(KgRPTv7Lp6Jg(pI-EtLmy)dflh6kU1*(lVXw8tZAUfw*yD~)
zl6a>LaG>p&b_~Qh>gUU$G{EDOqv^|ZMs$DFbod(#Q}k`oWiTcq{%An=wZuZvxZ!b_
z4kE?O5qdCLW%ibOSl{fA#*-Ct93ip*NH%k6055^CO{5CLMZlKj)ZCi9K;1s6NDVCb
zd)Z{FlT!dPhP!>fP&5?x`D)eF-Y3UQqb_E_EW^-cJQR!0lQJev2hk>h!E{ApLJu2q
zswylz&r{1)L4P2QPEt2Ak-Lo621Ad}R`rE%YDH>9pPGDTLso|_V}_ipn)+Pgt|2U*
zUb0p_(Pxx4N4?jlwEwv@%0<*;p&HQlj2ENjLV)Q2Ol_9*M;monfC>=laY!aC$B6sm
zx-Z%qF-$F#Nm@aw)h9wmd`;KLep8&&pfDc0T+@x*V>j&+Dy0|I8lz2f{_~vLFpUhv
zA}uo_>fpfbV?Z{IUQ1u#neb?)GaQN}*#c?Fbp)d7*q{=%^qxGGQ=LDa2WgVIV2zdQ
z(KCYw7#O30b4T`^#U9RF7|+mxj-z^XlpDSQ<t8PI#^;uyRRbteCyM;p#Z)}|ckH5m
z#iMJ#j_KJt@l}_^VCK2Su{s)Ny8JCIdQ_N>NGK|6x>}3-YFe%O{@_Uz{ZO;E(R7if
zw1(r-t6mvASzE6hLoVv3)kj?{HnVQ$OJ9>e?(;{5PdAM9T`faiu_tx|O>Qt27P|0}
z?3?9AjAHUZxMU#`>l;E*gQRvb)zy3QC8n3n0NkHt%v!`MncYOPj5SJsL`w8oOu+0#
z`r8tY#Z4z#<r36oJB_B;7LS`0_Sxa3<DnX+D-dg??Q3d2c^w9yD?8K0xQ$!8`j$MU
zO=huZvW_~Es<N_R9_W(k48`?G#&DMC50tZXK5e0erS<dE^58w07xf8$+;3SB(-{wu
zTpB=yPQc6!#G(NbQ<OetHY>fMZ!+W>^=@gID~NuJmAh=@fMw|2%kXd3eIX(D(~~Ec
z);POt<f2>Alnul_wDWf04QRS7gr;VgSwYhHuULTIG`&G9)1!fQpEUfZ)#<VoT5*?a
zq%-S+2Qm8%1kT=vlBIwsL5yZ@V>seV+Gy+~e_P0qrD}$!YGR7N$=#ORpr%`dB2pKv
zf%GnA?tP+;c@_^I2-G785Vl^yP@@s*(5a{e<Y1LO_EA>}+J~s;$F>Y+a|!-QYweb?
zAo}GYwQOAf+yI(b4u@6CxJijPJ)2Ic`@^O~kiNJ!(}+MbZwoWm2!vYNCjf^bPh%u2
zmy1Yfl{)d9<U{*N?gbK8TiD;AhfP<z_3dmU3&NTw@4(Qd04%oSC_!rSK!3pH0I-O%
z5THg)9zV?mODI>U*6||-ZAI}C0;c9S-7NvXZ?hgS<df>5@gud<>SyE6FXSy^pdX|{
z8#JeU*pQ^Yo<a@pN=u7;gyGGux#gQR7d#wE=AS1_9bZaC+H*7Um9%Co+myr@1q-D{
zjjWhAkT+{Di<b+5hgrWyZLcWKufsF8s3)5=@t%rlL*apC2LaP=_~|2h?iKY##UQOh
z4XCWo-;F_b)+70#+Eh7q{)HGh7hn;<<p9?J+z9Xhz&ilH21rxOC9GOj6PT{vuN<sZ
zcYRTLnZ}e`I*U%J)v8EUN$>P9*zJ|-aMhFxYk=i804R>>CiBE*LmoyQG>EsY>1vIJ
zHnr*!Iz|h}`sD|zaMF2;-lw9ObX-A6ZkIeXtG?-=JxzK>M$%{PNSWd~5%vL4{2qW;
z@-CHTs7;f~wRgL2pEOk~SOL7+1WV{;>WZ4Guk%VNj2W9ts5#Zu+91_ly>=B#(_u;T
zLtxFL{U7&5$R~uNjY%Gv+BsqaUEoqP?SBIym&@iE$mY>>M*U>W4ybXHOD?<<=(0l5
zV9YR`FnH`T<x~u3uJS6Gjfc=Rha6NlGmfb5PhPZm7F9!&Frp6N3xK~9Fx_Asanq20
zLmkiM5Uu12h^QXw!2G+Ax~%CMF~zP`q!cA-$>mKLeltL-DXsBREuUJOV6(9hU6}SG
zR5J7ZG<{>#XZaD!eAl5N%P*zB-C7`fEG*EIWi7%>Wy=m{2gBqi;X=}a?;$3vJO=jW
z*la3)TCG;ob;-1T+IW$3w6HVln1)kIywB0Tom!`DW52-^y=oN4b~>!5pq(PwT(z~P
z`n#w_hkCnaL81&K&I9-@0kRCgqp!w8x8U!@Om@{BePa|P%FX1>E!{2Iy3NzShecp<
zQXobC(G{R*t~FwbkEQ+x;68T;NobzD9U#G*4m*jxp=CH>edS#M5P&DoV&-a8tRi3*
zG8rI7#w~}(d+;N4_!xls{zsJl1b~eJ%@lfDrqpv5QQH=9eyqx`^=JiM)wQR!`Pe2B
ztVuEiUZ~kqXf2ENnCa#P=0F2^1*|x{Ug(TmsLE#)4`hX}!1Gf8Zt9H{sr56;wCU=q
z8S{oQwH;{8P{~KBZeQ20W{lDj#b{=#hM>fx<y+)%D8ZDZiD2zz5_yy7kN`P~ru~EC
zk$en*g~Ieutb~e(0l-R%CG<#3k65DYTxMnx`UH(RqaLW6yY5pe)l=;C8=(46fU^MD
zk#ZIQ_${$8i9tOO@!EO&enJhaArjgQHn)~D(Cu~)a*vXqbS<7WS#w#eP}j_!S;|EW
zgnTk(YubYoUsb=Iy|ls_(LEiBbTCAAF0-k3)oZ98s;Q4=4|i0Nzr1!%T~(@AcPOZ}
zn_7?z<q&ipmO{joucJABuD+NvNjs?~&aIgGIOa^76OLDcs4E0ubCF9W(?w@CogxfD
z4;0ijs+;F}7WHuR=#y=Q+7z^tGl_-MYdD7%kYN@g((^UOV!{__@>7ugoXVd!NgJbP
z&5MoVEoHX^K8<>3@PS43zWVvRk<-3ThADPf!+6ymV#0g$Vc$T9T92xlKW@%gsv**N
z^)H?~4<oTCm3q8YMty!GKt~etK)>q7W3y8h;OPA67cd^{bS{C1PA7ESh(#zYYibP!
z!xSNM(A0Ass`r9n1AhxtE*fF8=V*DFoIua><+ELL7dW!qBdD4NcBv(cR_2$`L-J4}
zPpXrPJlYUtE*j$+hWaw)UR>4t3>8f~ofKyf7CBR0yf{?NJF$lfrbUuwj#<F|l`<F<
zGlukFg=(uGsdDS5e#_C;D5MX7r;0}F{T~5Y!-)<vE8^c6lcbp^)yMT?v|?4f<gIBO
zpM8KntYC1T9M8&UQDS+F$6%+S>q3@eJ#0d1)Yzrt!gc7x)?^_{uwG_XytO5)%cH1Z
zm4*(R`IZ$<h7Oj!W(x$vZf5(#VW!!Tf5y0_0574%9yH|3czQwoXz7Z|=TYY{qK4T^
z=mte0y6Fn*jie6pIC^may&olxpFdM`myOhBsL-+@=S{GN>JigT(RoZZkT1?xhnAHT
zB4Cm$0eIvp_0qCY+BEgivLda&%2~c(D2s6p^Dr6<89F6?{6@392;8_>MVF7sO>6pi
zb!>S*?PK-G@?F~AuIVerYu$OP$jTAtGn3OBofq%&%f^VmO@hFrUsI2)TsHCqp4l7G
z+1H<W7Lzy7?5ow-ORKdTx;9+;rB(zppFVXZz49$HWwB~oy&}OoBgc~Rxdg}yFmN9h
zleRe$br6uJK<-D=0f6fXP+-cL805PQsZi55K(Uohkzs>wkdH;P+)iUj7us8NI!ubt
zf*^Ii5jaxIl5b#+S*Yv5J`C%PB$(72DYIaOaTG@ypwebn2i8o|4tG7b#-(XxYWOWh
z>g0_h_ykECw2pdwtI+KZHN{)Lq6ppX@(L`WAQ>|H;xgn98|0K&&3StOZ93GK8$fXl
z>cB@&Z`qEDyHT-2J?gDq287nZr-LfZ-mz%7-AB1WUs#Ve#+&4AXy1<hUH}e9GC072
zwv2^+!EkKL4eG+n$GMKvsDMLN&g`$Qx_p~<;QF;g%hIFLCAP9hsPn_bH91t(t&J?2
zh<RXv%p9=RI)*2+m{tOg7fNv=7{cq}iVUojBe$JWpRN5c!3gPz){QZ4%TjQp26JFK
zIzla+b(KZbZiNxdtsHB)V2s5AW=5K0tQ`k|BAc=(NnJXF-aTn?)ALb6{D_LJKe(J#
zlPS6jlYt_6CjnXtn;6K77RxJQF|k-&&0A>WePvr6qk&JVOMO)dps>^)Xi2>w0a_}v
zX9pS?8CWgB8Ky{Lyfsx%g2DOZwuuCdhoH}SM3qf*+ucsL%bnF<llRlefCIJk=eDV7
zGjmm+4T&*KM;ckWl#=vwmZTrs5MKhFW_D8!b~MTh0eI<kC}oJ_9|-#ewbkF1xCcw(
zNRSFh?lR((u;smI&p`aG<t$HN^j9QIbrA4Nv@-QzL)*Lu(2HsKhmw`|p>GD8AePtF
z#=s;?tU&kwEYYM-OD46RJx8bxcb3~pLVb`XN(mJr?;0Y0tJRl+k#;jD>~AqB4HptR
zA6CA)B@aPpSPN0hipT3pN&an1@+*TIregzqb^N#4((kAbgO6BJPVMWjY-vl%<sj$(
zi>U3k_tmJTYu7Ucti{AR>+B*QqN2qgyJdLUF}X5!G#hfuxH6Kj(STp8g3#3Yj1CfD
zoP0d8PD}GqS3r|BUdn$zpGn?i0V-=j!At_Cle|(p38PkB7aH8-^v_$2r3EuVFNK7#
z0;zh8&p}-;e_K<SLKZ9fS%&%`z-Iu6SLmq6GVZ`^ccb>}jLRxC@EUfELr2U^$+L!Y
z)C3x0=BCxlcuOtY=+Vw}ZQXcJZdFpD>ZsGV@|S<LtCgD<OobsJqdJg=B{Q$G3uVL>
zaxh9m8Kap+xv((#{iuv>G32xAXPc)4Dlu4-bLh*oZ;3U&`U?Ik>j)w0n7hfK<*b!T
zkT+q7-HK_b2=y+;EIbtl)R+$<@wZdIS{2!Hse2g)X_O`Xmn{=>Yk+w=K#6j-jh@AF
zd>GFdPHShw1HsB408{`B2U2aUhu$*M%P}U{74gev3L9&c|6qT0W!td6Dea<(*Mg{t
z>Qvi+j8x$=^^3OQWhXI}ofl{+5>$*x<hgNnrDH?+&7<f=m73W;(hc@f>jKr>UL0CV
z%^lgJN&+#8umW*2cd3ju_!~muP`q8vr$)2a3M)#I3xU;6@2K=(Tq5kJIWGcb7h4EZ
z<@r`d-q{d~H&GvuC+pQ8+M6#~1r$qqgg$9Ye;B>kw>w2B5Vy<$hqE<kSqSg}0Ow<#
zRmVHZhjLb?f-22ir0<gF5_tgKzEHpK7?{A>Laqhiqs2Ou)&uwe*n4M~92hiu85E;K
ztrJBz?I@$gAE1K(z1P%bi*@eH3x?=O(V)k-=z7$2#-j8(lcVh%id*UU7!T1IWkMy>
zW;tfkSmN<%%*z`R88+Dhun7QKD1)kNYpK?#ezA3=7E*uRT2kIb)n;#cE9z^HhXT!<
z<u`b{OM#{wG_|40kj?6%&I=YF$LcxIW(%+%kG%QETbrNv{Yuns#gyC*xtMu|-r~nQ
zI=tUbhK}mR&XQ4x&yz2Me3Uu1BJ6{}fYXDeW~M6GRyuM!J)3SEwg`nGufzmb0qg{T
zXQlM)%57Ccp`m83Kityfr>J$a9_BQ?WV?EaI=*cT>pyAf3mC)3X1MZg9&Z24-u3rw
zPOan%)Wa?Wa%oQdK@0QgGlxD#TUTk~#;pAP@w`q-SJK6sUhA&vdTwW{Q~R5me$Csi
z&aqQ$+Xk!QyU*A1)V|$y>gwI^Ydclw+GX0C>d9-5Yjw(Z-T7KvU47j+?Q_+2-IUxn
zlec2>w{)maskTFn*kcTRud_xBF0u(nvvYre#M1_P8f-nOQ+w{BH<qjSuB*6yfA39g
zlpv*^MYfHDJLjU;KwIaGX50SWvaPd#pJsNJ)zb5(HucNB=P!Pxb6)4{W*4>pk?LOQ
ztm`c7oYzbj27sf-=vO*tbrxjN^(ZlT3f(@awRO%u>gk+^v@FK01}y1Xz3&T6TdOWT
zFi8umYYq(6R;haqJg42LBG>;(+pGLaTF0DZmV{Gr|KOd>?aFy@t+s}2R#hTt^08UF
zEf#A>%8sQI+b0|(Aqp&zy8jl|UOoa@{}TX`26*elf219cqVh2S`{SDWcys@Rd;(8B
ze1UwDD*or3fqaV45!GzyVbYW2{B$P!mPxnV^5RVSL(uRvz(1#Zyf;(+2v5&^I%e&@
zAOC8m{6<HTJasNevb2s?`>nj_J88bIB<8R;)hHX2#JBLf;`=O(kI?;?h(8OvLgiXD
z<fck(g<5>mXfYS-{5CF6K1)5x#Xg4t?2?~H37kwa(OG8|Hms>57PHr6K95mv-Bdkd
zL3eLyLsa^;PHt*795`_+B_Br9ao$p+J<t@3=yd!J<OExz*4wb2t}q|2F!TXkW1zVg
zs2|;2Q33xaUj(pxpL_}BZ*+ac=&9fF9nT*1y<0At%YsSnH80W6k(LEG5?@Dkdz`6F
zBVQ&mh)Bz^qO8eRRNq6(U9jWwdbROTU4rAqw4+I9UNT$`d*Fqbl>IEb*ake-0l*<!
zFB>^4VfkLvy@r`y2ly$#8vrSHEZ@`;^lbTI^!^!|{{<k!A>WTWxLhlh$7eo1?qzu8
z0I+ZZ>s54qnjR_KrQ>x>#;w#6hYAg7WGA0Ud-<QE(n$apFUSWQ{JvTKE2@71@D{*}
znBm`0f*Y46N_|o42f%LoZK@+*%-Ja<o_zGm&)_t<{10?*1b7GFU0(ltD24#uCxG|?
z_oO6TI=^|5{~IL%fK0P&2KY6=Zvf!T<bR@c9Hn!;xkOUY%q9waxGrMwbpW{!9iIWX
zs2BDK`_t#~`~te?85Cy5>D|j_IuLp;&O8i!O;8}usy`f_xM4r2xB>v`iIgZOn%K4P
zK(&cE{(#bl03QK-4Dd&Q?P%p_^iT9?_D%)}-Pbkbl`3)DQ0w4)0R3-3e@<I<p!9px
z-3E}#-dLylov7oS%qJ+l0B{LL{5wjY0{j=ip8;Sm<X-@oD|b_gvi6|}M?0Ts6g$Yz
z(fU__>(!VeJFEXfxJ);4-O}IvuJmshmjk{dSfU&}Ga5TMxnb4&M+RuOslOcgJB1d1
zzWvu)se0$=z0*fP9>vs+j8n#*OWI`)%oq`G`Xe;atF-AU#rC6PbGolQ@%fGSc8@7B
z`R&yGb>+F^;}Q2@0X!WQY@VF2q*K)N<Ac}2zstD*kE2&^4}Nb4KIS)zoM+46JPUIa
zmXtvTk}7Yi)5lk{AtY<hqj_m8%lj~6#+YMN_wQ=zo#i1=O@{3_I@N%z@0JZCJIFIJ
z|MCgfD=?Y`gvh4l+1(Czw%gHPqpV4O*K2n++3ha7O4N*0mG}1^TZA=(E}EW?Ti+DC
z{9fH~kIu(1nho{_1620CgA*(oURug%f7rulQ`gjc9cEtz)zo|XSax>?T7k}VHtX$%
zB{Q}}ac%)%>$?}Fr2uShm!PE6IGJF(lCl!rXZ>*y67OH^J!@nR>ay_Y#c<xn+fjmu
zXj3E)d`WgN23{Q2g*zG%Qhh!p7{9~<kv6mSQl%%(*Cupba^eAPz-UZxKEN0N55Rwi
z*?m;#{>JXmvYXj7nH)<PiScRD&gFZ*tdz5?r#l=_C7I%|91n6DP>1Z9#a1@4I16?2
ztsI>?n6cKz#541qxyL+X(|0HOu-Rfu^z}@m7jH?XvC+r^eFJDVTXuog-(FxZw)dqQ
zxOT0t=E6<f3Ysk7LOcC+z4pKg`~MhAJ?h}YP5%v+j!Ux?N7ZEJ4BLkG0rWZaIfYJ=
z-%?td-S=zxJ{n^a8JeAz!#ApThISWbI<K^zr-j5St{|VZxY)Oh=1lK)kOGyLEbOmJ
zPL-VMuMJbxr>1L_YSSsP9K@SBd>xst&rre!Q4txrH-KaoCgPn9zmaL=3P_~$!9zDS
z+O!-pL&cgbQ0JcxT+Te{jp4l6k<p+_JA+Ws51<qHc&jqCtm`R-H*;|#+1I3pDf6%x
zGqB<i{xO=YQolO=K^+^nGD2eZLV71mmqv@NIXG&91F-Tf9PsF9a$CNZ^;)SuXUdA1
zg>apwBdTxdxja>WW`MR_>1RCJl&*be&T0uxprvMIeMyG+bPD3BlW%kuV(Rq#Ob&>k
z{ARSBBe8QR8H7#9mYE6R%j|p!ic)9HpX<S{aU8Ox7$kFPl&*3ft<ZA2Y93vZB#&uK
z6G{V9@0l{FVr6XrrOclIdgiA2W9H$Vm)_=&(0gg4`}HM-40OMR2~Ke?X>&(mfqaJy
zu}Ald8$^`naJBsL5^cMxe|(fSqbvUSWm*+8lqpV1?_R8ymFHX%Q`D(1pQySpb9yYK
z7L2@+nT5v`ugro*$|nUwk~9h_oP#R*<S1>nI{4)46VlUzS8@=6v``<Wx^z~Xu5_vc
zPgNwCX=%507sjROW#7eym%UaHtuVn>x+>eu24DFNuNCdOgY?S-X4A@4SC%?!X%*_A
z7HK99urfylxRj0MSjg^HYkx38o2Rb$!TaNRP0Y>YTPHqaETF~Ejf71UtY;HK(-9Cs
zS)s1|;jk$j^Sp)1tQOy9e>gpvHP?CZEit;P1Xb|SmCYDsCFz~@qw2jM7UeG?rk}7|
zuL6p?Tu)zZPjGshiC~^$kMjOVF*fy@y~DM&X2PW4Hu8GE!1T_z-WE@`($zQ@1PZrH
z));5sj=VX^u$At|mGGIBPp^$1EPHI&uHT({x^eif-z`f#Htge>oRs$TOD0x90dK<^
zlpY6Qi7{ia0_z;d!eV)bX*L-FV8i=4Mr249LZvt~L7mL(?r_?m{;Q}jBuLD}?@_F;
zD#=Wq4MbQ9UrNVoe@RzzGWp#v<yg}{kNL5B|HtLps;=J8{FgSJeE>UZ))HtsA~qoB
zjlrs)(trAu2C~omD~-2R-E$KXJf4R_OvK%j7z%`R1CN8-4^WY2Spy#VT+4jpb8{W)
z_$<8;CARW@u<R%vnaONblJELC?6@0EojiORiipL#Z3Sn@nf$K{LHYypE`aX=oV!9Z
ze~YA_4bb%m{HntEH?{criP~i<@%+Hto6!4*08gu9&yUsSs#l&b9^$IM?9%ee$@Csz
z^&DSi)g&`#!J@^!6?AZ#I7Ll;VM5iXpl%+(Gc=Z%o3gA_rW)(gNOwK@s8?Sor?lZ;
zUwA-$@RPa)(6nxud{fxH#FdKZ`1&uMqS#cN1-ZRREjy~zjF$!#;}=J`>682wK}a`z
zTes2w2`x0$RIl_?R}$j4nQ~2f$`H~P9+%|&3QpcVuJd$)f^C|5=tYnE^d}XG?=nSG
zdWw$Z+g^y4=%rU!4IQNs5>d#yHpllLjSBiLn00A^lMfSGOxox4=7#R`;yyCnJS=v%
z=zeL?Rd0I8h2<krku_?gScBdzCFT(IMt=g0^y{;AdE4mMXq{S-?SOrgbA_cV#2z@%
zZi!%jX6O)0;L9}a6ak$<FG$49qn~_`ZL(4=-A~E)77rsYayHMxj=gRD00wzBniXcZ
zBfG357A%Q`0`zksh&7*|yDwaK^zXh$gBL!?d`xKJo8pj+nafTWhm&ZW@Fh2F38r(+
zr`vW^xSXwB{d`Z>N-nMB+WFKQL2JgBedT=O!crn5*;?;;fYVE+ovxMl(xaP0G8U!R
zkgbuHb~##o^wRLr@OGxFrM=0BMELa%&Qtz|n2btTO1fM|r@LoV^~;O3*VK-eOSN0o
z-7jCF8S3MghfK|(n`S{7i^SXF62}a>Tukn`g$!Rj*1{6D(nYTlrXFv|-c)N|8I|}m
z&1~hB0%+!hmK{^}#`P|~{^qk@6d?eSBLT_*DgiD4s0Mf!;1hr`z~ccJ2QUF(F97a_
z$-@B00od9lIEX)k=Pv+C!8RPi<pgFRKqUYoVOa$*1pv{RoB=Q!U>*SbnN=uV2JjmY
zGZ&>zzzm~e=7nQ^+|i`pBL-t;UNkn}x?@VeGfBQ5!!eJ-tk#B&^kUPnUWxq;m|Q>t
zLu?kd#OPu%C-ewP>h(+Y(W?WOe2K<9{z*Xh6|whGm2LIKRsSKt-WNZ0j2XU$b`C9k
zxNMHJQ;#5UqamlORj=vVpVU*Y?Q{00Vd~n~Lo+#~vEI+)Y;6{><uY3IKqwfZaJC(k
zvLoQs3wxluKsblEgH+v5rzDC&gL6}BC@$Fwaqep}RmwHoB1kW#j^{Wph_YR<1U^^=
z8ASa)%)-w5C6rzUcmd#ffOKrewl?iZxHKMC!gSXX;jLkP2GM3Syvp%LRkfqQo@aN`
z`C|yBT}tQ&vLyt^I&z#PnmeZ;VN2Mu+(iY^oWTWIbd9A~PS8ErU8cfs)bw)G0GHjZ
zo_wQLeep(aw7^B@Vf!FLW~bxQFghLEwLyf|r4{B>YA)9xx=UNCIq5lv?$x^aAHMnH
o7(jjVXtF{(;nhatpPfjtmpJn5<Lv{eUmpD{pnq<yP;;;UKi5qu_W%F@
diff --git a/__pycache__/simplex_plot.cpython-37.pyc b/__pycache__/simplex_plot.cpython-37.pyc
index a8537b6f9da6be3317317d69450bbe8640da3809..d1a178aa1d78d9cd467d1b86bd96b0a477756845 100644
GIT binary patch
delta 803
zcmY*XJ8#rL5RQGZeUBV)$teUPh(btV1PLTah=%YIP=qKFDbN<I_Iw*>!}`H^FS3u~
zB>sY_k*H7_BorE|@CWz-Y=M;WDntb(V_y)k@_rue%zQKR?Y^7&G}FA>Y}UtceIM+1
zzkF!A=G6Ssa?%R|1{vKJeh+xEd~n13Y#Jp0HU4x!Y;-;601jqB4Z(lee*f^O|MAd&
zbks-l#~P_b+uL9kzR}UwHGj2S2W!I_Ja!E|o`^mtmJ<Gs5|@p!u`Rbk8l6s*QW1JP
zJRy~CEJ#BE3zL`u50yKsjlm|dM{Q|io!ILl<SX~J8OBKX0u-tE9wSxO?NW&@Vx=q<
z327x<(SYqnKxO1H6@fM%Z;{4EPs;Vt0kN_oRG)|;+U1^lB7#`*R98?8;v|@)IFRuN
z0cRm`P?a4Pb8W>SLPe?te=0>VvNE)wkQF-#!$>j3$wWCRl`M;TisY~Do0oJwVSAh^
zp(KaX7_KMwI4GWak&GbG)s&)+NUigq_JXmH&pNM+#e+A_p<&qhv)c3(3zLR60#p$;
z+$yvR8mrGCO0;#{BP!@~n9SeS+SkvZpV$gStWaUmO|^kvM$Gsc?iJjtuuu?|aIvu3
z`LO01?R=`ff1xPjn%eNne%(3kP}+?=O0^aCf;fX#{;fVYf3m333Ydyzq@?I9DdFLp
z!NTPqrt<TRbJnRMD!<*B9{d+Onx5FlDj4_3<UdR?Q|0S${ut)w9at>;M9@`ALl$sK
si7i=&%Y;}AfMukvc%tG6Zz}<}h^$B+!qp=2Y`eg}Q;$Z^hmBjm0sEc<E&u=k
delta 804
zcmY*X%}x|S5Z<0)W|mnt24e^T4Jz?4MiKA;m?!~E&;$r39CR;}-q~hpoj<F4$t)KV
zUtv!802(ev!^wESv(X2zZXUUD_hR+#8nH56U8(x2y6T(yJpO*XzFM#Q6_RhT+WPXT
z9$4(s$y4jA#jy09=ZEMf99{{4Dq+AN+i#DI9G|v-d^d8wjCQuW3zmsJ`E_vR!}`zV
z0qAJ?Y@FNy!yzRP1Ug!__q4pjDwUnU##*zP6`@R{J)zKU<q~Tec&c&^BGthpvKQOv
z?#gxt1hNfuH)Ja2AETA#RPZ%uwA&HKDsQ#Im_q0jypRe<mC#|rUuF=-S;Tb)Xl?Fb
zZL1T<8v{rAz@2=Co}#c8fJ2GkDS9B(9?z*z4l>nPHHC^;CVi_*7Sw1UdNNHj&9%Uh
z-bi80yIDtL`OaOs-S8FP7onCq7BFMrBl`lhD55ORKsBmGNQ6=8Mmg+FvxV}m_l7N<
zeDqEk^U8d9?9ip8GoxUb)K*Z12~);75RfT!20ap0Z2^-eVn^h+h(KG(l;E@EY`7X^
zu8CRMDv;6FCCD6|`G%F>hXum))uOScNv+UPS1YCCPqFFpmVbQHydPM|qTV3E6nbIU
z%AzntC+#G87bePcf9l2<!EiL;I*&6Q%jS%h5p4~+{Tr^Av$bo7W*4x20Tm@^MM=Y2
z+x#!N&l9?Xo2}vH3rRWGy@m1q8rJp-EK%v+2v`lnlqZ5ZbYtEWG4aI#c(>u-(?aPy
cqaUl_p=rgAMGEUCbIz!Lg!HV%@^|ggAEibUf&c&j
diff --git a/analysis_pmf.py b/analysis_pmf.py
index 03e55c1d..8be3eecf 100644
--- a/analysis_pmf.py
+++ b/analysis_pmf.py
@@ -50,32 +50,32 @@ if __name__ == "__main__":
state_types_interpolation = state_types_interpolation / state_types_interpolation.max() * 100
fs = 18
- # plt.plot(np.linspace(0, 1, 150), state_types_interpolation[0], color=type2color[0])
- # plt.ylabel("% of type across population", size=fs)
- # plt.xlabel("Interpolated session time", size=fs)
- # plt.ylim(0, 100)
- # sns.despine()
- # plt.tight_layout()
- # plt.savefig("type hist 1")
- # plt.show()
- #
- # plt.plot(np.linspace(0, 1, 150), state_types_interpolation[1], color=type2color[1])
- # plt.ylabel("% of type across population", size=fs)
- # plt.xlabel("Interpolated session time", size=fs)
- # plt.ylim(0, 100)
- # sns.despine()
- # plt.tight_layout()
- # plt.savefig("type hist 2")
- # plt.show()
- #
- # plt.plot(np.linspace(0, 1, 150), state_types_interpolation[2], color=type2color[2])
- # plt.ylabel("% of type across population", size=fs)
- # plt.xlabel("Interpolated session time", size=fs)
- # plt.ylim(0, 100)
- # sns.despine()
- # plt.tight_layout()
- # plt.savefig("type hist 3")
- # plt.show()
+ plt.plot(np.linspace(0, 1, 150), state_types_interpolation[0], color=type2color[0])
+ plt.ylabel("% of type across population", size=fs)
+ plt.xlabel("Interpolated session time", size=fs)
+ plt.ylim(0, 100)
+ sns.despine()
+ plt.tight_layout()
+ plt.savefig("./summary_figures/type hist 1")
+ plt.close()
+
+ plt.plot(np.linspace(0, 1, 150), state_types_interpolation[1], color=type2color[1])
+ plt.ylabel("% of type across population", size=fs)
+ plt.xlabel("Interpolated session time", size=fs)
+ plt.ylim(0, 100)
+ sns.despine()
+ plt.tight_layout()
+ plt.savefig("./summary_figures/type hist 2")
+ plt.close()
+
+ plt.plot(np.linspace(0, 1, 150), state_types_interpolation[2], color=type2color[2])
+ plt.ylabel("% of type across population", size=fs)
+ plt.xlabel("Interpolated session time", size=fs)
+ plt.ylim(0, 100)
+ sns.despine()
+ plt.tight_layout()
+ plt.savefig("./summary_figures/type hist 3")
+ plt.close()
all_first_pmfs_typeless = pickle.load(open("all_first_pmfs_typeless.p", 'rb'))
all_pmfs = pickle.load(open("all_pmfs.p", 'rb'))
@@ -111,27 +111,27 @@ if __name__ == "__main__":
# plt.xlabel("Bias flips")
# sns.despine()
# plt.tight_layout()
- # plt.savefig("./meeting_figures/bias_flips.png")
+ # plt.savefig("./summary_figures/bias_flips.png")
# plt.show()
- # fewer_states_side = []
- # for key in all_first_pmfs_typeless:
- # animal_biases = np.zeros(2)
- # for defined_points, pmf in all_first_pmfs_typeless[key]:
- # bias = np.mean(pmf[defined_points])
- # if bias > 0.55:
- # animal_biases[0] += 1
- # elif bias < 0.45:
- # animal_biases[1] += 1
- # fewer_states_side.append(np.min(animal_biases / animal_biases.sum()))
- # plt.hist(fewer_states_side)
- # plt.title("Mixed biases")
- # plt.ylabel("# of mice")
- # plt.xlabel("min(% left biased states, % right biased states)")
- # sns.despine()
- # plt.tight_layout()
- # plt.savefig("./meeting_figures/proportion_other_bias")
- # plt.show()
+ fewer_states_side = []
+ for key in all_first_pmfs_typeless:
+ animal_biases = np.zeros(2)
+ for defined_points, pmf in all_first_pmfs_typeless[key]:
+ bias = np.mean(pmf[defined_points])
+ if bias > 0.55:
+ animal_biases[0] += 1
+ elif bias < 0.45:
+ animal_biases[1] += 1
+ fewer_states_side.append(np.min(animal_biases / animal_biases.sum()))
+ plt.hist(fewer_states_side)
+ plt.title("Mixed biases")
+ plt.ylabel("# of mice")
+ plt.xlabel("min(% left biased states, % right biased states)")
+ sns.despine()
+ plt.tight_layout()
+ plt.savefig("./summary_figures/proportion_other_bias")
+ plt.close()
# total_counter = 0
# bias_counter = 0
@@ -287,7 +287,7 @@ if __name__ == "__main__":
plt.gca().spines['bottom'].set_linewidth(4)
plt.tight_layout()
plt.savefig("single exam 5")
- plt.show()
+ plt.close()
@@ -301,14 +301,14 @@ if __name__ == "__main__":
plt.ylim(0, 1)
plt.xlim(0, 10)
- plt.ylabel("P(rightwards)", size=32)
- plt.xlabel("Contrast", size=32)
+ # plt.ylabel("P(rightwards)", size=32)
+ # plt.xlabel("Contrast", size=32)
plt.yticks([0, 1], size=27)
plt.gca().set_xticks([0, 5, 10], [-1, 0, 1], size=27)
sns.despine()
plt.tight_layout()
- plt.savefig("example type 1")
- plt.show()
+ plt.savefig("example type 1", transparent=True)
+ plt.close()
state_num = 1
defined_points, pmf = all_first_pmfs_typeless['CSHL_018'][state_num][0], all_first_pmfs_typeless['CSHL_018'][state_num][1]
@@ -319,14 +319,14 @@ if __name__ == "__main__":
plt.ylim(0, 1)
plt.xlim(0, 10)
- plt.ylabel("P(rightwards)", size=32)
- plt.xlabel("Contrast", size=32)
+ # plt.ylabel("P(rightwards)", size=32)
+ # plt.xlabel("Contrast", size=32)
plt.yticks([0, 1], size=27)
plt.gca().set_xticks([0, 5, 10], [-1, 0, 1], size=27)
sns.despine()
plt.tight_layout()
- plt.savefig("example type 2")
- plt.show()
+ plt.savefig("example type 2", transparent=True)
+ plt.close()
state_num = 4
defined_points, pmf = all_first_pmfs_typeless['ibl_witten_17'][state_num][0], all_first_pmfs_typeless['ibl_witten_17'][state_num][1]
@@ -334,14 +334,14 @@ if __name__ == "__main__":
plt.ylim(0, 1)
plt.xlim(0, 10)
- plt.ylabel("P(rightwards)", size=32)
- plt.xlabel("Contrast", size=32)
+ # plt.ylabel("P(rightwards)", size=32)
+ # plt.xlabel("Contrast", size=32)
plt.yticks([0, 1], size=27)
plt.gca().set_xticks([0, 5, 10], [-1, 0, 1], size=27)
sns.despine()
plt.tight_layout()
- plt.savefig("example type 3")
- plt.show()
+ plt.savefig("example type 3", transparent=True)
+ plt.close()
n_rows, n_cols = 5, 6
@@ -364,8 +364,8 @@ if __name__ == "__main__":
a.set_xticks(np.arange(11), all_conts, size=tick_size, rotation=45)
a.set_yticks([0, 0.25, 0.5, 0.75, 1], [0, 0.25, 0.5, 0.75, 1], size=tick_size)
plt.tight_layout()
- plt.savefig("animals 1")
- plt.show()
+ plt.savefig("./summary_figures/animals 1")
+ plt.close()
n_rows, n_cols = 5, 6
_, axs = plt.subplots(n_rows, n_cols, figsize=(16, 9))
@@ -393,8 +393,8 @@ if __name__ == "__main__":
a.set_xticks(np.arange(11), all_conts, size=tick_size, rotation=45)
a.set_yticks([0, 0.25, 0.5, 0.75, 1], [0, 0.25, 0.5, 0.75, 1], size=tick_size)
plt.tight_layout()
- plt.savefig("animals 2")
- plt.show()
+ plt.savefig("./summary_figures/animals 2")
+ plt.close()
n_rows, n_cols = 5, 6
_, axs = plt.subplots(n_rows, n_cols, figsize=(16, 9))
@@ -422,8 +422,8 @@ if __name__ == "__main__":
a.set_xticks(np.arange(11), all_conts, size=tick_size, rotation=45)
a.set_yticks([0, 0.25, 0.5, 0.75, 1], [0, 0.25, 0.5, 0.75, 1], size=tick_size)
plt.tight_layout()
- plt.savefig("animals 3")
- plt.show()
+ plt.savefig("./summary_figures/animals 3")
+ plt.close()
n_rows, n_cols = 5, 6
_, axs = plt.subplots(n_rows, n_cols, figsize=(16, 9))
@@ -451,8 +451,8 @@ if __name__ == "__main__":
a.set_xticks(np.arange(11), all_conts, size=tick_size, rotation=45)
a.set_yticks([0, 0.25, 0.5, 0.75, 1], [0, 0.25, 0.5, 0.75, 1], size=tick_size)
plt.tight_layout()
- plt.savefig("animals 4")
- plt.show()
+ plt.savefig("./summary_figures/animals 4")
+ plt.close()
# Collection of PMFs which change state type
all_changing_pmfs = pickle.load(open("changing_pmfs.p", 'rb'))
@@ -477,11 +477,49 @@ if __name__ == "__main__":
plt.gca().set_xticks([0, 5, 10], [-1, 0, 1], size=16)
plt.gca().set_yticks([0, 0.25, 0.5, 0.75, 1], [0, 0.25, 0.5, 0.75, 1], size=16)
- if i + 1 == 30:
+ if i + 1 == 35:
break
plt.tight_layout()
- plt.savefig("changing pmfs")
+ plt.savefig("./summary_figures/changing pmfs")
+ plt.close()
+
+ # The biases of all first PMFs as a histogram
+ biases = []
+ mouse_counter = 0
+ consistent_bias = 0
+ hm = 0
+ for key in all_first_pmfs_typeless:
+ left_1, right_1, left_2, right_2 = False, False, False, False
+ x = all_first_pmfs_typeless[key]
+ for pmf in x:
+ defined_points, pmf = pmf
+ if pmf_type(pmf) == 0:
+ biases.append(np.mean(pmf[[0, 1, -2, -1]]))
+ if np.mean(pmf[[0, 1, -2, -1]]) >= 0.5:
+ right_1 = True
+ if np.mean(pmf[[0, 1, -2, -1]]) <= 0.5:
+ left_1 = True
+ if pmf_type(pmf) == 1:
+ if np.mean(pmf[[0, 1, -2, -1]]) >= 0.5:
+ right_2 = True
+ if np.mean(pmf[[0, 1, -2, -1]]) <= 0.5:
+ left_2 = True
+ mouse_counter += 1
+ consistent_bias += (left_1 and left_2) or (right_1 and right_2)
+ hm += left_1 and left_2 and right_1 and right_2
+
+
+ print("number of mice is {}".format(mouse_counter))
+ print("of which {} have a previously expressed bias in type 2".format(consistent_bias))
+ print(hm)
+ # number of mice is 113
+ # of which 77 have a previously expressed bias in type 2
+ # 17
+
+ plt.hist(biases)
+ plt.axvline(np.mean(biases))
+ plt.xlim(0, 1)
plt.show()
# All first PMFs
@@ -497,7 +535,13 @@ if __name__ == "__main__":
for key in all_first_pmfs_typeless:
x = all_first_pmfs_typeless[key]
for pmf in x:
- axs[pmf_type(pmf[1])].plot(np.where(pmf[0])[0], pmf[1][pmf[0]], c=type2color[pmf_type(pmf[1])])
+
+ defined_points, pmf = pmf
+ pmf_min = min(pmf[0], pmf[1])
+ pmf_max = max(pmf[-2], pmf[-1])
+ defined_points = np.logical_and(np.logical_and(defined_points, ~ (pmf > pmf_max)), ~ (pmf < pmf_min))
+ axs[pmf_type(pmf)].plot(np.where(defined_points)[0], pmf[defined_points], c=type2color[pmf_type(pmf)])
+ # axs[pmf_type(pmf[1])].plot(np.where(pmf[0])[0], pmf[1][pmf[0]], c=type2color[pmf_type(pmf[1])])
axs[0].set_ylim(0, 1)
axs[1].set_ylim(0, 1)
axs[2].set_ylim(0, 1)
@@ -517,13 +561,14 @@ if __name__ == "__main__":
axs[0].set_xlabel("Contrasts", size=label_size)
plt.tight_layout()
- plt.savefig(save_title)
+ plt.savefig("./summary_figures/" + save_title)
plt.show()
if save_title == "KS014 types":
quit()
# Type 1 PMF specifications
counter = 0
+ counter_l, counter_r = 0, 0
fig, ax = plt.subplots(1, 3, figsize=(16, 9))
for key in all_first_pmfs_typeless:
for defined_points, pmf in all_first_pmfs_typeless[key]:
@@ -537,6 +582,8 @@ if __name__ == "__main__":
use_ax = 2
else:
use_ax = int(pmf[0] > 1 - pmf[-1])
+ counter_l += use_ax == 0
+ counter_r += use_ax == 1
pmf_min = min(pmf[0], pmf[1])
pmf_max = max(pmf[-2], pmf[-1])
@@ -561,9 +608,9 @@ if __name__ == "__main__":
ax[2].set_yticks([])
ax[2].spines[['right', 'top']].set_visible(False)
ax[2].set_xticks(np.arange(11), all_conts, size=tick_size, rotation=45)
- print(counter)
+ print("left biased pmfs {}, right biased pmfs {}, neutral pmfs {}".format(counter_l, counter_r, counter))
plt.tight_layout()
- plt.savefig("differentiate type 2")
+ plt.savefig("./summary_figures/differentiate type 2")
plt.show()
diff --git a/pmf_weight_analysis.py b/analysis_pmf_weights.py
similarity index 61%
rename from pmf_weight_analysis.py
rename to analysis_pmf_weights.py
index 05985150..a2b06196 100644
--- a/pmf_weight_analysis.py
+++ b/analysis_pmf_weights.py
@@ -32,6 +32,115 @@ def pmf_type_rew(weights):
else:
return 2
+all_weight_trajectories = pickle.load(open("multi_chain_saves/all_weight_trajectories.p", 'rb'))
+
+
+# for weight_traj in all_weight_trajectories:
+# if len(weight_traj) == 1:
+# continue
+# state_type = pmf_type(weights_to_pmf(weight_traj[0]))
+
+# if state_type == 2 and weight_traj[0][0] > 0.8:
+# plt.plot(weights_to_pmf(weight_traj[0]))
+# plt.ylim(0, 1)
+# plt.show()
+
+n_weights = all_weight_trajectories[0][0].shape[0]
+ylabels = ["Cont left", "Cont right", "Persevere", "Bias left", "Bias right"]
+
+for min_dur in [2, 3, 4, 5, 7, 9, 11, 15]:
+ f, axs = plt.subplots(n_weights + 1, 3, figsize=(12, 9)) # add another space to split bias into stating left versus starting right
+ average = np.zeros((n_weights + 1, 3, 2))
+ counter = np.zeros((n_weights + 1, 3))
+ all_datapoints = [[[[], []], [[], []], [[], []]], [[[], []], [[], []], [[], []]], [[[], []], [[], []], [[], []]], [[[], []], [[], []], [[], []]], [[[], []], [[], []], [[], []]]]
+ for weight_traj in all_weight_trajectories:
+
+ if len(weight_traj) < min_dur:
+ continue
+
+ state_type = pmf_type(weights_to_pmf(weight_traj[0]))
+
+ for i in range(n_weights):
+ if i == n_weights - 1:
+ axs[i + (weight_traj[0][i] > 0), state_type].plot([0, 1], [weight_traj[0][i], weight_traj[-1][i]], marker="o") # plot how the weight evolves from first to last appearance
+ average[i + (weight_traj[0][i] > 0), state_type] += np.array([weight_traj[0][i], weight_traj[-1][i]])
+ counter[i + (weight_traj[0][i] > 0), state_type] += 1
+ all_datapoints[i + (weight_traj[0][i] > 0)][state_type][0].append(weight_traj[0][i])
+ all_datapoints[i + (weight_traj[0][i] > 0)][state_type][1].append(weight_traj[-1][i])
+ else:
+ axs[i, state_type].plot([0, 1], [weight_traj[0][i], weight_traj[-1][i]], marker="o") # plot how the weight evolves from first to last appearance
+ average[i, state_type] += np.array([weight_traj[0][i], weight_traj[-1][i]])
+ counter[i, state_type] += 1
+ all_datapoints[i][state_type][0].append(weight_traj[0][i])
+ all_datapoints[i][state_type][1].append(weight_traj[-1][i])
+ for i in range(3):
+ for j in range(n_weights + 1):
+ axs[j, i].set_ylim(-9, 9)
+ axs[j, i].spines['top'].set_visible(False)
+ axs[j, i].spines['right'].set_visible(False)
+ axs[j, i].set_xticks([])
+ if i == 0:
+ axs[j, i].set_ylabel(ylabels[j])
+ else:
+ axs[j, i].set_yticks([])
+ if j == 0:
+ axs[j, i].set_title("Type {}".format(i + 1))
+ if j == n_weights:
+ axs[j, i].set_xlabel("Lifetime weight change")
+
+ plt.tight_layout()
+ plt.savefig("./summary_figures/weight_changes/all weight changes min dur {}".format(min_dur))
+ plt.close()
+
+ f, axs = plt.subplots(n_weights + 1, 3, figsize=(12, 9))
+ for i in range(3):
+ for j in range(n_weights + 1):
+ axs[j, i].plot([0, 1], average[j, i] / counter[j, i], marker="o")
+ if j < 2:
+ axs[j, i].set_ylim(-4.5, 4.5)
+ else:
+ axs[j, i].set_ylim(-2, 2)
+ axs[j, i].spines['top'].set_visible(False)
+ axs[j, i].spines['right'].set_visible(False)
+ axs[j, i].set_xticks([])
+ if i == 0:
+ axs[j, i].set_ylabel(ylabels[j])
+ else:
+ axs[j, i].set_yticks([])
+ if j == 0:
+ axs[j, i].set_title("Type {}".format(i + 1))
+ if j == n_weights:
+ axs[j, i].set_xlabel("Lifetime weight change")
+ plt.savefig("./summary_figures/weight_changes/weight changes min dur {}".format(min_dur))
+ plt.close()
+
+ f, axs = plt.subplots(n_weights + 1, 3, figsize=(12, 9))
+ for i in range(3):
+ for j in range(n_weights + 1):
+ axs[j, i].plot([0, 1], average[j, i] / counter[j, i])
+ axs[j, i].boxplot(all_datapoints[j][i], positions=[0, 1])
+ if j < 2:
+ axs[j, i].set_ylim(-8, 8)
+ else:
+ axs[j, i].set_ylim(-4, 4)
+ axs[j, i].spines['top'].set_visible(False)
+ axs[j, i].spines['right'].set_visible(False)
+ axs[j, i].set_xticks([])
+ if i == 0:
+ axs[j, i].set_ylabel(ylabels[j])
+ else:
+ axs[j, i].set_yticks([])
+ if j == 0:
+ axs[j, i].set_title("Type {}".format(i + 1))
+ if j == n_weights:
+ axs[j, i].set_xlabel("Lifetime weight change")
+ plt.savefig("./summary_figures/weight_changes/weight change boxplot min dur {}".format(min_dur))
+ plt.close()
+
+quit()
+
+
+
# all pmf weights
apw = np.array(pickle.load(open("all_pmf_weights.p", 'rb')))
diff --git a/analysis_regression.py b/analysis_regression.py
index 575f234c..8e09fa87 100644
--- a/analysis_regression.py
+++ b/analysis_regression.py
@@ -11,11 +11,13 @@ if __name__ == "__main__":
regressions = np.array(pickle.load(open("regressions.p", 'rb')))
regression_diffs = np.array(pickle.load(open("regression_diffs.p", 'rb')))
- assert (regressions[:, 0] == np.sum(regressions[:, 2:], 1)).all()
+ assert (regressions[:, 0] == np.sum(regressions[:, 2:], 1)).all() # total # of regressions must be sum of # of regressions per type
print(pearsonr(regressions[:, 0], regressions[:, 1]))
+ # (0.6202414016960471, 2.3666819600330215e-13)
offset = 0.125
+ plt.figure(figsize=(16 * 0.9, 9 * 0.9))
# which x values exist
for x in np.unique(regressions[:, 0]):
# which and how many ys are associated with this x
@@ -28,7 +30,7 @@ if __name__ == "__main__":
sns.despine()
plt.tight_layout()
- plt.savefig("Regression vs session length")
+ plt.savefig("./summary_figures/Regression vs session length")
plt.show()
# histogram of regressions per type
@@ -38,7 +40,7 @@ if __name__ == "__main__":
sns.despine()
plt.tight_layout()
- plt.savefig("# of regressions per type")
+ plt.savefig("./summary_figures/# of regressions per type")
plt.show()
# histogram of number of mice with regressions in the different types
@@ -49,7 +51,7 @@ if __name__ == "__main__":
sns.despine()
plt.tight_layout()
- plt.savefig("# of mice with regressions per type")
+ plt.savefig("./summary_figures/# of mice with regressions per type")
plt.show()
# histogram of regression diffs
@@ -59,5 +61,5 @@ if __name__ == "__main__":
sns.despine()
plt.tight_layout()
- plt.savefig("Regression diffs")
+ plt.savefig("./summary_figures/Regression diffs")
plt.show()
diff --git a/analysis_state_intros.py b/analysis_state_intros.py
index 91b3250a..b8435d95 100644
--- a/analysis_state_intros.py
+++ b/analysis_state_intros.py
@@ -44,14 +44,14 @@ def type_hist(data, title=''):
plt.ylabel("Type 3", size=fontsize)
ax2.set_yticks([])
- plt.savefig(title)
+ plt.savefig("./summary_figures/" + title)
plt.show()
if __name__ == "__main__":
- all_intros = pickle.load(open("all_intros.p", 'rb'))
- all_intros_div = pickle.load(open("all_intros_div.p", 'rb'))
- all_states_per_type = pickle.load(open("all_states_per_type.p", 'rb'))
+ all_intros = np.array(pickle.load(open("all_intros.p", 'rb')))
+ all_intros_div = np.array(pickle.load(open("all_intros_div.p", 'rb')))
+ all_states_per_type = np.array(pickle.load(open("all_states_per_type.p", 'rb')))
# There are 5 mice with 0 type 2 intros, but only 3 mice with no type 2 stats.
# That is because they come up, but don't explain the necessary 50% to start phase 2.
diff --git a/behavioral_state_data_easier.py b/behavioral_state_data_easier.py
index 18a9de78..aed7c5dc 100644
--- a/behavioral_state_data_easier.py
+++ b/behavioral_state_data_easier.py
@@ -5,6 +5,8 @@ import pandas as pd
import seaborn as sns
import pickle
import json
+import os
+import re
one = ONE()
@@ -34,20 +36,48 @@ 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 = ['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'))
elif fit_type == 'prebias':
loading_info = json.load(open("canonical_infos.json", 'r'))
+bwm = ['NYU-11', 'NYU-12', 'NYU-21', 'NYU-27', 'NYU-30', 'NYU-37',
+ 'NYU-39', 'NYU-40', 'NYU-45', 'NYU-46', 'NYU-47', 'NYU-48',
+ 'CSHL045', 'CSHL047', 'CSHL049', 'CSHL051', 'CSHL052', 'CSHL053',
+ 'CSHL054', 'CSHL055', 'CSHL058', 'CSHL059', 'CSHL060', 'UCLA005',
+ 'UCLA006', 'UCLA011', 'UCLA012', 'UCLA014', 'UCLA015', 'UCLA017',
+ 'UCLA033', 'UCLA034', 'UCLA035', 'UCLA036', 'UCLA037', 'KS014',
+ 'KS016', 'KS022', 'KS023', 'KS042', 'KS043', 'KS044', 'KS045',
+ 'KS046', 'KS051', 'KS052', 'KS055', 'KS084', 'KS086', 'KS091',
+ 'KS094', 'KS096', 'DY_008', 'DY_009', 'DY_010', 'DY_011', 'DY_013',
+ 'DY_014', 'DY_016', 'DY_018', 'DY_020', 'PL015', 'PL016', 'PL017',
+ 'PL024', 'SWC_042', 'SWC_043', 'SWC_060', 'SWC_061', 'SWC_066',
+ 'ZFM-01576', 'ZFM-01577', 'ZFM-01592', 'ZFM-01935', 'ZFM-01936',
+ 'ZFM-01937', 'ZFM-02368', 'ZFM-02369', 'ZFM-02370', 'ZFM-02372',
+ 'ZFM-02373', 'ZM_1897', 'ZM_1898', 'ZM_2240', 'ZM_2241', 'ZM_2245',
+ 'ZM_3003', 'SWC_038', 'SWC_039', 'SWC_052', 'SWC_053', 'SWC_054',
+ 'SWC_058', 'SWC_065', 'NR_0017', 'NR_0019', 'NR_0020', 'NR_0021',
+ 'NR_0027', 'ibl_witten_13', 'ibl_witten_17', 'ibl_witten_18',
+ 'ibl_witten_19', 'ibl_witten_20', 'ibl_witten_25', 'ibl_witten_26',
+ 'ibl_witten_27', 'ibl_witten_29', 'CSH_ZAD_001', 'CSH_ZAD_011',
+ 'CSH_ZAD_019', 'CSH_ZAD_022', 'CSH_ZAD_024', 'CSH_ZAD_025',
+ 'CSH_ZAD_026', 'CSH_ZAD_029']
+regexp = re.compile(r'canonical_result_((\w|-)+)_prebias.p')
+subjects = []
+for filename in os.listdir("./multi_chain_saves/"):
+ if not (filename.startswith('canonical_result_') and filename.endswith('.p')):
+ continue
+ result = regexp.search(filename)
+ if result is None:
+ continue
+ subject = result.group(1)
+ subjects.append(subject)
already_fit = list(loading_info.keys())
-remaining_subs = [s for s in subjects if s not in amiss and s not in already_fit]
-print(remaining_subs)
-
+# remaining_subs = [s for s in subjects if s not in amiss and s not in already_fit]
+# print(remaining_subs)
-
-data_folder = 'session_data_test'
+data_folder = 'session_data'
old_style = False
if old_style:
@@ -64,22 +94,30 @@ names = []
pre_bias = []
entire_training = []
+training_status_reached = []
+actually_existing = []
for subject in subjects:
+ if subject in bwm:
+ continue
print('_____________________')
print(subject)
# if subject in already_fit or subject in amiss:
# continue
-
- trials = one.load_aggregate('subjects', subject, '_ibl_subjectTrials.table')
+ try:
+ 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'))
- .sort_values(by='session_start_time', kind='stable'))
+
+ training = one.load_aggregate('subjects', subject, '_ibl_subjectTraining.table')
+
+ trials = (trials
+ .set_index('session')
+ .join(training.set_index('session'))
+ .sort_values(by='session_start_time', kind='stable'))
+ actually_existing.append(subject)
+ except:
+ continue
start_times, indices = np.unique(trials.session_start_time, return_index=True)
start_times = [trials.session_start_time[index] for index in sorted(indices)]
@@ -99,16 +137,19 @@ for subject in subjects:
easy_per = np.zeros(len(eids))
hard_per = np.zeros(len(eids))
bias_start = 0
+ ephys_start = 0
info_dict = {'subject': subject, 'dates': [st.to_pydatetime() for st in start_times], 'eids': eids}
contrast_set = {0, 1, 9, 10}
rel_count = -1
- quit()
+
for i, start_time in enumerate(start_times):
rel_count += 1
+ assert rel_count == i
+
df = trials[trials.session_start_time == start_time]
df.loc[:, 'contrastRight'] = df.loc[:, 'contrastRight'].fillna(0)
df.loc[:, 'contrastLeft'] = df.loc[:, 'contrastLeft'].fillna(0)
@@ -137,9 +178,15 @@ for subject in subjects:
bias_start = i
print("bias start {}".format(rel_count))
info_dict['bias_start'] = rel_count
+ training_status_reached.append(set(df.training_status))
if bias_start < 33:
short_subjs.append(subject)
+ if ephys_start == 0 and df.task_protocol[0].startswith('_iblrig_tasks_ephysChoiceWorld'):
+ ephys_start = i
+ print("ephys start {}".format(rel_count))
+ info_dict['ephys_start'] = rel_count
+
pickle.dump(df, open("./{}/{}_df_{}.p".format(data_folder, subject, rel_count), "wb"))
side_info = np.zeros((len(df), 2))
diff --git a/canonical_infos.json b/canonical_infos.json
index 99ae626e..9e26dfee 100644
--- a/canonical_infos.json
+++ b/canonical_infos.json
@@ -1 +1 @@
-{"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
+{}
\ No newline at end of file
diff --git a/dyn_glm_chain_analysis.py b/dyn_glm_chain_analysis.py
index 1b806222..37bb5dad 100644
--- a/dyn_glm_chain_analysis.py
+++ b/dyn_glm_chain_analysis.py
@@ -22,7 +22,24 @@ import multiprocessing as mp
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
+import re
+performance_points = np.array([-1, -1, 0, 0])
+
+def pmf_to_perf(pmf):
+ # determine performance of a pmf, but only on the omnipresent strongest contrasts
+ return np.mean(np.abs(performance_points + pmf[[0, 1, -2, -1]]))
+
+# def pmf_type_temp(pmf):
+# rew = pmf_to_perf(pmf)
+# if rew < 0.6:
+# return 0
+# elif rew < 0.7827:
+# if np.abs(pmf[0] + pmf[-1] - 1) <= 0.1:
+# return 3
+# return 1
+# else:
+# return 2
colors = np.genfromtxt('colors.csv', delimiter=',')
@@ -45,15 +62,6 @@ def weights_to_pmf(weights, with_bias=1):
psi = weights[0] * contrasts_R + weights[1] * contrasts_L + with_bias * weights[-1]
return 1 / (1 + np.exp(psi)) # we somehow got the answers twisted, so we drop the minus here to get the opposite response probability for plotting
-performance_points = np.array([-1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0])
-def pmf_to_perf(pmf, def_points):
- # determine performance of a pmf
- # we use this to determine regressions in the behaviour of animals
- # therefore, we exclude 0 as performance on it is 0.5 regardless of PMF, but it might
- # overall lower performance. The removal of 0.5 later might also be a problem, let's see
- relevant_points = def_points
- relevant_points[5] = False
- return np.mean(np.abs(performance_points[relevant_points] + pmf[relevant_points]))
class MCMC_result_list:
@@ -305,13 +313,13 @@ class MCMC_result_list:
state_appear = []
state_durs = []
- state_appear_dist = np.zeros(10)
+ state_appear_dist = np.zeros(11)
for res in self.results:
for j in range(take_n):
sa, sd = state_appear_and_dur(res.models[self.m // take_n * j], self)
state_appear += sa
state_durs += sd
- state_appear_dist[[int(s * 10) for s in sa]] += 1
+ state_appear_dist[[int(s * 11) for s in sa]] += 1
return state_appear, state_durs, state_appear_dist / self.m / take_n
@@ -392,7 +400,6 @@ class MCMC_result:
state_start_save[20 * index // len(seq)] += 1
state_starts[session_bounds[i] + index] += 1
observed_states.append(s)
- print('observed_states {}'.format(observed_states))
plt.plot(state_starts / self.n_samples)
return session_bounds, state_start_save / self.n_samples
@@ -625,6 +632,7 @@ def bias_flips(states_by_session, pmfs, durs):
def pmf_regressions(states_by_session, pmfs, durs):
# find out whether pmfs regressed
+ # return: [total # of regressions, # of sessions, # of regressions during type 1, type 2, type 3], the reward differences
state_perfs = {}
state_counter = {}
current_best_state = -1
@@ -1077,7 +1085,7 @@ def state_development(test, state_sets, indices, save=True, save_append='', show
ax2.plot(np.where(defined_points)[0] / (len(defined_points)-1), pmf[defined_points] - 0.5 + test.state_mapping[state], color=cmap(0.2 + 0.8 * j / test.results[0].n_sessions))
ax2.plot(np.where(defined_points)[0] / (len(defined_points)-1), pmf[defined_points] - 0.5 + test.state_mapping[state], ls='', ms=7, marker='*', color=cmap(j / test.results[0].n_sessions))
all_pmfs.append((defined_points, pmfs))
- all_pmf_weights += pmf_weights
+ all_pmf_weights.append(pmf_weights)
else:
temp = np.percentile(pmfs, [2.5, 97.5], axis=0)
if not test.state_mapping[state] in dont_plot:
@@ -1382,7 +1390,7 @@ def state_type_durs(states, pmfs):
# Takes states and pmfs, first creates an array of when which type is how active, then computes the number of sessions each type lasts.
# A type lasts until a more advanced type takes up more than 50% of a session (and cannot return)
# Returns the durations for all the different state types, and an array which holds the state percentages
- state_types = np.zeros((3, states.shape[1]))
+ state_types = np.zeros((4, states.shape[1]))
for s, pmf in zip(states, pmfs):
pmf_counter = -1
for i in range(states.shape[1]):
@@ -1427,9 +1435,9 @@ def state_cluster_interpolation(states, pmfs):
def get_first_pmfs(states, pmfs):
# get the first pmf of every type, also where they are defined, and whether they are the first pmf of that state
- earliest_sessions = [1000, 1000, 1000]
- first_pmfs = [0, 0, 0, 0, 0, 0, 0, 0, 0]
- changing_pmfs = [[0, 0], [0, 0]]
+ earliest_sessions = [1000, 1000, 1000] # high values to compare with min
+ first_pmfs = [0, 0, 0, 0, 0, 0, 0, 0, 0] # for every type (3), we save: the pmf, the defined points, the session at which they appear
+ changing_pmfs = [[0, 0], [0, 0]] # if the new type appears first through a slow transition, we remember its defined points (pmf[1]) and its pmf trajectory pmf[1]
for state, pmf in zip(states, pmfs):
sessions = np.where(state)[0]
for i, (sess_pmf, sess) in enumerate(zip(pmf[1], sessions)):
@@ -1532,26 +1540,25 @@ if __name__ == "__main__":
loading_info = json.load(open("canonical_infos_bias.json", 'r'))
elif fit_type == 'prebias':
loading_info = json.load(open("canonical_infos.json", 'r'))
- no_good_pcas = ['NYU-06', 'SWC_023']
- subjects = list(loading_info.keys())
-
- 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()
+ subjects = []
+ regexp = re.compile(r'canonical_result_((\w|-)+)_prebias.p')
+ for filename in os.listdir("./multi_chain_saves/"):
+ if not (filename.startswith('canonical_result_') and filename.endswith('.p')):
+ continue
+ result = regexp.search(filename)
+ if result is None:
+ continue
+ subject = result.group(1)
+ subjects.append(subject)
- print(subjects)
+ print(len(subjects))
fit_variance = [0.03, 0.002, 0.0005, 'uniform', 0, 0.008][0]
dur = 'yes'
# fig, ax = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
pop_state_starts = np.zeros(20)
- state_appear_dist = np.zeros(10)
+ state_appear_dist = np.zeros(11)
state_appear_mode = []
num_states = []
num_sessions = []
@@ -1581,6 +1588,9 @@ if __name__ == "__main__":
regression_diffs = []
all_bias_flips = []
all_pmf_weights = []
+ all_weight_trajectories = []
+ bias_sessions = []
+ first_and_last_pmf = []
new_counter, transform_counter = 0, 0
state_types_interpolation = np.zeros((3, 150))
@@ -1589,17 +1599,29 @@ if __name__ == "__main__":
state_nums_5 = []
state_nums_10 = []
- for subject in subjects:
+ ultimate_counter = 0
- # NYU-11 is quite weird, super errrativ behaviour, has all contrasts introduced at once, no good session at end
- if subject.startswith('GLM_Sim_'):
+ for subject in subjects:
+ if subject.startswith('GLM_Sim_') or subject in ['SWC_065', 'ZFM-05245', 'ZFM-04019', 'ibl_witten_18']:
+ # ibl_witten_18 is a weird one, super good session in the middle, ending phase 1, never to re-appear, bad at the end
+ # ZFM-05245 is neuromodulator mouse, never reaches ephys it seems... same for ZFM-04019
+ # SWC_065 never reaches type 3
continue
print()
print(subject)
+ continue
test = pickle.load(open("multi_chain_saves/canonical_result_{}_{}.p".format(subject, fit_type), 'rb'))
+ # 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
+
mode_specifier = 'first'
try:
mode_indices = pickle.load(open("multi_chain_saves/{}_mode_indices_{}_{}.p".format(mode_specifier, subject, fit_type), 'rb'))
@@ -1609,6 +1631,9 @@ if __name__ == "__main__":
mode_indices = pickle.load(open("multi_chain_saves/mode_indices_{}_{}.p".format(subject, fit_type), 'rb'))
state_sets = pickle.load(open("multi_chain_saves/state_sets_{}_{}.p".format(subject, fit_type), 'rb'))
except:
+ print("____________________________________")
+ print("Something quite wrong with {}".format(subject))
+ print("____________________________________")
continue
# lapse differential
@@ -1619,18 +1644,37 @@ if __name__ == "__main__":
# _ = 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()
+ 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, save=False, show=0, separate_pmf=1, type_coloring=True)
+
+ first_and_last_pmf.append((pmf_weights[np.argmax(states[:, 0])][0], pmf_weights[np.argmax(states[:, -1])][-1]))
+
+ continue
+ # all_weight_trajectories += pmf_weights
+ abs_state_durs.append(durs)
+ ultimate_counter += 1
+
+ info_dict = pickle.load(open("./{}/{}_info_dict.p".format('session_data', subject), "rb"))
+ if 'ephys_start' in info_dict:
+ bias_sessions.append(info_dict['ephys_start'] - info_dict['bias_start'])
+ else:
+ bias_sessions.append(info_dict['n_sessions'] - info_dict['bias_start'])
+ print(subject, info_dict['n_sessions'], info_dict['bias_start'])
+
+ continue
+
+ 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])
# 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()
+ # 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)
+ # quit()
+
new = type_2_appearance(states, pmfs)
if new == 2:
@@ -1644,60 +1688,50 @@ if __name__ == "__main__":
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)
+ # 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'))
+
+ all_pmf_weights += [item for sublist in pmf_weights for item in sublist]
+ all_state_types.append(state_types)
+
+ 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
+ quit()
# 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
+ regression, diffs = pmf_regressions(states, pmfs, durs)
+ regression_diffs += diffs
+ regressions.append(regression)
+
# 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)
+
+ 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)
@@ -1745,13 +1779,7 @@ if __name__ == "__main__":
# 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)
@@ -1759,7 +1787,6 @@ if __name__ == "__main__":
# plt.savefig("temp")
# plt.close()
-
# pickle.dump(all_first_pmfs, open("all_first_pmfs.p", 'wb'))
# pickle.dump(all_changing_pmfs, open("changing_pmfs.p", 'wb'))
# pickle.dump(all_changing_pmf_names, open("changing_pmf_names.p", 'wb'))
@@ -1774,24 +1801,50 @@ if __name__ == "__main__":
# pickle.dump(all_state_types, open("all_state_types.p", 'wb'))
# pickle.dump(all_pmf_weights, open("all_pmf_weights.p", 'wb'))
# pickle.dump(state_types_interpolation, open("state_types_interpolation.p", 'wb'))
+ # pickle.dump(state_types_interpolation, open("state_types_interpolation_4_states.p", 'wb')) # special version, might not want to use
# abs_state_durs = np.array(abs_state_durs)
# pickle.dump(abs_state_durs, open("multi_chain_saves/abs_state_durs.p", 'wb'))
+ # pickle.dump(pop_state_starts, open("multi_chain_saves/pop_state_starts.p", 'wb'))
+ # pickle.dump(state_appear, open("multi_chain_saves/state_appear.p", 'wb'))
+ # pickle.dump(state_dur, open("multi_chain_saves/state_dur.p", 'wb'))
+ # pickle.dump(state_appear_dist, open("multi_chain_saves/state_appear_dist.p", 'wb'))
+ # pickle.dump(all_weight_trajectories, open("multi_chain_saves/all_weight_trajectories.p", 'wb'))
+ # pickle.dump(bias_sessions, open("multi_chain_saves/bias_sessions.p", 'wb'))
+ # pickle.dump(first_and_last_pmf, open("multi_chain_saves/first_and_last_pmf.p", 'wb'))
+
+ abs_state_durs = pickle.load(open("multi_chain_saves/abs_state_durs.p", 'rb'))
+ bias_sessions = pickle.load(open("multi_chain_saves/bias_sessions.p", 'rb'))
+
+ quit()
+ print("Ultimate count is {}".format(ultimate_counter))
- if True:
+ if False:
abs_state_durs = pickle.load(open("multi_chain_saves/abs_state_durs.p", 'rb'))
print("Correlations")
from scipy.stats import pearsonr
print(pearsonr(abs_state_durs[:, 0], abs_state_durs[:, 1]))
+ # (0.30202072153268694, 0.0011496023265844494)
print(pearsonr(abs_state_durs[:, 0], abs_state_durs[:, 2]))
+ # (-0.008129156810930915, 0.9318998938147801)
print(pearsonr(abs_state_durs[:, 1], abs_state_durs[:, 2]))
+ # (0.20293276001562754, 0.03110687851587209)
print(pearsonr(abs_state_durs.sum(1), abs_state_durs[:, 0]))
- # (0.7338297529946006, 2.6332570579118393e-06)
+ # (0.3930315541710777, 1.6612856471614082e-05)
print(pearsonr(abs_state_durs.sum(1), abs_state_durs[:, 1]))
- # (0.35094585023228597, 0.052897046343413114)
+ # (0.49183970714755426, 3.16067121387928e-08)
print(pearsonr(abs_state_durs.sum(1), abs_state_durs[:, 2]))
- # (0.7210260323745921, 4.747833912452452e-06)
+ # (0.8936241158982767, 2.0220064236645623e-40)
+
+ print(pearsonr(abs_state_durs[:, 0], bias_sessions))
+ # (-0.048942358403448086, 0.6099754179426583)
+ print(pearsonr(abs_state_durs[:, 1], bias_sessions))
+ # (-0.07904742327893216, 0.4095531301939973)
+ print(pearsonr(abs_state_durs[:, 2], bias_sessions))
+ # (-0.0697415785143183, 0.4670166281153171)
+ print(pearsonr(abs_state_durs.sum(1), bias_sessions))
+ # (-0.09311057488563208, 0.3310557265430209)
from simplex_plot import plotSimplex
@@ -1804,9 +1857,11 @@ if __name__ == "__main__":
plt.ylabel("# of mice", size=40)
plt.xlabel('# of sessions', size=40)
plt.tight_layout()
- plt.savefig("session_num_hist.png", dpi=300, transparent=True)
+ plt.savefig("./summary_figures/session_num_hist.png", dpi=300, transparent=True)
plt.show()
+ quit()
+
if False:
ax[0].set_ylim(0, 1)
ax[1].set_ylim(0, 1)
@@ -1850,11 +1905,12 @@ if __name__ == "__main__":
# ax[1].legend()
# ax[2].legend()
plt.tight_layout()
- plt.savefig("first_knowledge_state", dpi=300)
+ plt.savefig("./summary_figures/first_knowledge_state", dpi=300)
plt.show()
quit()
if False:
+ pop_state_starts = pickle.load(open("multi_chain_saves/pop_state_starts.p", 'rb'))
f, (ax, ax2) = plt.subplots(2, 1, sharex=True, figsize=(9, 9))
# plot the same data on both axes
@@ -1864,14 +1920,14 @@ if __name__ == "__main__":
ax2.bar(np.linspace(0, 1, 20), pop_state_starts, align='edge', width=1/19, color='grey')
# zoom-in / limit the view to different portions of the data
- ax.set_ylim(24.75, 28)
- ax2.set_ylim(0, 3.25)
+ ax.set_ylim(245, 300)
+ ax2.set_ylim(0, 25)
- ax.set_yticks([25, 27])
- ax.set_yticklabels([25, 27])
- ax2.set_yticks([0, 1, 2, 3])
+ ax.set_yticks([250, 275])
+ ax.set_yticklabels([250, 275])
+ ax2.set_yticks([0, 8, 16, 24])
ax2.set_xticks([0, .25, .5, .75, 1])
- ax2.set_yticklabels([0, 1, 2, 3])
+ ax2.set_yticklabels([0, 8, 16, 24])
ax2.set_xticklabels([0, .25, .5, .75, 1])
# hide the spines between ax and ax2
@@ -1896,13 +1952,13 @@ if __name__ == "__main__":
ax2.plot((-d, +d), (1 - d, 1 + d), **kwargs) # bottom-left diagonal
plt.tight_layout()
- plt.savefig('states within session', dpi=300)
+ plt.savefig('./summary_figures/states within session', dpi=300)
plt.close()
if False:
f, (a1, a2) = plt.subplots(1, 2, figsize=(18, 9), gridspec_kw={'width_ratios': [1.4, 1]})
- # a1.bar(np.linspace(0, 1, 11), state_appear_dist, align='edge', width=1/10, color='grey')
- a1.hist(state_appear_mode, color='grey')
+ a1.bar(np.linspace(0, 1, 11), state_appear_dist, align='edge', width=1/10, color='grey')
+ # a1.hist(state_appear_mode, color='grey')
a1.set_xlim(left=0)
# plt.title('First appearence of ', fontsize=22)
a1.set_ylabel('# of states', fontsize=35)
@@ -1932,5 +1988,5 @@ if __name__ == "__main__":
a2.tick_params(axis='both', labelsize=20)
sns.despine()
plt.tight_layout()
- plt.savefig('states in sessions', dpi=300)
+ plt.savefig('./summary_figures/states in sessions', dpi=300)
plt.close()
diff --git a/index_mice.py b/index_mice.py
index 96d35de5..6c42c4d2 100644
--- a/index_mice.py
+++ b/index_mice.py
@@ -9,14 +9,14 @@ test = ['CSHL059_fittype_prebias_var_0.03_303_45_1.p', 'ibl_witten_16_fittype_pr
prebias_subinfo = {}
bias_subinfo = {}
-for filename in os.listdir("./dynamic_GLMiHMM_crossvals/"):
+for filename in test:#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)
- if subject == 'ibl_witten_26':
- print('here')
+ print(subject)
+ continue
fit_type = result.group(3)
seed = result.group(4)
fit_num = result.group(5)
diff --git a/simplex_animation.py b/simplex_animation.py
index ddf6df79..dce03864 100644
--- a/simplex_animation.py
+++ b/simplex_animation.py
@@ -5,6 +5,7 @@ from analysis_pmf import type2color
import math
import matplotlib.pyplot as plt
import copy
+import string
all_state_types = pickle.load(open("all_state_types.p", 'rb'))
@@ -39,7 +40,15 @@ assert (test_count == 1).all()
# quit()
session_counter = -1
-alph = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'za', 'zb', 'zc', 'zd', 'ze', 'zf', 'zg', 'zh', 'zi', 'zj']
+string_list = []
+alphabet = string.ascii_lowercase
+
+for char1 in alphabet:
+ for char2 in alphabet:
+ string_list.append(char1 + char2)
+ if len(string_list) == 100:
+ break
+
# do as many sessions as it takes
while True:
@@ -58,7 +67,7 @@ while True:
else:
temp_counter = -1
- assert np.sum(sts[:, temp_counter]) <= 1
+ assert np.sum(sts[:, temp_counter]) <= 1.000000000000001
if (sts[:, temp_counter] == 1).any():
x_offset[i] = offsets[i][0] * 0.01
@@ -66,8 +75,8 @@ while True:
type_proportions[:, i] = sts[:, temp_counter]
- plotSimplex(type_proportions.T, x_offset=x_offset, y_offset=y_offset, c=np.arange(len(all_state_types)), show=False, title="Session {}".format(1 + session_counter),
- vertexcolors=[type2color[i] for i in range(3)], vertexlabels=['Type 1', 'Type 2', 'Type 3'], save_title="simplex_{}.png".format(alph[session_counter]))
+ plotSimplex(type_proportions.T * 10, x_offset=x_offset, y_offset=y_offset, c=np.arange(len(all_state_types)), show=False, title="Session {}".format(1 + session_counter),
+ vertexcolors=[type2color[i] for i in range(3)], vertexlabels=['Type 1', 'Type 2', 'Type 3'], save_title="simplex_{}.png".format(string_list[session_counter]))
if not_ended == 0:
break
diff --git a/simplex_plot.py b/simplex_plot.py
index 420dd26c..f9f8a070 100644
--- a/simplex_plot.py
+++ b/simplex_plot.py
@@ -15,7 +15,7 @@ import matplotlib.patches as PA
def plotSimplex(points, fig=None,
vertexlabels=['1: initial flat PMFs', '2: intermediate unilateral PMFs', '3: final bilateral PMFs'], title='',
- save_title="dur_simplex.png", show=False, vertexcolors=['k', 'k', 'k'], x_offset=0, y_offset=0, **kwargs):
+ save_title="./summary_figures/dur_simplex.png", show=False, vertexcolors=['k', 'k', 'k'], x_offset=0, y_offset=0, **kwargs):
"""
Plot Nx3 points array on the 3-simplex
(with optionally labeled vertices)
@@ -47,12 +47,12 @@ def plotSimplex(points, fig=None,
P.scatter(projected[:, 0], projected[:, 1], marker='*', color='r', s=np.mean(points.sum(1)) * 3.5)#s=50
# Leave some buffer around the triangle for vertex labels
- fig.gca().set_xlim(-0.05, 1.05)
- fig.gca().set_ylim(-0.05, 1.05)
+ fig.gca().set_xlim(-0.08, 1.08)
+ fig.gca().set_ylim(-0.08, 1.08)
P.axis('off')
if title != '':
- P.annotate(title, (0.395, np.sqrt(3) / 2 + 0.025), size=24)
+ P.annotate(title, (0.395, np.sqrt(3) / 2 + 0.075), size=24)
P.tight_layout()
P.savefig(save_title, bbox_inches='tight', dpi=300, transparent=True)
--
GitLab