Added initial metabolic polytope simplification with HiGHS

[DimitriosPavlouGR]

Summary

This PR adds simplification support for metabolic networks as a preprocessing step for sampling and volume approximation.

It includes:

• A new class for metabolic network polytopes, MetabolicPolytope
• A BiGG JSON model parser for metabolic networks
• An LP based simplification routine built with HiGHS that relaxes redundant flux bounds and fixes degenerate dimensions.

Metabolic Polytope

convex_bodies/metabolic_polytope.h defines MetabolicPolytope, a metabolic network instance.

Parser

io/bigg_parser.hpp loads a standard BiGG JSON model and stores it in a MetabolicPolytope.

Simplification

simplification/warm_start.hpp simplifies the MetabolicPolytope by converting degenerate dimensions to equalities and redundant inequalities to relaxed +inf/-inf bounds.

The simplification happens in three key steps:

• build_lp - constructs a HiGHS LP problem that describes the feasible region of the polytope.
• build_simplified_polytope - constructs the output MetabolicPolytope directly from the HiGHS model, once at the end of the simplification.
• simplify - iterates over all the variables until no simplifications can happen. For each variable it solves four LP problems to detect whether a bound is redundant, in which case it is relaxed, and whether the variable is pinned to a single value, in which case it is turned into an equality.

Tests

test_metabolic_polytope.cpp covers construction, test_highs covers general highs functionalities, and test_simplification covers simplification on general shapes such as a cube or a simplex.

Benchmarks

benchmark_simplification.cpp runs simplification across selected models present in test/bigg.

[vissarion]

Thanks for this PR!

I wrote a couple of comments inline.

Also, json.h is an external file and should be placed in external directory together with some attribution note, license etc

[vissarion]

this is not used and can simply computed as getNumEqualities()>0 if needed. So should we remove it?

[vissarion]

should be in above line for readability. What is RHS?

[DimitriosPavlouGR]

I will change the placement, RHS stands for right hand side of the equation A_eq x = b_eq. Do you want me to update the comment?

[vissarion]

volesti does not uses namespaces this way, it is useful to use them but need to be part of a larger plan

[vissarion]

similar to other places where namespaces are introduced in this PR

[DimitriosPavlouGR]

Should I remove them entirely?

[vissarion]

So you end up with an empty polytope here? Could you add a not that extreme example? e.g. a 3d cube in a 10d space, a square in 3d?

[DimitriosPavlouGR]

No, actually this test creates a 2d dimensional cube with d dimensions that are degenerate (we have 1 <= x_k <= 1+e^-12, for d variables). So after simplification the derived polytope should be a d-dimensional cube.

Notice that:

• result.dims_fixed counts the number of added equalities, which should be d after simplification.
• result.bounds_relaxed counts the number of relaxed inequalities, which should be 2d after simplification , because every degenerate dimension double inequality of the form 1 <= xk <= 1+e^-12, became the double inequality -inf <= xk <= inf (so we have 2 relaxations per degenerate dimension).

So this test tries to create the scenario you described.

[vissarion]

I see, thanks. But should hold xk=1 not -inf <= xk <= inf after simplification

[DimitriosPavlouGR]

Sorry, I wasn't clear. After simplification we will have xk=1, and this information will be encoded in A_eq and b_eq. However the flux bounds will be changed to bu(k)=inf, bl(k)=-inf.

The idea is that once we get to the transformation step, we can easily filter the non-redundant constraints by checking if the inequalities are of the form: xk <= inf or xk >= -inf. Once we discard these bounds, we will derive a simplified polytope A x <= b. And then we can proceed with solving the system A_eq x = b_eq (where xk= 1 is encoded) and deriving a particular solution x0 which can then be used to drop the equalities and form the final h-polytope A' x <= b'.

This is what I had in mind, let me know if you have any other ideas.

[vissarion]

OK this is how polyround solves the problem, right? We can either follow this or merge simplification and transformation in one function. But at this point this PR is OK.

[DimitriosPavlouGR]

Yes it's essentially the same thing. The difference is that in their code they explicitly remove the inequalities A_i x<=b_i that are redundant at the simplification phase. That's natural for their implementation because they use the representation: A_eq x = b_eq, A x <= b.

However, due to the formulation we use i.e.: A_eq x = b_eq, bl <= x <= bu, I think it's more natural for us to remove them implicitly (i.e. by marking them as redundant) at the simplification phase by setting bl(k)/bu(k) to -inf/+inf. Then once we get to the transformation phase of the algorithm, we can remove them explicitly by not adding them to the h-polytope Ax <= b. (*)

(*) The important thing (which both us and polyround do) is that we drop the redundant constraints bl(k) <= xk <= bu(k) immediately from the lp model, so that subsequent lp calls take less time. Keeping them temporarily in the metabolic polytope as bl(k)=-inf and bu(k)=inf does not affect the running time of the method and it's more or less a design choice.

[vissarion]

Also it seems that some tests are failing. Could you please check them out?