From 8dbd7441b6abff015e2bf849090e6b3de15a69e2 Mon Sep 17 00:00:00 2001 From: Arthur Goldberg Date: Wed, 2 Dec 2020 15:35:39 -0500 Subject: [PATCH] revise for models with only irreversible reactions, as made when reactions were split by SplitReversibleReactionsTransform in PrepForWcSimTransform; enables great simplification of determine_bounds(); simplify by not bounding individual reactions to stop species populations from going negative in the next time step -- just use separate constraints that prevent this; ensure that all reactions are irreversible; unify treatment of pseudo-reactions; set correct reversible attributes for rxns in dfba_test_model.xlsx; skip exception from Validator errors until https://github.com/KarrLab/wc_lang/issues/142 is resolved for fixtures/dfba_test_model.xlsx; simplify by making/using dict. of all wc_lang.DfbaObjReactions; stop warning when exchange reactions don't have the form 's ->' --- tests/perf_results/wc_sim_performance_log.txt | 12 + ...Solutions to test dFBA models, by hand.txt | 2 +- tests/submodels/fixtures/dfba_test_model.xlsx | Bin 60862 -> 60981 bytes tests/submodels/test_dfba.py | 347 +++++++++--------- wc_sim/multialgorithm_simulation.py | 5 + wc_sim/submodels/dfba.py | 263 +++++-------- wc_sim/submodels/dynamic_submodel.py | 2 +- 7 files changed, 275 insertions(+), 356 deletions(-) diff --git a/tests/perf_results/wc_sim_performance_log.txt b/tests/perf_results/wc_sim_performance_log.txt index 4efd27f3..01478e8e 100644 --- a/tests/perf_results/wc_sim_performance_log.txt +++ b/tests/perf_results/wc_sim_performance_log.txt @@ -322,3 +322,15 @@ Performance summary for Next Reaction Method on 2020-10-27: 8 2664 6.114 435.739 32 10643 25.212 422.138 +Performance summary for Stochastic Simulation Algorithm on 2020-11-28: +# SSA submodels # events run time (s) reactions/s +2 698 2.830 246.618 +8 2658 10.589 251.008 +32 10766 50.124 214.786 + +Performance summary for Next Reaction Method on 2020-11-28: +# NRM submodels # events run time (s) reactions/s +2 648 1.842 351.708 +8 2622 6.307 415.705 +32 10671 27.047 394.532 + diff --git a/tests/submodels/fixtures/Solutions to test dFBA models, by hand.txt b/tests/submodels/fixtures/Solutions to test dFBA models, by hand.txt index 1755defb..b454a4a6 100644 --- a/tests/submodels/fixtures/Solutions to test dFBA models, by hand.txt +++ b/tests/submodels/fixtures/Solutions to test dFBA models, by hand.txt @@ -155,7 +155,7 @@ VI: Modify V by scaling bounds by 5 Scale the stoichiometric coefficients of the objective function's terms by 10 and scale the bounds by 5 This should produce the same solution as in II. - +OTHER IDEAS: 0 <= ex_m3 <= 0 ex_m3 = 0 diff --git a/tests/submodels/fixtures/dfba_test_model.xlsx b/tests/submodels/fixtures/dfba_test_model.xlsx index 5180b36c7bea07e9c860455722bd4abbf9765748..4b861182ff34f5be8a812b2ae244b15de4a7af35 100644 GIT binary patch delta 13149 zcmZv@18`=+)94-BwzYAdB)hR~+qScDp4c`vcCs-xwr$(CdH28f-uM09Ti>ZVRln)( z>h7tjIW^POr*jRwaUHzo4HgJtkMd*)2?DZ(4g!J<0s`W0!{lb~WNl<`Z_Vg#Yx7$} z+qvL7hVPZk7q}r$y*v_BpuMrzQ0<=PpOp=^)QJLdsFhOqFe$R1#dBYqsE0TZO zm(d(1-59qxxtKrP;lM9&S>-I)<|f9nt9FKqh{BOrKW+>+2v6{^X>NfE0^~e%`pirG zlfVyC&%*P!nAs~lny_b?r{FN^kRDy)&@la39B~J8dCf2go0XnzmcHRVzCkhsgo-76 z98$9QY|FHor{ml=UdNk=R>?5(_-}?B0$}4ueBmBd3@%#*j3Vb=Xl}ZAgSkfHU}Snq zMXgw^i4&yx2p3=uZze#l7Bh`SmC%0Eq`8h}2vspNlL`HL^E9)=TklqC(9e`$U17&} z$Y3dj&wdh2HxG>63^ z0fkUM2g(im>$_8+k>DQwX?{0Ht2Q#fJ#NgmRD3xhVZW{Q?GeU8BBbV8YU+59+&l`f zh5C15CI59cPbdg_XDhzsv!$L)dCI9SRPTM$(gQVa-uU5yq4Z7m&?d;Dy9E6j=o;=~ zx12C)Z<@-y4LV?zPw)PAWICUVn1ii64)#+WqZ!vULQ4S89(04WTbzP$kNg<)xx`l^ zP0_@_09H0qjo4F0CMCfsiH{LD4~st=Iv9M@jR7KDMQ@I&ck&U#frrLYHc)^XPF4H( zif*c}_3_>HUIr)Hq=!~1VYEfHrb=?-C|uoygd-T?w>1zTtvMtkLBS59AQp`C`%Ho5 zC+uN{dw5=3attCSr$!fLiG)BseVl=A;t-pzQb?|X=f^+%2dzlF!bJ-;lM}_i!E8>MV zbzBS}$eOeOU-{))qpFrJ+rWjMPgROG@(!9MceUP|5FAKa^bSd2}<*5S@;d&@PTxW|cCx)^# zXqdX$HLNM41gEaC#}A#s!8QmJjIrWP)-(Xee!}P7h94Pn%AMRH%TB!XiJH%x$a9Uw zu&3d{P09z3UaKlj?MUVQmOLspZVZBOtRiEmCpnS;HUvXFyue9VUU|R}bIVElF5E)Z zD+4F52ri^cwvaN;4OEXdt@J%KfYykuSWAb+so-Cmcn;B=)vG=u3r+btW0sB(hZ>nN z8**ZKOzP;bf%Tm(^K5H^p@({+>SjT+3n^NDJW_Zw5zZ8$I)4`6_iLZ1D5l2l%Gk?$ zd-Q_o7U;JAi-- z!zI0{s*8%+IsJR_uOwoy^PD%bIeh%R_rURy2I`j|49o-}fFuYb^n7r=0RyEiJfymo z4G?F#CN8zFRn3RCl=@3I>`B*ZYou2-H-l=1Aw8|_R@G&^>`im4Fr98k$JacYBz)X# zD)__K%;YAdY0%CNlRzi4_;=Yp3&?%n@h`~wpU4=Z!09gul^bXf5dLosi{i@_ru=yxI7Hv4qC}5F>Or9URKspU)pV+dv=JZI8{D z-5gF`5t)SBB+DNknB&u%QyP5$jq*+8o4jmZ_xtEM<*G zwZqoX1VpQ)BQl$(3y>}2jXV}x=D5PIR3o=g(lI539sqrG|6O*xXV0WhG^3vc6zF_P ztKaWQ-0(&p%}Sf37kF3&z;aM1$EvlSmx=+HR`q5vOYVAbw)`x7K4n;D!8;WvSz^zTh^ z_K5?=&LV){RLrj=qwFHN%1kQ>us{KX@h8Ko#c!R5Hi*?&%Jx3H2k<=k5qMLK5(PYR zuM}aWNprr2!#)XqLFL5Pa@-!SC8WM;4sdwCnPH-V&EBztzSwNAix$rbw#``V8Ko%I zV%E;S;sr!nY)33}!x71}>EE#$w%LN7`cv1942gJ5&sZPL3>G5tCg`UlDS?TWkUL&B zf$t@&!v2Phycqngb)nK3-c1&h4_I8amsrn>JdGi(yEc=%QJ~p}VA^ihlXlAbL*_d? z71TOy@fm%Cl`i*j@@(Flg<#fK?tMM?;NuTi_uv%8k!-tw@6n{~;2699T~W8%KtWO! z*fe=`ybZ}1f%_LDYt9udC!lTD8YUcAj04BST!*#?yX7(1ov^_Dx;)HJ6ozy&1iKb$kOhCx zlS^X;4SXy?hg1yAInz&6kw|gJ-VjkHB6JOgN_@>s{TgiC7@crX4Cvb|%8VMC3$E!g zk`ntYK=vDA+jzwnbDG~T3V`)1(1*QrCWek1@JUB`+YJ}F#5P)MnOX)~m)J+6JqIa# zQ}0G#u7hDOXf9^}vCK9X2D!>LjJQ|9z6oYTS6kQJbZ}6+Kc+jb&>kn55?Xcz!#5^s zd}nlU-A@oRbIGnl3v?O6vBfL+;BsIr(!^1lP4;Q9K`N_)9Q=6_%?Al21xn)!* z&{VYSZAyLcn<;_F;c8+aUZbJ@J{lf;`%-Rg+yte1k(AkAt#+I;U!zS)P5`k@;LyA% zKwTJk>Zn;*fQp5>3SmdtFtSDxW$QV}rg1?vfzF7z7S+3c1pMaai<{Uo^=bWPu>j6h z7rLiiY9ndg`e@UomP&KF_JcCg?otoRaoH8xL?e{Wy-q8r-F;T;9*nV*adbZ4B7#lV z6cVK<5mJRnm>_crH;e&Bm1i5b5u=G`c-NH^O)AJ@qll40Q`zYXftg}b(fu%MMthB; z?V&c2(4Y3a1=v9~G>H!tUx7L73Z8Ni?=Ip?c9SMkvQ;^6-(QDht6yD_N&@DO#FXOA z|0zg(#AmsbboNolfbHE$KD+8ApSP;5mJVlnT5ceoqM`l0E`yvg_owxUM!!W(ZqE@u ztpe%*`{xRB>Fodv+MH9KJ|xnZr*Et-%4 z8>5wO)%}NTfZ6;ww=WOGsXAj=R*_8$#PZ?;I2}x=!xpHE6x8MJWI)OEae}Yz*w7A1 zlBBBw2oSLXo90mkjXf*BVrDwt!CWDOBmw9M+e9Z^<|)M9y}Cw5c(ad^&VUMD-O82W zbN=bWMwpec_l6~Jo(#Xy6?}?`PA1Aupocbo%ehcCjf?m0JsaVz;tnZH0~8=0m6n9IfASNH#p-7?5a<7&I1N@*Yc=*OJ?u z&3SZri~FsSpNa3ln-GrCnu3!J-`c8lRanE<60dIMCV)yy9I6O7~7hO4{eE8Nx&aQS^CymiC6x!aRq!p^do6!rDje za6TTerTjcBebf>*g3hR;Y|eu~Zl?be$f(A)qkpc&6C&R`Z7#7X-sCrl+u37?b2_U1 zblLPSP$U>|>zrJr7kW8m@Wjj%68@VnPl!kLPOUVQUO;O#?Vx6Z`xh2*;Yvo%M3igZ zWR&*+g;B0VI$AN5>XwB8l|H!JnmGwt55LwH#med=>_f z|IP;U>e}*a%ox5o4PPiP9TW_~jjCuO@X2w@K_b)_3zEO)dsBxU^GsKPTeVqK2THkn zbqq{5j^}$<-h7{O1Vjv+cU9r#r5LEEmypBcWXZUgqCuTg#aW=Yw9%p==@=A3bU>xN zg6F%EnprwA#-kt$WSuZtA?M1m_$&33Uk*AKBMC6i3DN_#P|BFTsXCJcQhG+R=%b~2n)<_gC&#~0ly-tZ=S}~_&NbRw^EgadmYYyI zcR*iyrysFVA}b<&j7!8g~Usdt0v$I=t%AR=e@ze&J-OrHawTo1m&?k zdI9R zhvz$ZML>gYSPkoLg@x7qk(Ig1T-JCp$+>lp|DL{Qx>;3&h)t6vX~>Gwe~wgT#>FVz z1bDTSbjBLDtTYyb?+M382_;FHygyW$(N42)S+t=8B+JarFLyU{s{zNnxoz>UAHTrb z@q@>3HCwwU_;`acpS6`?@xP2{ui}n0e=;S_Q!G^N_HAAX9Q*W}?u{A)Isi^k+u87@ zfj_4^k*6Yn89(W8kkMWyJVG-bv)fA?v|!&$*STZ4S-$uQBhJ*E|A0 zQ!B2QQ|A9D9EjI9(<(U`aF;ip$)8pZd?)tY9gX zD6Fod*+kne`Xl3x+ytCx-mwVTzR`^w2OwuEV&7|7oss<=ufJ=!AXg#;H0P*aXCm3# z=UD{#X2*0P$jg;7a zrG@%|uPUxRP1t`^(4v?wz6=O-c+EaAZ}2vSV2d$!%+WSxj)Z zxIblWlv5j$&3Ea#uDzS`-4?{jpRTWM&}Z>C)2jNZ7)Q#}ZJ z9g+h%9uS#=q(fEnDvg?sH3zh$&>PyJ(1fuJf%l~>Zee3WR7ABb;4dF#?03d?F||nC z&|dt9^@M?3tOhJ~Hj%v0GMM`~uLiqH$gwsSPs*x&V&mNI@ypyPJBm3mnz+DaRUyHH z?|*13;T%@CSQ#upuYT^wJ8GY(JhX)b(7t&ILi^MY7vy+~A;j8+$s&l!OeY55WM##( zoT07R(eaWv1V_<)^hhIt6bKQro}=A+vj&1^ahx#*Q@-$VC*^Ms)k6}dxfoQSCdf>q)ht9FjP-%!RaxGZH zW(xQa{`&hG>z{%b&t;YK?m)2sw$X{=;_udy?6pBGfi)J0I3cDxq_jZzTRDa~3eqEm zVVM&t4dApO@A&m+wm2~=0z~OnsBeg&>Ou>1Fye$QZbKOBi34L^Z38Cp_>N)M%z6cn zA6^uUQ=EW+qF`Q?$j%^UT{w$RxD$@z2(CmDdbks<;v~hWLjEu5B=~TrZ;DkE zqm3dh@ic$X4j5(2{G@>B9RapB0@c5ywf%SS>GlQkT_~mG^Y1A~^F%lkXog6OVCyA~ zRNFKw@PgttDB}_r;^r3OSS4ExV_fv>X0wj2F-fp1NK-3F)e~Hn&<@1%*QJu$*eH$V zysLcP*m_#Jf2Gc)W|vz2r^_(%J5`VWTY4^5@|`G`TtJDt;T-cZ^MyMvuzF2EItq- zVOduWNIG-rM3TR8ARy9~Sk5>4jbc?7Nu;3-r><$E2tMp*s&bBn8z+h*S<3&9SBy*6 zkdPkUN`ngASM5A*|Dz}(t?%S|`=*P=Omn20K@O=C;!s)6Ii*lULPCf;L+n{sy?%|? zh8lA(P+?0g1`vD^~TF12qRpg zYf&r?({ViQ>u^~X;k9Shy5{!DsMmfGmnfUw+(Lb!%XQoyAb0=?j04NgOfJz)SN?^f zLTM41aRjeJ0B3}!!fa?YL8mPQOe!F*n?_LJx% z5DBt}(7|yg3?!!p_P{h3)i5@wxfD+#9#w%Rbo8Juz>J9PR}mKz#XF%V)8dzy7lWGQ zvX^3WLBUBD5DJG2*<-?1XCpM6%_F9)@dImI8np$gkD8lak6VUh_0UkP=j3S6jeFoS zc&RJL7+;8|NtPJd=KND)aB>WKp7W;!pis{+PKTCnR)>!d1{(oPqjb@tFeJtz)G#Zb zvMO(goi1|zXYNUjOYuCEQW55#>=2tVb4cZ7F$@O-pkk&ReV`k`_^;qrKiu&iZoc%M zAc#Dn{)6566=Gcx-9Zl?^0X2H;*T?5xs@fw-X|{# zJAKk9I-h<91(n))TE#WuPR@8#SN4uf7Bb2%V`+~p720#YgFW2cgZq`OcHBv_ssXAH zapqe?U^IhR)J2Oa7Vu{n>n({313BD;vc;Kp^j#kt0I-@}{8T5FNp1M8`bb*IJoWbV z34HYmH4n*f^KQ^=T5@E%De0O4KHVJlEKtn`|(Ii)4Vc%{^oN9_36OGIzq z!$%w{t(!;41F>c|?8i9+cs@~1ErIPb>fD3{Vvp>Oke8He6c|z=Ki*cCQYek&6%(7O z+)5u|%ZBEZm-|?r3pJb6KE2~;r;k~ABR%fG(IbmiG9Whc#~$RnI)Xj3*of+x9O zgRF&AhKFWVPL2H{9F~T*dHGH~56e_9#BFJ4RN1bSA(4I-fqHT!F-Izqu=juo5i;VCQQ494|$)`+c-6oMQ#q5G_5eR=k4QN+BM* z)Z@xR9h{zycdzgqbh@NX3ET*XE9y>Z75>ii821|$$NqjgO0ic>z_sa4|DX1=JDiR- zGC8&n6V?d#D*>^Z;m#*0kV$S-`vD{pwEYgOFEG@6f9LX(l9dJi?k}BBAmf)eaOlDL z(T{)RSB^DJFbOT43t1`khHhfd@iPC^5}A{n!a(AJpgL0zl7w%l8NwtCQ05lnSdx z%@?9;Dq6Iy8My6hOMQ*a86!@kF`Oz+9P$|& zO&G74B=y_%3EnfJU@4*uOuZa2biU6*)+hIo>0X_Q58T3$6Rum{gbR>LLp<3DcIrmyfy#P_G7&h z{*vJGs%PR`$96lJ4PAhA^lM_h9`8Qvq-&ua+o5Bu98c%dp*P4kcKGk6HHB{u={J-o zvJz3mBZahcNZK@~3OT^1^@1w^%7B`S z4^~pD%`fKss7Xj4Y>?gxirON+TXfzLaBLyjFIEVL zzR(HC=uREPvXqK>58*L^VXZ17@A)$0HP!dT5M#J+xe8&`a#n*?rQnpuwCGwd`MG!k zfX*yng50nt?FFug)?iiECDbf}T`l)~qaRc`xOdeA+RH|1*GiA?^#;p{Ma1WjSRYwfLlXQ(PHHtZF|6M)nrq$geD1qaB(ura7;ak<)aQ$(#B%xL-l72yI6> z%sZd?W8jcxcuaD>dD9@25CsFoC4p+}s9Jaj;9%W?DKPpVgX^eN3K?=5^^=z893CniKYX%79^64m-ux8dt9D z`H3O*!Xu~(!rC~eK?`3?daIA42ZbhTk5GZ?B?gQL5n)1QN1|kDR3nV}v2c?YuwkHp zffdpz7|{nw1WCFi?vDHbjZefIwF6ReEl0Cvo*_ye=?D_hCBTb1l}S)}Nc#KMO8UM? zQ7V|c7L}fua^ZeZR${!4#o-UQ6VMJ9)Tzl@iX}gEM$fo2@T$iRMeGa3>(C}h>M7+C zx9t!eW77aWW)(cs<-~xJ6(JZafdCWESIm7yO*+d)F1i^O)4zjSvK9pjU}}U(LMq1b z(2NoCJv>IbuVuyYe@8nwR*w@mtt>RX_$?xyr~f|_{I&Eb&?OLmY(HWREKZ!!pDib4 z&sj;Ll;h;3V3I zqUp24-%G465ya`UjgCkEjJA->Y_fw1#0Y=yAcUf6uz(mym_9_13=l6ko|{GTg$)J% zv0f8?z3@2iJX-0IO=zcI`QXr~eDR_o8=49d_Rqu~B=m7=Ikw=D=DxIG_LdEg6-kr% zR;#!UYbf;aD!AnO#`Qr<3!4cOsu&a;1Xr&>N%CNp$RW)wVitQllL&Iw-Jv^IE~Z8b z$G~H*LUwG?Cvg29UAZz(B2(douoMndT3D@*R*pI^OySL%kKp^z^Lekon5(`t9SN4q z(H{#>%a%+oHP_N8?VK0R#=jxFWAa@1jn{0bbhRqT4ObjAF~V=R6!pc?n1w(7q(yYiXtPandHg@}wKE=q#0w4Ul@;@a(qY>djkc<{pbv~bX^YI5PV zAJV7s6kp08ANk4txPovdP1A@Sd8+dJs_t+3axC%e%d**021(XhBWeWkB{|Ek41-7U z8+jMm-ADc~<3`C7+C(38J9IcuiHg{D0FnW4vEgB z$F^>(y$c0#^c@e1;|d}di`oe~2uoh74|+}CA4rG!l{lVejVfG?*E$0%|0o&FvTUrv zpD5p~I`t%1?z67@hVgdPx48_+IGj-oUY{NiYUk0~5#a0tawMG4j7XIfuODZ(ad@Zp zNs;fbH%#UblUA>lg}McRnm*33rw-TngVnT|B=GYuY=`EzBuQ%sqtCytd~nLfhKebb zn(SL>?lrgFt!s}vy%n3$zbT72?ZgSMiS$j{wol(VQ=dn+$)tBMGsM4^kNZ-G(=^q& z(jRDADF>2Nb=1ABRaw)mzC)vN*&xX9KuO^L2^9YF2;Yg`yK}h-?=r9PIoQ5=SwbjILn9sn0tTxIGh-2 zkiCy`wF+chFtE85S`W&%)xGKYb65h+sg1PxKK+qrxYIX{KQuw!?h@&G4LOC0H;8`R zV6WxVJ{f?kbo+~D^kB5Pao-K29W{Q$oFKLFoNw#o>p#_SSiT;RO%iwzkj=#EX)>Vg ze@befe@kk_Y*m%gb?DvSLm=c0=-Mer+t24seqQx0&ZjVgwHs^qfA?6aXxs8-7{Sp7 z1dYd)d%d9T4Q*7a@iS;z(P8uRz7Tnm81B#Jmf@P}3R5^7&m$hnUZhT`7daz}vLF>N z7l2w?O;ODHT75h)VaKw)<1xU5qXIQx@{60o_KI}NL?#2CIcVo7VYf%3Lil2KzW-9s zHp?wHxVpMz(-TM~fl~(b(J^iM;h&@Xw3A+XgKhw<@SN08o3qJK z1kh&38O&F*`P?M;1*!xcMyZ}A@`yOVfbmy#g5WIyInNJ@&oeaY6~JOW_VI?~#z1SI zJlMl#Djjs$0=23ohK@-#>`ND^AfFwpTGTM?A7D#h#TC-XE63W6Hr#_#JDp7SQL41B@L$%jqS zE)dp0%$tgp3@-&HnStb~<8XfB?);bpr>%c>x@d%%L@V);Kf&zH9$pDelwTmA70}wDx3V|eh^a=3%nF(F|OXsBF@u1fK`zue2b4 z@a-lNE+r7?atT={cMTG<7j7j{`tT4(e+k|TMMe*Uc0v-5G|XBEN$OE6+0kqjsw^7S z$_$LBDq178BB1~3*r-fDJ;Owy=%wv$kNJglq;*P^Q-R@P{Bqezo0f=EzNHD{|_?D_Y z9?>7$4r);X_{LwJ(kmZtc=GW#F2H?d;^|&A+=J?|S}+tXn@V0h z+)Y4n7!}>g%|r*yd3ZUT-oX9ngAp1P`?Qr7_=6+|^=}Rw($Ff6hY3G+wFtNu>p3$o zQ!_7enVBFz!kOS(#yr??Vp!O?o{nKonv*=Imfv}%FodzJJ?PzELidaB#B5g9QX=;W zuP;TW`n?i1x;0EpbjbJ&G!BJbOQkDT-EiOiWpRM~FOD~rdk}oDK^(;j`z9+-9{zzj zeVZR}Ue`0-cgMzSeLYJq*X_O{(Y5lkEDwmdBke}Q3z??=B_xH z9&sQE@yCU)=pM2gSRKA4Bj3+l(AcP(m6>M^fY8U}kNLr}@N-3fpu zIn!&>eji?d1fEOi*M|XKX6W4h&a7&u8aDDTrdfpK z*o0psU$HxqN-fqa8~e&;J7!k0ikQe(GL+i#|8)gO9c1d2;4=;^j4)r8epP z(<2Abex~n+VfyOJk!T9pv~`GMFEs@R+@Hy8EMHEEXUIh<6IhYy;d;icCF}eb?5v?0 z5FMO!98@+RiRntdhHl4cNk^tG7^tC8a|OyS`{){d-N-O6Gd%sDGb4<}+Ztd$#?*k) z0i@unaWzdlcDO( z@tZ$Kw_YRCJTTN61&j^vm2EvrX{;oFRh>P$^6g2BU3KKys{zFDc6?FYW)tluge4CKa3I| zET<>Q63XQsN7q$ppBSOZb;=6VIL}!>*%qkra&wF7k9bc`1~TM@u>_mf6fN%b2R1C| zrh?|W-PU|mpQYd9U@xZY zJUUbnFR~BD{xFWY7@P)wv*Cu_(qX6_4_PS9NsrdW%#gi=IeILAvdGI;ac#Qk=_==b`qCAUa?*XxIC+rVNcw8H;hRRTz;DV z46RG%DnD12-bbApHyt~0ay2V~ozp&vvF*aMZ0(q>JZo(1XQS7{)cDNJ6$i&;?B7Xj zyOecYkc+Rzb?I~OkC3kPJ&P_TPEpu*NOft;dUH`K&>xNk6;r9Zd*qu)M_z4Ky5?iE zIve}&N>h0Kq%RkITf-z>ouAL~_FCi`eRjIP3>_mb|NhS6a@NlH6KdAEgd#i?L)~3{P+6KJqP!OEbQ8oL zX96>e^0!ru3Zp2*H33=KZH!*_B_-hB18F5=O}>pZe77Qw-Yy>f!MxHss1G~Z4Cu&8 zm_xf+h2zLufHu4Ef3_t zypR(G&oL7x517DCG5#UOAu|{u!9SEIfJ($XM@e)%!~?4*|0l}-tA+pe7qJpg51GKW zS^i0xkC;J`6Ah2Vz?=jVV+Ao2JCA6=2qpgMzK)o{zRCYX(PK>f|F^^VZwJx;dAt%` zkE!wh-}a*as|5kcPW*Gs1*T!}FF5Q8BiM=AKa@XV0@JYlhbez?#pPeD^}k5#^AE{R zF$w;=1p)*F=f8darI(O@Mi!?`U$-z8=c5ZwO}Fn8#mZT2}P!GGNzS^i%mh?##z u2Y(%EuP0(%@PO_nDqm0${M#*{ARtL#ARws!!>!eBV)Dh`{dDjj>Hh+l-4(U~ delta 13090 zcmZ9zW0WRAvn|}1wr$(CZQFKF+j`ozF>Tv6r)}G|r}@o0-?{hPbLvN}osp3{s#dMa z6&bm+c@wl}1GMfP3V@J?^k57I1hj_&1cU$t1mt1M;BN0^V{C74L+|m^wo3D7QUNF8 zm!H}T0p@tc=OH0vtoC_op-HuSW2f^&!^2!8B3Y|o8)>w>-WH%29`c7#%klGTLB#y$ z{#9N1WSRXTzrkK{2IpqVgLdd>wwB9n)Bzu<%$70 z#CVTA)n`pWV7P*-LAZx4=H`=Q_$yhAUR$wPG&*R&oojOp`^2SG_tNz$L!!Atf=E%8`_U$y(*-*yz<|!8Z z=s+i=x8o(Q<5+^ID~8bhGQ7Z+*}n+>aa+sA>P@1@VlGC1oxT}x^dM-6P85K_ii}vM z68>$ksJsP}D-7M_{h=X@89F(!I!S&TvB|W4&qL}%X9Vb1<;damup}h^2>uzr>VQmS zd>QeI?d#VqR-boA=c;6gyBQN^zz+L6c><0ipQz8??$#7X83AH*HK9-vXjTz9;9Tnq zpVvs8y^$2~x zE=n!aW9xi9jTrp0{TLPUEIkaV<_hUD5)!sRzfgs4U9#Dg%p8cBp5>TySEhfhPFWR3 z#+`I#0FPOYk&1``8Il8xcCmCUG*9NKV;mXDOp$7Ax*TbH*fZ$XrDDHwhB3R+u#Vmg zTL<71rIbxELG6z*J5mXCWxtgQhYnf%hxNX5w1R|SIisUguIw=@iER7fmj$~ucUg`$b^c)q zxqF;@dVji-QLWk;w_J24#XxCFo;G#vVF|#lAg;g{P2wmaAbVH;A%pvnrg|FikR#ZeKHJ-XXFWPq#$}GDU$GVIn{}*Voc|nt-1BK_0NFtr}Vz=>@F;Z5G4RidRWS8t#l{l>o*i@m5GjS~kN%SloR=Pzr$4QloJXP1dh9oAOe2|04^)^_s5^YCwL)`Kk!J zMH*aFwQ*&_L2z=nV6qWzb8Anw;!m%vA97jzxqTE2vXS^Q7Rfa+;&z-4_=k6!datQQ4-qhUwsAzv@)#a{n8>8SYQOlZYO9204Me3qZURWhye;|l;AEQfPiSIfPhf{|LNG- z(#XloRK>-~%Fe?1KXY=EwvPPTKNQSDV#j5z#%y%)F45b zrYW@o)a{A^s01KUbKVl(l`|9oB*r@1wh$J_5?PSX{6b{bSq9g$*YM#>dVtxp5@G%ncqsrJqM00^dw$ z`X@_qV>*>`4)gU+sF3V%fx^7kdC7hOf-ay7`|3$blnv+Q6Axg)hz%#+n>2Uf3xN1~ zhY<1gU}WuFp|Ir=+_9*WbNm?gNravXR`)HK@|A}n|LvQxJ>V%Dl4!9W-cUqwhO1en z(ROoS4zktM27%3|PY`qbHJOB!bWBZ1jc&bE!l^)l-U?CrNRM{ak5a$Sj0Po{#|LoD zovLR;b7AqAkwU48oyGJuW*kcU{Gu-JF({Nnvc_B+-HNuv~qY^C5U}_1s zPZ~CH7Gd{L@Nqv|;}g{nnYvG?6Sg5fR6Y{}6}AMqr9 zBaJRsKJYUd^(7V*UPACF&F<=GgcP9R1cOt@3?32C_Z2j;CqvsqhBK%~<0ueELEB%l z+{)QuDhorG^BjO}FfOJxj|QXTjN0Y9<7(ARhfr91i|&TfY&N%Kfq66@GX;nXz2a{W zewRNh8fyMT7>=vC+>4B4K%Z*g7N%J42JT@~z!A}NWH*x<0a9>|sP;CQwo};|vDo9S zq|oh5$Q&H0c6p3ffb-cdLAJRu9qxY=nS4Tj1SKtLWId*~!Vz&npd1=B#6NNGj}|UN zZt#c3-l2&VeEeg);n2{0$PW0l)1ZXp9H1JC0(Yy=Yk0vt@W3vG!5<)gu$=oiLO$I4 z<@;MXzOL){1m!C(+;$BV9pV_5$B(cCeCJ}a+@leUeb?yy{tIG|RHR@J8vP+`%-^<{ z7a4}blCZQgjxm@}8yms2u@C9vhA=jU%Dn61b>3Y+-GCYLe5)9ZY;ggiIzLdHA-UsI z(rm%P5l9E5fq+x5IDkJvM=4GEsa7?Du(KS^Bt`($X877hccCL(>{rZef^JTz7Bbu> zlzz|hiP01Lbm@~&(Lh<|*p0KA5bb6@^+P&>wqzLSt1Vr;^T9tEE`SSXKyuGFQs3Ff^>B`ALhk=mYz*0gR*X<_L>N@h(v=&vN{e$s_Hpex zR0l{&LM$CLNd6|IC2oS+l(tNAk7BGp!&)$GZnesBC?7IMz48U*IGv<1Ww<^YZ2T7I z%OmmTY!4bwZ78vsdA4zlDCqH3X*fN;HbV^4bC(QtZ40ER6dR3-iV~+_>2?1=l63K? zKa}{UaQ^jP)ZL=HGSQ-Y`zLnP@&+`-+Fr|yI@zis5Kqha=h~=sY~mKj`g4rV#YdP+7XuHx9Woe&#-%g%GF70 zEuVGG&tC-TWl^*CkUzKO-srzFKk=FO`Co{(NPb`QJ6+n{_O@(zn|5|DyCJj-y1m$P zxbfyHoz`F;=&?a*SDJ6%9S|w@Itg8pnmus!wUiTR3ULO$!{@0mJ7Bn*^M;70L7!{$ z!$aiDS_gsvEH%3Y>q`?YH>f(W;qai`L5iUYUNzDTJDHO!C{#hKy6NVT6yOuqkeO~P zq7KGNuXiLC=E}MBTRpfdb<@A-RBO|dsVer;y$~5Ks&6>)#B^Zk%iL#;AOX}wkp#M{ zmz#B9CGgxoI^em3OSW)E*;t*kHO#aqrA>)71$V1+Y&%Q| z+7n|sZwonA1c_iMCGEf~BpI^#HW7R6d6hIz^3*ooo1OWIBHaGAgx)uhuRB6SDhyHs z#IYB?G4peN@ovme*?hM~Zd_F-dxbxGJOYw!MHW_Mb&~a55_yzgM~oQedCNC8F1Do> zZm+`uK{8>@Tpfoa+K$JsmnUm37b_7iTMd`tGgK?7$*qUTSxko5Zj##u(Zmc#ddH2J z6$i^Xa!e2g7<+0fx3aVZIGd|Dnn!8OrJjczJnxrAjqgg4 zEz1^W@_1AapdoQxvN=)dyhrnFuG58-=FesT3@0J5Z|jwPC?Didevy0rdMuaeTe#SQ zQPu4Z_tql_YB=QN+)9@xZdm|$XP6a}tm+dK+J1Rvt1zD_@c8@)ZP?ZJUezn%)bfX| zc@>PDIjUs|BB)gN0d zM*}m~;2A8rCv<)!lWvi~ujPr*+xmk&KoSuI8X%u)d0iS52#5t12ng=KvI7@STeE*@ z!KJ#M{W>G6Z*J2!;-74AN=8Ac3J|f1MY)vmqHgX8l&UupCs~;6=NGn5m;PA`nNv59LJHH*@vgoK*_L}c7i z{C<$q-A0GzQsK`)!pTJ1)_T#ty|Ssn2lt}rqJ5I&Xuq*#Q{K(;A?0Nn!&!Pa>=1ln z%LU<$CX>?q9jfd!#Ap(_`Z7n7mAdOSM-RqcnK6|2(t(#;OT&`rTG(>>oHHr4#$|JX zelavglOIw|BnaFCQ$aIRxtx&!oXA|l<$>2dFn6Fm8IH{%ku|3x1MmAg2VuRqM;{6a z6J)6sM6CTYE|AkvI4mNuV$hmVSSG19+P-S8>Y=L5ouK#D^mpPNlM4vh>0|FaYADSi zKCXy*nahUh1K7=~rVsl_;kXNc+*DD>jtxE3yz8k5O+erc&|0!zCm(eIdvc{40>)xg zSaNW2KGyEiM~>mg5BZ<8^13mK35je{BCMxeAK+({w|I>fuXwlY%{y71-eHi_SMEMP z*nc|E8<#617GB}8tc%&4gR8M zL3%WbQ}4#b0NcQnYS2EQ)og5>#+m}7uoFi0|KeQGkrpE`%uk?Cax}$$ z67iHJZxy%l-*nDroriCD6u6plwnaM_6Y_iN8gzb0^;FnxLH2kj*!o=v|IzEgRUIGb z!v4ccA;%{BNwy!gfT1KWc)b|C(?cwyWM5+*Ht)=<^!Xa!N(FyCEbkc!?&#VLaF8VY(dY0UI;Fd=oKs%=(erQv_$k1A(1jVO_Ca#u<7ljRSo6BC))h>*PFXw#pUy;>*0vn85-bm1CyU?>wqX6kARJ zr+19ShjZmuVI9hKgx4c^i8t|>J_9agN-+of<<|0!uHDy}?cp-3>{j`v zr$>~MeN`9Kqxm?ez=yKE!#C3NfFSXiS`awA>Oke$gC%mF+E5q}qlkLE1YYy2%@4JQ%b9hsViLU2!p8uB5;m+A$_qjw|VUS==+ zbQ<79hq`6T{-dA4k=AZdG5c;&+1kS+8mv4G2TtziVxn~kvOJckn2jN3%CLwo`h9Gb zU#|Es5+asj8l_Jo2-oPo)wSnUk^rVcKLS8@M848eJ7Wmk)K4t!p^Qr>dLo9V)v##- zxAT~*anMD>k@)T(AYvgQ>&V3q)M@5;}N1NlF2d= zPGm|oWaBZS?J1P*R?=Y9a*xm(ngRM;X{mus+-mDNqE<4=v7%?v$-WUzgi7!v(x3np zabqlJJvR>GmS67ntQqhpYtyFHsK~@0i-e&x|TAQ5)=@`tG*l%#E)k^ zs!u}lM2@6~fUPk8t(c0w2Tnaf%X2At-+?;rsQ}2qb6M;2c2PikJY^d4 zk#vFX|MPZRMmudFod>!NZmndK#NJS(^_{umBE!|Z{ z=iSoi3D@=HegESy6ruSa*8fe~wu6BM{?as1ZGLDkb{)(Fwc29_u^7cSJOwB(Jyo=q zY;Qe5Ns?<*Pz?%;S0s5gk~`^StV8GjwdB$Ryze@V+0PJy@Gx9Q_i@7SI*`HhYt5e2 z1AkBr_Eg0|ddrAW__FUI-5rJ3IqlEft!^;VYzvY%uhZ|O-RQR-rXh-+R`B|XIV)q$ zk7a0&yF$Oonqn8dZI2*X-3x%9FjCxYUz(=FIEanmo(dx~NQ0dZDSFxp7EA6k6b@n* zIOSPdIFF;32O*@XzYIYf`9s-vM%BxJuv|(y$zV|Lq|&5_6oCPk&lW5aJq)2xZf2^< zDj9Qug7`g1m%VlZtGVtia@i7)~(85g2jhks3hW?A6@oMitqN^x?N|6 zy>~iTQ!n7S81Wrxrja3P;GpM7r1W8?Z?u87dGAdSF*KQK%jrDwATSlrr)>Cq0p?W`Y&5#+JTT$U#VNayXyy4vwM>m-fz(=@Q z*QfCI5!;*Ks*RsImZ>0j{JbJ<3J03(<WjH~j7*=EiE|w2-c`A-FBRzgEx8CRex(s1kHh?ulRz<4V2H z`H8So*EPU4?fZ1`yRmkw>ZdQ;CBJHd<0=B%C*mH_fNkA;(@ElOg{g4tr66%*o| zU(eG$PzcnXzQJesRg6AboOe?D#XZs5Lr@^iXX;Gw{- zrF{cJVO|G$)=$#Tx3JgTiSL=C+?zH(`tRJ*{@(;!)?$r<+xLKyGmBdNFO7yNzY8#a ztyY`c3y69HMojTWl-Fsz~J;JASkVpVeV8kkfzlTVyW}ntSz5|d`@-)|BhB4>9R<&AI!}u2P#bC+&W^jwGP|&Lsqc5+Fo|h zxA*>K>ZzBrqTIX;o56c$T~z>@oY5wx?E7UqTJvZnb2<)$uKYgbQe=IS?Fy<_;T~`= zMHFrLPHF~EWd{>{xz7hPL)`sD1XG0Jao?>8QOp0Hx0I>f zSg&R?&1270>1Sp^A9)^w=2UscIhU%v&BW=< z?j8{@U9+=mad&@V=XnQ|T+e>?!FAkj?@qJZ$)ACppw(AHGLWAuOglZzQUHEkNN!0_ z9^mJC_*_i=K~lCxXE)QF_tlm3k@7-`K)O79r5XtVr(SB1HFu`}E}JmtnxK#%xUBMh zWu$3(240(*l8t2f@l&IXnuSbcW1DN~Q3dx1gwiX{e?OqJuE`KGTzgQoaA#JLW-MhX z`}^-o^ZV~jQ3mt}DiAmjWYSDNSR+rOGaGa|u6h<3)ZZOL$|PV5YyiuD%F>(KKb_Yl zk^gQLz7vpFR693>2?_hhEPF(0;ACBv>|b6>q#=S@(plB9zLIIU_kQLi1jm~PugE)4 zyBYwo@y$H)ZavIbbNxB$e;LJtBS{0pn<_iU&JP$gx2rp~>;+)s;oFezyD>KMq0h!X zOOjlrEf9P;I?W7Z0W77cC~Q+(7As>IO~>>;7!Ip89jy#&}#bwrKns>t%;Fe-`!Z>l(G2;L>cQy zS1dOQ?jvcK$eQVM?cL9>>TM~$rY|RpdbCWmh2=;Vc>RyEARtzwPI2`L)oTa-;_$>U zhpPSHw$AAgBDYdL8nQ^Dp`VC=ol-Z}%4&5H`=v(xr1>!nVECl-F=6~+u$(Bwbj=Fb z34bOYl0sHgL=cc78bxE;V9DS=&Pmea9KgOwf#Y_7l&2mts+`Fa=zGlsUCKlmG3GHD z=*Otfy9*8K0OAD(F*$~Xu~P=C;#q^xP0Eo2MWI-C-?f|)*1~bGv@?zP-;syH3A)sY zQu@ky1f6^Mrx=v?4#oM8%-PXjiRy4v4G?CdPLDLBg$?TS)~@OWCgXn_CSp$clJsa} zqznM6*|@^t|BPb(oL8TlaAAIv(WZJ#Bbu12C5#)_mYP`Ciij2N4dhiP#fUsqPLRD% z1o`Bwy$bLPAfDSvT=|lvDC{j#6VYw}C*zjO@B{NyHIM7cPCKRiA91K=8SbxY4c@^b zWqX^Jp8o6UtU>_z$^SeHnY~GW=cN3V;>71li}60t;q6yi!X^P5539a2jX_zYo>~p>j&<4(FN&gGxlmys91PJs zZ--8Ea~`Q~n{$S*naBhYj2X(tg*7Mx!RJpwH5bpWZ)&Q@>`3s1kT6ha=J|qx=Tqd4 zsczv#_$h@Hz{>&Noh529Rl8e)_3#^^ZDFmKHNF9!jYKsxBRR8?jH2^CXZS7mCrzRx0*8m=l*dDsn=N)Q z@b0#UwZ-S2-}s;btZnofIvqK!2Mmjyg*Hlirk&XDwxLgC$?H*KuGGu!YoB9fM`KUl ztedRlVI^#|BZr|MdszS_Rwcy5?F{;| z!Ma@+;nb#0@cV>=Gdm?na5qDpP`o~ZAMdJ{!$0h;N0%bG>0=J_#@*?|mi4{27N6y0 ze!V*kB+84|{=UH^bGe?@W4zlJ8{BFkvlDc%D>U+dkO#7MBz70Hyene$sry8GZB@Qg zFdfm2!t`32jZwdIR^NtP6n(qv~k6?$Lq-|E5B#Kb!&_J2?`jO8OJAQhr_BsZS_y z#JH3UJ5^@c;sJ1=EKR#UHNd3ARchlmV)HL<{LqtHzMgOO&>;aH>^ol?md#IrO+@S4 zZnB@*s3042qhH_KVF?{Rw%WvI@2OG3xnQxa=uy`b7$H0f%wnV%@LWxDte2Z3Fnj{; z13hgM$!jnR$mO-;m+pz-#xzl_^`8{fFc6du0AyulGTpQOfsg0NC9Im6qR{8hw zJoffa22Rs{-CKZrlI#V=s+jk544t~(Z1g{Vy*}!^TSn(W5nB4~Qk5*V+1);(144=* zdpe;24246wYWj*p!UA|0a|kpiM1q}LP^|ubkka;`DL+N+R3xc9|E=S6@bZKobZhx; zHBz2Q?1zASx%TMCKO;r4Tr6~`ANOOYZa)+VmWwlCugw8fim=}HVK5kATnO$eLtAi* z8iR-T5+69N+QwA++@$WB>;n zw1@du1qF4%A1ump8#qQ!Rx*6pGK7%w)l80abBN; zBadyqt9XD%?XV!X-k;B%TRB{qpA_zEr!+O#!x&Za3J%5n>uQ~o-1oUF2)+)MUvKJI zZ>^aItYaPS*Hm}rnL9l@U4PPBB68~Gf6f+TN6YDlXDK#S)9}A8e>@UeO!Q|FbPBN3 z`e1q;w0Q=n+111M?6wOJdVpj<%T>m$XemnPpopQh>teO0e#8cH5HtUYK4?hQJL35y z_)n1=S{-G-5EbUHB8ll=Zo&MY+~QE@!g-$)`TMU7=|82nG*b$H1Hs!LEdwM0x2xr2 zm)l|0<#a|rJbt%?|GZzV?Cky2PGxE_h0GvQdwE4Z$nDgr{%STYuh#u~JZwJ=_>nL2 zSyGUq28Og#!t&dIHN}A6AdMasReUwqDDADNoHt4FPx;eIuP*;QVArONlvLQTA0ca; zIZ|D{76Sg-JF7`dZF4b8Mx#6dueFq||mVteB!|x=3Z8kWoCbF=g&JJ()*iN1v#kQAeX}#oBZC=4v((Yo;GH~Ur z5fuX@Q;*b!j=5F}Q@^*0SFmIwLY5~bAOY3Jm9&V4h;OaI&f0TgC2=_>_>s(%kf;@zIXONdIGW5f zQH&b#gmsmXIu_B^yto`lBTKM_tcZk259&;6v^jXw;RiCW$fX$nhc7)SR=1A1ga2gp z;+vu{KpdGlHC_sREQH`Vi6l5z%Zhqda$T>RsHlV1xXd#lmbM_myaSB6AX`5nnv*}$ z^s!5TLyS38&?O-{BMobmnazl}jV?SEu`i!i84VX7EgbSI)(lA{&q$e2S;jkA1|QR2 zLK5FRgh6CJV!jM|Eor<6rMOBBZGSsO%tZtfP}W^FnLIN;Agg`2&5nOpynX24Q?InU5tsDk(gI%sfWS zmAaSW^&>R{Ey3*0HGRB@c0wC@q9NQ4(xfykh045rlrQ1voJLdL48gFy41#z<*)kLv zz?>LgA1vD-8EfW|ix*&KE@rBQ>2#$^Zl<-IsF^C};<$fqX2v-&u;`|Dm@~?pw)X+` z_Q^K`4u+3S9HqPxhWvjC#55go1H}FkDEZ`;GE+r0qz)Me7)1cXM_9W8f6>D5SG)Gs(hlP5*$PeL=PUnDF{Wsbgji3d7u9q}DW`$}4T zK=7vYsG{)y`N(&J)&jp*vshS&MfNjU>}b%}0oT)Jd=Z7!>uxs|&}TanP-u;+rV<5G z&-Dx zGSzliHBE{s(|YIVS~SdCTJIHgLlUQIS6z3gV6MyXA~d10_ROda|)2|Fy@e< z0lPwiL7k@=De44{yCQVB%*Ij#5MS1KF;bYqbSWBvK$lEsw3w9uH{ufhbh6`2=^2ac zA(5pY?BK{+PZXta`jM#N0wM&w(?z(Xf`m z?)sz_KY+3t+INmFVcu>V1e*#gnPHfFewTu!bb*AhPtmacAuj5pV!N*TEVX~EM8>Y; zE`Z1j=PM4o8IFwC>kqcqMWU~dR;Q+oPyoqix6qS1sCpe#)!^$`T$DW?#qAaapbsg8 zm|~)5M&|Mx_wjmVH3R{oMuuR6{<*g#=;VpSI*>0hw%mSH@)!-rnjeQSmKH-(iG|e^ zJ;U#wuw{l|)gTy)nt+T+=!;-5smN)isbkEdYhb1NB6uzq<4VaqLV@%PPy&~0SNvk1 zrJUAx3ox{A^u3|asmu*x+v{Bo?g!|ZaL5WeR*bUK9uMZE{z=;!`k0U-SGvwX*P%5h zhYdN)`Ta#zo5A%;d^)Aok7RU_JDGeB(+0Uh{R ze)^dSi~_~D(YARAYhDoy;Ao?BU2-`^a+l(0HMqj@mhAA#ZGlW)D|ex-GS}CXr<#1^ zypdqP%c8Z>eI$s(({oVWon0H=E@Pwy!-_60u!J5KS(G|-R${(&$kJ2vv~6?lMpjrn zx1cjn7=*#~e$b+~M{5`w#>@9bY!RK?Ls5Zqg;CI9FG;4Cq#REo(u0iDR6|&`x^wwzaW>6mlTAUrqKsRB! z%!Eje&Za`W@aCKfd zQ;oTrI(~cK8 z5k3?r!wnoJ;IDNa@~!V`%MbneAgy;XPN1p5ydSGJhvJ> zCprRdQ-dv2i;6`Xf!){Wyh_2RW9 zyMfWwCe3ypz{Bxw_YcMM=~?HDs82;qbD;n)reVhL)Iafun%BU1$$ zCiD{ z^6c)0lLtRH-Cl3tOxb2dBDRU>rp|3NNv`4KE~<(mS8CfU#jiET%SAj-6PE!cGg5Oe z28B#V+vmlu5Tz$QL+2t(i+HWkL?4cOlOAENrR3{S>tR`HBrVBZK=;#7^HxG-D#{wK z=Wdvp$^rP;DbkU?mAhnV&JYJi8Fx3)SczX9BRp)KuZZ*;hr*of^CWt?_n(He(P|R` z91L>*KEtfhcb^m(5S=(-H%N#m_LB03fbNYW@+O&ZRfbTURZ_92THR5|C{?v;2wa}j z%g}=I+n2T@C}U!8JHWJoE@~Ggda79(DZ&`lA6fg@I-?-%4=ZkAA9oEYLsKB?RbPRW zn^$6ZB-GE{u~X>iiGJ{$B3cx+NcCwoSQ8p`MJF83ko&0+06UqxpG|7pzVTF=zZI^* zE;4S~rsF2iUI=)=$|vkfJ_CL8x&pS>nJrZHSU40KO#;s98p*-Vw{y75tjA z%ane1vnNB<1$Tu+HQY_^ZtL2hGY;FN8Ve0DrA(0_y*rpG6hHl^H0G{+R?VjU2kS~* zWpe_%_hXQ}(QVy4Z0GLs)JG?+Q#?oK<$BxF$d8wfcml)J>aVaN9ZM8a1ARx^D=KPE5LjO7K{PiD8_y@@@7(l|){}I~$g~s{+tPOu5 zYUMxcpTBUr{vQ;)WB^u7a<~)+S?&90n7%~A{qGtx5D?b?2*&&SpA*^@8t#8D7XPw^D_&shr12|q+'", str(w[-1].message)) + expected_exchange_rxns = set([self.model.reactions.get_one(id=rxn_id) + for rxn_id in ('ex_m3', 'ex_m1_forward', 'ex_m1_backward', + 'ex_m2_forward', 'ex_m2_backward')]) + self.assertEqual(dfba_submodel.exchange_rxns, expected_exchange_rxns) # test dfba_obj_expr model = self.get_model() @@ -218,15 +234,25 @@ def test_dfba_submodel_init(self): # test exception in dfba_obj.expression that contains an exchange reaction with self.assertRaisesRegexp(MultialgorithmError, f"the dfba objective 'dfba-obj-{dfba_submodel.id}' uses exchange reactions:"): - self.make_dfba_submodel(model, obj_expression_spec=(('biomass_reaction', ), ('ex_m1', ))) + self.make_dfba_submodel(model, obj_expression_spec=(('biomass_reaction', ), ('ex_m1_forward', + 'ex_m3'))) def test_set_up_optimizations(self): dfba_submodel = self.dfba_submodel_1 self.assertTrue(set(dfba_submodel.species_ids) == dfba_submodel.species_ids_set \ == set(dfba_submodel.adjustments.keys())) self.assertEqual(dfba_submodel.populations.shape, ((len(dfba_submodel.species_ids), ))) - self.assertEqual(dfba_submodel.reaction_fluxes, - {'ex_m1': None, 'ex_m2': None, 'ex_m3': None, 'r1': None, 'r2': None, 'r3': None, 'r4': None}) + self.assertEqual(dfba_submodel.reaction_fluxes, {'ex_m1_backward': None, + 'ex_m1_forward': None, + 'ex_m2_backward': None, + 'ex_m2_forward': None, + 'ex_m3': None, + 'r1_backward': None, + 'r1_forward': None, + 'r2_backward': None, + 'r2_forward': None, + 'r3': None, + 'r4': None}) def test_set_up_dfba_submodel(self): # test exception when the ids in DfbaObjReactions and Reactions intersect @@ -244,15 +270,32 @@ def test_set_up_dfba_submodel(self): for k,v in self.dfba_submodel_1._conv_variables.items(): self.assertEqual(k, v.name) - expected_stoch_consts = { - 'm1[c]': {'ex_m1': -1, 'r1': -1, 'r2': -1, 'r3': -2}, - 'm2[c]': {'ex_m2': -1, 'r1': -1, 'r2': 1, 'r4': -2}, - 'm3[c]': {'ex_m3': -1, 'r1': 1, 'r3': 1, 'r4': 1, 'biomass_reaction': -1}, - } + # test the stoichiometric matrix + expected_stoch_coefficients = { + 'm1[c]': {'ex_m1_backward': 1, + 'ex_m1_forward': -1, + 'r1_backward': 1, + 'r1_forward': -1, + 'r2_backward': 1, + 'r2_forward': -1, + 'r3': -2}, + 'm2[c]': {'ex_m2_backward': 1, + 'ex_m2_forward': -1, + 'r1_backward': 1, + 'r1_forward': -1, + 'r2_backward': -1, + 'r2_forward': 1, + 'r4': -2}, + 'm3[c]': {'biomass_reaction': -1, + 'ex_m3': -1, + 'r1_backward': -1, + 'r1_forward': 1, + 'r3': 1, + 'r4': 1}} test_results = {species_id: {lt.variable.name: lt.coefficient for lt in linear_terms} for species_id, linear_terms in self.dfba_submodel_1._conv_metabolite_matrices.items()} - self.assertEqual(test_results, expected_stoch_consts) + self.assertEqual(test_results, expected_stoch_coefficients) self.assertEqual(len(self.dfba_submodel_1.get_conv_model().variables), len(self.dfba_submodel_1._conv_variables.values())) @@ -260,14 +303,13 @@ def test_set_up_dfba_submodel(self): self.assertEqual(conv_opt_var in self.dfba_submodel_1.get_conv_model().variables, True) for constr in self.dfba_submodel_1.get_conv_model().constraints: - if constr.name in expected_stoch_consts: + if constr.name in expected_stoch_coefficients: self.assertEqual({linear_term.variable.name:linear_term.coefficient for linear_term in constr.terms}, - expected_stoch_consts[constr.name]) + expected_stoch_coefficients[constr.name]) self.assertEqual(constr.upper_bound, 0.) self.assertEqual(constr.lower_bound, 0.) - self.assertEqual(self.dfba_submodel_1._dfba_obj_rxn_ids, ['biomass_reaction']) self.assertEqual({i.variable.name:i.coefficient for i in self.dfba_submodel_1.get_conv_model().objective_terms}, {'biomass_reaction': 1}) @@ -289,7 +331,7 @@ def test_set_up_dfba_submodel(self): test_results = {species_id: {lt.variable.name: lt.coefficient for lt in linear_terms} for species_id, linear_terms in dfba_submodel_2._conv_metabolite_matrices.items()} - self.assertEqual(test_results, expected_stoch_consts) + self.assertEqual(test_results, expected_stoch_coefficients) self.assertEqual(len(dfba_submodel_2.get_conv_model().variables), len(dfba_submodel_2._conv_variables.values())) @@ -298,17 +340,16 @@ def test_set_up_dfba_submodel(self): self.assertEqual(conv_opt_var in dfba_submodel_2.get_conv_model().variables, True) for constr in dfba_submodel_2.get_conv_model().constraints: - if constr.name in expected_stoch_consts: + if constr.name in expected_stoch_coefficients: self.assertEqual({linear_term.variable.name:linear_term.coefficient for linear_term in constr.terms}, - expected_stoch_consts[constr.name]) + expected_stoch_coefficients[constr.name]) self.assertEqual(constr.upper_bound, 0.) self.assertEqual(constr.lower_bound, 0.) - self.assertEqual(dfba_submodel_2._dfba_obj_rxn_ids, ['biomass_reaction']) self.assertEqual({objective_term.variable.name: objective_term.coefficient for objective_term in dfba_submodel_2.get_conv_model().objective_terms}, - {'biomass_reaction': 1, 'r2': 1}) + {'biomass_reaction': 1, 'r2_forward': 1}) self.assertEqual(dfba_submodel_2.get_conv_model().objective_direction, conv_opt.ObjectiveDirection.minimize) @@ -320,8 +361,8 @@ def get_species(id): dfba_submodel_2 = self.make_dfba_submodel(self.model, submodel_options=self.dfba_submodel_options, dfba_obj_with_regular_rxn=True) - expected = dict(r2={get_species('m1[c]'): -1., - get_species('m2[c]'): 1.}, + expected = dict(r2_forward={get_species('m1[c]'): -1., + get_species('m2[c]'): 1.}, biomass_reaction={get_species('m3[c]'): -1.}) for rxn_class in dfba_submodel_2.dfba_obj_expr.related_objects: for rxn_id, rxn in dfba_submodel_2.dfba_obj_expr.related_objects[rxn_class].items(): @@ -329,12 +370,12 @@ def get_species(id): expected[rxn_id]) # also test reactions with the same species on both sides - r2 = self.model.reactions.get_one(id='r2') - r2.participants.append(get_species(id='m1[c]').species_coefficients.get_or_create(coefficient=1)) + r2_forward = self.model.reactions.get_one(id='r2_forward') + r2_forward.participants.append(get_species(id='m1[c]').species_coefficients.get_or_create(coefficient=1)) dfba_submodel_2 = self.make_dfba_submodel(self.model, submodel_options=self.dfba_submodel_options, dfba_obj_with_regular_rxn=True) - expected = dict(r2={get_species('m2[c]'): 1.}, + expected = dict(r2_forward={get_species('m2[c]'): 1.}, biomass_reaction={get_species('m3[c]'): -1.}) for rxn_class in dfba_submodel_2.dfba_obj_expr.related_objects: for rxn_id, rxn in dfba_submodel_2.dfba_obj_expr.related_objects[rxn_class].items(): @@ -386,31 +427,35 @@ def check_neg_species_pop_constraints(test_case, constraints, expected_constrs): test_case.assertEqual(set_of_terms, expected_constrs[id]) # generate constraint involving 2 pseudo-reactions that use m3[c], ex_m3 and biomass_reaction - dfba_submodel_2 = self.make_dfba_submodel(self.model, + dfba_submodel_2 = self.make_dfba_submodel(self.get_model(), submodel_options=self.dfba_submodel_options) constraints = dfba_submodel_2.initialize_neg_species_pop_constraints() - print(ShowConvOptElements.show_conv_opt_constraints(constraints)) - # for each species, set of expected (rxn, coef) pairs contributing to species' consumption - expected_constrs = {get_const_name('m3[c]'): {('ex_m3', 1.0), ('biomass_reaction', 1.0)}} + # for each species, set of expected (rxn, coef) pairs contributing to the species' consumption + raw_exp_constrs = dict(m1={('ex_m1_backward', -1), ('ex_m1_forward', 1)}, + m2={('ex_m2_backward', -1), ('ex_m2_forward', 1)}, + m3={('ex_m3', 1), ('biomass_reaction', 1)}) + + expected_constrs = {get_const_name(f'{s_id}[c]'): constrs for s_id, constrs in raw_exp_constrs.items()} check_neg_species_pop_constraints(self, constraints, expected_constrs) - self.assertEqual(get_rxn_set(['ex_m3']), dfba_submodel_2._constrained_exchange_rxns) - # TODO (APG): later: test with dFBA objective that uses multiple dFBA objective reactions - # avoid the need for any constraints by remove biomass_reaction from the dfba objective - dfba_submodel_2 = self.make_dfba_submodel(self.model, - submodel_options=self.dfba_submodel_options, - obj_expression_spec=(tuple(), ('r2',))) + # add another dfba obj. species to biomass_reaction to cover the remaining branches + model = self.get_model() + biomass_reaction = model.dfba_obj_reactions.get_one(id='biomass_reaction') + # m1[n] won't need a constraint because it is not consumed + dfba_obj_species = wc_lang.DfbaObjSpecies(value=2, species=model.species.get_one(id='m1[n]')) + biomass_reaction.dfba_obj_species.add(dfba_obj_species) + dfba_submodel_2 = self.make_dfba_submodel(model, + submodel_options=self.dfba_submodel_options) constraints = dfba_submodel_2.initialize_neg_species_pop_constraints() - print(ShowConvOptElements.show_conv_opt_constraints(constraints)) - check_neg_species_pop_constraints(self, constraints, {}) + check_neg_species_pop_constraints(self, constraints, expected_constrs) def test_bound_neg_species_pop_constraints(self): dfba_submodel_2 = self.make_dfba_submodel(self.model, submodel_options=self.dfba_submodel_options, dfba_obj_with_regular_rxn=True) dfba_submodel_2.bound_neg_species_pop_constraints() - expected_constraint_bounds = {'neg_pop_constr__m1[c]': (0, 100), - 'neg_pop_constr__m3[c]': (0, 100)} + expected_constraint_bounds = {'neg_pop_constr__m1[c]': (None, 100), + 'neg_pop_constr__m3[c]': (None, 100)} conf_opt_model = dfba_submodel_2.get_conv_model() for constraint in conf_opt_model.constraints: if constraint.name in expected_constraint_bounds: @@ -422,111 +467,40 @@ def bounds_test(test_case, dfba_submodel, expected_bounds): """ Test whether the reaction bounds specified by `expected_bounds` are set in `dfba_submodel` """ for rxn_id, bounds in dfba_submodel._reaction_bounds.items(): - print('rxn_id',rxn_id) - for ind, bound in enumerate(bounds): - print('bound, expected', bound, expected_bounds[rxn_id][ind], 'bound == expected_bounds[rxn_id][ind]', - bound == expected_bounds[rxn_id][ind]) - test_case.assertAlmostEqual(bound, expected_bounds[rxn_id][ind], delta=1e-09) + if rxn_id in expected_bounds: + for ind, bound in enumerate(bounds): + test_case.assertAlmostEqual(bound, expected_bounds[rxn_id][ind], delta=1e-09) def test_determine_bounds(self): self.dfba_submodel_1.determine_bounds() self.bounds_test(self.dfba_submodel_1, self.expected_rxn_flux_bounds) - # test changing flux_bounds_volumetric_compartment_id + ### test changing flux_bounds_volumetric_compartment_id new_options = copy.deepcopy(self.dfba_submodel_options) new_options['flux_bounds_volumetric_compartment_id'] = 'c' new_submodel = self.make_dfba_submodel(self.model, submodel_options=new_options) - self.expected_rxn_flux_bounds['ex_m1'] = (-120. * 5e-13, -100. * 5e-13) - self.expected_rxn_flux_bounds['ex_m2'] = (-120. * 5e-13, -100. * 5e-13) + self.expected_rxn_flux_bounds['ex_m1_backward'] = (0., -(-120. * 5e-13)) + self.expected_rxn_flux_bounds['ex_m1_forward'] = (0., -100. * 5e-13) + self.expected_rxn_flux_bounds['ex_m2_backward'] = (0., -(-120. * 5e-13)) + self.expected_rxn_flux_bounds['ex_m2_forward'] = (0., -100. * 5e-13) new_submodel.determine_bounds() self.bounds_test(new_submodel, self.expected_rxn_flux_bounds) - - # test additional branches in determine_bounds WITHOUT negative species populations bounds - model = self.get_model() - # reaction missing a forward rate law - r1 = model.reactions.get_one(id='r1') - self.assertEqual(r1.rate_laws[1].direction, wc_lang.RateLawDirection.forward) - del r1.rate_laws[1] - # reaction missing a backward (reverse) rate law - r2 = model.reactions.get_one(id='r2') - self.assertEqual(r2.rate_laws[0].direction, wc_lang.RateLawDirection.backward) - del r2.rate_laws[0] - # no bounds in an irreversible rxn - model.reactions.get_one(id='ex_m1').flux_bounds = None - # no bounds in a reversible rxn - model.reactions.get_one(id='ex_m2').reversible = True - model.reactions.get_one(id='ex_m2').flux_bounds = None - # bounds of NaN - model.reactions.get_one(id='ex_m3').flux_bounds.min = numpy.nan - model.reactions.get_one(id='ex_m3').flux_bounds.max = numpy.nan - self.dfba_submodel_options['negative_pop_constraints'] = False - dfba_submodel_2 = self.make_dfba_submodel(model, - submodel_options=self.dfba_submodel_options) - dfba_submodel_2.determine_bounds() - expected_results_2 = { - 'ex_m1': (0, None), # no flux bounds irreversible rxn => (lower, upper) == (0, None) - 'ex_m2': (None, None), # no flux bounds reversible rxn => (lower, upper) == (None, None) - 'ex_m3': (None, None), # bounds of NaN => (lower, upper) == (None, None) - 'r1': (-1. * 1, None), # missing forward rate law => upper == None - 'r2': (None, 1. * 1 + 2. * 1), # missing backward rate law => lower == None - 'r3': (0., 5. * 1), - 'r4': (0., 6. * 1), - } - self.bounds_test(dfba_submodel_2, expected_results_2) - - # test additional branches in determine_bounds WITH negative species populations bounds - print('\n# test additional branches in determine_bounds():') - - def test_determine_exchange_bounds(self): - # test bounds on exchange reactions that avoid negative species populations - - # test exchange rxn with reactant and no max bound - # add a metabolite - met_st_new = self.model.species_types.create(id='m_new', - type=wc_ontology['WC:metabolite'], - structure=wc_lang.ChemicalStructure(molecular_weight=1.)) - comp_c = self.model.compartments.get_one(id='c') - met_species_new = self.model.species.create(species_type=met_st_new, compartment=comp_c) - met_species_new.id = met_species_new.gen_id() - conc_model = self.model.distribution_init_concentrations.create(species=met_species_new, - mean=100., - std=0., - units=unit_registry.parse_units('molecule')) - conc_model.id = conc_model.gen_id() - - # add an exchange reaction that consumes new species - ex4 = self.submodel.reactions.create(id='ex_m4', reversible=False, model=self.model) - ex4.flux_bounds = wc_lang.FluxBounds(min=0., max=None) - ex4.participants.append(met_species_new.species_coefficients.get_or_create(coefficient=-2)) - - # test exchange rxn with reactant used by another reaction and large max bound - ex5 = self.submodel.reactions.create(id='ex_m5', reversible=False, model=self.model) - ex5.flux_bounds = wc_lang.FluxBounds(min=0., max=None) - m3 = self.model.species.get_one(id='m3[c]') - ex5.participants.append(m3.species_coefficients.get_or_create(coefficient=-1)) - self.expected_rxn_flux_bounds['ex_m5'] = (0, None) - - dfba_submodel_3 = self.make_dfba_submodel(self.model, - submodel_options=self.dfba_submodel_options) - rxn_id, species_id, coefficient = ('ex_m4', 'm_new[c]', -2) - species_pop = dfba_submodel_3.local_species_population.read_one(0, species_id) - expected_max_constr = -species_pop / (coefficient * dfba_submodel_3.time_step) - self.expected_rxn_flux_bounds[rxn_id] = (0, expected_max_constr) - - dfba_submodel_3.determine_bounds() - self.bounds_test(dfba_submodel_3, self.expected_rxn_flux_bounds) + # TODO (APG): to cover other branches in determine_bounds(), test without using PrepForWcSimTransform def test_update_bounds(self): new_bounds = { - 'ex_m1': (12., 20.), - 'ex_m2': (27.4, 33.8), + 'ex_m1_backward': (0., 0.), + 'ex_m1_forward': (12., 20.), + 'ex_m2_backward': (0., 0.), + 'ex_m2_forward': (27.4, 33.8), 'ex_m3': (0., 0.), - 'r1': (-2.5, 1.8), - 'r2': (-3, None), + 'r1_backward': (0., 2.5), + 'r1_forward': (0., 1.8), + 'r2_backward': (0., 3.), + 'r2_forward': (0., None), 'r3': (0., 30.), 'r4': (None, None), - 'biomass_reaction': (0., None), - } + 'biomass_reaction': (0., None)} self.dfba_submodel_1._reaction_bounds = new_bounds self.dfba_submodel_1.update_bounds() for var_id, variable in self.dfba_submodel_1._conv_variables.items(): @@ -543,27 +517,32 @@ def do_test_compute_population_change_rates_control_caching(self, caching_settin with EnvironUtils.make_temp_environ(**config_env_dict.get_env_dict()): # use fluxes and expected population change rates from - # test I: 'No scaling (scaling factors equal 1) and no negative species population checks' - fluxes = dict(ex_m1=-12, - ex_m2=-12, - ex_m3=0, - r1=1, - r2=1, - r3=5, - r4=6, - biomass_reaction=12) + # test I in "Solutions to test dFBA models, by hand.txt" + fluxes = dict(ex_m1_backward=12., + ex_m1_forward=0., + ex_m2_backward=12., + ex_m2_forward=0., + ex_m3=0., + r1_backward=0., + r1_forward=1., + r2_backward=0., + r2_forward=1., + r3=5., + r4=6., + biomass_reaction=12.) + self.dfba_submodel_1.reaction_fluxes = fluxes self.dfba_submodel_1.compute_population_change_rates() - expected_rates = {'m1[c]': -12, - 'm2[c]': -12, - 'm3[c]': 12} + expected_rates = {'m1[c]': -12., + 'm2[c]': -12., + 'm3[c]': 12.} self.assertEqual(self.dfba_submodel_1.adjustments, expected_rates) # test repeated calculations self.dfba_submodel_1.compute_population_change_rates() self.assertEqual(self.dfba_submodel_1.adjustments, expected_rates) def test_compute_population_change_rates(self): - ### test all 3 caching combinations ### + ### test all 3 caching combinations: # NO CACHING # EVENT_BASED invalidation # REACTION_DEPENDENCY_BASED invalidation @@ -661,17 +640,23 @@ def check_expected_solution(test_case, dfba_submodel, obj_func_value, expected_a def test_run_fba_solver(self): # Algebraic solutions to these tests are documented in the # file "tests/submodels/fixtures/Solutions to test dFBA models, by hand.txt" - test_name = 'I: No scaling (scaling factors equal 1) and no negative species population checks' - # TODO (APG): later: report strange behavior, presumably caused by obj_table's shared related attributes - # since ex_m1 flux_bounds.min == ex_m2 flux_bounds.min, multiplying one of them by 1E11 increases both of them by 1E11 - self.model.reactions.get_one(id='ex_m1').flux_bounds.min *= 1e11 - # self.model.reactions.get_one(id='ex_m2').flux_bounds.min *= 1e11 + # TODO (APG): update this file for split, irreversible reactions + TEST_NAME = 'I: No scaling (scaling factors equal 1) and no negative species population checks' + ''' + print('flux_bounds before') + for r_id in ['ex_m1_backward', 'ex_m1_forward', 'ex_m2_backward', 'ex_m2_forward']: + r = self.model.reactions.get_one(id=r_id) + print(f"{r_id}: ({r.flux_bounds.min}, {r.flux_bounds.max})") + ''' + self.model.reactions.get_one(id='ex_m1_backward').flux_bounds.max *= 1e11 + self.model.reactions.get_one(id='ex_m2_backward').flux_bounds.max *= 1e11 + self.dfba_submodel_options['dfba_bound_scale_factor'] = 1. self.dfba_submodel_options['dfba_coef_scale_factor'] = 1. self.dfba_submodel_options['negative_pop_constraints'] = False dfba_submodel_2 = self.make_dfba_submodel(self.model, submodel_options=self.dfba_submodel_options) - dfba_submodel_2.get_conv_model().name = test_name + dfba_submodel_2.get_conv_model().name = TEST_NAME dfba_submodel_2.time_step = 1. dfba_submodel_2.run_fba_solver() @@ -681,27 +666,27 @@ def test_run_fba_solver(self): 'm3[c]': 12} self.check_expected_solution(dfba_submodel_2, 12, expected_adjustments_I) - test_name = 'II: Add negative species population constraints to I' + TEST_NAME = 'II: Add negative species population constraints to I' self.dfba_submodel_options['negative_pop_constraints'] = True dfba_submodel_2 = self.make_dfba_submodel(self.model, submodel_options=self.dfba_submodel_options) - dfba_submodel_2.get_conv_model().name = test_name + dfba_submodel_2.get_conv_model().name = TEST_NAME dfba_submodel_2.run_fba_solver() self.check_expected_solution(dfba_submodel_2, 12, expected_adjustments_I) - test_name = 'III: Modify II by scaling bounds by 10' + TEST_NAME = 'III: Modify II by scaling bounds by 10' self.dfba_submodel_options['dfba_bound_scale_factor'] = 10. dfba_submodel_2 = self.make_dfba_submodel(self.model, submodel_options=self.dfba_submodel_options) - dfba_submodel_2.get_conv_model().name = test_name + dfba_submodel_2.get_conv_model().name = TEST_NAME dfba_submodel_2.run_fba_solver() self.check_expected_solution(dfba_submodel_2, 12, expected_adjustments_I) self.dfba_submodel_options['dfba_bound_scale_factor'] = 1. - test_name = 'IV: Alter II so that a negative species population constraints change the solution' + TEST_NAME = 'IV: Alter II so that a negative species population constraints change the solution' dfba_submodel_2 = self.make_dfba_submodel(self.model, submodel_options=self.dfba_submodel_options) - dfba_submodel_2.get_conv_model().name = test_name + dfba_submodel_2.get_conv_model().name = TEST_NAME dfba_submodel_2.time_step = 10 dfba_submodel_2.run_fba_solver() expected_adjustments_IV = {'m1[c]': -8, @@ -710,29 +695,29 @@ def test_run_fba_solver(self): self.check_expected_solution(dfba_submodel_2, 10, expected_adjustments_IV) dfba_submodel_2.time_step = 1 - test_name = "V: Modify II by scaling the coefficients of objective function's terms by 10" + TEST_NAME = "V: Modify II by scaling the coefficients of objective function's terms by 10" self.dfba_submodel_options['dfba_coef_scale_factor'] = 10. dfba_submodel_2 = self.make_dfba_submodel(self.model, submodel_options=self.dfba_submodel_options) - dfba_submodel_2.get_conv_model().name = test_name + dfba_submodel_2.get_conv_model().name = TEST_NAME dfba_submodel_2.run_fba_solver() self.check_expected_solution(dfba_submodel_2, 12, expected_adjustments_I) self.dfba_submodel_options['dfba_coef_scale_factor'] = 1. - test_name = 'VI: Modify V by scaling bounds by 5' + TEST_NAME = 'VI: Modify V by scaling bounds by 5' self.dfba_submodel_options['dfba_bound_scale_factor'] = 5. self.dfba_submodel_options['dfba_coef_scale_factor'] = 10. dfba_submodel_2 = self.make_dfba_submodel(self.model, submodel_options=self.dfba_submodel_options) - dfba_submodel_2.get_conv_model().name = test_name + dfba_submodel_2.get_conv_model().name = TEST_NAME dfba_submodel_2.run_fba_solver() self.check_expected_solution(dfba_submodel_2, 12, expected_adjustments_I) self.dfba_submodel_options['dfba_bound_scale_factor'] = 1. self.dfba_submodel_options['dfba_coef_scale_factor'] = 1. # test raise DynamicMultialgorithmError - self.model.reactions.get_one(id='ex_m1').flux_bounds.min = 1000. - self.model.reactions.get_one(id='ex_m1').flux_bounds.max = 1000. + self.model.reactions.get_one(id='ex_m1_forward').flux_bounds.min = 1000. + self.model.reactions.get_one(id='ex_m1_forward').flux_bounds.max = 1000. dfba_submodel_3 = self.make_dfba_submodel(self.model, submodel_options=self.dfba_submodel_options) dfba_submodel_3.time = 0.1 @@ -743,24 +728,21 @@ def test_run_fba_solver(self): "'infeasible' for time step [0.1, 1.1]")): dfba_submodel_3.run_fba_solver() - ''' species = ['m1[c]', 'm2[c]', 'm3[c]'] expected_population = dict(zip(species, [ - 100 - 120. * self.cell_volume * dfba_submodel_2.time_step * 1e11, - 100 - 120. * self.cell_volume * dfba_submodel_2.time_step * 1e11, - 100 + 120. * self.cell_volume * dfba_submodel_2.time_step * 1e11])) + 100 - 120. * self.cell_volume * dfba_submodel_2.time_step * 1E11, + 100 - 120. * self.cell_volume * dfba_submodel_2.time_step * 1E11, + 100 + 120. * self.cell_volume * dfba_submodel_2.time_step * 1E11])) population = dfba_submodel_2.local_species_population.read(1., set(species)) self.assertEqual(population, expected_population) - ''' # TODO (APG): OPTIMIZE DFBA CACHING: test that expressions that depend on exchange and biomass reactions # have been flushed, and that expressions which don't depend on them have not been flushed # test flush expression - # TODO (APG): fix # self.assertTrue(dfba_submodel_2.dynamic_model.cache_manager.empty()) - return # TODO (APG): later: test using a different solver + ''' self.dfba_submodel_options['solver'] = 'glpk' del self.dfba_submodel_options['solver_options'] dfba_submodel_2 = self.make_dfba_submodel(self.model, @@ -780,6 +762,7 @@ def test_run_fba_solver(self): 100 + 120 * self.cell_volume * dfba_submodel_2.time_step * 1e11])) population = dfba_submodel_2.local_species_population.read(1., set(species)) self.assertEqual(population, expected_population) + ''' def test_schedule_next_fba_event(self): # check that the next event is a RunFba message at time expected_time diff --git a/wc_sim/multialgorithm_simulation.py b/wc_sim/multialgorithm_simulation.py index 0e16804e..d92513a0 100644 --- a/wc_sim/multialgorithm_simulation.py +++ b/wc_sim/multialgorithm_simulation.py @@ -179,6 +179,11 @@ def prepare(self): self.model._transformed = True errors = Validator().run(self.model) if errors: + # TODO (APG): skip exception on Validator errors until github.com/KarrLab/wc_lang/issues/142 is fixed + if self.model.id == 'dfba_test': + warnings.warn(f"cannot validate model '{self.model.id}' in fixtures/dfba_test_model.xlsx", + MultialgorithmWarning) + return raise MultialgorithmError(indent_forest(['The model is invalid:', [errors]])) def build_simulation(self, prepare_model=True): diff --git a/wc_sim/submodels/dfba.py b/wc_sim/submodels/dfba.py index 717ff0a1..681a5020 100644 --- a/wc_sim/submodels/dfba.py +++ b/wc_sim/submodels/dfba.py @@ -11,9 +11,9 @@ import conv_opt import copy import enum +import itertools import math import scipy.constants -import warnings import wc_lang from wc_sim import message_types @@ -50,8 +50,6 @@ class DfbaSubmodel(ContinuousTimeSubmodel): _multi_reaction_constraints (:obj:`dict` of :obj:`str`: :obj:`conv_opt.Constraint`): a map from constraint id to constraints that avoid negative species populations in `self._conv_model.constraints` - _constrained_exchange_rxns (:obj:`set` of :obj:`wc_lang.Reactions`): exchange and demand reactions - constrained by `_multi_reaction_constraints` _conv_model (:obj:`conv_opt.Model`): linear programming model in `conv_opt` format _conv_variables (:obj:`dict` of :obj:`str`: :obj:`conv_opt.Variable`): a dictionary mapping reaction IDs to their associated `conv_opt.Variable` objects @@ -59,11 +57,12 @@ class DfbaSubmodel(ContinuousTimeSubmodel): species IDs to lists of :obj:`conv_opt.LinearTerm` objects; each :obj:`conv_opt.LinearTerm` associates a reaction that the species participates in with the species' stoichiometry in the reaction - _dfba_obj_rxn_ids (:obj:`list` of :obj:`str`): a list of IDs of dFBA objective reactions + _dfba_obj_reactions (:obj:`dict` of :obj:`str`: :obj:`wc_lang.DfbaObjReaction`): all + :obj:`wc_lang.DfbaObjReaction`\ s used by the :obj:`self.dfba_obj_expr` _dfba_obj_species (:obj:`list` of :obj:`wc_lang.DfbaObjSpecies:`): all species in - :obj:`DfbaObjReaction`\ s used by `dfba_obj_expr` + :obj:`DfbaObjReaction`\ s used by `dfba_obj_expr`, keyed by their IDs _reaction_bounds (:obj:`dict` of :obj:`str`: :obj:`tuple`): a dictionary that maps reaction IDs - to tuples of scaled minimum and maximum bounds + to (minimum bound, maximum bound) tuples _optimal_obj_func_value (:obj:`float`): the value of objective function returned by the solver """ # default options @@ -115,7 +114,7 @@ def __init__(self, id, dynamic_model, reactions, species, dynamic_compartments, Raises: :obj:`MultiAlgorithmError`: if the `dynamic_dfba_objective` cannot be found, - or if the exchange reactions are not all of the form “s ->”, + or if some reactions are reversible, or if the provided 'dfba_bound_scale_factor' in options does not have a positive value, or if the provided 'dfba_coef_scale_factor' in options does not have a positive value, or if the provided 'solver' in options is not a valid value, @@ -184,27 +183,36 @@ def __init__(self, id, dynamic_model, reactions, species, dynamic_compartments, if 'negative_pop_constraints' in options: self.dfba_solver_options['negative_pop_constraints'] = options['negative_pop_constraints'] - # determine set of exchange reactions, which must have the form “s ->” + # ensure that all reactions are irreversible errors = [] + for rxn in self.reactions: + if rxn.reversible: + errors.append(rxn.id) + if errors: + rxn_ids = ', '.join(errors) + raise MultialgorithmError(f"DfbaSubmodel {self.id}: reactions are reversible: {rxn_ids}") + + # determine set of exchange reactions, each of which has only one participant + # TODO (APG): should we use the exchange reaction pattern from the config instead? self.exchange_rxns = set() for rxn in self.reactions: if len(rxn.participants) == 1: - part = rxn.participants[0] - if part.coefficient < 0: - self.exchange_rxns.add(rxn) - else: - errors.append(rxn.id) - if errors: - rxns = ', '.join(errors) - warnings.warn(f"exchange reaction(s) don't have the form 's ->': {rxns}") + self.exchange_rxns.add(rxn) # get the dfba objective's expression dfba_objective_id = f'dfba-obj-{id}' if dfba_objective_id not in dynamic_model.dynamic_dfba_objectives: # pragma: no cover - raise MultialgorithmError(f"DfbaSubmodel {self.id}: cannot find dynamic_dfba_objective " + raise MultialgorithmError(f"DfbaSubmodel '{self.id}': cannot find dynamic_dfba_objective " f"{dfba_objective_id}") self.dfba_obj_expr = dynamic_model.dynamic_dfba_objectives[dfba_objective_id].wc_lang_expression + # collect all wc_lang.DfbaObjReactions used by dfba_obj_expr + self._dfba_obj_reactions = {} + for rxn_cls in self.dfba_obj_expr.related_objects: + if issubclass(rxn_cls, wc_lang.DfbaObjReaction): + for rxn in self.dfba_obj_expr.related_objects[rxn_cls].values(): + self._dfba_obj_reactions[rxn.id] = rxn + # ensure that the dfba objective doesn't contain exchange rxns errors = [] for rxn_cls in self.dfba_obj_expr.related_objects: @@ -228,7 +236,7 @@ def set_up_dfba_submodel(self): """ Set up a dFBA submodel, by converting to a linear programming matrix Raises: - :obj:`MultiAlgorithmError`: if the ids in DfbaObjReactions and Reactions intersect + :obj:`MultiAlgorithmError`: if the ids in :obj:`DfbaObjReaction`\ s and :obj:`Reactions`\ s intersect """ self.set_up_continuous_time_submodel() @@ -236,11 +244,7 @@ def set_up_dfba_submodel(self): # raise an error if the ids in DfbaObjReactions and Reactions intersect reaction_ids = set([rxn.id for rxn in self.reactions]) - dfba_obj_reaction_ids = set() - for rxn_cls in self.dfba_obj_expr.related_objects: - for rxn in self.dfba_obj_expr.related_objects[rxn_cls].values(): - if isinstance(rxn, wc_lang.DfbaObjReaction): - dfba_obj_reaction_ids.add(rxn.id) + dfba_obj_reaction_ids = set(self._dfba_obj_reactions) if reaction_ids & dfba_obj_reaction_ids: raise MultialgorithmError(f"in model {self.dynamic_model.id} the ids in DfbaObjReactions " f"and Reactions intersect: {reaction_ids & dfba_obj_reaction_ids}") @@ -258,19 +262,15 @@ def set_up_dfba_submodel(self): self._conv_metabolite_matrices[part.species.id].append( conv_opt.LinearTerm(self._conv_variables[rxn.id], part.coefficient)) - self._dfba_obj_rxn_ids = [] self._dfba_obj_species = [] - for rxn_cls in self.dfba_obj_expr.related_objects: - for rxn_id, rxn in self.dfba_obj_expr.related_objects[rxn_cls].items(): - if rxn_id not in self._conv_variables: - self._dfba_obj_rxn_ids.append(rxn_id) - self._conv_variables[rxn.id] = conv_opt.Variable( - name=rxn.id, type=conv_opt.VariableType.continuous, lower_bound=0) - self._conv_model.variables.append(self._conv_variables[rxn.id]) - for part in rxn.dfba_obj_species: - self._dfba_obj_species.append(part) - self._conv_metabolite_matrices[part.species.id].append( - conv_opt.LinearTerm(self._conv_variables[rxn.id], part.value)) + for rxn in self._dfba_obj_reactions.values(): + self._conv_variables[rxn.id] = conv_opt.Variable( + name=rxn.id, type=conv_opt.VariableType.continuous, lower_bound=0) + self._conv_model.variables.append(self._conv_variables[rxn.id]) + for part in rxn.dfba_obj_species: + self._dfba_obj_species.append(part) + self._conv_metabolite_matrices[part.species.id].append( + conv_opt.LinearTerm(self._conv_variables[rxn.id], part.value)) # Set up the objective function dfba_obj_expr_objs = self.dfba_obj_expr.lin_coeffs @@ -400,12 +400,10 @@ def parse_neg_species_pop_constraint_id(neg_species_pop_constraint_id): return neg_species_pop_constraint_id[loc:] def initialize_neg_species_pop_constraints(self): - """ Make constraints that prevent species populations from going negative and span multiple reactions + """ Make constraints that prevent species populations from going negative - A separate constraint is made for each species that participates in more than - one exchange and/or dfba objective reactions. These constraints prevent the species population + A separate constraint is made for each species. These constraints prevent the species population from declining so quickly that it becomes negative in the next time step. - Also initializes `self._constrained_exchange_rxns`. Call this when a dFBA submodel is initialized. Do nothing if negative species population constraints are not being used. @@ -414,62 +412,47 @@ def initialize_neg_species_pop_constraints(self): :obj:`dict` of :obj:`str`: :obj:`conv_opt.Constraint`: a map from constraint id to constraints stored in `self._conv_model.constraints` """ - # TODO (APG): AVOID NEG. POPS.: either check that all reactions are irreversible or handle reversible reactions - self._constrained_exchange_rxns = set() if not self.dfba_solver_options['negative_pop_constraints']: return {} # make map from species to pseudo-reactions that use the species reactions_using_species = collections.defaultdict(set) - # add species from dFBA objective reactions - for rxn_cls in self.dfba_obj_expr.related_objects: - for rxn in self.dfba_obj_expr.related_objects[rxn_cls].values(): - # ignore wc_lang.Reactions, which cannot have a net species flux - if isinstance(rxn, wc_lang.DfbaObjReaction): - for species in self._get_species_and_stoichiometry(rxn): - reactions_using_species[species].add(rxn) - # add species in exchange reactions + # add dFBA objective reactions and their species + for rxn in self._dfba_obj_reactions.values(): + for species in self._get_species_and_stoichiometry(rxn): + reactions_using_species[species].add(rxn) + # add exchange reactions and their species for rxn in self.exchange_rxns: for species in self._get_species_and_stoichiometry(rxn): reactions_using_species[species].add(rxn) multi_reaction_constraints = {} + # TODO (APG): consider expressing the constraint in terms of ds/dt, instead of -ds/dt (consumption) for species, rxns in reactions_using_species.items(): - if 1 < len(rxns): - # create an expression for the species' net consumption rate, in molecules / sec - # net consumption = sum(-coef(species) * flux) for rxns that use species as a reactant - - # sum(coef(species) * flux) for rxns that use species as a product - # which can be simplified as -sum(coef * flux) for all rxns that use species - exchange_rxns_in_constr_expr = set() - constr_expr = [] - for rxn in rxns: - for rxn_species, net_coef in self._get_species_and_stoichiometry(rxn).items(): - if rxn_species == species: - net_consumption_coef = -net_coef - scaled_net_consumption_coef = net_consumption_coef - if rxn not in self.exchange_rxns: - scaled_net_consumption_coef = net_consumption_coef - constr_expr.append(conv_opt.LinearTerm(self._conv_variables[rxn.id], - scaled_net_consumption_coef)) - if rxn in self.exchange_rxns: - exchange_rxns_in_constr_expr.add(rxn) - - # optimization: only create a Constraint when the species can be consumed, - # which occurs when at least one of its net consumption coefficients is positive - species_can_be_consumed = any([0 < linear_term.coefficient for linear_term in constr_expr]) - if species_can_be_consumed: - # upper_bound will be set by bound_neg_species_pop_constraints() before solving FBA - constraint_id = DfbaSubmodel.gen_neg_species_pop_constraint_id(species.id) - constraint = conv_opt.Constraint(constr_expr, - name=constraint_id, - lower_bound=0.0, - upper_bound=None) - self._conv_model.constraints.append(constraint) - multi_reaction_constraints[constraint_id] = constraint - - # record the constrained exchange reactions - self._constrained_exchange_rxns |= exchange_rxns_in_constr_expr + # create an expression for the species' net consumption rate, in molecules / sec + # net consumption = sum(-coef(species) * flux) for rxns that use species as a reactant - + # sum(coef(species) * flux) for rxns that use species as a product + # which can be simplified as -sum(coef * flux) for all rxns that use species + constr_expr = [] + for rxn in rxns: + for rxn_species, net_coef in self._get_species_and_stoichiometry(rxn).items(): + if rxn_species == species: + scaled_net_consumption_coef = -net_coef + constr_expr.append(conv_opt.LinearTerm(self._conv_variables[rxn.id], + scaled_net_consumption_coef)) + + # optimization: only create a Constraint when the species can be consumed, + # which can only occur when some of its net consumption coefficients are positive + if any([0 < linear_term.coefficient for linear_term in constr_expr]): + # upper_bound will be set by bound_neg_species_pop_constraints() before solving FBA + constraint_id = DfbaSubmodel.gen_neg_species_pop_constraint_id(species.id) + constraint = conv_opt.Constraint(constr_expr, + name=constraint_id, + lower_bound=None, + upper_bound=None) + self._conv_model.constraints.append(constraint) + multi_reaction_constraints[constraint_id] = constraint return multi_reaction_constraints @@ -480,7 +463,6 @@ def set_up_optimizations(self): # pre-allocate dict of reaction fluxes self.reaction_fluxes = {rxn.id: None for rxn in self.reactions} - # TODO (APG): later: ensure that ids of Species and dFBA objective species cannot collide # initialize adjustments, the dict that will hold the species population change rates self.adjustments = {} for obj_species in self._dfba_obj_species: @@ -493,112 +475,48 @@ def set_up_optimizations(self): upper_bound=0.0, lower_bound=0.0)) def determine_bounds(self): - """ Determine the minimum and maximum bounds for each reaction - - Two actions: + """ Determine the minimum and maximum flux bounds for each reaction - 1. Bounds provided by the model and written to `self._reaction_bounds` - - 2. The bounds in exchange and demand reactions that do not participate in a negative - species population constraint are further bounded so they do not let a reaction's - species population decline so quickly that it becomes negative in the next time step. + Bounds provided by rate laws or flux bound constants in the model are written + to `self._reaction_bounds`. """ - # Determine all the reaction bounds flux_comp_id = self.dfba_solver_options['flux_bounds_volumetric_compartment_id'] - # TODO (APG): since 'flux_comp_volume' is used, cached volumes must be invalidated before running dFBA + # TODO (APG): since volumes are used, cached volumes must be invalidated before running dFBA if flux_comp_id == self.FLUX_BOUNDS_VOLUMETRIC_COMPARTMENT_ID: flux_comp_volume = self.dynamic_model.cell_volume() else: flux_comp_volume = self.dynamic_model.dynamic_compartments[flux_comp_id].volume() + self._reaction_bounds = {} for reaction in self.reactions: - # Set the bounds of reactions with measured kinetic constants + # defaults + # the default minimum constraint of an irreversible reaction is 0 + min_constr = 0. + # None indicates no maximum constraint + max_constr = None + + # if a rate law is available, use it to compute a max bound if reaction.rate_laws: - for_ratelaw = reaction.rate_laws.get_one(direction=wc_lang.RateLawDirection.forward) - if for_ratelaw: - # TODO (APG): if reversible rxns are split, call self.calc_reaction_rate(reaction) - rate = self.dynamic_model.dynamic_rate_laws[for_ratelaw.id].eval(self.time) - max_constr = rate - else: - max_constr = None - - rev_ratelaw = reaction.rate_laws.get_one(direction=wc_lang.RateLawDirection.backward) - if rev_ratelaw: - rate = self.dynamic_model.dynamic_rate_laws[rev_ratelaw.id].eval(self.time) - min_constr = -rate - elif reaction.reversible: - min_constr = None - else: - min_constr = 0. - - # Set the bounds of exchange and demand reactions + max_constr = self.calc_reaction_rate(reaction) + + # otherwise use the fixed bounds elif reaction.flux_bounds: rxn_bounds = reaction.flux_bounds - if rxn_bounds.min: - if math.isnan(rxn_bounds.min): - min_constr = None - else: - min_constr = rxn_bounds.min * flux_comp_volume * scipy.constants.Avogadro - else: - min_constr = rxn_bounds.min - if rxn_bounds.max: - if math.isnan(rxn_bounds.max): - max_constr = None - else: - max_constr = rxn_bounds.max * flux_comp_volume * scipy.constants.Avogadro - else: - max_constr = rxn_bounds.max - - # Set other reactions to be unbounded - else: - max_constr = None - if reaction.reversible: - min_constr = None - else: - min_constr = 0. + if isinstance(rxn_bounds.min, (int, float)) and 0 <= rxn_bounds.min: + min_constr = rxn_bounds.min * flux_comp_volume * scipy.constants.Avogadro - self._reaction_bounds[reaction.id] = (min_constr, max_constr) + if isinstance(rxn_bounds.max, (int, float)) and not math.isnan(rxn_bounds.max): + max_constr = rxn_bounds.max * flux_comp_volume * scipy.constants.Avogadro - if self.dfba_solver_options['negative_pop_constraints']: - # bound each exchange reaction to prevent it from consuming its participant species - # so quickly that its population becomes negative in the next time step - # if exchange reaction r transports species x[c] out of c - # then bound the maximum flux of r, flux(r) (1/s) <= p(x, t) / (-s(r, x) * dt) - # where t = self.time, p(x, t) = population of x at time t, - # dt = self.time_step, and s(r, x) = stoichiometry of x in r - for reaction in self.exchange_rxns: - # to avoid over-constraining them, do not further bound exchange reactions that - # participate in constraints in `_multi_reaction_constraints` - if reaction in self._constrained_exchange_rxns: - continue - - max_flux_rxn = float('inf') - for participant in reaction.participants: - # 'participant.coefficient < 0' determines whether the participant is a reactant - if participant.coefficient < 0: - species_id = participant.species.gen_id() - species_pop = self.local_species_population.read_one(self.time, species_id) - max_flux_species = species_pop / (-participant.coefficient * self.time_step) - max_flux_rxn = min(max_flux_rxn, max_flux_species) - - if max_flux_rxn != float('inf'): - max_flux_rxn = max_flux_rxn - - min_constr, max_constr = self._reaction_bounds[reaction.id] - if max_constr is None: - max_constr = max_flux_rxn - else: - max_constr = min(max_constr, max_flux_rxn) - - self._reaction_bounds[reaction.id] = (min_constr, max_constr) + self._reaction_bounds[reaction.id] = (min_constr, max_constr) def bound_neg_species_pop_constraints(self): """ Update bounds in the negative species population constraints that span multiple reactions - Update the bounds in each of the constraints in `self._multi_reaction_constraints` that - prevent some species from having a negative species population in the next time step. + Update the bounds in each constraint in `self._multi_reaction_constraints` that + prevents some species from having a negative species population in the next time step. Call this before each run of the FBA solver. """ @@ -689,6 +607,7 @@ def scale_conv_opt_model(self, conv_opt_model, copy_model=True, # scale bounds in constraints; bound values of 0 are unchanged for constraint in conv_opt_model.constraints: + # this 'if' will always be true for properly formed constraints if isinstance(constraint.lower_bound, (int, float)): constraint.lower_bound *= dfba_bound_scale_factor if isinstance(constraint.upper_bound, (int, float)): @@ -784,7 +703,7 @@ def handle_RunFba_msg(self, event): class ObjToRow(object): - # TODO (APG): later: cleanup; unittests; docstrings; etc. + # TODO (APG): later: cleanup; unittests; docstrings; move to conv_opt, etc. def __init__(self, col_widths, headers, attrs): self.col_widths = col_widths @@ -816,7 +735,7 @@ def obj_as_row(self, obj): class ShowConvOptElements(object): - # TODO (APG): later: cleanup; unittests; docstrings; etc. + # TODO (APG): later: cleanup; unittests; docstrings; move to conv_opt, etc. @staticmethod def show_conv_opt_variable(header=False, variable=None): diff --git a/wc_sim/submodels/dynamic_submodel.py b/wc_sim/submodels/dynamic_submodel.py index 48aacc95..e73bf178 100644 --- a/wc_sim/submodels/dynamic_submodel.py +++ b/wc_sim/submodels/dynamic_submodel.py @@ -111,7 +111,7 @@ def get_num_submodels(self): return self.dynamic_model.get_num_submodels() def calc_reaction_rate(self, reaction): - """ Calculate a reaction's current rate + """ Calculate an irreversible reaction's current rate The rate is computed by eval'ing the reaction's rate law, with species populations obtained from the simulations's :obj:`LocalSpeciesPopulation`.