KarrLab / wc_sim

A multi-algorithmic simulator for whole-cell models
MIT License
6 stars 2 forks source link

Double-check coefficient scaling in unscale_conv_opt_solution() in wc_sim/submodels/dfba.py #113

Closed artgoldberg closed 3 years ago

artgoldberg commented 3 years ago

Yin Hoon writes:

It was really nice to see the things we discussed implemented in neat code. Below are my comments:

There is a bug in line 641 under the method unscale_conv_opt_solution. Specifically, in contrast to bound scaling where we do the opposite (division vs multiplication) to un-scale, for coefficient scaling we do the same operation. This may seem unintuitive so I’ll use a simple example to illustrate that:

R1: 5A + B -> C                 Bounds: (0, 2)
Obj_rxn: -> 5A                 Bounds: (0, None)

Obj function: Maximize Obj_rxn

Let’s assume there are other unbounded reactions that will balance B and C. To balance A:
-5* v_R1 + 5*v_Obj_rxn = 0

Therefore the true answer here is:
v_R1 = v_Obj_rxn = 2

Case 1: Scaling bounds by a factor of 10

R1: 5A + B -> C                 Bounds: (0, 2*10)
Obj_rxn: -> 5A                 Bounds: (0, None)

The scaled solution becomes:
v_R1 = v_Obj_rxn = 20

To unscale and get the true answer, we do the opposite:
v_R1 = v_Obj_rxn = 20/10

Case 2: Scaling coefficient of objective reaction by 10

R1: 5A + B -> C                 Bounds: (0, 2)
Obj_rxn: -> 5*10A           Bounds: (0, None)

The steady state for A becomes:
-5* v_R1 + 5*10*v_Obj_rxn = 0

The scaled solution becomes:
v_R1 = 2 (same as the true answer)
v_Obj_rxn = 0.2 (reduced by a magnitude of 10)

Therefore, to unscale and get the true answer, we have to do the same operation of multiplication:
v_Obj_rxn = 0.2 * 10
artgoldberg commented 3 years ago

Not a problem.

From my email to Yin Hoon:

I believe that the method unscale_conv_opt_solution() is correct as it stands. My analysis differs from your example below because the terms in the dfba_obj_expr (to use the WC Sim terminology) are reactions (wc_lang.Reaction and wc_lang.DfbaObjReaction instances) not species.

Take the example from tests/submodels/test_dfba.py, the model in tests/submodels/fixtures/dfba_test_model.xlsx. Its objective is

dfba_obj_expr = 1 * biomass_reaction

where biomass_reaction is 'm3[c] ==> '. The dFBA problem is solved in "tests/submodels/fixtures/Solutions to test dFBA models, by hand.txt”. In lines 31 - 50 I obtain the reaction fluxes.

If the problem is changed by scaling the dFBA objective coefficients by s, then the entire solution remains the same, except that the objective functions given to the solver — let’s call it objective_scaled — is

s * biomass_reaction.

Since the solution to the remainder of the FBA conv. opt. model remains unchanged, the solver reports that

objective_scaled = s * biomass_reaction = s * dfba_obj_expr

Therefore, the unscaled dfba_obj_expr is given by

dfba_obj_expr = objective_scaled / s

My unit tests confirm this in test "V: Modify II by scaling the coefficients of objective function's terms by 10”, following line 709, which uses the same inputs as test II, but scales the coefficients in dFBA objectives by 10. The test confirms that the same value for the objective is obtained as in test II.