From b9bd83620d3a2000fd38529c5c7c5791733bce06 Mon Sep 17 00:00:00 2001 From: ScottCampit Date: Tue, 2 Jun 2020 11:05:53 -0400 Subject: [PATCH] Added quadratic programming example for Gurobi in Python --- .../inspectionProfiles/profiles_settings.xml | 6 + .idea/metabolicmodeling.iml | 8 + .idea/misc.xml | 6 + .idea/modules.xml | 8 + .idea/rSettings.xml | 6 + .idea/workspace.xml | 4 + MATLAB/DFA/src/DFA.mlx | Bin 5822 -> 5843 bytes MATLAB/DFA/tutorial/dfa_tutorial.mlx | Bin 7897 -> 20088 bytes MATLAB/expression-constrainers/CFR.m | 245 + MATLAB/expression-constrainers/CFR.mlx | Bin 7125 -> 7254 bytes MATLAB/expression-constrainers/CFR2.m | 136 + .../addMediumConstraints.m | 105 + MATLAB/expression-constrainers/knockOut.mlx | Bin 4712 -> 4629 bytes .../expression-constrainers/perturbMedium.mlx | Bin 5437 -> 5138 bytes .../QP_example-checkpoint.ipynb | 314 + .../Untitled-checkpoint.ipynb | 307 + PythonCOBRA/QP_example.html | 16694 ++++++++++++++++ PythonCOBRA/QP_example.ipynb | 314 + PythonCOBRA/Untitled.ipynb | 307 + PythonCOBRA/constrain_flux_regulation.py | 46 +- PythonCOBRA/knockOut.py | 59 + 21 files changed, 18556 insertions(+), 9 deletions(-) create mode 100644 .idea/inspectionProfiles/profiles_settings.xml create mode 100644 .idea/metabolicmodeling.iml create mode 100644 .idea/misc.xml create mode 100644 .idea/modules.xml create mode 100644 .idea/rSettings.xml create mode 100644 .idea/workspace.xml create mode 100644 MATLAB/expression-constrainers/CFR.m create mode 100644 MATLAB/expression-constrainers/CFR2.m create mode 100644 MATLAB/expression-constrainers/addMediumConstraints.m create mode 100644 PythonCOBRA/.ipynb_checkpoints/QP_example-checkpoint.ipynb create mode 100644 PythonCOBRA/.ipynb_checkpoints/Untitled-checkpoint.ipynb create mode 100644 PythonCOBRA/QP_example.html create mode 100644 PythonCOBRA/QP_example.ipynb create mode 100644 PythonCOBRA/Untitled.ipynb create mode 100644 PythonCOBRA/knockOut.py diff --git a/.idea/inspectionProfiles/profiles_settings.xml b/.idea/inspectionProfiles/profiles_settings.xml new file mode 100644 index 0000000..105ce2d --- /dev/null +++ b/.idea/inspectionProfiles/profiles_settings.xml @@ -0,0 +1,6 @@ + + + + \ No newline at end of file diff --git a/.idea/metabolicmodeling.iml b/.idea/metabolicmodeling.iml new file mode 100644 index 0000000..d0876a7 --- /dev/null +++ b/.idea/metabolicmodeling.iml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 0000000..28a804d --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 0000000..ad10015 --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/rSettings.xml b/.idea/rSettings.xml new file mode 100644 index 0000000..97e6f1f --- /dev/null +++ b/.idea/rSettings.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/.idea/workspace.xml b/.idea/workspace.xml new file mode 100644 index 0000000..db2a81c --- /dev/null +++ b/.idea/workspace.xml @@ -0,0 +1,4 @@ + + + + \ No newline at end of file diff --git a/MATLAB/DFA/src/DFA.mlx b/MATLAB/DFA/src/DFA.mlx index e28863d7c76a615a5b108e519adca40ff1f1c057..cd5073996ddb4c06ae893edb60b6dde3dd844e2b 100644 GIT binary patch delta 2768 zcmZ8jX*3iJ8y!Pqodq*wFT1htWXryL$&xLy?-VHwnIcOQGL1dUjD08B%GfD}1|dt9 zl*3Nr{iv%flnX)+y5{H3&b4H$9y$7$U%+VhL6`1qUU^s<0ADo>2EP zCG3r4u4;>ffKZQcZ1|z~KaPCSV!k~+;v7_3!jAjHWJ+uj9XT?4d6|YO&q(SfUQ<}i zR$;6z#-h6@(7xHP^h$-5FFB)6-3f*ib@eYbmO#i#P)Vm{zpH8h}VQKT+cb)+}sqK<+L6+NHYIVWcs@NPpsuq%Ba+26LVDQ=` z$-|N*Of^8zwp-gcEQY&;0jV}4VrPQ0Pa1*C4D58{25?Ic(cNN>dMI*4(Wq+TLGkk= zD`oX8Vt-^nhvP+{wulM0o>^J$r$Sjb_R`liPS%A^tFBI%IUjJmhv8C3&SX`sfktQa zu&S(;Sxs>Y!WMJw#E0)*I&00Od!e;fnsE-ts25O@MUq7>Icp)`*3y1v~GUIa(3uczfmIQW@q$|h{$j2>)jWQrTSu|ED-q)r32@okYOGHrQF3j zRVCNKUC_}8)e#3sQ8jX2Q$s@|l4kCi?nZ{$`+S5b5l$gzh1-!ob^_*4*(SMCIaWv!^7_bPA~+Q znnM85-hkTy(!LGSq8-}4D-lg5A5Q(de+IL@4N{hg$?~3qvt|iPYtGp1_DA=&E`GSp zG0nY?)}sE-V?JJI3Hs3?HWOn@Rb~eol$N;L4`0v9RMdg!wF=7J7QE`H|2Em3^=hL~ zv&3bDUZ0V2VUGZTF{8e7h$tn3UTIEka3u1IAEo(R_b9nT_wQYvYb#kL-0S`s`8#{QBG)> zfuz3$QG(VMxc4PU-U;USK76hteLJ0BoE-^hYEYGzf9gNL|8lnArfFHT0s+2egPA062J2L2E*TR#UxPv~(k4;{L0gUn`gQA?akU)GbdmK1Ub&le124 z>%^SC(^Ca|gZF2YEZTv_>rsaKfaBFwpUm-qulEM!uV7tC6~KHv7Rw;@maT};g6%V0 zTya3j&&h6Z5b1x<0^PkAaifOq-a8|oX*nTWn$D3sEV^|UeT`C4RyAni@^pTz8D~=_ z$NMq~viAt`f|2jL?c&5cZG|?YbQ^wN^2qM>>D3^rY&q@lGn>aG;?h&(zOoGdFIEK_ z`g5!k)qma9M^c^6E~L%X_{7YGnnm%wbR)9dD9k9jzQecC_ud7`6Bzn=T1)4P)G9RP zUC9QpaWdz_yMoH7iNg*whbZAbI4l75;B$od%D{t}!dQ&Gf9>3>iJMRhBuj2%TyI)k zeY|0iHL>!QUbk3!NJHX{TpAi3;VFli>nC!hF8zB7&s>|GPHu`_M0G~CP2UnbtWI-Z zC9v#%eC(Z%xN?xS6kPkNyuDSftQw^_-M0p7p%wV9SVgVfoD6u`#YQztL1%p(A!eAz0Kga zU2u-ExD}^yCD$Pt@!YVvVr)0t$As5VCgl!oxc&<0VROJCUR(FgD}_rMI5qWe=r7S~ zo8IGR`P<}0)4K}yE7nzS|huz7O6e9>pk zA;H}3Rc03(?V1Ong>5Cjs?2ovTP<#x8co~tz}k85T@`5~yZ89r=#DY78DO%C?wsnA z)njars*4u#y;-Hq_Ox!ux?O)uOKA*~oqpA(y;UN7^4Wu278hRn&OvEB^880xOEizf zL1g*q*XwlDdeybmwdA44;NeMEQ2UOToQVCXgPsYlK6|1q5E=C?x^0u)g-#0!6-iN> z(f$fBYJ7m!-1Jg$f?Q8A{Gy=tq+3~fPb@=J} z$f2bb`djwR7q=Hb#|zh(uuB*o4Uef~CN1s^Dy!6?tod`8)9z#M6`45|PNDJ^3j1i? zIL8h^CXKXDp4!?P4SDT;oH<}~w)}__mCHAEDh#5_@Ywpo*7(&y(2PmBeVo7VUDGt0 z;aJSED=AgzJxos^ri<3i-$&pQl1o*5J9E=OTR`>}zMJNWpTXXvN#vHK^LAE`{H6Qj zzmV#}@-M&nm&Usd#78ZEI)!&?X>=r2m6EtEgs5!$t=^VUGvkd5@1H@GghGxfy6pxC zfuv&3-{F#f9oDnHU)04;5h0{>&XmK^9P8K2={gCylH?>D9{IIE#7^slBr^mE8B*HzOffHesOKj+c8k!@`m-r~~au7J7% zd;FYuSouCRqE7(GGGjy#x}Xpg9Ng?YIIP(Mmu>3Q@LT(Q>f|bT$zOH8@v!H?A zVQaw{CGsnwoJs0w^N#Ab`4j6xNtLFVIeCFz-m28~v+J)d>GRhjYlh}e(W5xEA_sBq zZiu26@lxRkJ6D}}fQh+JZc{&L#&qXkNF)F_AU%a2?pVV0Lq_A?^4hfi8}X(cLL5bz zHE?&HUm^Qkisi-V5a_1qv0+wak6*45)llG>{NfhH;=6ebAi{D#n3z}OAVTun)y?jP*x*zz+f``s@PTb5d=u}@IE>-C8Eg9B#w z17Ai8j@s7jTafEU#=ciTrF9SwqfE-zW3fet7EzuAw6$=f5?`z7&TxR3{PXDR9Wg1n z)RH4+_@S6+{99dU5wtJL#;(u9`QtWVs>XNgtBx7F^y^@@Sp+x&lHT?;Dd*OCDXZdB zQ5-!3M@At$(;qM6`oBP3o{Gr%*%$a?ea9^m=v~GiMv29OpZ&O^=aT%iAWPY@p7f~c zL;ja_jpUE3p1g7cL%nDF#!hZULJ*xTwxASf@#4X``NEpJZ9lp}RQ54&G)8Q}N~#hxFDHq~?IWn|i~9*hieV!+ z&JYx*Mk~mHB{sd~NPmKPij_jivXRY=(fz-S9brapwK7C>ee*%v!{CU4pG+d4j%8f)B(de`-R{h-qWhe%H7)d8|7T+yLrca{Da zFJt!@EWT;mguKda_;Nr0o1>MWTWH@JoM_$dG^!p~oW=#$Xdb+onGe zV{>Xn?{kwf*c`jfAoeK1s`7;5O9LW2IMB3Vca^p61nqIq*N948noQSLKSv1nD`<}QAF@C-+IWg4J5k$g(q@`O{B6}e@0m=P(^vt>CX28y1Te#cuo>T(#R>^Pw7in zs_%M_i;TkaB}8tdJ+yo$^Pa-@Ehuetc%io$ZWS0e>b#iu8F8)j@a?nMZ!cczor-CI zN!~Y(VVM@ic+dpI5CQW@=TzbOMm1cmqkP1905oxQ*~~(#(m_jmD}E z6A!|x+2UGVF)g0qRdW?ndEs$v6lr!-;ae{tsqY=ovef++>ne>K@=$C4_@I7?I&%Fa z3J+tC@YwAPe2tcAAnz<}71R?lFV$wt(Y0zf{cS99412!$V;z2<^lp6a+7ug3UnsS< zwJi)7Py>AQ8VWb0JW{#}$6a4Fge|{xdh7{!Nt!fnpVhWppEIUUj|-T*13(sLoLr(H z5QrD_%fVuW9c}G-K%lkiN*+@#Kr>{B^*nojDe{wGT9u?U`zDZIZwQf)C_l;GxybkG z=9y$4q3tjYz>dp>l8(XR@yl)*gSULh9hz3jj-4goy!_;3{+HeGf=6}9ttK-|DiIb6DK5T4P|+v2 z5G6Fn?)Y@-xfjkIn!Wrb z*@$o2DV7aYb@au0|4Iya0)rPvzRH8(fX+rMlbo zKfS!DB>m;`L!hp&Yr1M0Ur!aEFJZ_2skP8KshgA=B#8EXr|3bz*dDdUvsaqT)TS@( z$xiv_u6ZKI2wAH+|1pDIew|PG*kW&pLLR@gpT<(nouh zvrxQdE5sF;HaF<+J?HOo8?Q@QJSM#+Srfi8T=nfjkZu>}AKv(h3E056z5pESs}Ndv zrdr~}jr%0)f*H-;R7~P^#Jb?>z3wBto}+97YD{9FdjR=(Q1H@0=icp+9>#&*#9pFf z=8l#?WCMSy$@aEyM3-0k>w(opd>1miyH;*Tr7{C>V5nxEkkE5G|2jPV!kM97Ic71% zo^#PC$ZT&pGiJdz_kDU&MWwDZS#P`S((T-~JA<>7Ls`X1;`H4c<#(N9;p8pCo@(Fx zm2$MaW5(-{xG{i4hiRTDFDCoit^SDz>!6qzkqFuh#zvjxoWHMYyDDa}IY^Ei)?*Ea zgyx@}x*F@MN*$y3WfPCym_$D%|8OZH^_fEgVZtUNu1YA{hK^JoBJyy!bTLvA zyP8_6k%%9~GB^3EqRj9G7vmG15ewY$yv815qGdeZM&>3nPJqI~kIH{*Kdi%ge#KKyNsHzvsi#7utG>&P6oLEPW!#s{U$t!5mh))F=SjvtF+G>e&$B0$v`VL`AUFp@o5|I>j~IY1yq5KSL0 zNh^axfWSmE2t{4!&id}Xw9rSQSaqev_8x%^UwPOS8K!w$`$M%f6K3#4mbYtA$mcct z;70?SeXpM_TS!)@#Z;}WSBC2{kI<@PJKDmWRN)QVv!INZ+RYd zbBsR3d&j&dVW}%^@@w=Vme@IV?wd!6%bNOy@n;*Jc*=SZP)j|N4z*6GRl{M(mp^9L z%D1g_=-3T@WiIPHk#>#Z;tJx|YG}l)W&L+|_5Plk=R51m5SpQA%>O-v!=k5=|1wx< z5D5Bn=wJ9b_fMqI&@frfpBx$vCdna1y9JX+{MsJq&$JL21d{xl^*{C;gmzd+iuMME m2CE4Ftr-Na%ArF04X%Uux1j%&{nJAQ`d@jG+$z6JAM_tZj2F-V diff --git a/MATLAB/DFA/tutorial/dfa_tutorial.mlx b/MATLAB/DFA/tutorial/dfa_tutorial.mlx index 1675d6c03fc7080b5301b35d309410e72edc6f4a..3736e8acf11c9dd65d85b9fdafad346f154170e5 100644 GIT binary patch delta 14920 zcmZ9zV{oTW)b|-XGnqJGR^8 z>aOnU{(M(?pEJN~8X-_GLn`I1$6xRrP{6>Z6~MsIz~UM4G3qxUV2}ZcC;Z=@6=IXt z8YN|z(!VH;FE-S-b+@Fwf_^G!4^#eOP$LVN_X2Rlyl38@lCeh<3+VD~r+rXPS<#?K zs>!mGx3hu~8GW%V)lL$Zm8a)1=c`Pb^_bPeBz6x?OS{+iEl&n_m9NdmwUXNQ=eJ+RMz1S?Nr2V6T~3*iT1JRRX8 z_eAet9QOFV{oi@a&mtjHn2N#vjvq*h|9S2fF5KdP(0yZ5)b>UKx;>`aBnH_EYEY^7QXqbBc>+cbPUg~=vYGtQItLYYiHE68bOnZ~8*mKPh zTM?vVW)l{vi27#Whu(MH!KfcMnG`!+@@&*< z6bphPVDAOCAPK^)Yp|Kn9i2Cie~kfNYMWORXI4e72n9OS5G!=?8j0t-lS?1OidPNa zBp+nkKP-pcyH1?}0=Q~%m>jS#yzr8c?SF>;kM~0nn`ezn2%<+1*5S}{OFy#{NsAa< zXyn5$VTsU_Y^fS8T2(8XbQfVWAN5OBEu3*P6^nIjlDL{1TdyQJ!?7&wB5fztDjkvv zsgFmwQ)&ebSokK{@>vgJP`7=z-n>&;%3upcJ2hbTZA)PRQTtx911vo#?>CW6{mTmN zP5!>oDEfR7U8vaC@do25UF&!25LR|vB*Cxg+56e(g^j zJKAt?SCB`0i}3_OC^v`QUq4pdX^G%5XACh*h@ju#g*4jo=-gEmNVHufr;1}%%X$QP zqfuQcj6ACWBag3KG&5Vl`a+t0QX~`292#(^qi$y~FUgAsLYj+`GiVxpN~K(wNyA7(y=^@ydB69zEWq5W zn<`5ACvlk@>_1c6J{dW)MT+}5-1V+Y)PB4ELAUG#N?LNeUn)zwVK3XE!z|jD{0e#h zC8YgV%j#aYV5!2W)U4t+{L`4qHRVCxT17cw@B*jfBWtxLV}IM)KaJ^^WUd8-y%WU{ z50H$$3%U7QCkeQXL$5R7`iGz~bvav-a>@PZ_T&b$>h3?aOmFr>%U_W4Qjk#SU|?Wy zV3#34rThb*QGhNB7?=qjIM|Q>%vSo2mih+tHcpPVPL6c0R+blD9uAu9Nypna#VvYM}GWy8+-(7dD??L9SguC zRdCO5(1fr9!^6v9C)oICW7ys>o0;eJn~Q^%Z944ZO}&HmcD3rS0*}C{ zot4FKYW)~f;a}^+%Wn%;mzE8i5)L7T=k&?(&5N%WB6p)7p?4J(~mi1_SZ=CN{W!@%AnIoWh@zayZqdl!^4~0M2`NrjK8}S|g)Ah#}(pS_J zF7DC1>D%J00VB6gN2w@AwVRi=3&H%C1#lHW*GB@#)<+O4s|7uNeJoB_{_GIiRTwAr z4d;u>d}l?uI(~eT{m1nAxvKWb`Z@c-|M8Z*`br;vYPs6+e(<^K?saoBH3jVT)j!|Z zGkLRWXKNEg7l&uc>kOxEeI=`%b))yaUb>KV3ADenzr3^(72Q~s`_~^Gh{d&9yw>IU zmYtNP&AH2(Fl%hwxSoeH+LQ%X7IqB=8v@W}I}UeJ%9K7@tR-}_w`W^8xbbYuk~NtJo=vg@=;>w6 zjZ16Vc=&WPbf(f&9&EitA+tR=Gu~H}<4=q0Pbyw!T!OqVva^D8OzIidla;oeZ|r@{ zt8;Ej+&D=k2|Ok9wUiaYcxTE$YG)#i?jL3TJQ}h` zbbUtWg+zTycScv6R!}OfNrORwg05V0c4OnxCZ30+7hdG2DzCTw z^I*8YO_t~5N4U*nbk}kbbj3ACnnWJ#N5 zMQbFe`jTtjkI=DG=7?p=^mv>W^nP<<9lh)(p<6|^MNs{tF<`$PL|L=UR-fp`Wt`FdmuLnP0 zAMWj*RaCwR0I6NIhlhBZOD~4D#B2{jVUQW?11y?G`m#)M%f|?Y$1uHxoZlOf_>o7G zF8x#}j}?~C#hJ_d?ejF9_c`t482r(}QPaIb>D#n*>@;}b^KK>ma z0|K|}>@6u3s-^zF!)%&tUTj#&sg-B=zErl~^qg;q;2bO&{_4ocZf!w*t;g->LsMRR z?(ihCq0dr7x*hxd9*f|y&Jf1IYR8>0rjC_w|eA-QS%OZgGdv;J_b&&yxc6>^lvuC+J-JF!&=1sS^U zn88$7a&MI~Vjyb`+d5HF{4B*?dw7w{g{jGqo~pGc%X!)Fj^@A1MSF z)wTF-oAAIz7L7^rvBT?#kQ4NFuYx$+ z^NrA%{KkXy?77OiH_y<46x~rZ&(d!Mq`&u$c54%Tj*1t5{B54Te`{e9@>%xj(j)B& z-o`)YuRn}o|EGd1^vzFa`2yzpZo0L-YW?@j-}X(gEw#BsAxxsT^S9|hGhXT|b6&0xyXca^>~~kT*Td%1axOQR#>&|;Y;^a5@j_&dvC(s!(QvE2 zFLppj(KMk zmTqt}^KP4n#b9852Eqje7*rM<{w1X{?GX0>`;BjU;)Hg-q*0j zK&>Nz@YJVuAv?60J@4oOH%G`BdwuTbcietHICn^B$_EBt*xyoC*J{Dr+Qlxxf5O|n zURUFpl7IUpJ-$5@SICn=n3Y<)ezu|(ez?ytfw(z5?MKc~i8ll>6?e;l;&TmXPB0D0 zfSnHQek7q7W8q)RGP?V5-vlt@c4nKUF<+y=8*0|+v#AZEt-&-V^uBUmoKEe)ZpjRW z!tp*a8*28w+#;7*kTfI%hjn?l->X9Y5*g)6aC_Q2X8JROfx^z!$`+7sV%V>(wOU>F zw`Q+k?%83NPjAX&23Poi)}bRL!7UOtoF(}F7!e&eqwAj`GOa-Ivu(G3OkwQk59<18 z@H&or-sjsh`;j@F?!!}<7-J|``{$4S;%(loDJWV!obR3=rIeEtt}}vG*b|$qxF&FQ zEMf0)ce`K~D0LB4e7cUATeyBSEPJ}l-tzx)e;gpAgLemO zn6xYPHd%a4%M~jOMx_{lbLCza1n0uxBQv8LJ^gLAZ*R4KH~7Emg^R1`Ui?41(Sa9) z$dA#1eT0bXfDX|+yoYJauG|X$f#Qa5>IYi>+jP0xmxREVS&a)=AU8uJpg6+~_SQbX z6D_jJF92Z7Y^MXdfD6PsygRwxS5M9jg6TI;CeRcXHDGMHZ>pDaSu6`go!V(LPtp$O zl{o3<@jd!Iu@eb{`d*z zNi8B%+E9ZUOTjwKeATcD=b2jlkobJxc}X~xV!)hP$VfsMg*DzM3r!**!FCfxKUN&s zIAa=WLX)9isrl2=LN_LW%NE(rbk0(_lo?a8XV?!B9B}~O)EuLwo^DstKOc$N;1o(C zzr_2`Wr3ivQ9zJHO9>XyRsqOOmsF!F9Je{X1Z+4Fl2wCyndweEI#7~JHInm83S8bc z5fB&8E7D5yt&cr{oMXW^Mk6uA8OK4!6LUbxF-lyAGL~YE}{=>R_Mt|lYl0h@4lZYMdri|z3y7W?)@`} zE#>epeNFajAg&4G9=Se~@$VuF?A;%&UO=EfJ_)i>6x|6ySQPaN^iP7(FvkdEB7LUu z*JWA?FQpnCx?I0)sX+Kbp9VUEhDs(&p8Z^|BBLP5&#MT#9^;G(Bt0o8NOz{a(W#s! zBLPGb?plk7noy=?3}{n0t@?+a1N)xKp+&qoLRYViBSuF zMZMxcOX<@7ce|mZJ4yR;ed7ogFo{p>G+-{<|22$98*zoHI7I?_awfnbU6h0umYAy0tqyujlPHC0nsX7(v4JDHh`vmiQvYEj zcv@202R*KVTw+i1ukG~d)JCJo@SIs?@sAQRN$A#>#vr~H>3LFq9DwCFya2L4 z_q~pUsbG$P#Sqz9rLt*ozpE3bNz78y$N>rM9~S$05erGKO`keaIW|iAx(8GU=D(Kq zP*Cim$SH&Mf!%D~9Tr&C%;82NEPb;*EUmIWf(jbR1wPNxi%ac=O)^y_`OvIsx!|5 z_0{lXdW|AdIeZB#yac4+q^M;#F&d|VA&qr(I|ex&VO>q^fm_+vwKfey1(N8XFm!tg znk6!7XOr6gl;5Bb(;wuhXak|etAPfz`dp3UOGE7}9~^d%CpMLw!2l6wVH{+M;fVB5 z8}7+;s~xCPiYkyLA56VYmFOm^vb}~mPn?bxNocA_R{7A7Nf1?1k+I}H);=xc8UYoa zm^ChABx#~;7Xl=M5XD@)5-SkOvjg>4f+=nwcGx&qcoz#+9C~a$b**fHQsm~_@=9Y+ zHErAgILrYnsK=+75s;OF))0vjiFJ|=kHG(F?nDL3OiHo+r>w7xHAfOH3BPKpWrl>S zP=LA$2>Yxe@fzqHSebHVg+nKmy`pfUiubrnFb%^0l-NZm8I8sDSKUQ-$4NCUr>&h(ebt?i}NCB^J zY2yLKX;)OeYnfuY5k-B4)>ilZ^T}|9Q^FZ0p=$lY+cUEN?9sg5;+Ok~n~Y7-{(VfE z>UKGgUP>Cc4uDBcc@ts{t{|brp;PE7#yu(Gye5`8O|R)!1Q8J@LMjnoMBkAe!!kq2 za@a~y90kXn9ZDHE)l$jU66n{0zk87M=~9j*&?E201@t#938of zrGRnBuYpfdMnvX@^IKktSr?W=ScjeycA<qytVAFCkvB)iP4A^TF8t)sG_s)|`T;2hEf*3KvjYB?J z5<%oJNvWQPJYu;~5LI#7M8DAsd`YNY#84BX0w|&#k@QZTuqI+H+ea|Bz>W$O^*RfB z?#ZD8N{$-9sT?$yS14hf$#qcwuor6+{z4@m1^7xS+q83=5R1Y83LwH(7`jon#5y2Z z7}}pJB#Wyf*bfEz$IVMd>=laGA&i5_j1MLdn(D zMSPJNYJNR@o29cOniJuYwlOEtMEI6h0Q~VabwjV-^nrdS`?2O`JVQHM_4FzfYwRY5 ziZ@h6`4WPWvT#M_hnc9zgtutI-4td})=+b>TzioEikGpZegJ_rz%Zc)Z!kH)5YfW_ zAeBY#&T-J-sCl1Du~s>hWM*}~lBAWn;B$s%lY1~be1kD&VnsS5d=K9tue@d?3dB1k zSv0b>aE^Pnp}7+cqs{m89mDja-sBn@@2Z0-1oTfdfQ_MD^>5XZ;6s)7r9zPrvm-ln z6sXos%INz>w*06;6#u~_*&8L*Fxjy`In!~zcxjw!Cz zNIQYViChr6FfyFCS(#W~OwF<8fH+(7{+*P|3wUT~rrNe_l&DEdWGGdur84T?K5W#p z0o5S$mB%{|ikcI15(|0X^cnAZY>udSe+ce{0)!Z=k#xnHvoVtMam1=WZGlE1QK~xx z@r8T#Jv82L|I%m6M=kB^ic;*z6HJoH?_E}7)61C9N#ge|{d8HAaKluLf%G5c-17ZO zrE%!JXlC3hY^u0aZVj+a^h;Ir?&^tplx8ko`~;BHxrI2?|HLB-#cPN;O_+WrIFuY3 zlIJW5<7q7`tpqrGDGk|eS!P@prds?8`?mSm`h%~D;YOk2!ugSO6k55gc8}vmWA|6m z19E?gw?zk`$eps(2LgBx5inPg-tzPBKMRr}+?bN{mL>*=vZ>O69BMfhu*cOTOBE4W z_+Gs_<%88dj`ZJyK8+J5glXa>3F(Rp}3kwpxuJe_&2mQ6U^(A#~IROgCN zoeF$mJ=sQ0D#%wG@yq{&R#A?k-(n(4&$tXgI5g6FB_T+Xu<=k!&;X3_krWO~)8^BT zLB7Jf4LH1@0tq!rCUt~;R_S%faV<-q-vcXP%`J>_g&88eR`d2~=E#z)$Or)@tHPZl z;g%*jJPRr7%3|e+WkS$Vr;=Ps$UJJuWVdSU0aGT02IEPl;u|0;5ms14WEvRmF>t}lsJ_!XKVo{Jo7Om7w-A@oKL}; zeeQ&7*0X{MUN9AZUVKFas=8@h)V&4@N56>GSWy1p;of9q-XLH8f^2aRvm^mPF&EdO ze4DNfRyK_LN1Fl6$zEsF4aIhBE|{cB*BpDmM$DF@VxNfl2{1RqKvMn*m&ynmQ1b12 ztvd#sjMuO$g1A_#ytW*Y*db1ZE@ix^O6-?1ES1%g6>?dkc+oFD^i_BU1I#KWC^UqE zS7xS^9g};HaI!iv!DZZ%hx|SUDue9*av|1;RAQ~g>;o6MG`HWRsHxQs`c+4&m;2m! zK>jI@YJ2|;5z$6(M#!rY6=KH}HBy%<>J@$)g5*LaBM$T#oD<+Qr8!6E zsO3vuzVc}}@9^qz0ApXpe(;KQ{2lD=Kujvc(yxfh9Yx12W<}(ZhNbr9Zq(EFE^DgNf5Z{_|@EG`=xMa_9!3{1Sumhtz@Z>fWaj%^&;P-h%sXRei? z^Jy4olDMBHvpb4hzJ-ycZ4)MbvX&w;LVKwdM#)EnjfyjEoJL%EJb@z5;xYkTU5*KdW z(ZH_J=>&~^BfTCi>6qr>|@TDu6-DbA2%1ufOXUU)RqE2hsOK!5e!EhKuP zCh`UtXK-TweLnX{>BHl9Vsn!XB7BQR1OYzZW+E|~_RNw!#pA`(;A5synDFp6jehLU z9_Sh;X+!+HvAh705fHFC4%||{=+xsJUHVVV#HPL#z~DuK*K>d}X_AKHonD|Lr&Lw& z6_Ci3g_9WN8GBCagOP~!2}V~-ZS$xr)f)w}c=S&@jY_RsIV!Zky#LYOs~n zjOY^AZi%C+5y&1vDO=@ivRDh01$2(4gdd;v=DEx=Cnw>2I?=6yVwN=|IAY+GSr}i% zk>baZP1HMY1C+C)hVrNH`xp93&tM#MVa^z&4)2aqwoL3sK5O)+&Up|->7%p`WqgxG zk&wj+y~tMwvownH2WJLa@X?W@LOp2M6x(PUPo^5EG>oZ7P{7%bLXlu70pSgJ4?dNR zDRFq6#7gk>8R!vao_(DE-cAbaQP*E9L{Epw|G0|&RF7TSJkNxRokhM~s7$w7R~eAO z&K77(`pxb{SSAx&`CDuh@mNg1)sSTMsjof`ZMDq^kG%dODM~hU-=~L56DPj#OGf#+ ziW-{F4$|aiHu-H^S%^Q-0|>vd4ITVR+W42=x6yPVC178=aY)E}4UQwg*kUuA@e!6& zSvS#q1Jc@+gYl zjYU6en>UYrT-qC@v(je6N!;{J)FKi|HNchiVJJsOt!!e2U>@G#t3JPkb1Ah0Pe2id z^B0N4fBUfKDI05A ziGyc4Gc^{$z$(<&0z_@D`{Ln26~_7R*Oc+QED|v2jX(_qgKFD(2%@pYlU9PC=uHON zE-yVPK=4E&NMnCc-zj-{mih~VG-N%44V0K$Ag_=%TH)@7(1Iq!{pXj0K`+XLYmw$QZQrsLXB{g%`v}0B=u>JJN=eR(V+mABWo8 z(vmTRb8edosT`$_qf`qfVNrG*){_+hYElVFwh>A=rE~9*Y$LN(&PeVj-~5}EBuyDX zE}yMRJ+@mQcz-cRkP8fTj2|-xjTACJaXvgo}1(aE&s4gnv>x8o7k%i*2 z8()~OEjTI)V2`NdN{07ncvrYao^F$(VidaOL?fBrp2J~Mt1|BXIo{-CO+}97N+)k8w1BgQEU9}^Im?9PJkGX^sxOs z>HE5fe7TIakY*N$%UVR0a=x&qhp=-fM0g^izM?)@Z3J)UMjgPx)|KDGJ5WHiLqV(- z11kp3S1Rva0%OaTAexyCt`5%Lm|`+Oy)i#h1-&ZTUB0amE!$3y@?#`OXLbk2o2NzS zzr!1uNnr1V6F#fLafb8+K~cts{gMhF_J$sj_|vt-P_RueN_iWsHBT{zqfRa};_xr@ z0IFqSiS?DM;Ivod-o-gLZh?ZrW0PJ9U33b_$po7p|7)m^^Zo#Lo}i50##F>xgFi{V zrF8q9 z?EUmk>y9HBGR}1t;qNgJOwO|P8LurcN@3=Y9ks{8VLRa2DL~*?FtI8D=$qeBPnv)aPxWAy(FE+;HxC})GZ#DMASjcjcPpIL_C0wF(+cdZyR?4MOg03?1l?uYaVp+xN52# zS|B_|DZHq%_HVg9C2am~5Pxw(TtQWPCrrmmJtzFq_D)cBPKkPqvO6 z*}FV!G)+FK(VSXEr6cJKLzVBvAlBV>*LG!`Dri zD9@V7+8oPyLEPh_h1{{fe zCe*O87j8h8upbInL_J1?kqB5_P@kZ-&E%}d%y1^dE;>nxU4=FcBb}FyJligVA8eDg zQ3~WfPfC_j5Jp81H0$s@ku`D0dhqG0T(bf$tIe8cQ-+j@%NcdA&D+V?ZA)@+lqC>D9LLu7mRgw=GI?%?XM2TbpG#=dhZdlgrt9 za>i+v1tRJ4V8#WLkad7++&GJm3>J;Z7n?{*B_q}R6e`omr>I?<3I){BQKlipc7uG( zKY4V+{votc2Cbh72y~xs7w|?ZRtj10?5@KB976qnTQ`86dc*;R>u;2Upe)bpGn*3ZWkn_0E?5LV>RWF-3ThtH`ScI$;ofaEs5k2 zFe*~FH1ojPbMgzbbu&&26)8##2?E{g&031mt!RX_$bY!u)7HsDQTwH zD7&Y58MCE~B<}uZ#EeZcKO}PIuv7n@i)r@(G_iAyBOexfNG=0FuK>cg4 zf`6~bfxiDLu1~e{W=b=yNlv*pi*#Hy_8^V8`Kfxo30G@hOWNvY_U90rWa`#{i*W$- zBn{9q5=C$ydO%m{zcF4gs;$d`^t*1jp`|BTCBPS%b*bA1T?(}mW+{JdaS3#btHn}> zj#?nnM86{}OkfrPAB$QM_oC5-P0$^Q)&L62ix zk!q02m)e`0^2g{f-Y-7I$J*Wy6mgQ^wVN@4(t^~EUqus&$*vZw6Cc}izEuqGAOI<9~F^m>)sWNE0sVjLYQo?;o0Y+ zj_~M}=FX;2)h;i(P1z6qG!DLyywgVxQMCaZtAB(hr$gqXwVJJ5#Qf8?2@$oq1{jK3 zk7M@ta(5anS#y;}Us=Z9Ur^Npj>|wZlKck0u5oOHHqB}*M0Zla@~WCQ1Nn_$@fRPQ zR(2aR*{3t1*1LSqCGBHo{<2KbThw;|)R${k;L(q-bef<8@f+24tm{JAK1z{z-H;mp z5=jT^?sUDUmWnLK;B!-`pm@<~iihlHaBn6vleJ!i@g|{>-D)5Q?6R2#1Sb$e5-NAy ziC!f#lwIW2Uew<-iRWm3q~lTOp@Q$=7i3xoTu%0BFU#s!Ed2Rw<&Kz}oqgE$n)i<0@V3Sin1Gzz2yfU($!`y>b zyJHSgn`|1EQn*<#x2D+#JHYkA_bo;1NhD9wAIukS6S88rECu&T1^LpajmFUS+%OY+ zmTi2R85y)D>7a!>Actz$Z(A4by| zdQ8vQ^Q4g8q9uW0th(Vj$<(K)V8z{ zsQX(rM4G2il<9v<2)KDB=zmXWAM|EtFiLOtGq)4PW7li<$mMOSk^l6o7^@j)v$^%~ z{X7iGSYcZy+$Z{5W=qO)E40!zZtbnUqX=OB#$t@G%J%B`BT5ilYZqF{=N_}1*LW1c-$h$vyNl33Shg!7j`lyuAr(e^{`^u z@OuQUtkj#Vdhl;L{DS|x!x@{2Fq)Lt?wH8-S0A1=YpMjwy1;q(8<*@%x(%tZNBHWc z)p6H-`;;rMsIbUPwp9-#R$iiyM>}^ldHuMFjmyVl-m0ctlAUb0@1q)T6WNw5&3RN0}lJTu#QxS9%4PGOL zh9kJS&=Ds1@XE}UWf%%;NAvBI02f3)D=VhV8Ze_^Ym_xldRSM?F?b1wu8r`jQep6LCNKLXe4SzK`*do`|kosznGH z2tG?+s0pO=1`QRXR7{N0!%$CHVN5m94A4Qhi%eBESHB?(niho_Rim8ab|#(u2ZgC6 z2zJKHAL8ur-^RFLe(-b9dc!Kd)9H01FBG9Q@GgRi^nhZ{q+OMiEjq{hF|8y`lk7gq z4Qrcb?thyD7f+f2?Il|PbJoU;P@5%C{a$4@tBX+7-f1U(Z^gl>!H#aH!Ju zAC+XCTx&bT=fWO4!yft7XS0Z@=cFeRBSylV!CoE>>9We&a_GM_EaB1Vxw2Ipc@mS& zbidwh;S*X5A=m!<6SD|tPdjgpGJZk5p`Aw~K59n35d2=QB=}B4EF2rIh90H0 z52s@X3cK_z&C|1Z(_c4e?%r7n3#M9YMK7axrdfxk<(X4M2TR1aq(m$ls{)0_aV9~p zSiTLG*h?mC0u3Te%8rydbQ=bg8AiTu`-JE1yfD=59%`8RnASPTnKS_Hf7!R22q=nu{6zPKznyJ|mH_Q3t*QlFrkK9<8Oo;*)K+{*)16;E13zrGu^wseFvP8WaW zY^{nt$(h72hPm@-YXJn)xVXpT|3lbyXJx$2(qopNxVNS7M{O!NS$-1E&p+TaCSP{* zyJvfMZ}M~Oj9ES&pXx4bB059&YIN#G{~z2feG_ye<=!5aRo|V%?^M20y54{Ja=cOc z)}3EgC*%IU-ou*bHXztpTs4 z=xQl%Ud2lf>^VQNfo?kfx6}5k?RfRBlMBYh|Ha^%hm>oZ*A2_)CCIs5M1s5M4I4H$ z^P-C3ul^4U*Sx(Cmh|p*>@~2h7Wh)sJZvSG>M6%lav1)kLa$4JR~s|7odM~_-ybz~ z3DghDK~8=5T61>$4}2P}PGtu#e0R8YuJZJ5;?S)~J6sG~JRFU)VB^!y?&QtR(h2;0 zoMq2fi(fr3lCg3YH2*Jd#k;&JzIiFBYH6=;(bCnb%17hP|Mov5 zrFByIm>OUY)M@?4r*5xXnM}Gd>s~Q7x6glU0wZU+-&bdX@%T1;cOPMV_B{Z#YVvVP zKDWTDC zmylqEn@Ara!g7T0O4cj;K^lNoO(tJyM0mk=H{YK;y)TwoHOBXM*j<*&{nO&gXh}M0 z(TMF0Q30A-qbB4H!3CkKk1*k#5a;X}L9$b7`-J3J zSn(Eb3yVvoh@2cwUpLxls2SW(4qui9+q_<#UauCeucuS8K$m)!n)eBzXJ-og_bnN) z`b;qQ6tLH+lMv@mEc$42{m@J+Fu2zRQxhk`}JYX?C{O)8lO#BaiW{~aW zuK&ZKyTH`B<-19f``u2y1rI^+@;xCB|KBuLuj`$i>s3?r8&1Z#01E=?5=vE`e4TX) z7m;j_4SrAZ)J;{Uas6nKBs(0a?CurAR%8iF{I_TEJMzWEM4V#jp67P90l0;Yy=NQ!@d4IziPHAupEyQzIRv?ihGk>jwJ8;CQXGYl zo)vzdG`;UK&jsQ+1l#-JMcZ!8?~w z5LWeV>q2ncmS*QRSvb@#>Rh-^hvi0N+UZ1+R^QRjYuswDVv*K;(QrN;aBuD#IC+sO z)jYnw|G&B?CI3A6v&!r!bZ|H@Fcb(du%BS@wfOizWbC+XKLd&gKzhqpyC0%hL89%` z1_2-X3x#_;s|d5!a5bm_i2j+@K&4%Sx5+GaN1?{ONxHq$aVisd8<7JQ@vo4{Do`ig zsZBWa`-b?|0}NFIY;>H8Ho~~@)z0bzrk4CpJOxu0yk@QyT$29zASTNbh(sERy>p&t zV6G8ZvA-GB>`91$UP5U(@ZI-kXsrK%pQtz?`F7Z*@C!QV8%xTLUDo0AfCa`QW|gz_ zE)IT=Ie}vF_dYv->(c}?|I>IPJm@H5gUQ{+Bm>GR^oytV6TChg zCkPVx|0CjzC%|VR{2$BzUp)4oehdu!e;P2kcsqO&y#G_>|EI|R>-fL@n_~Z!TnvIr iBK{R$0)k5NzwT$e9sxT9mDK+Z8wdztsHFcV;{O6X)@#QA delta 2608 zcmV-03eWZUoB`QAu(bjb9E7bQb`Ssn*HHif69ALYOck@-0tO9#p;3v`en$u& z>r%5D7pl4i6G=*1J2%?*mp4)ChJbgTY7;w6@1zfUdgB15%qhUzGPx2Gkp3@GsW*-Ac#RI4}Jd$sm{qgJiAORdG*_W^_HN<4W47{pgGWJG)?3b-$n z&B_%DwNM!VKPg;R(uX)a0Cpm{*Zv-QEn(nQ_dN4@yD6Rmsn?1Cy33?iC3|?ob^nOe zJFSCKz0+<#+wq=%uezlN!Yh&iD8N6*A}ZoH``Dm8y)k3=dabpBk$Jt_V!@O3#%E=T zsmvV<)S~-=2!$n_O08P1Rjch+^+Bh7SgMtpokpkBsn?sO&TVb(J5$c$1zTJrT=P}y zoMq|CdGZcstPmNEPqocx4LPL7ccAyaUnY&}L8(!%mcclGTk8rZ%ODV}0LZtpOZJG0z8^0rB&CLdb(EoERxnPdOyCA1D7{Ua|bs0m3hyRM; z?mZ587XrfjfVp@DrH3IzeEIKx1h!s_NP-3Qm<6gx!G6M!maeFao>de{*PW8{DP*`R zuu#vV1Re^18Q(A2-hb#s8l2E1-E9&oq`)+}v+OD-= zX%9N3=51~6>u}!TQmNyn=jpzOI&7;pS8;Ew)DuwvY2CTG>y>TUAY3d^?WVg59_dFg z^Ql1nKujQY=%{x(rADo>K6F4>d6sh46CBB^68YoqB`x0TFH)?`zNGQ>{{T=+0|XQR z000OR033v^Pze+DQ@;cN0Kp#s01*I_5ls_+&6#0un>ZN9-`}UOBAp;o3P~CuSt(Pm zTcxU7b836jN!3(2;G|rQjg0BiZ+|wBjHEu?Ra`_~#g8%Y{P^SNZ{vSlCm1{c z5=e50eS1IdWB5XIUN#~c<_RL)x9d@#7N?huJPR0t(Ti{I#QZG9G3SU?15!p)%H&0Y zp*8e|v+-y)b*3Z79ge3H->&uh_TE_gXAb2x&KQF)fWM^ZvVXo=HhD3sUvL)`XLz*&yO1?wSyACYLjye&h)65>o8u@kGVb<3H6AA}e|1|Y^K zzP-~?9d@0IXTpNvrnhSQc6B=(GX5ndC>V{zoork3Hn!o-{|dk*iO?El`v~l9rcEhV z88zfwLTN}jCNV#!3kL5f#S6MV0pTMgD^$d8DHbv{jaAO*A~yQXSVnp*$7(Ep$L(Ti zK?#_`$TZF=Ltn-Dv6VosOMk3Sf2vM@+BsE(BJqb6A}FNE0-RQwm@od3M?WKTr+)y9 zMNAWzWe?Zqm+nvp>R1bEJXC^`;g0Dfg$yQ{^#s!t0s(5&occtY`td}Wx(pvyBfpmW zSmyru$ZbyF(WmdI(|0UAnlU~U}PFR=aj_|TWg%HjRCMM>omD2CQMZp2w{=#OidpiIQ&z4se2PqmmQ>Owp@xXRvt!WriY$w2WF zm1rNk^{4QgNcc^ZS>K>QF7Owqnx#0?Gn%Ox&AP`ZrCBV)S;?r7;JCOC^A#`igV5e) zz$Gc{jnUUxJH6#~)HT~_xn{S43r=R!he)hK-iRxtgu|F*W~D^|@JyUSP}5=4TVhAO zQg9BwQh19eUdPFET#|c#AGR;^zc0UjyL$6{eFB=f z=*_dEwopz7E!0L)|AnT`r5TIfSUYNC<#e#*y_U;nvf5O^?F@EOMw;R1ZMUPgT~5ok z>wnWsjn65l&r=%N=>52(Ub8y~U$ehQHat#HN@Kz+V{ktl0OxyushQuwo#+Kcw{%qu zWho_c4-)){`6ECud<5mp-ONyL>K(O{at`jKv}+orxM~!)2%C}U4Y;dzOK$geIPMEs z$pSx+vME{_#pHV_sZSSYW<@Iw6KN&~;FS=Ik>z?(SNHY4+*Q{%_u%zStET&noJx9| zR!L8dP9?y>+t7`F^&Z_-*Fd-B8rXvBwlcHU<}xPH9w;m9H!p$5dgtw_UqP>f3zoWa z-U~!;s6F*k#OvTvB#M{Ihb7L};x)?h#H`C@^Z_HX;*a};*qUVhJ~#q@6aWAK2mn)qqEJ!xW>XTA@JAUPMVll zqwz-ITSqLE0NU-}OY7Dxx`RjZzTXqzt2PJSCtu*a8cdDyjGzRJ7^E^?jb7Sqx<-i) z+(_OJ25NK!kCIJO3CqOb;K>aZTu|VX2pP|WtG^PfGHs*))2Vlf#V5G?INNRDaAtCPuI_Y6c~EYV1mbamg7rTt8lb%bK)FcwF5Zq zZ{!Xe#_?iOn~&DoKU|4IG56 zP=IM|!gde<0N0bfOf(-Hgso5s6ZKQS1ONcR9{>Om00000000010000008*0>O(Yyq z_GVKO0RR9H0ssIT00000000010000006-X%P)$byiy4!dO-TZa8k3n#8) and the PROM algorithm (). An informal description of the algorithm is provided +% below: +%% +% * Reactions associated with upregulated genes have their flux values maximized. +% * Reactions associated with downregulated genes have their flux values minimized. +% * If the |pfba| flag in the |hyperparams| structure is set to |true|, the +% sum of fluxes through the entire metabolic network is minimized. This provides +% a unique flux distribution. +%% +% A more formal description of the algorithm can be found in . +% +% *INPUTS* +% +% |model|: A structure of the genome-scale metabolic model in COBRA format. +% +% |hyperparams|: A structure containing parameters for constraining the flux +% solution. You can read for more details about individual parameters. +% +% |onreactions|: A cell array containing gene symbols or BiGG reaction names +% that are upregulated. +% +% |offreactions|: A cell array containing gene symbols or BiGG reaction names +% that are downregulated. +% +% *OUTPUTS* +% +% |ConstrainedModel:| A structure of the transcriptomics-constrained genome-scale +% metabolic model. +% +% |solution:| A Gurobi structure containing the linear programming solution. + + % Unpack params + epsilon = hyperparams.eps; + kappa = hyperparams.kap; + rho = hyperparams.rho; + mode = hyperparams.mode; + epsilon2 = hyperparams.eps2; + minfluxflag = hyperparams.pfba; + kappa2 = hyperparams.kap2; + + % Set upregulated / downregulated to empty cell array if empty + if (~exist('onreactions')) || (isempty(onreactions)) + onreactions = {}; + end + if (~exist('offreactions')) || (isempty(offreactions)) + offreactions = {}; + end + + % Find reactions associated with gene symbols + if (~exist('mode','var')) || (isempty(mode)) + mode = 1; + elseif mode == 1 + [~, ~, onreactions, ~] = deleteModelGenes(model, cellstr(onreactions)); + [~, ~, offreactions, ~] = deleteModelGenes(model, cellstr(offreactions)); + end + onreactions = unique(onreactions); + offreactions = unique(offreactions); + + % Set penalty for downregulated reaction + if (~exist('kappa','var')) || (isempty(kappa)) + kappa = ones(size(offreactions)); + end + if numel(kappa) == 1 + kappa = ones(size(offreactions)); + end + + % Set weight for upregulated reaction + if (~exist('rho','var')) || (isempty(rho)) + rho = ones(size(onreactions)); + end + if numel(rho) == 1 + rho = ones(size(onreactions)); + end + + % Set the minimum flux for upregulated reactions + if (~exist('epsilon','var')) || (isempty(epsilon)) + epsilon = ones(size(onreactions)) .* 1E-3; + end + if numel(epsilon) == 1 + epsilon = repmat(epsilon, size(onreactions)); + end + + % Set the right hand side for downregulated reactions + if (~exist('epsilon2','var')) || (isempty(epsilon2)) + epsilon2 = zeros(size(offreactions)); + end + + % Set Parsimonious Flux Balance Analysis + if isempty(minfluxflag) + minfluxflag = true; + end + if minfluxflag == true + kappa = [kappa(:); ones(size(setdiff(model.rxns, offreactions))) * 1E-4]; + epsilon2 = [epsilon2(:); zeros(size(setdiff(model.rxns, offreactions)))]; + offreactions = [offreactions(:); setdiff(model.rxns, offreactions)]; + end + + % Setup Gurobi model + tmp = model; + tmp.A = tmp.S; + tmp.obj = tmp.c; + tmp.rhs = tmp.b; + if (exist('model.csense','var')) && (~isempty(model.csense)) + tmp.sense = model.csense; + tmp.sense(ismember(model.sense,'E')) = '='; + tmp.sense(ismember(model.sense,'L')) = '<'; + tmp.sense(ismember(model.sense,'G')) = '>'; + else + tmp.sense = repmat('=', [size(model.S, 1), 1]); + end + + tmp.lb = model.lb; + tmp.ub = model.ub; + tmp.vtype = repmat('C', size(model.S, 2), 1); + tmp.modelsense = 'max'; + M = 10000; + + % RESET OBJECTIVE TO BIOMASS REACTION + %tmp.c = zeros(size(tmp.c)); + %biomassPos = contains(tmp.rxns, 'biomass'); + %tmp.c = biomassPos; + + % Maximize the flux through reactions that have upregulated genes + function tmp1 = max_flux_through_oneactions(onreactions, tmp, epsilon, rho, M) + tmp1 = tmp; + for i = 1:length(onreactions) + onpos = find(ismember(tmp1.rxns, onreactions(i))); + + % Set reaction in forward reaction to be active + newMet = size(tmp1.A, 1) + 1; newRxn = size(tmp1.A, 2) + 1; + + tmp1.A(newMet, onpos) = 1; + tmp1.A(newMet, newRxn) = -(1 * epsilon(i) + M); + tmp1.rhs(newMet) = -M; + tmp1.sense(newMet) = '>'; + tmp1.vtype(newRxn) = 'B'; + tmp1.obj(newRxn) = 1 * rho(i); + tmp1.lb(newRxn) = 0; + tmp1.ub(newRxn) = 1; + + % Set reaction in backward reaction to be active + newMet = size(tmp1.A, 1) + 1; newRxn = size(tmp1.A, 2) + 1; + + tmp1.A(newMet, onpos) = 1; + tmp1.A(newMet, newRxn) = (1*epsilon(i) + M); + tmp1.rhs(newMet) = M; + tmp1.sense(newMet) = '<'; + tmp1.vtype(newRxn) = 'B'; + tmp1.obj(newRxn) = 1 * rho(i); + tmp1.lb(newRxn) = 0; + tmp1.ub(newRxn) = 1; + end + end + + function tmp2 = min_flux_through_offeactions(offreactions, tmp, epsilon2, kappa) + + tmp2 = tmp; + % Minimize flux through reactions that have downregulated genes + for j = 1:length(offreactions) + offpos = find(ismember(tmp.rxns, offreactions(j))); + + % xi + si >= -eps2 + % si >= 0 + % rho(ri + si) + newMet = size(tmp2.A, 1) + 1; + newRxn = size(tmp2.A, 2) + 1; + + tmp2.A(newMet, offpos) = 1; + tmp2.A(newMet, newRxn) = 1; + tmp2.rhs(newMet) = -1 * epsilon2(j); + tmp2.sense(newMet) = '>'; + tmp2.vtype(newRxn) = 'C'; + tmp2.lb(newRxn) = 0; + tmp2.ub(newRxn) = 1000; + tmp2.obj(newRxn) = -1 * kappa(j); + + % xi - ri <= eps2 + % ri >= 0 + newMet = size(tmp2.A, 1) + 1; + newRxn = size(tmp2.A, 2) + 1; + + tmp2.A(newMet, offpos) = 1; + tmp2.A(newMet, newRxn) = -1; + tmp2.rhs(newMet) = epsilon2(j); + tmp2.sense(newMet) = '<'; + tmp2.vtype(newRxn) = 'C'; + tmp2.lb(newRxn) = 0; + tmp2.ub(newRxn) = 1000; + tmp2.obj(newRxn) = -1 * kappa(j); + end + end + + % Set params for Gurobi + params.outputflag = 0; + params.Threads = 2; + params.Seed = 314; + params.NumericFocus = 3; + + up = max_flux_through_oneactions(onreactions, tmp, epsilon, rho, M); + up_and_down = min_flux_through_offeactions(offreactions, up, epsilon2, kappa); + + % Solve using the Gurobi solver + ConstrainedModel = up_and_down; + solution = gurobi(ConstrainedModel, params); + grate = solution.x(logical(ConstrainedModel.c)); + + % Check and see if a biomass constraint was added onto the model. If + % too constrained, square it until a feasible flux solution is + % obtained. + %if model.ub(logical(ConstrainedModel.c)) < 1000 + % while grate == 0 + % tmp.ub(logical(ConstrainedModel.c)) = tmp.ub(logical(ConstrainedModel.c))^2; + % up = max_flux_through_oneactions(onreactions, tmp, epsilon, rho, M); + % up_and_down = min_flux_through_offeactions(offreactions, up, epsilon2, kappa); + % ConstrainedModel = up_and_down; + % solution = gurobi(ConstrainedModel, params); + % grate = solution.x(logical(ConstrainedModel.c)); + % end + %end + + k = 1; + while grate == 0 + epsilon2 = epsilon2 .* 0.1; + kappa = kappa .* 0.1; + + up = max_flux_through_oneactions(onreactions, tmp, epsilon, rho, M); + up_and_down = min_flux_through_offeactions(offreactions, up, epsilon2, kappa); + ConstrainedModel = up_and_down; + solution = gurobi(ConstrainedModel, params); + grate = solution.x(logical(ConstrainedModel.c)); + + if k > 100 + break + else + k = k + 1; + end + end +end \ No newline at end of file diff --git a/MATLAB/expression-constrainers/CFR.mlx b/MATLAB/expression-constrainers/CFR.mlx index 6c44e0d48251d1702cf2c9701c048bf6bf616b3e..b3d402db2664fe2d74035632b78cac32a8711a3d 100644 GIT binary patch delta 4732 zcmY+IWmMD+)5ezu0YREwLY8hMgrz&ByStQbSK&`AAP5pmE)9~>Ard#;B^}b;f;0m1 zdOznqZ(QHznlm5fcjm*)%m~klRJWt!^0zl>-y1E(1B&Ymq^ZeCs9*V)5Qh& zEU_%vxjE9!ekcje|I{nOoAg1NoGjb`^A&4mE^d0Y4a#&hL#`{s642m-;ZRni;q;AU zG6wP26MrWNc$sU3VnzEGQvFi{%RGK@F z{A31usyP>D2vDDLld%)<5n>9&wRIV>%JDhpH*1Xu$8OPFX)*K)01^m~WT8>U3+ciTB9&y3$9ktKarX#)hY^od*Fzx>{rmF>;etfwSd{4wys89s*nrzL#j zRt29gK+nO9Vy(_}%~@LoW5VQj43%KjGce8PsG%t-OnE&mLJa=Tvux=@Wwul#*-#M{ zLoNp?r?n*{b91JCzB%>;wigo}hZ)PP!ncYY3V6x8P z;+=qv=kxv+-H*tW|8_Z?q*of!(ptYtoc5CAQoXEX^rcYIhGLs|Cj-&m8A%q}VYapixk2u?sYqPbW#0MpNH9IztgFU z5IqF~-G;ts+$8#SooM{rVb#w%P3%1`$dLXBSI+ynSrI5EL>>|wqz8{>-?$j;kA zAK#cRcu1tOG%?rLVy9k9=}U^XsLvD)ft`eRqa$0`2M5*F;)s7?@-r5C|8<-wrfs^+H|w zAD)0fQ+yy0>0jH`+Q-G(hS$#B*4NeE&4(wz)n&%?tu6#c@hFuetES)QH+!=epN?l$ zM1XJXC!5=`FdAXQ?%U|X>C-5h2$u7EnG=}vS#mlu3U82XABAQLhr|0%Xmyast8QH& zb5b6*kBV=)SF8vV25cpPScig!Mr7Eff`!ZP5j9y770io-d9#?=9i6_6hqlul zrk?zM8Ob-q(gkV6gtR|>U`7vJ1eTGNbiv}-l%BU8Ihoh@Rh*R`MF`bGf zZYe^ak;ZA1M#@5KhS6^X2@@euGR;z$clRUy23lQmVa5m^$FJ`}R2*~VD(!<=5)^n^ z;fg&SQgNcA$QZ`vVNdUULmuCD?!CIOCq?f@@Ox~y(`vI9#e8^&V1&xVSWv>k z)^iW0%%T&6M9vy`Sd06Wh_tP1j4n2CA@^@u6>6nu?k9m&s;kKM`{rNyl+{iS!!5#* z|G+4=o?v*L@WP2=J30GQ(axee+r0D-vSlHsw1`ciOEg^`ym_1R>|j7WTpaC%BTc_| z*gYE_A@h5Y3xO`{Cl&Ag!a*=*BRK*P3+6dn%JUs`OeK{O!W(NVMJRNf#*(ytjO?=v*cHR|X0i$?{vY}juQ zHiLfgex#Z#nI!?wh~u1(ftoVFT)_SN_l& zt_~OIb3mJDuu;w9zj;^}Yv#Ny-n0Az4Y*z}e7v5DyO+NIb9dMIOU-kGzV*aXx3-$J zS}(<$xw;wICrOr{B4t?1PLat_a}GP zxiXU29#`BO^6{H|^WO*LL*Mdn$@CW>a-;wbOU>^TrGno+zjya64qX1+yS%;kJNqzM z=$j^{!|q0;1G3n{$cXP6ut7RpPaG+mw*+R;y#j$w|*?#3n5#Y2|~j-^IuCrxDmKcP2OF&E!C zSaPv&bH@&y-WGnUHAb+tOLHX?w(4&?E!xZIM&j`2a6B!C+gG zbr7f+mBUX_rywNdG$Je}t7)9gCnJfDC2ft?-y$}HP+x1+2u+H;nl-Gg6*MCj%}}ez zrG=}_gk?^a)4Vb&EA7!sEc)~TWf3iB{N4nqmTkquR!@<{+_e;XvH?<}ld;9^Fh#>jdDSKUgPP7} zg=VuVGnGHNn$iR*wjk`8D0ax?BvJ-b5?yB4LCp+!=3TlL@an}MhNiJvZ>CKGM-P2Q zuOwz$Rp*ZCqiYvS07dZGb1llOH!f6?j^=_#9gmdiGM}k1-THWh71t=}ln?gT7+r%t+RBi;U!jeUmSQ3Q{$ z=2e?%*8R2_ie|NiuyZTQOz!#(|(qkVO~EZ>nM zuMh(3#nYc)b201$Q~zhEW6@}jeOUS}1w#{01Gl zD7}7ru~Zi3Wmjl;k15NvN$4JklDNfMVa)TxMK9SGOBJL3Jhv>a8fnY? z0$CYNncR!i0aut`+7XR#seC4+HM;9gCD|v$Lf}RFJ&!8MNFAyV8o9{jFBBggJs4Z| z+&1zA^b5KA)9PNw!2DvIKl*^>Mx1eL9GcPcTB zqU)IS&ge?yO^FutWV|s;UmNxhie_GnU*U3x;kgG7ZuWQmF0;CMLL+FmK>q|MHvQ#4 zwUen|L_O1ynb~SpX2#G>`JL;pO zaGVh}go1s;jK|V@4@hgWS94DSQ0jI25dCUW%T`Cz!N5M`U{B#H8~QFTFQTcav`v=H zxHoa+1L_}(60nx;cfoVNheAY+bJ=l>T*YgnD&s9)b+;dA8QMfRLVmq@>@%M38m*iU zK)+`wv9Cz~0<0w93-<&%J|gha(B7lkW+z|0PA+mPDZ!4wxT`+bbs7uqG?@4=|wB25uo;8wpqc}F6tdL11n-%KJv8SH6){esQiF1v)Dl2VND>rl7`#L$-gIM1O-l!ox9GbwRvbM#(Ik9cPRaWm7n9|ugsSg}i^w0Wa%II% z*^Flk>&>1%&0ezV#KMw{f^Ig|lm9?i!5A#ur15hUyLvn%M7U+S%$3T4c=~*dL~Rpd zLYaCfZY`ya7-T|bVaV87%3KQ1L?su8kR_l?M^U5Pz}j)$rvHU_0fxviZ62{GZ?j=x8T5E=-?213d}AjoP65i**J5O78lVrJjoiN${=k%(PJob5cZ*Pc4W zXwf}aW{I!tUUFYUaF?J9J#PnKcp#LEHX^2^FGtkorK>KizGVj znbMt&18X+@k&|pUNwarBs_3aLwMvH2`!nxuU-MzHUe;C_@xqrQ9>tx3uZDtL=i;ii zNEut)8M3RGt6{=xrrbJomLNPkrhxET%?QiMyVGzyT<@Lk!BDyb6a!`|@c9h7 z`Gc+me~4WCT~Nuz3te07&L_V2xCF4&QBwLZ&58eK-?+h8rb#M z?OKq6F4N22%Nk8Bf2!Jd+<5YpMz$=2wt3Y{1OykU1=RIO(8$>hlW%#6*0h^qNlZ(2 z`RnmOx~v>f|k3~d*UL567ME4TTUf5`2i62JeW@?VT?t}5tbt>2cnY0nN zKgV|s+Zel&(jZjHDIv7?P|Xrbpa;o-hjR?i$1B65q;nhuz;kcy`ykTxSP!f6$9Dxl zIw@Wee!`1eb+s+CR;cJmIARg1y zGeg`ChAGS@cv1YJ4p`49NU8ZL_t;hpPBD65d>ST>KLba4e^2Y4AqvP$X<)b_4I>J1f7o+c!HLN zSr2E}pASPoHta?zLR7yRG6YJxpIr1rHBzA_aHjKn4DL%UQ+URIGf&x4?b*k1X5vP| zpkf>PTpv^G^BY{vD4c_g6?fl zW!Qp}QwtXn?^BxdO9(?3h_=_(dv^k}CDV9Us;Yci+Cs%;N>4<}r>#t$X7e3Ay?~2Y z^<=*pD=9L{cGxpOD+juBNgj&RH3>!`l1T3Orz@T)NRyWNJGI!FJ2E3}>dRn~UGY85 zIMI(w<=6xTYgqavK9u3;8yTWLu$2x%xH43BD+c2&zYB1@t^y@du zuUe?m0DSW;ZvDD%G);|_lyx&?ws9MZQGK;QEzsga&$liYJ}$v|nWF*bZKaPea-R39 zrAIQI3@*!^F3zZ^PgU<-t_?B#56+s%R?>c`FWlm4t6^bN0RR9z07WaTLi=Q(6Bpt^p3pfG@C{5IQNo=4;D2|SLw}m$jTfOnZx}m?^m0BP<%G*li&qo4U9Nm2%L9=NdmQDxK~gOb%-x%|U7)WG~Wxsdb8&fV48ObcsX^Pi747BqlL_ziW7#02Y^mOvpuD%n_`IHri9rFfJ=RW~UXAq81*@>R0(7Pn#2 zCs_{|Q7_b7H>ZQkg*o6+c8e3W)9xc;CZy{i$R=l8`sd2P`8w9iE<)L>l$N-%JhJ?s zpCXnq_6IbnJlOUWWu@hDuhGjPnmVw~>=U)Vjtaxo2+Oc~-&#EpSr#g%VN-osHirta zDG4mfdRz0=&-uTD?rA7vSZD!VrHpDs7oB#nxD2+O;K~fQ1gux~YACwY6a)VyCHUON zS@IceG@n$ye}AUDy&EplVkfZg2kwKL;rHr8^JbXzLku(ABG6AnDu8uj1t}J4s%pTA zp(O17^tITK@h?HP50p!PqO6=2$QU}IZJA~t&X-+q>~mWB$_GzL(m@=F>x`)}*EVV} z@dc5*ikTrabLK&gmZ(d9HS>}{AC;n(ZGGf^`btQdB6d8uK+XMy>&DE5XBt(ZWj4R^ znHEjC`xdD zn#p2*Lh~s`C+5=*hn$WNa2;oEfD#An{p*E3IpyepQRKS-=)>y|0%)gu-!D?_nbK^n z5E2ir&nMoy3KoNzg}>3suQLe`U(XI=HYaYTn+4K)*iHL$CXoXMSW`#aiJdG1EJ*BE zu4TCxKdSE!*j{QG;=x&XOBD+R93*j(Ozoo$^72Ndxb@T#c z!Vxl|de$vOWs63rA+jJ`cl4!bdg?H z+yri2gqd5hGq_}RKmP|ziIEqU=rEYNyX!eTN(WSUI6nk~gH83F3)=FR!$R^_=&DfJ z*lT#LJilThl@TeAUkpF_@_yfRY?O1DzLyFZfe+%K@eK(AN#h~85p6e<3P%>?JIpp3 zT{K;x(CV^{LaC7=YPM6%Mm?gV<_c7F{{V>4iVvqcRXl)+H#N6kjR&T@q25moi2!#o z0W@6#qm0yP^-ppIEpqTgZEfUfudx{bkwAvS120M!%q_299FBTlpf8ht{*0RJ7uF$3 zFPz_Ia#L;h5VUIMgqK=sP?J-47DVJ5l1^GTgWtm?1bbha?ssKmmeahH4?1V!L`0J@=xZ{W{%nIDHkj z!Y{32@$Hs#6KpP~V7)t+sa0oT(x2boVJ7^cdM`?u_Xw?LMnm*%LWIn)`1R|Xry()M z5?M=blbYOIBHf=kJx^^N*+kp-e|!DS`JS;j@h%30QctHw8EYAp6f7yk6jexZD1@zhH@yYbAa zUe2Tj>2=vqZ5A_h91wEjZN}3)OOH{?)g1}wV$^CC{S^HB^TsNJx+$0F=X7p@Ivi1| z=lW*e!fz>1eklaAokwLe*oYn+YvcL~lm2s!-ZCERcqxa6Cf=!~X12%Ox0Mh| z;iarxPRI!e0u3B;c8+S%3NqFBb|te?=2i~de>6;~c#bwj)j%D}Th`&?pOAvFb%VQU z8%1MIJQmD^(^B!Z1j4k~zMY;Llcy+exm_shVNFPTfAFn= zcTr^HSaZ>9&NOx8h?L(N#m?J;Gf@P?rfSh$L@^n4uduV@%|$VTEJwTch~%u#3kg^e zyorz`tI4iiC!u}X#7FNn-`)sm!|rcRfj_gK{lPDxLhxbVzQlli8HL^sIp9EjCHk!~ z=PgB~XoPDtmM1gCh{Mn=gKay3-X}pBp$Cc2q?cVNaQbxqLF@SON*CEgnih{6(h zA`A2`TPongP1?rY^MRuq|LRTP8OxuP{m?s$@xwp=!(=NEU#R+l>#zJVX%RuxTrd<$ zvzLVmXU*_UbFpl>uFA5RgLfeuE($utlivFnP9fXY-d<6)DxBi0a{urFL$h{em;V}P zOW#$sPbZ2(u+BB{joO(>^%0-y2eKk^GuPbbl1JRe2H_2_E=aP?T)TllC_Mi7)e?oX zd`%{qfLLz0*6ZFrw^d0#d8yW9BN(fUA|V**DF?wuePkJ8>hq>S1v{;)FA2=Y$g<;m z=5}Li8I26pmg_GVH;W5}bL|5(Z$E7arDd2_Z_Wu4?Zn5m`@&v)&tb9{D7Ln}u;9SZ z?+@Ob3dWlRldpBsW_)hCyU#Pv<3*xLL?T2?>Vfu#nNx8gnd8iIj6U@SbPpft`_Dl2 zUyGi_IRtvJY0B!y%cxZZ0_rSVvbGJAQT10;&Tb(B`ejI3_iy4sgb za$fF`Q%JX5?5wSJOhoj{7aTSV`7?N7>fUhI~LOhI&8xig!i zImhPsidmx-KgX{8Ea_?F+ClY=&a-^wV}d8xy4e&1LM0_N4rle0fom^YFW(5F#v7V% z54T7W{qzz4RZ$SSAW4#cB$fYvV-T4j5$1m&hft}`q`yPn|8?@e{*Ph&8wC&tATYzf z$MP2-!vp{x{fqx!{*(&A|40c_fd;Wm_Xwf%hze2k2!x@4=zjE!;lHc<+sXYuUC-(M Qmh&QH8OU&+)BltC4;G286#xJL diff --git a/MATLAB/expression-constrainers/CFR2.m b/MATLAB/expression-constrainers/CFR2.m new file mode 100644 index 0000000..2ce2525 --- /dev/null +++ b/MATLAB/expression-constrainers/CFR2.m @@ -0,0 +1,136 @@ +function [ConstrainedModel, solution] = CFR2(model, hyperparams, reactions) +%% CFR2 Constrain metabolic fluxes using the modified iMAT algorithm and PROM lowerbound constraint +% + + % Unpack params + epsilon = hyperparams.epsilon; + kappa = hyperparams.kappa; + rho = hyperparams.rho; + mode = hyperparams.mode; + epsilon2 = hyperparams.epsilon2; + minfluxflag = hyperparams.pfba; + + + % Set the right hand side for downregulated reactions + if (~exist('epsilon2','var')) || (isempty(epsilon2)) + epsilon2 = zeros(size(reactions)); + end + + % Set Parsimonious Flux Balance Analysis + if isempty(minfluxflag) + minfluxflag = true; + end + if numel(minfluxflag) == 1 + kappa = [kappa(:); ... + ones(size(setdiff(model.rxns, reactions))) * 1E-6]; + epsilon2 = [epsilon2(:); ... + zeros(size(setdiff(model.rxns, reactions)))]; + offreactions = [reactions(:); ... + setdiff(model.rxns, reactions)]; + end + + % Setup Gurobi model + tmp = model; + tmp.A = tmp.S; + tmp.obj = tmp.c; + tmp.rhs = tmp.b; + if (exist('tmp.csense','var')) && (~isempty(tmp.csense)) + tmp.sense = tmp.csense; + tmp.sense(ismember(tmp.sense,'E')) = '='; + tmp.sense(ismember(tmp.sense,'L')) = '<'; + tmp.sense(ismember(tmp.sense,'G')) = '>'; + else + tmp.sense = repmat('=', [size(tmp.S, 1), 1]); + end + + tmp.lb = tmp.lb; + tmp.ub = tmp.ub; + tmp.vtype = repmat('C', size(tmp.S, 2), 1); + tmp.modelsense = 'max'; + M = 10000; + + % RESET OBJECTIVE TO BIOMASS REACTION + %tmp.c = zeros(size(tmp.c)); + %biomassPos = contains(tmp.rxns, 'biomass'); + %tmp.c = biomassPos; + + % Maximize the flux through reactions that have upregulated genes + for i = 1:length(reactions) + onpos = find(ismember(tmp.rxns, reactions(i))); + + % Set reaction in forward reaction to be active + newMet = size(tmp.A, 1) + 1; newRxn = size(tmp.A, 2) + 1; + + tmp.A(newMet, onpos) = 1; + tmp.A(newMet, newRxn) = -(1 * epsilon(i) + M); + tmp.rhs(newMet) = -M; + tmp.sense(newMet) = '>'; + tmp.vtype(newRxn) = 'B'; + tmp.obj(newRxn) = 1 * rho(i); + tmp.lb(newRxn) = min(tmp.lb); + tmp.ub(newRxn) = max(tmp.ub); + + % Set reaction in backward reaction to be active + newMet = size(tmp.A, 1) + 1; newRxn = size(tmp.A, 2) + 1; + + tmp.A(newMet, onpos) = 1; + tmp.A(newMet, newRxn) = (1*epsilon(i) + M); + tmp.rhs(newMet) = M; + tmp.sense(newMet) = '<'; + tmp.vtype(newRxn) = 'B'; + tmp.obj(newRxn) = 1 * rho(i); + tmp.lb(newRxn) = min(tmp.lb); + tmp.ub(newRxn) = max(tmp.ub); + end + + % Minimize flux through reactions that have downregulated genes + for j = 1:length(reactions) + offpos = find(ismember(tmp.rxns, ... + reactions(j))); + + % xi + si >= -eps2 + % si >= 0 + % rho(ri + si) + newMet = size(tmp.A, 1) + 1; + newRxn = size(tmp.A, 2) + 1; + + tmp.A(newMet, offpos) = 1; + tmp.A(newMet, newRxn) = 1; + tmp.rhs(newMet) = -1 * epsilon2(j); + tmp.sense(newMet) = '>'; + tmp.vtype(newRxn) = 'C'; + tmp.lb(newRxn) = min(tmp.lb); + tmp.ub(newRxn) = max(tmp.ub); + tmp.obj(newRxn) = -1 * kappa(j); + + % xi - ri <= eps2 + % ri >= 0 + newMet = size(tmp.A, 1) + 1; + newRxn = size(tmp.A, 2) + 1; + + tmp.A(newMet, offpos) = 1; + tmp.A(newMet, newRxn) = -1; + tmp.rhs(newMet) = epsilon2(j); + tmp.sense(newMet) = '<'; + tmp.vtype(newRxn) = 'C'; + tmp.lb(newRxn) = min(tmp.lb); + tmp.ub(newRxn) = max(tmp.ub); + tmp.obj(newRxn) = -1 * kappa(j); + end + + % Set params for Gurobi + params.outputflag = 0; + params.Threads = 2; + params.Seed = 314; + params.NumericFocus = 3; + + % Solve using the Gurobi solver + ConstrainedModel = tmp; + solution = gurobi(ConstrainedModel, params); + %ConstrainedModel.vbasis = solution.vbasis; + %ConstrainedModel.cbasis = solution.cbasis; + %ConstrainedModel.slack = solution.slack; + %ConstrainedModel.pi = solution.pi; + %ConstrainedModel.rc = solution.rc; + %solution = gurobi(ConstrainedModel, params); +end \ No newline at end of file diff --git a/MATLAB/expression-constrainers/addMediumConstraints.m b/MATLAB/expression-constrainers/addMediumConstraints.m new file mode 100644 index 0000000..d26f204 --- /dev/null +++ b/MATLAB/expression-constrainers/addMediumConstraints.m @@ -0,0 +1,105 @@ +function medium_model = addMediumConstraints(model, queried_medium, path_to_file) +%% ADDMEDIUMCONSTRAINTS Add medium constraints to metabolic model +% |addMediumConstraints| sets the lower bound of the metabolic model to reflect +% different medium conditions. +% +% *Note:* The medium condition file is assumed to be an Excel sheet, where each +% sheet name is the name of the medium component. The following details should +% be included in your medium map: +%% +% # The medium component common name +% # The default lowerbound you intend to set for the metabolic model (Default +% lowerbound = -1.0) +% # The scaling factor (alpha), which I have defined as the proportion of the +% new medium component with respect to the default medium component that the metabolic +% model's lower bounds are set to. (Default base medium for RECON1 = RPMI) +% # The adjusted lowerbound (default lb x alpha) +% # The BiGG Reaction ID +% # The BiGG Metabolite ID +%% +% *INPUTS* +%% +% * |model|: A MATLAB structure that contains the lower bounds for the constraint-based +% model. +% * |path_to_file|: A string denoting the path to the medium sheet used to set +% the lower bounds. +% * |queried_medium|: A string denoting the medium condition being simulated. +%% +% *OUTPUT* +%% +% * |medium_model|: A MATLAB structure that contains the constrained lower bounds. + + if nargin < 3 + path_to_file = 'FINAL_MEDIUM_MAP.xlsx'; + end + + medium_map = path_to_file; + medium_model = model; + + % Support for MATLAB versions less than 2018 + if verLessThan('matlab', '9.6.0.1072779') + [~, sheetNames] = xlsfinfo(medium_map); + for sheets = 1:length(sheetNames) + + % If the sheet name matches the queried medium, set the lower + % bounds using the BiGG Reaction IDs + if ismember(string(sheetNames(sheets)), queried_medium) + [adjustedLB, rxn_ids] = xlsread(medium_map, string(sheetNames{sheets})); + rxn_ids(1, :) = []; + rxn_ids(:, 1) = []; + adj_lb = adjustedLB(:, 2); + + for rxn=1:length(rxn_ids) + pos = find(ismember(string(medium_model.rxns), string(rxn_ids(rxn, 5)))); + medium_model.lb(pos) = adj_lb(rxn); + end + + % If the sheet name cannot be found, set the lower bounds to + % RPMI conditions. + elseif ismember({'nan'}, queried_medium) + [adjustedLB, rxn_ids] = xlsread(medium_map, 'RPMI'); + rxn_ids(1, :) = []; + rxn_ids(:, 1) = []; + adj_lb = adjustedLB(:, 2); + + for rxn=1:length(rxn_ids) + pos = find(ismember(string(medium_model.rxns), string(rxn_ids(rxn, 5)))); + medium_model.lb(pos) = adj_lb(rxn); + end + end + end + + % Support for newer versions of MATLAB + else + [~, sheetNames] = xlsfinfo(medium_map); + for sheets = 1:length(sheetNames) + + % If the sheet name matches the queried medium, set the lower + % bounds using the BiGG Reaction IDs + if ismember(string(sheetNames(sheets)), string(queried_medium)) + dataArray = readcell(medium_map, 'Sheet', string(sheetNames(sheets))); + dataArray(1,:) = []; + dataArray(:,1) = []; + adj_lb = cell2mat(dataArray(:, 2)); + + for rxn = 1:size(dataArray, 1) + pos = find(ismember(string(medium_model.rxns), string(dataArray(rxn, 4)))); + medium_model.lb(pos) = adj_lb(rxn); + end + + % If the sheet name cannot be found, set the lower bounds to + % RPMI conditions. + elseif ismember({'nan'}, queried_medium) + dataArray = readcell(medium_map, 'Sheet', 'RPMI'); + dataArray(1,:) = []; + dataArray(:,1) = []; + adj_lb = cell2mat(ataArray(:, 2)); + + for rxn = 1:size(dataArray, 1) + pos = find(ismember(string(medium_model.rxns), string(dataArray(rxn, 4)))); + medium_model.lb(pos) = adj_lb(rxn); + end + end + end + end +end \ No newline at end of file diff --git a/MATLAB/expression-constrainers/knockOut.mlx b/MATLAB/expression-constrainers/knockOut.mlx index 7d59d340d146d5a2f2ee9c191d9fc911511d1068..97989bb48a8383dd8702cdf8aa108ce904ffd79f 100644 GIT binary patch delta 2633 zcmZ8j2T&8r77e{gOF)DOC`EdcNDW;;6H4e!MUWOjKoJNL2oQ=AijYJB>AhDG6d_Uq zNJk!Be01qWqzim>=KmT0?(ELoJ?HGq-Pv>JE(sBZDoK>gE>CL^d?erQcq#xO58p#6 zfk}I-YX|pBpjon{E{ke8gsf3Xr8x<36eSSsPfy|q^PN>gELx$`$*F+EbivKIjyWG{ zZp_$4%8^25t0z3EtQ{1Y2IgJP>dc-O`CI*)qij-(bHx$1E=AmK9Gw?TclITd$Tq^D zsD8)pvQOh>=AYL;R&N)c#`VUbJqO=Mg$vN(TfRnIy&PKeq)Sq& z8i>p*l{3%9n??C^c}044Nu?go*cr}zc8lbHsc|W~!2u%}Mz=lTN+X<(@$k92kk1iw zEe@`ztUuANZxj>1%WtU9TL;35Xn|UC9&hT3ocx@2jF_et6hYeY213bMCn(sqi!C}R zn4&~W@#%Yr^FJd3g%^14d*3+QT@^marnF>_<#KjrmlQTGA8s!Xq35_dR1BD>cK5Dj zAmod8JhV_p6<50)q=kbqyBx>-`D5E#oF`4KHU6>a>UCMuh7)V?1ci0=pu;(r!HTGY z2EHSyj{b`7CRw$iY-s2+5*)ZHUPa75N~t9)KSd|kCTWMhVsG8bl^$zNzl-);9$(mA zGv&|zrLjcDGYMaJ3R{@=*!K9|G-XMeYR2ZFel{Ew^hlgYpGILeX1^=qg289GTsQ6A zFu_@674pG8`nUV*C+@YRJQZhV&`o$eZaLbgm2x$C#0@5FCS!SjNvu-cVlpFlH+%?W zWULe~av9@w^!0GjUak)t?h_bzM*k%giC}$aq<2Y?ReL1Uqu2S!naYBWFQ&DpO004TxZ;NVfVKQ3JEjq-4_on5+SWS9PDzOFi zTlFxC3Q(1=d(@Fz;_lQ>U84_c(_Mkm1iTcpwS*_wCKj&;pRg3+e0SO()s3;qIw!cg zpceIxf2}LHY!P9c3>(U>*aRv9i2Ww?XX-JT+^vyTH?-5M#3mFgHX4E|ZwL3=xTnQL z>FAC?_J%cu2y@?Honx5M1M$`wz0NBMVTe4$-JVENyq;HN760Ro!aLFTUDIErsl$#B za{XgkTz(C^rsoLu7Rp1w8MmuCmtb#d9mJeO)EXDV0z(6eo7a%Im)si+@=r_R-&I== zeK>XOGT{nykFb6IBh}hF3z#u*$1HAeNM^AR^+j42%~(PJ#Zq8G%}r?3KW2~`-e^fu zc4;zQ5il|})S5HV7E?eof9cpKa>xzPzZMQBn6UwccI7Qba&O5&NCRrh&Cb}3;q{JV zvP>O1b+wGE;Y2;MJt_3dbzv3aQnv;B^Vc?U?mU_XwImKn0UK}H{_wfoz$*BXiI9|T zsgbpu1~m}@EW`LMP6=5h$Y+{NoxN1tIVSJizHDYPx&t>uMyaW`-h)9#c(sUQON_>r zUFI|C1asjqIzPJ*Wc;m9BOW#d{o(Jv{t&U@Cm&x~_ro0-t)T7>_T39QooqcOx!}fi z#u5D0qI=C>(A2EamTMD?8#)T>jTTd;tsN4Y z3WXrr)%WjA)V6SOl9P3MZUy}Ghj~T!xeA#(h)nv_yh=oJUy?5#eM4YGT=hfHvttxo z2vO*DjGx%l8;DNF?ID`>>{_`Z!G(=n29Mj%LY`de0SkKdzJI<erm^~ zKdxw`kkQeB?1&F@*cN$Ge}+{#q5I3__(4W8H8M4%9G}7@hPi8umF0*$V_a5yOCvW1 zXc$<7fk1LRn<-D)d638n@kvr$YxAr3OftsbzC)Gy7akTIErCn%*m5>Q`e;pD7V$&L zfO}0G6NhrYZ}+j5PQ>0c;%zb9A?cDNbh0^Q7kG!`=7gMs_}krYvN366VYAW<;Vd&Q zmC?tD2^O$4Mm>N|NqGLE#MZ}!^f|1o0*;-RSY}OF?B2s?a$Z6sG8&@c+oau=XUzOh@=9I2IPrL!0&M3L~<-0meuNGp`!LC|&XElsV zF@Tzp3#^n3gQ|+vpU!LJ4g3ZR6oaI=AuP7A7|X=2mPlL7mW~`;BE9Ysy!Wv;SYEY3 z5?#v-&HG~UZ^XM*Hn7v{ecjIJqEs3($6<3O7UcHCxyXhf^3_0fY^?$hArQk8X5|Gef^^xZU$3<0pLejBmb}W zddQ~0|97be06^zS{vH&6c7HsHP2dm54GSd^&(BWJ7yBRVzZ~~Ob|LWZKg4;iVWR{9 zz<@dn8X8a3wFn*Jr_wxT9^jswRmmrws4~d19^&$dJL;D+JIG@kX KZw_bu9sLJL9?vKM delta 2715 zcmZ8jc{tPy7an7*4TemEu|zY*Hnxym_AxUumePVSZIg((SwenDC|it?eHR%TvM)DD zwg!hbEPbB3dA;b5;*uH=0q>io9}d4<_lGboKEtl?{+ChNydzVNUM# z7A|GW!hKV*DlUINu2oBRS{dkKdnAOnF-9;D1Kg9V7oiPFT0LjF+r4@nW9&z6S3U3WlnZ29J6vkTFF<{-oP=j@R*pWcBh0WpeV^uB=cCKreO+#) zeNQ+plB})cG2wue>XKW zKeO-7cm%SS`yAy$-TTPc_&oPC!s_2*D!J@P1Frikv&K<3%_SEwbfoz0K`I`LeZosA znsh6{$riI2oxcy&-s(GN-P%nFv8mnj5Lt2FGuh;KuX)9d*BzTR?f|UOBGzILyxijp z%`Gj=IyVv?=HDIT>C7SzP?%fRnB~oBYTUG#?Bk*(-?e%9#&63DFBl6YU%zvdv1Kh$ zFP7FbOi4(MXVX{L_g)d!I5^5=cvTqVGNwn^?HR|XKrSZ=hv2t*s{OpE_Tz*AW+#FxwLReiIvl=z(@ zyoUxoDUHeGg&J^Zb-M2qRD?Nlx`+nM`ZUj|=;!E_2OCt1=TsSmJGmaxh}XQ!WxI2D z!*v0hX3Hu&l9glF>@%gUjb0xLQeQYPkD=|bR=Qh@wtq^I2<#AH5D3H#`n`fq2qOtS z2nggTNK4=p15SHh&%nX<-kn#w+^wGC^kI^4a;%j;!6=Zs=UmiREGMO^l&bid+C-?_ zS-S2w>iP`jGq~h@9}2#dop#f2?##z8{u8K3W<*T1ZSh}P;=ZAzgztiUS$)e__kMWJ zzunZA+Ua2lnV@qG;2$Ozt`!nHYq8CIR!(Xc*U`^t+AzKG7}vK7AR;b9CkI z(v{7Q6O++?!XMMW>-k!LT{t-U;1v~du&eX@+cyhcxQ;?H0|DG@Z%c*P8B`H4w~_`G zUpA{$^q@-o?aq0Vx6DMI&@q;0JaG~&+RNp0vG_2#Y;Rqj0%GA~%Wt?Kq|R=>T|03i z^rLWiED$J!fT$Zo>et-9?_J{L`C^OS2}1SrkD{JxA67O?v_R#;eer zPHl|L#i7wmby9j%!=kZL)ao6FuDZ=`u#L+4Y4X!5;^cK4<4yKHb?g1Bv8#TL3Nu+} zGTA-}t+a-j5U>V`?5b7$K#@X!v9i<@`9UPi-}9E6!-+0C_~690JG$y(yIF}#CNGksUkWX{e-VFL@HtMk? zj|b-p*fuQ;^s)CneNVc1p9zGE6QSk++~kHHvNIx*O0GfsE5&Ejgj{6j4IYLX@$MwDPi82YSjeI(( z#cvAx<5;znM;BX(;Iq&6JuaN5m`h}=8~pEkPJLI3a$jnv6`S^6qZGbW)8c((_~q=k zmx(R5{s||38ZRJ+bk^d$Wt9YhTSt6Ta#!?68*22VX!@a+cY&0UN#Pv=BXX@Y(jni_ zQ77q!2ow8`c`jbo`$9@hnt{Bj5egnPTV28R>c@i(U3ajv@qEA)IdL1qUBY!o;n+wg{6k@)MW&w z_*7eW4!^HL?2p`ux4VUc=7a%u#7g{ zl8F0k|FRca0B<`4@?d0;W_&;%V4Y}Fg_#5B4a4;#)um-9=g>wl&)vf?JyzPoLzU#Q z$Z8+~nO>dmTp|AUsO1;Ak)XHesp;lL>$_JhNMMU^aq6CTU6!mHP3{Uj*T`SuOa26s zL@$-zPhq`lty;p?=RU*OUeJF~Yt{r)>U_hJxbnUVmnz~FRBRqyu!QGBeHJ(Dw5X@XRYv2&E!;mHzs(~#7&~*6aZ>wq{g~Pv#0^mKu5|$M5_NT$ z$k%qC?|WDKcli(?Fj&;lny!K8_9d?2-)-zf1>RP$5V7L1UMQ7~ zt=s|hHI92)Cj@|?7J}CGp>1JX~M~({JvN>QP z?fNDEn;R=<1M*&f?eDCeCm}0TBA=|RR^HNmh(Y8*d& zkrWphYiohN@FstVJA&97H^-t+Zr3D^Am{!DQHfacGG!B-;yGIShsTJwEVXFboU=A^&0hj~z|> v48ws(p#R9lDJV>f;e#Fc&G=Jngcd5OBmQq?|K$I{NBRG&vlQp3z;EC`?+fw+ diff --git a/MATLAB/expression-constrainers/perturbMedium.mlx b/MATLAB/expression-constrainers/perturbMedium.mlx index eea0ed4fd6f003b6725d49c91ec72830ce5ced71..d1bf74c2fad51251982832565f0eaff3ecdd3ef5 100644 GIT binary patch delta 3157 zcmZ8jcR1XO5?wWVT~T%sMB7B~J<1ZjFQWIhdWjZ4tM|5gh!VX-EP@~|tJfrpRU(OC zwIqUYBh7Q~eeaL^X8xQr-+Xhv`OcX^SiaEHA!0f`&@&U#GFFRPQV?iW3_V8-1lb`+x7~CEi=T1b+E{NMW|?F=|`#|j2D;Ml`k1P-Y2J`Q_!?9f7v zY8(3T{ItSNkCirkT65_%Lf;Yo(bJ&m_TiPnJD4{2Aa6c_Zw7RCC%dxLfxxs^&H6}W&oj^_sT6`hygFxA*6+FRX*x6^*U#Ix1wt$9sLiXuql)rBn}Mg?w^54eSO zW1SZ}Rf=|whmtVIc`#`m+J@d<#(bRu3jt3`1jpbYNv@u!fdvNcZXoj2URT0&t$$=A zmtv9ESWODkE!P>elRgf(?-yd1wf4bbZkAkrpIoVT_1`}3)t+Hq zPVOLNYR5aj)*_dwj3s>k_2Rw~PU_^RmQ4$W2rBJH++gw1QnRKdW!ZOf68d_FG(FxG zw#a+WOB<7>9nK;$xyKe<^jzLXZi{7c1XDqq=bxEyj{rHzCsYB!g|2~ZNZ+Cw86B!V z-%RlSw#DYz@cR5=jBdr=Ky41bVM@n5YW7c-=Qv?Gy55Om$RpbvdX9Ju7mKi5pVL>0 zzR;ZgqGxYK_@a8w-bIPmUF+#Z`CH5nEM1(nwS`A*UX()`y$N5H$m)p=4Qun+W|>K7 z4G>Q0%B*`|?$QI`0s&QymdFnPqR^(E%A3|MggiX^J&IhqRVcmG0A=xVDUqtKek=1( z@z6(K4oLj=BW5pBNceCSfHh(*=6LkXrrAY0oK z=K^0~)}1qxyhpJpgmH(UUae&xeFL}H{3Tan>dHncMcS_k!EiO-*c!?QAje~N?9{iP z3&#N#=B$*4!MSEZMytK|VtQj#G@B$d4Mv1eo?h%hse(6Ny%rjxvWZM#2V+jKTP+>; zBkni)r2*}U=Rf(h*YFX z#h?fGTc@)?Fl$Im7oF5Dez+PCR27{ofq)uborKQ(a{gX*N}v z$pY2z?Rw+u8uw>oVUqlM{q~yixsDb>^cBK@PM&OoCaqc@{8}eEozS zvkU3#+^oCvs6xgWaZTdFzuu7Fb!zFE;Tjomu zwCmoAg4guKO=)fZHw^?bL1Bzr>~3B{ZVyQU3e|X0ifou5J%;PrZ{^NP|5{ zVjsbkoPc|Gt1@~wBm5>-THQI0K)vA*ddu9^y zl?k)gERzw0O58gbq%zWy#VDqMp&EKU}`eVir_6*B5 zF|O*1yx#Y4-`;r?N9u1{YLQjKlT#&TJ@DFYeRaIebamQFcDyd~^aWNJFKz_i3ooz6 zK)vothIqh3N~EsI0GWs?-0E%Z8R0w{Su_w9HX4OTmn~G+|c7qT9mr(qN-Nu8r6rjkq$?$iIa0sPAG2QA2XBC zbaBZ(rb;zp5H@wi1)vLRCY=50~g&|Tr zULsW8gt*I`mBcKlt^y-91u{p34V?FRrdtw=5)zuUEtdZ9eyhZGQU zi*HbL+cu6*&T^a8WbkTq_<4$A6d&AA4@QoI^+WrVxVSYwWH%R`IZ>mMo9p#EGKM!g zQpH*zV@wgi1-zK7E5Lp>!#Fj_0Aqr54@-pM^0b3hcrUl48UK))qdMW-Udak#H!3Yz z5tkAa@slYmj&tSVucIlim5^fNELexcYBV9?yf8Z^{oU^+u22^9$4|2igO$tDp;Jln zGb^M2m1 zUY6VzrJ*)5)-A=oCi!X|{MAD|g5HXI%A(p?zJTAS>Qyvx#pTzC<2tpEYX0BHCN%OF ztM6i&L!BE+o1bs>Imjo3at^Qz(Wx`UM4k1X+b6*g?^wd-h#X9^n3#eiE^7@NFEyviPi1W;yE_h4Z#r@Wzh&RoC z(D9v({J1t78`{(MrW3{TC7+6oupTqqa|u?{!b4ItWAi~s0}1@wSxhDGh=-Hke8TQ+ zwe_^*`q0M~w#E^RO}XBQ*BK-2d2Iub0Z*s_#?`*(Vk67PykCwvZl>xi$#BYKyiD@& z*0{ijZWfds5ZZsYnP~kek9aW@QnQ4V<@vU%B)C3Ceeo@T^I6SWunVNeC;QC6@+g0S zYky$7YV-ni%Su*iy}PZ>_b|b|dnw!0V0N(*U6A2E`L2~W6LhJju;y*d7j-#y6bxD3oO(U7GIKV!6^ z${5%>=aXmb%eS2h#E&ceL+a%rV~raba1CR|R^RmQXqC?&r8$Pn)Rg$R!9pxH0d{yvRDN!?DN#B#XeA{>qp7&f$GW>PxxR%x?S&x-5irU|n30btXDY~$iR`E6 zQ{l+DqtA!5&1~*q9293cC$vx$6Ljj3>-fvFI(ulf8DJ|e=_nJWlXow7XC2oyVcA|_ zZVBjYr>7>Sd!;0Y`Qyq2WQaf@ArRUD3_-K8&;w!F`tZp0`bG1KbZYW^G2e*XT??51W_dq6A15${b>BiJfj)kZ?6`=G?YT_3R~P7VSwieQIH zIDyCR$S-jE^IMc)LHYJQ%gow37_0M)4=k;y9iy*dFg8B7uAMCWpyyB9QZmM^$#mSo zqW%sR);=Dq$dxpD*WW=VxaI}9RX*J6Qv-Gyh9F^N%TZW6 zWF$Q%_=>$~z0$&7CaXxwXOA$3y5*UY2h^8xlp3pPaRX0zZOt_o#|n^=z}bYVY`XQS>=}$ zR@+VrUzgun0X2+;&?!%9Vn#n=e?DDxvxLzW zU+W(g#KYjE&K2^^VRdej4rP`Ch)Gv+!h5yv<4NCQ@=%Gz8@E32dGv26R@cnyL&g3w z_b@-^sutjM9-*R+egI*c4BpHB8DzcVhbvk8`@Y(=Cj@Pcq0Ta>VdakKHImlx3&pbuP-|1s@%z8^3Ot}#@Y=hXeMS*N?XF5V}=&c^q|E6ZqD zATePtzX#L#I6^lz`$dIzcWQE`$?v8wd{3?6wo0qUKkw8vWhUg@0;0W*@ug4hy*Ui! z_6h2~U7>eWje8Hny>E^`K8-X@nYHIHXj$l+&}ZEM4v^dorRZx{Z1?V}KXqQ!u*^mZ zk)NjR2U#W#o9@{}5UoMyfur_E%gY>YNX4a`N8cvR4CsStsupAC?Si4NN(AD1~Jln4ZWvnEV zo6Axxs8an>D^PVpsX27ttTY=zvb(okh|e5?qd9ri+# zgYU&Aaj!JoGC_&GPV*x99a+euWaIro{Hc_mekm>k7j6+x5g zvlXwe{m9#fK(2QPfBD&D9F^yu)QGF`#SzV^0bWUI5n;Z>FnM(Sr=Rc-?8-^eG2mqs z?C14_n_0Q6aOKI)`lYH2)H^iX#Xgw&$vWjHZfH}pBmB)Vw@Rb=-s+Lbd5ihmrtC1B z+e1}*y8(Mc!JDcJLI;3PwAP|TO^iNXdWkit~ zUy^DZ2M7eY^t#*nxY^o?J0R_S-4Px>VrX}_DT@Upha5eyVkBb?Pw5`fn;FaikR zGUMTx&uc<)dXD~}BfqWGE!we|(=W(h3sL;LxL5GzF7WbUVxM$d#^1T<8D!dnYyNnd z*LJg;veREX{#x7o&z9NcGfXz#s1LXtx|zu>i3En{arRg$*_w*2Mgv+O&4cpfnw zcnF(xUb|an@r>6bMVP%{x=1Wmp$}Nteyi8n=wzM>VXye*|3gVE_4LHWRI;-G=v%Sv z@`TtI+S>%v40f=U^d|aFW9MV!STItMzto!} zHN3SGbNIyu>iRnAobAXq1|D8RnA#v%=8ZY<7SkG~8^swLGQ-oPnWQ>yz@+w+Rk{SE z)_JL9y2XJV==%#A6f08sjLbvMkA)HU+-<<^vovYGalKE`!_WbD8g77T;%0#P^%rhWTj}j#j%K}hq(T%_9h5@71 z>Fn{GJNHq%@Eg*u1NlH_fYG-+Km~Me^%bO+aAiyWL}I^j7S&xK?HE z(wg|YO1gm8B`wqH@o9&DKG!HJ8pg)cm}QB6mvIYdcR$n?8{%Z@iKD1X3a7vikbPRdeT*?QSK+7cNJL&7U=^y$J zxOD|vZdz$|ZxG{S zIf%+U7*IBenuw8%?wbWwX7px|yzPAiHx46C7U>n=W-R_b(sTAZC5*v%Eq` z@k((oZjX2xrF-LmYH&YO`zZQorG+u20v=4?6lO-6o2DI(9vAuQWI?O8IlZ0&jL;_r zT3^Tu7W64Lwl(hO%|tDOE7duaqirOTLcE+MKaugGYieat;P+!h_ZUp z^SgfKjx`NSpUtfU8AaZ-hQ`}DuAjFQ2lMlDyC<;7*9&|Y@+XHyR&cDM6o@MNE-T(` zQA26%bCiQ@*GG?slN3;x-2A5W>l49@8&)6dv8skd)Ys5PJa`~n=mS*7CTDSi$xDY8 zQnH0`1(3CEg_!2=$4<{)WtWap;w$mQA(Ry&P5e#vhz)kti?-s|fN-7hBd5hb$6mNy z{kdFcF5g4mi^ky3 z>?43G-$^>{@^WoZ-zEKL`)%W&+Wv!& zu>se5$XzM#wzFW>XvKk^J*D^A2Vv!W*3X+ssAGBP6$UeKEbr~-qHBR9=;vBt>nhF# z>HGVFrV0!BNY5vSxNebRwy@-oLS4(1nZWFM+qDRZgv~%w{>6wsCllsz)SLp^ua%sY z3C*O@>+P)5Lfe5-ZJ7^eQA&u6{fhH{9BJO1X3nduLBm(*qXNi~5nk~0tEA6#au3% zUK0p{t_1V^?FxQ{Za?H;_HKBAgBy_-Z=1lYDsS>*a0tte!p%5r`3jzjw z6LrAh67-Qf#JJR;d?qakp>MszG`n7W%%;38x-WZk#Z&LMpuc1w9sL9LwjZPO#}CK) z>NgCaRpQar+dJ<=C9gQEr&UzE(gOovM16v7}OI^Ix|-(j|3k}6Vn zIAMJ-b^aiVyq#6;S9X7X;)b+xQM|&3x9-YILwT)Jpb!JUp5Z9Z2VGK zsjx1Lyrh>PY$79+B#RdNmhlGXzn13`Ku-b!LH~vOfBaKAtS&7KO9>IaQi%l@*8!|UVIpaSA{{UN-gdhL_ diff --git a/PythonCOBRA/.ipynb_checkpoints/QP_example-checkpoint.ipynb b/PythonCOBRA/.ipynb_checkpoints/QP_example-checkpoint.ipynb new file mode 100644 index 0000000..60c5187 --- /dev/null +++ b/PythonCOBRA/.ipynb_checkpoints/QP_example-checkpoint.ipynb @@ -0,0 +1,314 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Quadratic programming example\n", + "**Author:** Scott Campit\n", + "\n", + "This tutorial solvers an arbitrary quadratic programming (QP) problem." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# Import Gurobi Solver\n", + "import gurobipy as gp\n", + "from gurobipy import GRB" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# Example taken from https://www.gurobi.com/documentation/9.0/examples/qp_py.html \n", + "# Copyright 2020, Gurobi Optimization, LLC\n", + "\n", + "# This example formulates and solves the following simple QP model:\n", + "# minimize\n", + "# x^2 + x*y + y^2 + y*z + z^2 + 2 x\n", + "# subject to\n", + "# x + 2 y + 3 z >= 4\n", + "# x + y >= 1\n", + "# x, y, z non-negative\n", + "#\n", + "# It solves it once as a continuous model, and once as an integer model." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solving a QP Problem\n", + "First, we need to create an object that will hold the variables and constraints. This is equivalent to a structure in MATLAB." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using license file /home/scampit/software/gurobi901/linux64/gurobi.lic\n", + "Academic license - for non-commercial use only\n" + ] + } + ], + "source": [ + "# Create model object\n", + "mdl = gp.Model(\"qp_example\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we need to create variables, which we can do using the `addVar()` method, and we can assign the upper and lower bounds as well. \n", + "\n", + "In a large (ie genome-scale) problem, we would of course have a more elegant solution to assigning individual variables. However, because there are only 3 variables in this example, we'll assign them individually." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "# Create variables\n", + "x = mdl.addVar(ub=1.0, lb=0.0, name='x')\n", + "y = mdl.addVar(ub=1.0, lb=0.0, name='y')\n", + "z = mdl.addVar(ub=1.0, lb=0.0, name='z')" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "Next, let's set the objective function, which happens to be quadratic this example. We can write out the expression into a variable `obj` and use the `setObjective` function to set the objective function." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "# Set objective: x^2 + x*y + y^2 + y*z + z^2 + 2 x\n", + "obj = x*x + x*y + y*y + y*z + z*z + 2*x\n", + "mdl.setObjective(obj)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's now add our linear equations that we are using to constrain the model" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "# Add constraint: x + 2 y + 3 z <= 4\n", + "c0 = mdl.addConstr(x + 2 * y + 3 * z >= 4, \"c0\")\n", + "\n", + "# Add constraint: x + y >= 1\n", + "c1 = mdl.addConstr(x + y >= 1, \"c1\")\n", + "\n", + "# Add non-negativity constraints\n", + "c2 = mdl.addConstr(x >= 0, \"c2\")\n", + "c3 = mdl.addConstr(y >= 0, \"c3\")\n", + "c4 = mdl.addConstr(z >= 0, \"c4\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, let's solve" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (linux64)\n", + "Optimize a model with 15 rows, 3 columns and 24 nonzeros\n", + "Model fingerprint: 0xc6c7063d\n", + "Model has 5 quadratic objective terms\n", + "Coefficient statistics:\n", + " Matrix range [1e+00, 3e+00]\n", + " Objective range [2e+00, 2e+00]\n", + " QObjective range [2e+00, 2e+00]\n", + " Bounds range [1e+00, 1e+00]\n", + " RHS range [1e+00, 4e+00]\n", + "Presolve removed 13 rows and 0 columns\n", + "Presolve time: 0.02s\n", + "Presolved: 2 rows, 3 columns, 5 nonzeros\n", + "Presolved model has 5 quadratic objective terms\n", + "Ordering time: 0.00s\n", + "\n", + "Barrier statistics:\n", + " Free vars : 2\n", + " AA' NZ : 6.000e+00\n", + " Factor NZ : 1.000e+01\n", + " Factor Ops : 3.000e+01 (less than 1 second per iteration)\n", + " Threads : 1\n", + "\n", + " Objective Residual\n", + "Iter Primal Dual Primal Dual Compl Time\n", + " 0 1.69015022e+05 -1.71012100e+05 1.50e+03 3.33e+02 1.00e+06 0s\n", + " 1 3.60255402e+04 -3.91306233e+04 2.28e+02 3.82e+01 1.20e+05 0s\n", + " 2 4.14685168e+00 -4.40925173e+03 1.80e+00 4.00e-01 1.83e+03 0s\n", + " 3 2.81937163e+00 -1.92736174e+03 1.80e-06 4.00e-07 2.41e+02 0s\n", + " 4 2.81628339e+00 -1.81287557e-01 8.60e-10 1.90e-10 3.75e-01 0s\n", + " 5 2.26977145e+00 2.06670895e+00 6.89e-12 1.53e-12 2.54e-02 0s\n", + " 6 2.11498124e+00 2.11029644e+00 2.22e-16 5.55e-16 5.86e-04 0s\n", + " 7 2.11111498e+00 2.11111030e+00 2.22e-16 2.23e-16 5.85e-07 0s\n", + " 8 2.11111111e+00 2.11111111e+00 2.22e-16 0.00e+00 5.86e-10 0s\n", + "\n", + "Barrier solved model in 8 iterations and 0.03 seconds\n", + "Optimal objective 2.11111111e+00\n", + "\n", + "Objective value: 2.1111111149799298\n", + "RHS [4.0, 1.0, 0.0, 0.0, 0.0, 4.0, 1.0, 0.0, 0.0, 0.0, 4.0, 1.0, 0.0, 0.0, 0.0]\n" + ] + } + ], + "source": [ + "# Optimize the model\n", + "mdl.optimize()\n", + "\n", + "# Get the objective value\n", + "print(\"Objective value: \", mdl.objVal)\n", + "\n", + "# Get the RHS\n", + "print(\"RHS\", mdl.RHS)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Sensitivity Analysis\n", + "First, let's get some metrics out of the model, like the reduced cost and shadow price." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.7777777774265019, 1.9944607818577893, 0.0, 0.0, 0.0]\n", + "[0.22776144408201618, -0.8833496720106538, 4.18484965407449e-11]\n" + ] + } + ], + "source": [ + "# Get the shadow price associated with each constraint\n", + "print(mdl.Pi)\n", + "\n", + "# Get the reduced cost associated with each variable\n", + "print(mdl.RC)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's now modify a single variable. For instance, let's see the impact of the RHS for the first constraint on the objective coefficient if we gradually increase the value by 1." + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n" + ] + } + ], + "source": [ + "# Create a new model and hide the output table\n", + "mdl2 = mdl\n", + "mdl2.setParam('OutputFlag', False)\n", + "iterations = 10\n", + "\n", + "# Divide the RHS of c1 by 0.1 (relax the constraints) for 10 iterations\n", + "for i in range(0, iterations):\n", + " c1.RHS = c1.RHS * 0.1\n", + " mdl2.update()\n", + " mdl2.optimize()\n", + " print(mdl2.objVal)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It appears that relaxing the second constraint in the model does not affect the objective function all that much if at all." + ] + } + ], + "metadata": { + "hide_input": false, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "latex_envs": { + "LaTeX_envs_menu_present": true, + "autoclose": false, + "autocomplete": true, + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 1, + "hotkeys": { + "equation": "Ctrl-E", + "itemize": "Ctrl-I" + }, + "labels_anchors": false, + "latex_user_defs": false, + "report_style_numbering": false, + "user_envs_cfg": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/PythonCOBRA/.ipynb_checkpoints/Untitled-checkpoint.ipynb b/PythonCOBRA/.ipynb_checkpoints/Untitled-checkpoint.ipynb new file mode 100644 index 0000000..3a0877f --- /dev/null +++ b/PythonCOBRA/.ipynb_checkpoints/Untitled-checkpoint.ipynb @@ -0,0 +1,307 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Quadratic programming example\n", + "**Author:** Scott Campit\n", + "\n", + "This tutorial solvers an arbitrary quadratic programming (QP) problem." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# Import Gurobi Solver\n", + "import gurobipy as gp\n", + "from gurobipy import GRB" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# Example taken from https://www.gurobi.com/documentation/9.0/examples/qp_py.html \n", + "# Copyright 2020, Gurobi Optimization, LLC\n", + "\n", + "# This example formulates and solves the following simple QP model:\n", + "# minimize\n", + "# x^2 + x*y + y^2 + y*z + z^2 + 2 x\n", + "# subject to\n", + "# x + 2 y + 3 z >= 4\n", + "# x + y >= 1\n", + "# x, y, z non-negative\n", + "#\n", + "# It solves it once as a continuous model, and once as an integer model." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solving a QP Problem\n", + "First, we need to create an object that will hold the variables and constraints. This is equivalent to a structure in MATLAB." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using license file /home/scampit/software/gurobi901/linux64/gurobi.lic\n", + "Academic license - for non-commercial use only\n" + ] + } + ], + "source": [ + "# Create model object\n", + "mdl = gp.Model(\"qp_example\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we need to create variables, which we can do using the `addVar()` method, and we can assign the upper and lower bounds as well. \n", + "\n", + "In a large (ie genome-scale) problem, we would of course have a more elegant solution to assigning individual variables. However, because there are only 3 variables in this example, we'll assign them individually." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "# Create variables\n", + "x = mdl.addVar(ub=1.0, lb=0.0, name='x')\n", + "y = mdl.addVar(ub=1.0, lb=0.0, name='y')\n", + "z = mdl.addVar(ub=1.0, lb=0.0, name='z')" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "Next, let's set the objective function, which happens to be quadratic this example. We can write out the expression into a variable `obj` and use the `setObjective` function to set the objective function." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "# Set objective: x^2 + x*y + y^2 + y*z + z^2 + 2 x\n", + "obj = x*x + x*y + y*y + y*z + z*z + 2*x\n", + "mdl.setObjective(obj)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's now add our linear equations that we are using to constrain the model" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "# Add constraint: x + 2 y + 3 z <= 4\n", + "c0 = mdl.addConstr(x + 2 * y + 3 * z >= 4, \"c0\")\n", + "\n", + "# Add constraint: x + y >= 1\n", + "c1 = mdl.addConstr(x + y >= 1, \"c1\")\n", + "\n", + "# Add non-negativity constraints\n", + "c2 = mdl.addConstr(x >= 0, \"c2\")\n", + "c3 = mdl.addConstr(y >= 0, \"c3\")\n", + "c4 = mdl.addConstr(z >= 0, \"c4\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, let's solve" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (linux64)\n", + "Optimize a model with 15 rows, 3 columns and 24 nonzeros\n", + "Model fingerprint: 0xc6c7063d\n", + "Model has 5 quadratic objective terms\n", + "Coefficient statistics:\n", + " Matrix range [1e+00, 3e+00]\n", + " Objective range [2e+00, 2e+00]\n", + " QObjective range [2e+00, 2e+00]\n", + " Bounds range [1e+00, 1e+00]\n", + " RHS range [1e+00, 4e+00]\n", + "Presolve removed 13 rows and 0 columns\n", + "Presolve time: 0.02s\n", + "Presolved: 2 rows, 3 columns, 5 nonzeros\n", + "Presolved model has 5 quadratic objective terms\n", + "Ordering time: 0.00s\n", + "\n", + "Barrier statistics:\n", + " Free vars : 2\n", + " AA' NZ : 6.000e+00\n", + " Factor NZ : 1.000e+01\n", + " Factor Ops : 3.000e+01 (less than 1 second per iteration)\n", + " Threads : 1\n", + "\n", + " Objective Residual\n", + "Iter Primal Dual Primal Dual Compl Time\n", + " 0 1.69015022e+05 -1.71012100e+05 1.50e+03 3.33e+02 1.00e+06 0s\n", + " 1 3.60255402e+04 -3.91306233e+04 2.28e+02 3.82e+01 1.20e+05 0s\n", + " 2 4.14685168e+00 -4.40925173e+03 1.80e+00 4.00e-01 1.83e+03 0s\n", + " 3 2.81937163e+00 -1.92736174e+03 1.80e-06 4.00e-07 2.41e+02 0s\n", + " 4 2.81628339e+00 -1.81287557e-01 8.60e-10 1.90e-10 3.75e-01 0s\n", + " 5 2.26977145e+00 2.06670895e+00 6.89e-12 1.53e-12 2.54e-02 0s\n", + " 6 2.11498124e+00 2.11029644e+00 2.22e-16 5.55e-16 5.86e-04 0s\n", + " 7 2.11111498e+00 2.11111030e+00 2.22e-16 2.23e-16 5.85e-07 0s\n", + " 8 2.11111111e+00 2.11111111e+00 2.22e-16 0.00e+00 5.86e-10 0s\n", + "\n", + "Barrier solved model in 8 iterations and 0.03 seconds\n", + "Optimal objective 2.11111111e+00\n", + "\n", + "Objective value: 2.1111111149799298\n", + "RHS [4.0, 1.0, 0.0, 0.0, 0.0, 4.0, 1.0, 0.0, 0.0, 0.0, 4.0, 1.0, 0.0, 0.0, 0.0]\n" + ] + } + ], + "source": [ + "# Optimize the model\n", + "mdl.optimize()\n", + "\n", + "# Get the objective value\n", + "print(\"Objective value: \", mdl.objVal)\n", + "\n", + "# Get the RHS\n", + "print(\"RHS\", mdl.RHS)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Sensitivity Analysis\n", + "First, let's get some metrics out of the model, like the reduced cost and shadow price." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.7777777774265019, 1.9944607818577893, 0.0, 0.0, 0.0]\n", + "[0.22776144408201618, -0.8833496720106538, 4.18484965407449e-11]\n" + ] + } + ], + "source": [ + "# Get the shadow price associated with each constraint\n", + "print(mdl.Pi)\n", + "\n", + "# Get the reduced cost associated with each variable\n", + "print(mdl.RC)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's now modify a single variable. For instance, let's see the impact of the RHS for the first constraint on the objective coefficient if we gradually increase the value by 1." + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n" + ] + } + ], + "source": [ + "# Create a new model and hide the output table\n", + "mdl2 = mdl\n", + "mdl2.setParam('OutputFlag', False)\n", + "iterations = 10\n", + "\n", + "# Divide the RHS of c1 by 0.1 (relax the constraints) for 10 iterations\n", + "for i in range(0, iterations):\n", + " c1.RHS = c1.RHS * 0.1\n", + " mdl2.update()\n", + " mdl2.optimize()\n", + " print(mdl2.objVal)" + ] + } + ], + "metadata": { + "hide_input": false, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "latex_envs": { + "LaTeX_envs_menu_present": true, + "autoclose": false, + "autocomplete": true, + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 1, + "hotkeys": { + "equation": "Ctrl-E", + "itemize": "Ctrl-I" + }, + "labels_anchors": false, + "latex_user_defs": false, + "report_style_numbering": false, + "user_envs_cfg": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/PythonCOBRA/QP_example.html b/PythonCOBRA/QP_example.html new file mode 100644 index 0000000..ded5031 --- /dev/null +++ b/PythonCOBRA/QP_example.html @@ -0,0 +1,16694 @@ + + + + +QP_example + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+
+
+

Quadratic programming example

Author: Scott Campit

+

This tutorial solvers an arbitrary quadratic programming (QP) problem.

+ +
+
+
+
+
+
In [5]:
+
+
+
# Import Gurobi Solver
+import gurobipy as gp
+from gurobipy import GRB
+
+ +
+
+
+ +
+
+
+
In [6]:
+
+
+
# Example taken from https://www.gurobi.com/documentation/9.0/examples/qp_py.html 
+# Copyright 2020, Gurobi Optimization, LLC
+
+# This example formulates and solves the following simple QP model:
+#  minimize
+#      x^2 + x*y + y^2 + y*z + z^2 + 2 x
+#  subject to
+#      x + 2 y + 3 z >= 4
+#      x +   y       >= 1
+#      x, y, z non-negative
+#
+# It solves it once as a continuous model, and once as an integer model.
+
+ +
+
+
+ +
+
+
+
+

Solving a QP Problem

First, we need to create an object that will hold the variables and constraints. This is equivalent to a structure in MATLAB.

+ +
+
+
+
+
+
In [9]:
+
+
+
# Create model object
+mdl = gp.Model("qp_example")
+
+ +
+
+
+ +
+
+ + +
+ +
+ + +
+
Using license file /home/scampit/software/gurobi901/linux64/gurobi.lic
+Academic license - for non-commercial use only
+
+
+
+ +
+
+ +
+
+
+
+

Next we need to create variables, which we can do using the addVar() method, and we can assign the upper and lower bounds as well.

+

In a large (ie genome-scale) problem, we would of course have a more elegant solution to assigning individual variables. However, because there are only 3 variables in this example, we'll assign them individually.

+ +
+
+
+
+
+
In [10]:
+
+
+
# Create variables
+x = mdl.addVar(ub=1.0, lb=0.0, name='x')
+y = mdl.addVar(ub=1.0, lb=0.0, name='y')
+z = mdl.addVar(ub=1.0, lb=0.0, name='z')
+
+ +
+
+
+ +
Next, let's set the objective function, which happens to be quadratic this example. We can write out the expression into a variable `obj` and use the `setObjective` function to set the objective function. +
+
+
In [11]:
+
+
+
# Set objective: x^2 + x*y + y^2 + y*z + z^2 + 2 x
+obj = x*x + x*y + y*y + y*z + z*z + 2*x
+mdl.setObjective(obj)
+
+ +
+
+
+ +
+
+
+
+

Let's now add our linear equations that we are using to constrain the model

+ +
+
+
+
+
+
In [22]:
+
+
+
# Add constraint: x + 2 y + 3 z <= 4
+c0 = mdl.addConstr(x + 2 * y + 3 * z >= 4, "c0")
+
+# Add constraint: x + y >= 1
+c1 = mdl.addConstr(x + y >= 1, "c1")
+
+# Add non-negativity constraints
+c2 = mdl.addConstr(x >= 0, "c2")
+c3 = mdl.addConstr(y >= 0, "c3")
+c4 = mdl.addConstr(z >= 0, "c4")
+
+ +
+
+
+ +
+
+
+
+

Finally, let's solve

+ +
+
+
+
+
+
In [23]:
+
+
+
# Optimize the model
+mdl.optimize()
+
+# Get the objective value
+print("Objective value: ", mdl.objVal)
+
+# Get the RHS
+print("RHS", mdl.RHS)
+
+ +
+
+
+ +
+
+ + +
+ +
+ + +
+
Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (linux64)
+Optimize a model with 15 rows, 3 columns and 24 nonzeros
+Model fingerprint: 0xc6c7063d
+Model has 5 quadratic objective terms
+Coefficient statistics:
+  Matrix range     [1e+00, 3e+00]
+  Objective range  [2e+00, 2e+00]
+  QObjective range [2e+00, 2e+00]
+  Bounds range     [1e+00, 1e+00]
+  RHS range        [1e+00, 4e+00]
+Presolve removed 13 rows and 0 columns
+Presolve time: 0.02s
+Presolved: 2 rows, 3 columns, 5 nonzeros
+Presolved model has 5 quadratic objective terms
+Ordering time: 0.00s
+
+Barrier statistics:
+ Free vars  : 2
+ AA' NZ     : 6.000e+00
+ Factor NZ  : 1.000e+01
+ Factor Ops : 3.000e+01 (less than 1 second per iteration)
+ Threads    : 1
+
+                  Objective                Residual
+Iter       Primal          Dual         Primal    Dual     Compl     Time
+   0   1.69015022e+05 -1.71012100e+05  1.50e+03 3.33e+02  1.00e+06     0s
+   1   3.60255402e+04 -3.91306233e+04  2.28e+02 3.82e+01  1.20e+05     0s
+   2   4.14685168e+00 -4.40925173e+03  1.80e+00 4.00e-01  1.83e+03     0s
+   3   2.81937163e+00 -1.92736174e+03  1.80e-06 4.00e-07  2.41e+02     0s
+   4   2.81628339e+00 -1.81287557e-01  8.60e-10 1.90e-10  3.75e-01     0s
+   5   2.26977145e+00  2.06670895e+00  6.89e-12 1.53e-12  2.54e-02     0s
+   6   2.11498124e+00  2.11029644e+00  2.22e-16 5.55e-16  5.86e-04     0s
+   7   2.11111498e+00  2.11111030e+00  2.22e-16 2.23e-16  5.85e-07     0s
+   8   2.11111111e+00  2.11111111e+00  2.22e-16 0.00e+00  5.86e-10     0s
+
+Barrier solved model in 8 iterations and 0.03 seconds
+Optimal objective 2.11111111e+00
+
+Objective value:  2.1111111149799298
+RHS [4.0, 1.0, 0.0, 0.0, 0.0, 4.0, 1.0, 0.0, 0.0, 0.0, 4.0, 1.0, 0.0, 0.0, 0.0]
+
+
+
+ +
+
+ +
+
+
+
+

Sensitivity Analysis

First, let's get some metrics out of the model, like the reduced cost and shadow price.

+ +
+
+
+
+
+
In [18]:
+
+
+
# Get the shadow price associated with each constraint
+print(mdl.Pi)
+
+# Get the reduced cost associated with each variable
+print(mdl.RC)
+
+ +
+
+
+ +
+
+ + +
+ +
+ + +
+
[0.7777777774265019, 1.9944607818577893, 0.0, 0.0, 0.0]
+[0.22776144408201618, -0.8833496720106538, 4.18484965407449e-11]
+
+
+
+ +
+
+ +
+
+
+
+

Let's now modify a single variable. For instance, let's see the impact of the RHS for the first constraint on the objective coefficient if we gradually increase the value by 1.

+ +
+
+
+
+
+
In [47]:
+
+
+
# Create a new model and hide the output table
+mdl2 = mdl
+mdl2.setParam('OutputFlag', False)
+iterations = 10
+
+# Divide the RHS of c1 by 0.1 (relax the constraints) for 10 iterations
+for i in range(0, iterations):
+    c1.RHS = c1.RHS * 0.1
+    mdl2.update()
+    mdl2.optimize()
+    print(mdl2.objVal)
+
+ +
+
+
+ +
+
+ + +
+ +
+ + +
+
2.1111111149799298
+2.1111111149799298
+2.1111111149799298
+2.1111111149799298
+2.1111111149799298
+2.1111111149799298
+2.1111111149799298
+2.1111111149799298
+2.1111111149799298
+2.1111111149799298
+
+
+
+ +
+
+ +
+
+
+
+

It appears that relaxing the second constraint in the model does not affect the objective function all that much if at all.

+ +
+
+
+
+
+ + + + + + diff --git a/PythonCOBRA/QP_example.ipynb b/PythonCOBRA/QP_example.ipynb new file mode 100644 index 0000000..60c5187 --- /dev/null +++ b/PythonCOBRA/QP_example.ipynb @@ -0,0 +1,314 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Quadratic programming example\n", + "**Author:** Scott Campit\n", + "\n", + "This tutorial solvers an arbitrary quadratic programming (QP) problem." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# Import Gurobi Solver\n", + "import gurobipy as gp\n", + "from gurobipy import GRB" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# Example taken from https://www.gurobi.com/documentation/9.0/examples/qp_py.html \n", + "# Copyright 2020, Gurobi Optimization, LLC\n", + "\n", + "# This example formulates and solves the following simple QP model:\n", + "# minimize\n", + "# x^2 + x*y + y^2 + y*z + z^2 + 2 x\n", + "# subject to\n", + "# x + 2 y + 3 z >= 4\n", + "# x + y >= 1\n", + "# x, y, z non-negative\n", + "#\n", + "# It solves it once as a continuous model, and once as an integer model." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solving a QP Problem\n", + "First, we need to create an object that will hold the variables and constraints. This is equivalent to a structure in MATLAB." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using license file /home/scampit/software/gurobi901/linux64/gurobi.lic\n", + "Academic license - for non-commercial use only\n" + ] + } + ], + "source": [ + "# Create model object\n", + "mdl = gp.Model(\"qp_example\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we need to create variables, which we can do using the `addVar()` method, and we can assign the upper and lower bounds as well. \n", + "\n", + "In a large (ie genome-scale) problem, we would of course have a more elegant solution to assigning individual variables. However, because there are only 3 variables in this example, we'll assign them individually." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "# Create variables\n", + "x = mdl.addVar(ub=1.0, lb=0.0, name='x')\n", + "y = mdl.addVar(ub=1.0, lb=0.0, name='y')\n", + "z = mdl.addVar(ub=1.0, lb=0.0, name='z')" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "Next, let's set the objective function, which happens to be quadratic this example. We can write out the expression into a variable `obj` and use the `setObjective` function to set the objective function." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "# Set objective: x^2 + x*y + y^2 + y*z + z^2 + 2 x\n", + "obj = x*x + x*y + y*y + y*z + z*z + 2*x\n", + "mdl.setObjective(obj)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's now add our linear equations that we are using to constrain the model" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "# Add constraint: x + 2 y + 3 z <= 4\n", + "c0 = mdl.addConstr(x + 2 * y + 3 * z >= 4, \"c0\")\n", + "\n", + "# Add constraint: x + y >= 1\n", + "c1 = mdl.addConstr(x + y >= 1, \"c1\")\n", + "\n", + "# Add non-negativity constraints\n", + "c2 = mdl.addConstr(x >= 0, \"c2\")\n", + "c3 = mdl.addConstr(y >= 0, \"c3\")\n", + "c4 = mdl.addConstr(z >= 0, \"c4\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, let's solve" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (linux64)\n", + "Optimize a model with 15 rows, 3 columns and 24 nonzeros\n", + "Model fingerprint: 0xc6c7063d\n", + "Model has 5 quadratic objective terms\n", + "Coefficient statistics:\n", + " Matrix range [1e+00, 3e+00]\n", + " Objective range [2e+00, 2e+00]\n", + " QObjective range [2e+00, 2e+00]\n", + " Bounds range [1e+00, 1e+00]\n", + " RHS range [1e+00, 4e+00]\n", + "Presolve removed 13 rows and 0 columns\n", + "Presolve time: 0.02s\n", + "Presolved: 2 rows, 3 columns, 5 nonzeros\n", + "Presolved model has 5 quadratic objective terms\n", + "Ordering time: 0.00s\n", + "\n", + "Barrier statistics:\n", + " Free vars : 2\n", + " AA' NZ : 6.000e+00\n", + " Factor NZ : 1.000e+01\n", + " Factor Ops : 3.000e+01 (less than 1 second per iteration)\n", + " Threads : 1\n", + "\n", + " Objective Residual\n", + "Iter Primal Dual Primal Dual Compl Time\n", + " 0 1.69015022e+05 -1.71012100e+05 1.50e+03 3.33e+02 1.00e+06 0s\n", + " 1 3.60255402e+04 -3.91306233e+04 2.28e+02 3.82e+01 1.20e+05 0s\n", + " 2 4.14685168e+00 -4.40925173e+03 1.80e+00 4.00e-01 1.83e+03 0s\n", + " 3 2.81937163e+00 -1.92736174e+03 1.80e-06 4.00e-07 2.41e+02 0s\n", + " 4 2.81628339e+00 -1.81287557e-01 8.60e-10 1.90e-10 3.75e-01 0s\n", + " 5 2.26977145e+00 2.06670895e+00 6.89e-12 1.53e-12 2.54e-02 0s\n", + " 6 2.11498124e+00 2.11029644e+00 2.22e-16 5.55e-16 5.86e-04 0s\n", + " 7 2.11111498e+00 2.11111030e+00 2.22e-16 2.23e-16 5.85e-07 0s\n", + " 8 2.11111111e+00 2.11111111e+00 2.22e-16 0.00e+00 5.86e-10 0s\n", + "\n", + "Barrier solved model in 8 iterations and 0.03 seconds\n", + "Optimal objective 2.11111111e+00\n", + "\n", + "Objective value: 2.1111111149799298\n", + "RHS [4.0, 1.0, 0.0, 0.0, 0.0, 4.0, 1.0, 0.0, 0.0, 0.0, 4.0, 1.0, 0.0, 0.0, 0.0]\n" + ] + } + ], + "source": [ + "# Optimize the model\n", + "mdl.optimize()\n", + "\n", + "# Get the objective value\n", + "print(\"Objective value: \", mdl.objVal)\n", + "\n", + "# Get the RHS\n", + "print(\"RHS\", mdl.RHS)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Sensitivity Analysis\n", + "First, let's get some metrics out of the model, like the reduced cost and shadow price." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.7777777774265019, 1.9944607818577893, 0.0, 0.0, 0.0]\n", + "[0.22776144408201618, -0.8833496720106538, 4.18484965407449e-11]\n" + ] + } + ], + "source": [ + "# Get the shadow price associated with each constraint\n", + "print(mdl.Pi)\n", + "\n", + "# Get the reduced cost associated with each variable\n", + "print(mdl.RC)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's now modify a single variable. For instance, let's see the impact of the RHS for the first constraint on the objective coefficient if we gradually increase the value by 1." + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n" + ] + } + ], + "source": [ + "# Create a new model and hide the output table\n", + "mdl2 = mdl\n", + "mdl2.setParam('OutputFlag', False)\n", + "iterations = 10\n", + "\n", + "# Divide the RHS of c1 by 0.1 (relax the constraints) for 10 iterations\n", + "for i in range(0, iterations):\n", + " c1.RHS = c1.RHS * 0.1\n", + " mdl2.update()\n", + " mdl2.optimize()\n", + " print(mdl2.objVal)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It appears that relaxing the second constraint in the model does not affect the objective function all that much if at all." + ] + } + ], + "metadata": { + "hide_input": false, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "latex_envs": { + "LaTeX_envs_menu_present": true, + "autoclose": false, + "autocomplete": true, + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 1, + "hotkeys": { + "equation": "Ctrl-E", + "itemize": "Ctrl-I" + }, + "labels_anchors": false, + "latex_user_defs": false, + "report_style_numbering": false, + "user_envs_cfg": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/PythonCOBRA/Untitled.ipynb b/PythonCOBRA/Untitled.ipynb new file mode 100644 index 0000000..3a0877f --- /dev/null +++ b/PythonCOBRA/Untitled.ipynb @@ -0,0 +1,307 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Quadratic programming example\n", + "**Author:** Scott Campit\n", + "\n", + "This tutorial solvers an arbitrary quadratic programming (QP) problem." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# Import Gurobi Solver\n", + "import gurobipy as gp\n", + "from gurobipy import GRB" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# Example taken from https://www.gurobi.com/documentation/9.0/examples/qp_py.html \n", + "# Copyright 2020, Gurobi Optimization, LLC\n", + "\n", + "# This example formulates and solves the following simple QP model:\n", + "# minimize\n", + "# x^2 + x*y + y^2 + y*z + z^2 + 2 x\n", + "# subject to\n", + "# x + 2 y + 3 z >= 4\n", + "# x + y >= 1\n", + "# x, y, z non-negative\n", + "#\n", + "# It solves it once as a continuous model, and once as an integer model." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solving a QP Problem\n", + "First, we need to create an object that will hold the variables and constraints. This is equivalent to a structure in MATLAB." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using license file /home/scampit/software/gurobi901/linux64/gurobi.lic\n", + "Academic license - for non-commercial use only\n" + ] + } + ], + "source": [ + "# Create model object\n", + "mdl = gp.Model(\"qp_example\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we need to create variables, which we can do using the `addVar()` method, and we can assign the upper and lower bounds as well. \n", + "\n", + "In a large (ie genome-scale) problem, we would of course have a more elegant solution to assigning individual variables. However, because there are only 3 variables in this example, we'll assign them individually." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "# Create variables\n", + "x = mdl.addVar(ub=1.0, lb=0.0, name='x')\n", + "y = mdl.addVar(ub=1.0, lb=0.0, name='y')\n", + "z = mdl.addVar(ub=1.0, lb=0.0, name='z')" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "Next, let's set the objective function, which happens to be quadratic this example. We can write out the expression into a variable `obj` and use the `setObjective` function to set the objective function." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "# Set objective: x^2 + x*y + y^2 + y*z + z^2 + 2 x\n", + "obj = x*x + x*y + y*y + y*z + z*z + 2*x\n", + "mdl.setObjective(obj)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's now add our linear equations that we are using to constrain the model" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "# Add constraint: x + 2 y + 3 z <= 4\n", + "c0 = mdl.addConstr(x + 2 * y + 3 * z >= 4, \"c0\")\n", + "\n", + "# Add constraint: x + y >= 1\n", + "c1 = mdl.addConstr(x + y >= 1, \"c1\")\n", + "\n", + "# Add non-negativity constraints\n", + "c2 = mdl.addConstr(x >= 0, \"c2\")\n", + "c3 = mdl.addConstr(y >= 0, \"c3\")\n", + "c4 = mdl.addConstr(z >= 0, \"c4\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, let's solve" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (linux64)\n", + "Optimize a model with 15 rows, 3 columns and 24 nonzeros\n", + "Model fingerprint: 0xc6c7063d\n", + "Model has 5 quadratic objective terms\n", + "Coefficient statistics:\n", + " Matrix range [1e+00, 3e+00]\n", + " Objective range [2e+00, 2e+00]\n", + " QObjective range [2e+00, 2e+00]\n", + " Bounds range [1e+00, 1e+00]\n", + " RHS range [1e+00, 4e+00]\n", + "Presolve removed 13 rows and 0 columns\n", + "Presolve time: 0.02s\n", + "Presolved: 2 rows, 3 columns, 5 nonzeros\n", + "Presolved model has 5 quadratic objective terms\n", + "Ordering time: 0.00s\n", + "\n", + "Barrier statistics:\n", + " Free vars : 2\n", + " AA' NZ : 6.000e+00\n", + " Factor NZ : 1.000e+01\n", + " Factor Ops : 3.000e+01 (less than 1 second per iteration)\n", + " Threads : 1\n", + "\n", + " Objective Residual\n", + "Iter Primal Dual Primal Dual Compl Time\n", + " 0 1.69015022e+05 -1.71012100e+05 1.50e+03 3.33e+02 1.00e+06 0s\n", + " 1 3.60255402e+04 -3.91306233e+04 2.28e+02 3.82e+01 1.20e+05 0s\n", + " 2 4.14685168e+00 -4.40925173e+03 1.80e+00 4.00e-01 1.83e+03 0s\n", + " 3 2.81937163e+00 -1.92736174e+03 1.80e-06 4.00e-07 2.41e+02 0s\n", + " 4 2.81628339e+00 -1.81287557e-01 8.60e-10 1.90e-10 3.75e-01 0s\n", + " 5 2.26977145e+00 2.06670895e+00 6.89e-12 1.53e-12 2.54e-02 0s\n", + " 6 2.11498124e+00 2.11029644e+00 2.22e-16 5.55e-16 5.86e-04 0s\n", + " 7 2.11111498e+00 2.11111030e+00 2.22e-16 2.23e-16 5.85e-07 0s\n", + " 8 2.11111111e+00 2.11111111e+00 2.22e-16 0.00e+00 5.86e-10 0s\n", + "\n", + "Barrier solved model in 8 iterations and 0.03 seconds\n", + "Optimal objective 2.11111111e+00\n", + "\n", + "Objective value: 2.1111111149799298\n", + "RHS [4.0, 1.0, 0.0, 0.0, 0.0, 4.0, 1.0, 0.0, 0.0, 0.0, 4.0, 1.0, 0.0, 0.0, 0.0]\n" + ] + } + ], + "source": [ + "# Optimize the model\n", + "mdl.optimize()\n", + "\n", + "# Get the objective value\n", + "print(\"Objective value: \", mdl.objVal)\n", + "\n", + "# Get the RHS\n", + "print(\"RHS\", mdl.RHS)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Sensitivity Analysis\n", + "First, let's get some metrics out of the model, like the reduced cost and shadow price." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.7777777774265019, 1.9944607818577893, 0.0, 0.0, 0.0]\n", + "[0.22776144408201618, -0.8833496720106538, 4.18484965407449e-11]\n" + ] + } + ], + "source": [ + "# Get the shadow price associated with each constraint\n", + "print(mdl.Pi)\n", + "\n", + "# Get the reduced cost associated with each variable\n", + "print(mdl.RC)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's now modify a single variable. For instance, let's see the impact of the RHS for the first constraint on the objective coefficient if we gradually increase the value by 1." + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n", + "2.1111111149799298\n" + ] + } + ], + "source": [ + "# Create a new model and hide the output table\n", + "mdl2 = mdl\n", + "mdl2.setParam('OutputFlag', False)\n", + "iterations = 10\n", + "\n", + "# Divide the RHS of c1 by 0.1 (relax the constraints) for 10 iterations\n", + "for i in range(0, iterations):\n", + " c1.RHS = c1.RHS * 0.1\n", + " mdl2.update()\n", + " mdl2.optimize()\n", + " print(mdl2.objVal)" + ] + } + ], + "metadata": { + "hide_input": false, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "latex_envs": { + "LaTeX_envs_menu_present": true, + "autoclose": false, + "autocomplete": true, + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 1, + "hotkeys": { + "equation": "Ctrl-E", + "itemize": "Ctrl-I" + }, + "labels_anchors": false, + "latex_user_defs": false, + "report_style_numbering": false, + "user_envs_cfg": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/PythonCOBRA/constrain_flux_regulation.py b/PythonCOBRA/constrain_flux_regulation.py index fedb7ac..e26c65a 100644 --- a/PythonCOBRA/constrain_flux_regulation.py +++ b/PythonCOBRA/constrain_flux_regulation.py @@ -1,16 +1,44 @@ """ - +constrain_flux_regulation integrates gene expression data into the metabolic model and outputs a flux distribution +that is constrained by gene expression data. """ import cobra -import gurobi - -def constrain_flux_regulation(model, geneExpressionObj): - lp = gurobipy.Model() - lp.Params.OutputFlag = 0 - lp.Params.FeasibilityTol = 1E-9 - lp.Params.OptimalityTol = 1E-9 - lp.addVar(lb=model.get('lb'), +import gurobipy as gp +import numpy as np +import scipy.sparse as sp + +def constrain_flux_regulation(modelFileName, geneExpressionObj): + """ + + :param modelFileName: + :param geneExpressionObj: + :return: + """ + + # Read in common file types for COBRA models + if modelFileName.split('.')[-1] is ('sbml' or 'xml'): + cobra_model = cobra.io.read_sbml_model(modelFileName) + elif modelFileName.split('.')[-1] is 'mat': + cobra_model = cobra.io.load_matlab_model(modelFileName) + + # Create a model structure + model = gp.Model("ConstrainedModel") + + # Build the A matrix from the stoichiometric matrix + S = cobra.util.create_stoichiometric_matrix(cobra_model) + A = sp.csr_matrix(S, (len(cobra_model.metabolites), len(cobra_model.reactions))) + lb = + ub = + obj = + # Build the rhs vector + rhs = np.zeros(np.shape(cobra_model.reactions)) + + # Specify additional model parameters + model.Params.OutputFlag = 0 + model.Params.FeasibilityTol = 1E-9 + model.Params.OptimalityTol = 1E-9 + model.addVar(lb=model.get('lb'), ub=model.get('ub'), obj=model.get('c'), rhs=model.get('b')) diff --git a/PythonCOBRA/knockOut.py b/PythonCOBRA/knockOut.py new file mode 100644 index 0000000..f0f87a4 --- /dev/null +++ b/PythonCOBRA/knockOut.py @@ -0,0 +1,59 @@ +""" +Perturbation analyses +This script performs several perturbation analyses using COBRA-based metabolic models. + +@author: Scott Campit +""" + +import cobra + +def geneKnockOut(path2model): + """ + geneKnockOut returns the objective values resulting from single gene knockouts in a COBRA model. + + :param path2model: A string denoting the path to a metabolic reconstruction. + :return geneKO: A list of single gene knockout objective values. + """ + if modelFileName.split('.')[-1] is ('sbml' or 'xml'): + cobra_model = cobra.io.read_sbml_model(path2model) + elif modelFileName.split('.')[-1] is 'mat': + cobra_model = cobra.io.load_matlab_model(path2model) + elif modelFileName.split('.')[-1] is 'json': + cobra_model = cobra.io.load_json_model(path2model) + + geneKO = [] + for gene in cobra_model.genes: + tmp = cobra_model.copy() + tmp2 = tmp.genes.get_by_id(gene.id).knockout + geneKO.append(tmp2.optimize().objective_value) + return geneKO + +def reactionKnockOut(path2model): + """ + reactionKnockOut returns the objective values resulting from single reaction knockouts in a COBRA model. + + :param path2model: A string denoting the path to a metabolic reconstruction. + :return rxnKO: A list of single reaction knockout objective values. + """ + if modelFileName.split('.')[-1] is ('sbml' or 'xml'): + cobra_model = cobra.io.read_sbml_model(path2model) + elif modelFileName.split('.')[-1] is 'mat': + cobra_model = cobra.io.load_matlab_model(path2model) + elif modelFileName.split('.')[-1] is 'json': + cobra_model = cobra.io.load_json_model(path2model) + + rxnKO = [] + for rxn in cobra_model.rxns: + tmp = cobra_model.copy() + tmp2 = tmp.reactions.rxn.bounds = 0, 0 + rxnKO.append(tmp2.optimize().objective_value) + return rxnKO + +def mediumPerturbationAnalyses(path2models, mediumComponents): + """ + + :param path2models: + :param mediumComponents: + :return: + """ + pass \ No newline at end of file