From 2ea2f8eef269cbe74ab53376fd2e1f6ca512fb66 Mon Sep 17 00:00:00 2001 From: kgostic Date: Fri, 29 Mar 2019 14:03:27 -0700 Subject: [PATCH] Re-ran all code to check outputs and updated the numerical algebra check script --- Numerical_algebra_check.R | 26 ++++++++++++++++-- code_outputs/sessionInfo.txt | 13 +++++---- manuscript_figures/Fig2_par_ests.pdf | Bin 13849 -> 13849 bytes manuscript_figures/Fig3_model_fits.pdf | Bin 31420 -> 31420 bytes .../Fig4_prob_inf_given_dose_and_route.pdf | Bin 10652 -> 10659 bytes scratch_figures/Profiles_basic_model.pdf | Bin 29126 -> 29126 bytes scratch_figures/pI_profiles.pdf | Bin 8034 -> 8034 bytes 7 files changed, 30 insertions(+), 9 deletions(-) diff --git a/Numerical_algebra_check.R b/Numerical_algebra_check.R index d0c6ad2..7a15d55 100644 --- a/Numerical_algebra_check.R +++ b/Numerical_algebra_check.R @@ -7,16 +7,36 @@ ## Let Dw represent the reaches the within-host environment ## Let DI represent the infecting dose ## P(DI|d) = P(D0|d)P(Dw|D0)P(DI|Dw) ## The probability of a given infectious dose, given the expected inoculum size -## P(I|d) = sum_DI=1:Inf sum_Dw=DI:Inf sum_D0=Dw:Inf (P(D0|d)P(Dw|D0)P(DI|Dw)) ## Long form probability of infection, , given the expected inoculum size +## P(Infection|d) = sum_DI=1:Inf sum_Dw=DI:Inf sum_D0=Dw:Inf P(D0|d)P(Dw|D0)P(DI|Dw)) ## Long form probability of infection, , given the expected inoculum size ## Calculate the summand after algebraic rearrangement into three different poisson densities compare = function(d, D0, Dw, DI, pp, pc){ - ## This is the long-form equation for P(DI|d), as given in ## + ## This is the long-form equation for P(DI|d), as given in equation 4 of the main text methods original = dpois(x = D0, lambda = d)*dbinom(x = Dw, size = D0, prob = pc)*dbinom(x = DI, size = Dw, prob = pp) - ## This is an alternate representation of the summand, as shown in ## + ## This is an alternate representation of the summand, as shown in equation 5 of the main text methods rearranged = dpois(x = DI, lambda = d*pc*pp)*dpois(x = D0-Dw, lambda = d*(1-pc))*dpois(x = Dw-DI, lambda = d*pc*(1-pp)) return(c(original, rearranged)) } + +## Test the comparison function, which will return the value of the summand using the original long-form, and the algebraic rearrangement given in equation 5 compare(1000, 998, 100, 20, .15, .001) + + +## Test for 1000 randomly chosen values of d, D0, Dw, DI, pp and pc +d.test = round(runif(n = 1000, 1, 10^9)) +D0.test = sapply(d.test, function(dd) rpois(n = 1, dd)) +Dw.test = round(sapply(D0.test, function(D0) runif(1, 1, D0))) +DI.test = round(sapply(Dw.test, function(Dw) runif(1, 1, Dw))) +pp.test = runif(1000, 0, 1) +pc.test = runif(1000, 0, 1) + +## Output comparison +test.outputs = mapply(FUN = compare, d= d.test, D0 = D0.test, Dw = Dw.test, DI = DI.test, pp = pp.test, pc = pc.test) + +## Verify that both versions of the summand give exactly the same answer +test.outputs[1,]-test.outputs[2,] + +## Q.E.D. + diff --git a/code_outputs/sessionInfo.txt b/code_outputs/sessionInfo.txt index e2b8149..034c2c6 100644 --- a/code_outputs/sessionInfo.txt +++ b/code_outputs/sessionInfo.txt @@ -19,10 +19,11 @@ other attached packages: [22] ggplot2_3.1.0 loaded via a namespace (and not attached): - [1] tidyselect_0.2.5 haven_2.0.0 lattice_0.20-38 oompaData_3.1.1 colorspace_1.4-0 generics_0.0.2 yaml_2.2.0 - [8] utf8_1.1.4 rlang_0.3.1 pillar_1.3.1 glue_1.3.0 withr_2.1.2 BiocGenerics_0.26.0 modelr_0.1.2 -[15] readxl_1.2.0 audio_0.1-5.1 bindr_0.1.1 plyr_1.8.4 munsell_0.5.0 gtable_0.2.0 cellranger_1.1.0 -[22] rvest_0.3.2 Biobase_2.40.0 labeling_0.3 fansi_0.4.0 broom_0.5.1 Rcpp_1.0.0 scales_1.0.0 -[29] backports_1.1.3 jsonlite_1.6 hms_0.4.2 digest_0.6.18 stringi_1.2.4 grid_3.5.0 cli_1.0.1 -[36] tools_3.5.0 beepr_1.3 lazyeval_0.2.1 cluster_2.0.7-1 crayon_1.3.4 pkgconfig_2.0.2 xml2_1.2.0 + [1] tidyselect_0.2.5 haven_2.0.0 lattice_0.20-38 oompaData_3.1.1 colorspace_1.4-0 generics_0.0.2 + [7] yaml_2.2.0 utf8_1.1.4 rlang_0.3.1 pillar_1.3.1 glue_1.3.0 withr_2.1.2 +[13] BiocGenerics_0.26.0 modelr_0.1.2 readxl_1.2.0 audio_0.1-5.1 bindr_0.1.1 plyr_1.8.4 +[19] munsell_0.5.0 gtable_0.2.0 cellranger_1.1.0 rvest_0.3.2 Biobase_2.40.0 labeling_0.3 +[25] fansi_0.4.0 broom_0.5.1 Rcpp_1.0.0 scales_1.0.0 backports_1.1.3 jsonlite_1.6 +[31] hms_0.4.2 digest_0.6.18 stringi_1.2.4 grid_3.5.0 cli_1.0.1 tools_3.5.0 +[37] beepr_1.3 lazyeval_0.2.1 cluster_2.0.7-1 crayon_1.3.4 pkgconfig_2.0.2 xml2_1.2.0 [43] assertthat_0.2.0 httr_1.4.0 rstudioapi_0.9.0 R6_2.3.0 nlme_3.1-137 compiler_3.5.0 diff --git a/manuscript_figures/Fig2_par_ests.pdf b/manuscript_figures/Fig2_par_ests.pdf index c19b7f87988667978430c91b0b2c3b69d55b0a82..6c2c5105f8ed2d524061a4a62d9d29c06e696776 100644 GIT binary patch delta 29 ecmbQ4Gc#v`1-qr8v7v#9=|p=`7;|IvK2rdN83>aA delta 29 ecmbQ4Gc#v`1-rSSxq+F9*+hF$7;|IvK2rdN+X$lo diff --git a/manuscript_figures/Fig3_model_fits.pdf b/manuscript_figures/Fig3_model_fits.pdf index 199f116f01b997865736f2ec9f4930d5f528d263..e615003f14452983fa117823316c3d9e2f11ef17 100644 GIT binary patch delta 31 gcmdn$5Li(78A|ri|C7e^cSN4W zOX9Ku8u4@uZp}mGg)B59wf7dwO;wfa5C0rY-jhtPTql9X2T7*z3&zTDf z>{EaD{wgA8J`OJ{;-%a>uw@SNnyE{h37o!ucf31eUPnz9xByo>)Aw9Oo(`-A0|E8( z*B<^p$69qa`vVg*E@md%VawQBM#mY^s~ac4e?JXqJRn5Z-ZWpZT#dH4+`m3MtaVv1 zZ`d|IThP)kySZ#N%X`;td#_P4xMI#^HFN?D_tzU?QY4tk3MJ;D{$lh-X3MPo9g%HE zEn{}t&xUv58xw3#ZVfp+%O#I{?S@dx&GC$-UNH~1kOWZncgxu!I4QKE&jEg&f zE_#y@vGSAGb<_7;S|?|&Dnz|aLYu>}J7LYin-z7|k3;`A{t>7Y^_mTBrU4%%PZ~=& zcg{yrUqo%kNLu|7sC~>YpZb;3U3M&%EcH%XsNBet)s^Y+hb{nmV40phnQ zFStJVX+CH6CuyDxJGLF-{M^3+$eN)g=l3+}G5f>>Gk<@*abow*>|c2H&TOx8W^?g8 zZ)JSj>wIxP;Lg$G>CL6;afESQ^0!B>>T~lr=&Tt$Y{*H6xeeoOGL1#^PTNm&n9G;2 z>Ftl#a{jdaG+SsKEWwFlMJrJL9|4bX);fl{pf9k2Kg}&%JOQyp3V9%?Is7j+tF>(& zF?>x#Pk6#TyFT;1Y)CZZ5}r1BdY0ww_qaH$ zUSD3Yfo3*mc&_$a{>r+TULeq(JlhO zw~M*HK>r3a!1O5bTk_@tTYWf|Y~p5rit+kr>vT8Jm(1G1-y3_eGjQnxyu04q#^{UI z`Ajf6qi_8D&-vzQ=e6RiYUovH^0hfuS9D3`bKWgo(NatNtsrc7se61vNE&4C<7huf zC;Olz5Vv#mi!ssbR@h$qku>#?xnOm)`Ryr=KwJ@UL9lXuFZ~Bwp^t4z;EO|->kpL| zT_dv9t!n~kod?u>~Z;t z=a2i2_LLqU4~}T@e^Twt-%p)gjDW1CJ(fw_7O1;`8nu5!A28bl1~P%7k;IDK_$)Rb zC?HD!1rsO+8KKy%z+$t^Z9}BRB;Nreq3M8jBehfILW4&rA?-vo4;7hha5(sE+FZ7b zx^ei8y59@-8g;ihdGp!KW{L3GG!fBKK^AvwkS}h`5|OyGA09&53O!}1AWx&eAwf)d zv20oGazET|!jW!*95u-Oz3S}XXn`>5d?Y};=tJYf1hQH>Wn&t3ffgWN)b7ZkLe$?o zj9r=slQs>_e?{)Ieh)=*v2KMTzp!3?MMlVeL^i|@>~3SX4ojq8`x;u>g``n}t;!A*kUkGXGJQRhn@TodZ4W(#s^5LhqqJrn-Etnk`$mZpRdzlYISq1V{K9i!9*1 zW9uSxohdxO;|vu+xB15KEC|;f=z!s4mrDDF5V+DQ8cu!R{N(y4-01!SRVQX@&kF+2 zti24{Y;Jpf9tFEn5x!<;I8+)3%nboMY{qfexVaJFh%G(FjGNN}&a=%Fc(1uj&P>4DtJCiW%z|Idjj}LedV_S>*4H7&0Rayn7M|^WIJe^_@z9?$af?K3+OI3F`Z}>#3O>b|^3;=m0U# z9BOUy8qDq+A}88{?>*(^eDGd2NRt12=&&Oc2Ado?(s@!R>@xS>9W)+eAZk5_^nyBh zS+jN)q(!mcr~nZ7w;r1Vv}qEzsIkU%b>G@jzF{c2=pY1M@0YYBf472^5fWs3hNR%_ z4C}X`;fH)zq3(p(cDp8lIXkQ4#5p@_NabjD64d2r$~y1rx)}EGCM+5D2qr8!^nT^< z+)QvY4dqFhu=Ul3Xxp42BW~M7Hwg&E_G$A#OQUN5$Z6CNz%XEAu<^|XqRH7gq^uIl z*OP1AJIQTdT#|U;j7eyCwd0VNJNa651BF14sS?G1v1HLq z>6UYGEhAr{bn>~RmrTkD5D>JZI>0I zu*Z{O`<^KNQTWn+)Uf{ZG0dvMHog>pjg6~GQhu|a#`55+y0;COHyVMc;MPa<=_!4e|0 z-&$;s&wvdia!D-?7Oxe*WE|%kADgrEX4$xFePAyDeEpObKt;={V6I}0m0rskN0eTf z3SSE{945DZRu5Hhk|a(1n#ENAA-+$(eL4w@26sMduobn2F+V^LYIrPSFT-pB>nwlO zl=6pn?t=n{wUvE6L=}7{N4HqdD5WpC0*JIYRgS?@AAXPmAA7CMT68ohy}Lz&3yuWk zTHI!O*1M_OU@b-9JGgt9W7LiTutD(nr2Rg1^koZ-`6x)Qbm!@qH+sNYsv~{J5`PGn zSH9cPi{lPqLck!Cx|g;g;Q$!OJt$dwY6s*P%|(1BZt7KJb`BhL``w)@_Brg4mplG) z_G_P?iD$y%`HfV1`48yWLP0vaPANkMfzyyxck6;pPDrh?nY5 zV#Vs;kUYxQKtHFO?_`w~X8=2J9^i(5Aq|G(TN!0oGJZ%w<_4%6_Ub+WptbOeLetwj z_13%;8+KjBE8a)x^}d=68@Kf$bX#9s`p9Fy-le3L5`DHO10=?-?@J0lXG1>|mdmTS zY9Gr`U9g8CRb?Du393`)S(e-mzq8!9LsnJsj<1&mOX=6S=`-5v7xX-nl&wj?2kTVR zCMn%+q(kb#e=DD*5X}Sijk0#k)DdNz?2#4?cGnO`~-rckL5*n2bW;shxc zWky3H+&`n3)$zrsx)k15BL)RQlL?~uv!0%z1Uedi(6%G^N$q=o+*DudHWiL&oHy)p zCN*`0Lf)l!-Ir5(9dTVJ6H{@TSfe=`#Kx5E@*mfHDF2dE-9HRO6w1DAR`kJ;EEV@=k;y!Hp1+wyaN~*H^@RdJF zSSD_r(pjCJjX$IXG08S#_-ZUdQA=eTkx1Dh zBf(90VsxwaJv;JC+RfOK)~f6Ni^jeWQwD;KLux?zrVTnR;_W<`hP}!UP+Jw}fPQHz zIA~&tE&9Znn&;o!smZ=9wQZ68_MkR>o!-y17f*|Qq}O^46%xTZ#H%iR@UvfGSa;GE zx^FzubkhG31R9{D5E|~z0fAzqB!!0iJ3%0o%s>SzmRPOZ>MOd40WTb@>`Kk2>cZm` zpCVWbSZ@W>?Q3@r1elKZ1U8^P88d@>q?B2cSl3lix6h`8C$REQ$(M5@rNKJo-{$G5 zyIdZoIEODti`|J|U39|l*iF1V*X^YwVzji$Qi~Ccn-oDOiX49B7{n)LR19QX>TC1D zg-(e)ia7SVc}eqb5T(`PtY}CDElP}mXylhM8`-|ufOG?#0PriibtI*037ycY_JeVw zfA!<_uM1;Z8cu*d_Uur=I%CTKC18$oB99^O0@2&T&MgJs`hSSfXumQ`4Na1}Cvf#d zkTBbb;0+=UHA_&x5iA{y*LM6eZ%2@Jzo%u?*I8XOoybFSRA+6dAvOVEEE^1AL$%Gx z=1T$$A^P8dgy;UtiOx!&`Nc{-rSAl^#}QGB@p+0oB*t4GGgDC@P1pBlLxEilx9B5~ zL7P_Ps=2=LFeWJw*sSMEgdVdJfj6jFK|A9p*VSlSFM>TSTq3owDuOrdPgR6q8mCZ% za2l*JIlMYTBCWf6r|TP9&Em9-+8_U}g3Xraz0DZ0uUz9j03Yh1Y&sd)oWwJm`f-P&V>~N7H;Z z%voQ@dY4L1#Wkzp-he4==karoIY{)Uk_<{jogH1fPA73y8VA#Clg7F0-;BDta+S+=rmiy8_VbWH? zFPj*q6|-~&q?SLCY>rb~eGkby@^E`F!7Ry2^N|>liNlH-eRH+gBPpe3qa zb;6gd{y2DIGnU`zz%!D;e9I0b@wCYdC}B;MVE)na#H%_&Ja}M*G5VMFMq`fiv-7jy zi9W4k?$gx~AJMe>q^ab+@p>iZi{LH`QT?^YInut|`FCNRkMALpUC;$M@ z(6R+ET$*j3rf8)vsLfHZ|AI5$SbfE)|3zY2_AXS$Kp2kKKBOwLl`1OaH3iQjS}=FJ zBuLgk>xJ4r(J!P?0C*a8`N40w(#_L_3PPqE?0LdVSAW7BTW!D8_N98Lp<3H6F*#tt zUPCUZRKQhMf^fS8boHyH+5@%-aGYdhP~gc{?rxn?UHzJ`&hc6t2H;ZdSf_2wpnc@m z(vpezH9l+Pm7e+YOIHoeP!|ulR<20~ogJ%|X~5e&gE9p6Ui!BSgRXLuZU5Ef9(zU= zyY4ZACDT|_C!cEcQI+I&j4w=7vWzC_S26>LLE-(2E(UyA89T;T1}8(Er3{7 zCP^P)_vz%L)|rdhTVc5jbU`*(oBt+m;mJA2?}&}X9t+=ZNZZs(@nPXzmyGWGxh4+t zB$B%k>gda=QZ+<4`5x}u#AhE zkDQM;x92Xm8}%6%2hDL~W&Ay`ma~R3GvjLG&&xS7*}GPvqeMWD8po5@XjTY?F$H_O zY71--KyLEnHzU3UzvN*H3JqFBqpRT0MZ2L=99Q@5- zBj;{5S^V?Sthw0_@ULrb`JoxLx9d|r-z2%;G>pqddYX5@J?0nSZhLP;tV;hx*$_!wDhhkf&RmZ8<|Wo6LCh za@|*GB^QB@3KgK;uc)l}yT zu+yC1y}KbzP9?=B;gd@gkGQRV5S~PeAqq{A%~q6dk|l3+9wGQ2mBDhu=WOSrx+t3E zSAI!U_NAw;N8?y!fT=$Tk&v85cYaSmY1iBHvgb#Kc1Qy&V6d$MrL#@BMQ@?HI?dyd zQAdFY!;nKq0b4lbU~zK-_*H$$Nt^|=p7+H4hw9R|$vzr^kWDKALOhbXF0;XUu30L7 zacfa9vB^Kw7XQ&dNJsZ}rWA*xe~Z3Sgf3OxLOTB-d)-H6QO)Tny?VYAcSRqoUf_&7 zoyva<)s2lMBv+8R{9QPz!~%@8cjsltNxAV?`CdtSKn~W$ene|Zv9CkOuGqwel#=Y^ zJSmeU=nYfzxNT}j4ZCwPK06k~Bb6lEr}&8{D(Iq0rvFm@9@ccNJ*4Sg*(NIT)(lr` z+^T!B{Y_GFxBn!TR)u5mrFQ@w0Q|1(zF0ogdTCZBbE9qJb*_{Val6G_qq;YbmB^k; zmZE)uqKi;+7hQq@P~^D0+pzqcPCSF0@bQC?!~NQAIAcSK}FKH`Q2z%2=j;~kD? z_<-PY-pJa@>n|A?=OtAv})R9%rv4=IV+0tLNZzIwDP>U86A*WgPDG_fN3fJOFm!V*?*}ofEPeV(L%a)yGGE~i=f4h0?`t^yag;Nry_mljWB+HXMwZ4vZZ~Mosyuk- z{jr@sM;L|jAMz}JJ*7iS*n(RB6jY+5=*=j~%JJ!!PPU(-@jPszz}EBXjfb0Sg;tpX zVomsl@(eAmIt-UUMGRTsQ!ivj8UI}jwbL|w+bDWltUOJqp|wl?mYeQy`}ToVTfdE>D1rR z6ykx`f1JFY31v+g0XZ@$DL>yj8KV2tVYT0yH>^`=io9Dl|Amvqqso3#J%7>*1G~^7 z`dEbizd?rA3Jl$RaUN@XxAmB~0Mc1<9W0-c;emGN^gwTMklD%z6!^crTSD6tyY3P- z>QC!9YD|kX`!#-dJ#Kz_aZ%(}{pPNz!KG1ma(ey&fEv@%^|RH~SpIQmMlRH)l3u@+5A`i%=?im z8dM*b&mJtkaI5;0%ce2HneeW?N`iuC^P)Tcb+(Ib+ zGn)#gko)3qm96Bd>IUm54k+hq>toLq*EtS7)jRZ)N_EQO zVUZdxVkNG?*e&64_5ApIiO!)%d|?EEzB~_gKSgp1xk1m7i=c?K_|?#S&XKznBVHC$ zh`)DfNP^KuWwV63>BqIvXtPH>KlsQQYu0+H%h-G-J-o(K`vepv8JwiQZ?7D8bOtso zdMmuyk0_9VoSa|*>cIV`@W_x%emnnT`zp)Wch0) z5W%ds`|j#lq$Z85LGH#r%%?n1>#XF~q@ypm4qRi_<4x^Za(SH_skO>HR<0y>xN(Cc zBzx3Z{hcP6Oa>zORPf(=>im@CKSJ=CsG!h)d17M!U3dHIBR+*w5)k{p)#75$ME=VY ue+KzK9#~NDnecykVDSIFI#>|=|6w6u@c**$P_F01&jd->*j{QWlKc<$1k)4% delta 6544 zcmZvfc|6qL_s2s;3__96d@R{x3_=*$$=LVEItJO7EE(RFtwHweOvtV*Q?|w)Gcl-S z-wBnyu`}aW-|yq``{Vol-apQHo##F0-2cwKx7ak#G+CZjPD&1PTShKL~MMFQ;Xfqx$NkQ%TOko=H4E*MtBo8f_@D-ajEpHM_S z=eo`!TIhGT&iu_MD{k0Kc8ip3Y?yC^r6U-*e1H)HUKAeFRaUFc!QrQy6EXqvS$zOl!haZ6;(lbJB+Hz|lt&L@wFC z5p`1A$Qwgp4I?F9JJwr|?!V0OUT&EDs{P|EE1l*0#=wu-Mqy%0mHW5rIgaDkUzFwV zbna_y*3gkP_thGW6n8l$7-pXAJJ3wz0*1US-ky@g%kLyrft1#ln9Nu0gxSF%J*l)- zssiBhCc>w@?DcbeQ)nn+VBZv3XmM-eX!#Ock;AtzTedM}3%SK6s;X&jd9?(&daP1giSl1+O zF_SOsvud442Xc>lKn5DoGy`&x?rr+q?rofKlGkGTgk4XlYZ$G>L&*hK4G)bwbsDEt z9c6hsRvT-pYMzocHQV?Tim6aIFW)lB45VJHxqR9P3G8z>jGOZuI%PgjJ+wJ3 zk*5&xYW$Qs#ud_Ni5n_48)duCXjPUq38%l+1QmGG@{)B8 z2m(TgugYqIiU4m~e)1;+!gq55q(Chu`pXq)LRqq7p-eQ0p+p3c$_5orFe99kupm5{ z8-J5bg|gzM=_t zcJW5PntBvGnPSY#=3sE2munPLZYX( z4}p}OJja}MB(5w(GziwAvg)v+=LQ8MU<)cj(sOrRqGty2G+>VvvCqg_lERYL0Lt*k zOAxrj3%ktFq>n{WFUi)RhP(T;W0%B!fxBIMc6jE1-;O~jU9#iRJB2o7Zm=bOtY;^5 z!u#S|Icl3F$$FB|#8ED45F83Dle2+?;92lv;3v7_pm7!)25g)*2OJoHLMm$MF9}KK z(fP5uut0p?ModS{Cokr@I`&}$i7hjdy~} zVPL?-vWj(#p&VdsOG!gi_+bZ@RdjDN%iV&an*350x(DJxVCY@MVe%!u?%0IOoG?TT z!S^MvYWKAi5%rk(1XWJE7erM~Z94kdJ9^ z3P=u`gb0B`ja&|s(SKXvM5KrMjU-K56ZW77@k0PmFn|d8)q;R-8pU8zQa5c!L;{MZ`|s&C*Q9HtB)$OI3NE>2q*5EYvxQ!E$s)X#(bLuCE99Jb=O!Yu zlnhNVUKUZ}E7Q8<&Sh!NsMq&kcT>gm4ZS74j)7kDTvLSZ!eI%gypSr@F@t9fh5Wj=g7x0>O# zqtXLJxsp8sw-{0-D4pqt#%Arzwjk~xVBFTdSRw}`9^XzZtJceb&N;Mkj}W!93}=ux zRQW!$mJPaQC5h`I+LlmO0{7S=OxwzH0FkS(PpQpDZB%XY;^3Q~0)_20KUE83G#8RS z1<$xZnedt7Uj#Ai_(t&^L9i=+jMx}6^kM0w`zZ3)siq(0bskA>794tM<7UH7@WfF! zYS5k%2j#r}egk6G7ZWp2Z=meWFxn7zl)%~W^VX+?LiAl!3}oQxK^e1Hui61OfJePd z56L3Lmqtc62C6fawGs`U8LCbgZm5TJ8Nlsxzw-U-XZ5?;)b}P*oY(N@24B>27r;*urSLc^u>i2mt<9-?* zbJ2ahI7Z7oFCTopv`*{B^)f#zZ<)(`nNjX!CCiB)pkA_fk+j6uYafL%HtyJ+F;+o( z)@}wi$H&!%g3sgx4*((fd44JG7&O$%JG&e~sqfP$xfDRhX1-I}r(r(&B^`avk}KaF z&GbjCc0*X^O&HmzkSlCHRz0Y_j4TmQZ!Vui?-X^ZqEj)p z5`ObxHPYCyX!P<_6w7>frSeCYKRd``>A!wK(+gQ3QmIw-mdm57T<4i~k@QC|D2tN> zVfc|UXU+9|bA5+TEe=;#crFA)rye@zH$PO(tCpr&iOQ&!=9frAU$znCvlV1xX;>p> z!3WdObU^MS4K^mCE>?1FqKxS$%mz#!|nZ$3P#Mib5t*J%2SNL#G=%1U=^ zSkvQ1q?z(8!$z?GveL|X7JoD|H%C}3^qRcHK9h@-uhlwFA7Z|CTrpjr<7RQ6cAaSP z(*wsE-6w|CeN`UsGBSSz3pZ}2Ki~i;X3Ihsu$7=M*!#=Rzf)NH}AWzso70TKBs=OG_j8-cF6q$jCAq| zFksOdPY+`S(&?`8aCihZn8I?-i0(7@YPFo~{e$idZ8`ed*SFrDgmJw%ahuT~xg?fl zD?=upI45d%c_mu*8D!}s)K+2B2%ZSBxNOO&0aVdCx>wwaA-d}CE3LG6kdGW9Ila9s}OiM_uW)%#VVFY(`30*xP z%CHDw5%b$U4{ZdacfTb|&w2?2Qxr%mMhRs{b;erHY9B#$A9R=n&1g8i+iVx^7FPKE zX>3NMzx(Aw^D=mRRTaT*qNOFlQ#!J{{TkpR9pA0O2>_OPD(|xUZT>z$U{DC<4rXW7 zzW2^VGoO9)x%JpYaGVhBAf&I`Vcg`I0lb2C&0T6&K4d}Du2?rn->#n44w@#>VJ@`J z%YE)C`%8jj&r0^N8?@~n!Q(= zh3!RN&bWX-+(E8R_xyWZ!;6_Cc4=NByN9h?0w8i5Sy$?W7fTA3!sdGMPys(FV3yEk$sjdaP9`%w^R>g)~3$~QNkK1-pzm7&pb*mou)q*XbLmFU9 z5~s35={@1!s==v{D1#^JTav9Rdz8@=XK==0?agYE*N>-zYb^ZYse7Mc^d)YgMoyyD zHf+yXEC~(lbYpy%>JvK82;;iL*RdlZ_m+^N; zlXs^ToYQG~zYZjv(Z)Qv;`u8uLB9b&(dR>Yl~YJm5PzdWyO2qZ1eJ7Z@Plyg604fQ-3yy7$F@B+c03g_e`+wg$R8tC3p0ZcHiEHp zKV5=nx%{TQgKVXh=x-^rD?L9K8K*4L1`uag2KL$@84mQs6KESC4%ATKTQN3}`ev@B z6Yc$hS%)y!_8)9~BE?7Ndk>PxdoqlzHm}NvAEZV}JiU%?ByJhl`KJ`(wtso`y%A?F zyY`dA#brDeUojI(n{_rn_XE+j2DVq4S!`220v}d+8rp@c81SPG$*fn+hRJqqz`-!k z9(J_zetN4X*RTlyv+oMD&BbrQl%g$~z8qi>*SEw=Y{;7jzGW91j|cWQ$@??$%)YBN zPqZ<(h&{(?Yn3Ujn;uqf62sr0%P;v?!Fvjazv3KDpKKmpIU=jR?;{-~fcHDqeT4A= z#rbl(Of_Ss0*z;nfpGG%BITKnotGrc8w{YrTf>uREdUq}yQ@n3k|8hrbG)3-oINmY z*Vbjwzr6*v$3+C=Nni76*}g(Odm{~x&L;y+p-*B9Vt?0?6k%R^<40*qP|ZgbJX6bX8$AVX3h`piYD~Sd+tOBH-E5 za-Ufu2#9CT@t-IBOkuI?nlhTmkIFYs9C=y(&O(Im`>U(5YcS`@22y=~6fV5F%0z>E z`#C+E)Dlevp8?;m=-e3#R`ZzzhFAA(A^nJ|plWMBoug{YGGyD8|@qcDt zKPP`=(uf_q=S6UQYib?z&jh^J+dtnhF4p;^GrZ0c>V5H`$L`a(K9PECzWz_TYUJ8O zN!KfkZ>w6xI~ZXWVG70W+jZk&oW_w;PBjQKM(;m;Qe>MX_wQ*1I}ILQw`C6-NFnS3 z_&2nCD#zc;KO|{UtZDhdj4!|QtnyXi&wdzaQMmD0J81m6ZST`|(CQIklK4Dn5n9IY zfwhPx2EiCY|NbUXtiRn{<)>s9ax$D#^c6ABlqNX`gA`~&aXjk75|{kE#>f`MGeP(t z`4bd({-22fDKsyZ43*d?>`o6&Q-I7WFd$ste3~|>VOIm3rV!TlT$%2Euu=I$F(Wbh z%$3vU2=Q1jFKaTV6eE?R+8!85@+sQxqFC;VSxAn9yJ`s(y|c`;aO9Q;A1)Y}~0fwFdLm2Su*R=K-E>MH6f0 zk#ZNL#t>mJyt0sgN6VQ-eg0W}McpP^FT1c4XH9zYOTEI}aq{fkkjl{1(#GD_(%jfy>{bmI`RB@Up>(W=;=qWnn%LwGLI$dIeboH971Z5y@t zKQ55iQdr=l!Z*y}K2g2Wfr+CSP<7!!hhv?kdT>V&RxEnO>-Y|6CRgEz8NncTYcyLC3(zPan;+);-LryTwcCXni#1T#aMH?SCw@ z-!Kv!kRCLwc@vYD9VD{f0F3IH(nla`iElE9X2@9bvYMxQ2E&&A{lqGe@g%qW8>p6h z@I0>XpOOu9va???N}Ows34>v;%i-U2!= z`FpGiWCo`pW9GVvkdXxzJG{dkrOTqsMUj$S@3fUHm$PI`74EuZ0DL#SQi`FraL#KFQ#tpv&KcrmGog%{pvqC8Rc+3nJ%^Q`vtf*7mCI0QdI}>ST#7esK!z+#b$^NKk38HT5NY$@%1z9$ZyIsh5%|Y&dlA<7ybx^1ekknEHG2HjU9lQvSW$;E_ z?r~PMbu2SmB9k{BTkWW0HcvbpPEK>1D-=c3rKkc8z4myn_O;HL+$Uq1yfY`Q^Uw`^ z*4siBi~G(Z>d3?LJJTI!4@!rOOD2cQ@nmkzN-3k$+i?aL=I+0~zT3<-9*FWC>(Em$7Zr@h8 z{a;L8?!WUCAh*+A2`d6p|Cy~IFDw6FOhH!mKbRB*A}jY_ObYU!r%OSk{zokoBJ=Os TwEH4=FUd;N-@K`5s73#ODmuwV diff --git a/scratch_figures/Profiles_basic_model.pdf b/scratch_figures/Profiles_basic_model.pdf index 15c55be1a3bfcb43f9084140ccda79f1d3d440df..e2a0577b3758c0828c92f68f261b6c8ee52a44a6 100644 GIT binary patch delta 31 gcmX^1nDN+S#t9bemWIX#hNh+y?L}eCjnUbK0Icc>$N&HU delta 31 gcmX^1nDN+S#t9be=7wgb7KRoR?L}eCjnUbK0InDb-~a#s diff --git a/scratch_figures/pI_profiles.pdf b/scratch_figures/pI_profiles.pdf index aa9093f48ac2ff6005a5d7ca803602af04fc80bd..b2f7e29065f9e821cc1aed7b6c7c40e77ff08d14 100644 GIT binary patch delta 29 ecmaE4_sDL71-qr8v4Nqf=|p=`7;|GZyF37i?g#z= delta 29 ecmaE4_sDL71-rSSnW=@L#YB5i7;|GZyF37k90(Qw