Merge pull request #60 from COBREXA/mk-highs

try switching from glpk to highs, see what happens

Project.toml

@@ -24,7 +24,7 @@ ConstraintTrees = "1.1"
 DistributedData = "0.2"
 DocStringExtensions = "0.8, 0.9"
 Downloads = "1"
-GLPK = "1"
+HiGHS = "1.9"
 JSONFBCModels = "0.1"
 JuMP = "1"
 SBMLFBCModels = "0.1"
@@ -38,12 +38,12 @@ julia = "1.5"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
-GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8"
 SBMLFBCModels = "3e8f9d1a-ffc1-486d-82d6-6c7276635980"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
 
 [targets]
-test = ["Aqua", "Clarabel", "Downloads", "GLPK", "JSONFBCModels", "SBMLFBCModels", "SHA", "Test", "Tulip"]
+test = ["Aqua", "Clarabel", "Downloads", "HiGHS", "JSONFBCModels", "SBMLFBCModels", "SHA", "Test", "Tulip"]

README.md

@@ -74,7 +74,7 @@ COBREXA.jl is best installed using Julia's packaging system.
    dependencies.
 4. You may want to install [additional constraint-solver software compatible
    with JuMP](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers).
-   Type e.g. `] add GLPK` to install the popular GLPK solver.
+   Type e.g. `] add HiGHS` to install the HiGHS solver.
 
 ### Guides
 

docs/Project.toml

@@ -4,7 +4,7 @@ COBREXA = "babc4406-5200-4a30-9033-bf5ae714c842"
 Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 ConstraintTrees = "5515826b-29c3-47a5-8849-8513ac836620"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8"

docs/src/distributed/2_parallel.md

@@ -20,9 +20,9 @@ case, it should return a vector of _worker IDs_, very likely equal to
 
 Each of the processes contains a self-sufficient image of Julia that can act
 independently; in turn the additional processes also consume some memory. Each
-process with loaded `COBREXA.jl` and a simple solver such as GLPK may consume
-around 500MB of RAM, which should be taken into account when planning the
-analysis scale.
+process with loaded `COBREXA.jl` and a solver such as HiGHS may consume around
+300MB of RAM, which should be taken into account when planning the analysis
+scale.
 
 !!! warning "Using Julia environments with Distributed"
     In certain conditions, the Distributed package does not properly forward
@@ -36,7 +36,7 @@ analysis scale.
 Packages (COBREXA and the selected solver) must be loaded at all processes,
 which may ensured using the "everywhere" macro (from `Distributed` package):
 ```julia
-@everywhere using COBREXA, GLPK
+@everywhere using COBREXA, HiGHS
 ```
 
 Utilizing the prepared worker processes is then straightforward: We pass the
@@ -47,7 +47,7 @@ argument, and the parallel processing is orchestrated automatically:
 model = load_model("e_coli_core.xml")
 result = flux_variability_analysis(
     model,
-    optimizer = GLPK.Optimizer,
+    optimizer = HiGHS.Optimizer,
     workers = workers()
 )
 ```

docs/src/distributed/3_slurm.md

@@ -37,7 +37,7 @@ The Julia script that does a parallel analysis in a Slurm cluster may look as
 follows:
 
 ```julia
-using COBREXA, Distributed, ClusterManagers, GLPK
+using COBREXA, Distributed, ClusterManagers, HiGHS
 
 available_workers = parse(Int, ENV["SLURM_NTASKS"])
 

docs/src/examples/02a-flux-balance-analysis.jl

@@ -29,10 +29,10 @@ download_model(
 )
 
 # Additionally to COBREXA and the model format package, we will need a solver
-# -- let's use GLPK here:
+# -- let's use HiGHS here:
 
 import JSONFBCModels
-import GLPK
+import HiGHS
 
 model = load_model("e_coli_core.json")
 
@@ -43,7 +43,7 @@ model = load_model("e_coli_core.json")
 # is captured in the default behavior of function
 # [`flux_balance_analysis`](@ref):
 
-solution = flux_balance_analysis(model, optimizer = GLPK.Optimizer)
+solution = flux_balance_analysis(model, optimizer = HiGHS.Optimizer)
 
 @test isapprox(solution.objective, 0.8739, atol = TEST_TOLERANCE) #src
 

docs/src/examples/02b-optimizer-parameters.jl

@@ -51,15 +51,14 @@ model = load_model("e_coli_core.json")
 solution = flux_balance_analysis(
     model,
     optimizer = Tulip.Optimizer,
-    settings = [silence, set_optimizer_attribute("IPM_IterationsLimit", 1000)],
+    settings = [set_optimizer_attribute("IPM_IterationsLimit", 1000)],
 )
 
 @test !isnothing(solution) #src
 
 # To see some of the effects of the configuration changes, we may e.g.
 # deliberately cripple the optimizer's possibilities to a few iterations and
-# only a little time, which will cause it to fail, return no solution, and
-# verbosely describe what happened:
+# only a little time, which will cause it to fail and return no solution:
 
 solution = flux_balance_analysis(
     model,
@@ -71,7 +70,24 @@ println(solution)
 
 @test isnothing(solution) #src
 
+# To see what failed, users may examine the solver output. Because all solver
+# output is silenced by default for efficiency reasons, we need to explicitly
+# pass in the [`unsilence`](@ref) setting:
+
+solution = flux_balance_analysis(
+    model,
+    optimizer = Tulip.Optimizer,
+    settings = [
+        set_optimizer_attribute("IPM_IterationsLimit", 2),
+        set_time_limit(0.1),
+        unsilence,
+    ],
+)
+
 # Applicable optimizer attributes are documented in the documentations of the
 # respective optimizers. To browse the possibilities, one might want to see the
 # [JuMP documentation page that summarizes the references to the available
 # optimizers](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers).
+#
+# Default solver settings can be examined and changed via
+# [`Configuration`](@ref).

docs/src/examples/02c-model-modifications.jl

@@ -81,9 +81,9 @@ model.reactions["CS"].stoichiometry
 # Let's first find a "original" solution, so that we have a base solution for
 # comparing:
 
-import GLPK
+import HiGHS
 
-base_solution = flux_balance_analysis(model, optimizer = GLPK.Optimizer)
+base_solution = flux_balance_analysis(model, optimizer = HiGHS.Optimizer)
 base_solution.objective
 
 # Now, for example, we can limit the intake of glucose by the model:
@@ -96,7 +96,7 @@ model.reactions["EX_glc__D_e"].lower_bound = -5.0
 
 # ...and solve the modified model:
 #
-low_glucose_solution = flux_balance_analysis(model, optimizer = GLPK.Optimizer)
+low_glucose_solution = flux_balance_analysis(model, optimizer = HiGHS.Optimizer)
 low_glucose_solution.objective
 
 @test isapprox(low_glucose_solution.objective, 0.41559777, atol = TEST_TOLERANCE) #src
@@ -185,7 +185,7 @@ model.couplings["total_energy_intake"] = CM.Coupling(
 # The values of any coupling constraints can be inspected directly in the
 # solved model:
 
-solution_with_coupling = flux_balance_analysis(model, optimizer = GLPK.Optimizer)
+solution_with_coupling = flux_balance_analysis(model, optimizer = HiGHS.Optimizer)
 
 solution_with_coupling.coupling.total_energy_intake
 

docs/src/examples/02d-constraint-modifications.jl

@@ -32,7 +32,7 @@ download_model(
 )
 
 import JSONFBCModels
-import GLPK
+import HiGHS
 
 model = load_model("e_coli_core.json") # flux balance type model
 
@@ -99,7 +99,7 @@ fermenting_ct = ct * :fermentation^fermentation_constraint
 solution = optimized_values(
     fermenting_ct,
     objective = fermenting_ct.objective.value,
-    optimizer = GLPK.Optimizer,
+    optimizer = HiGHS.Optimizer,
 )
 
 @test isapprox(solution.objective, 0.633738, atol = TEST_TOLERANCE) #src
@@ -110,7 +110,7 @@ solution = optimized_values(
 
 ct.fluxes.ATPM.bound = C.Between(1000.0, 10000.0)
 
-solution = optimized_values(ct, objective = ct.objective.value, optimizer = GLPK.Optimizer)
+solution = optimized_values(ct, objective = ct.objective.value, optimizer = HiGHS.Optimizer)
 
 print(solution)
 

docs/src/examples/03a-flux-variability-analysis.jl

@@ -30,13 +30,13 @@ download_model(
     "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
 )
 
-import JSONFBCModels, GLPK
+import JSONFBCModels, HiGHS
 
 model = load_model("e_coli_core.json")
 
 # The "usual" form of FBA is available via the eponymous function:
 
-solution = flux_variability_analysis(model, optimizer = GLPK.Optimizer)
+solution = flux_variability_analysis(model, optimizer = HiGHS.Optimizer)
 
 @test isapprox(solution.ACALD[1], -2.542370370370188, atol = TEST_TOLERANCE) #src
 @test isapprox(solution.ACALD[2], 0.0, atol = TEST_TOLERANCE) #src
@@ -52,15 +52,15 @@ solution = flux_variability_analysis(model, optimizer = GLPK.Optimizer)
 
 very_close = flux_variability_analysis(
     model,
-    optimizer = GLPK.Optimizer,
+    optimizer = HiGHS.Optimizer,
     objective_bound = absolute_tolerance_bound(1e-5),
 )
 
 # Here, we relax that to 1% of the optimal objective value:
 
 one_percent_close = flux_variability_analysis(
     model,
-    optimizer = GLPK.Optimizer,
+    optimizer = HiGHS.Optimizer,
     objective_bound = relative_tolerance_bound(0.99),
 )
 

docs/src/examples/03b-parsimonious-flux-balance.jl

@@ -40,11 +40,7 @@ model = load_model("e_coli_core.json")
 # The pFBA is, in its most default form, implemented in function
 # [`parsimonious_flux_balance_analysis`](@ref):
 
-solution = parsimonious_flux_balance_analysis(
-    model;
-    optimizer = Clarabel.Optimizer,
-    settings = [silence],
-)
+solution = parsimonious_flux_balance_analysis(model; optimizer = Clarabel.Optimizer)
 
 #
 
@@ -71,13 +67,12 @@ solution.fluxes
 # fluxes. We can set the optimizer and parsimonious optimizer separately using
 # keyword arguments:
 
-import GLPK
+import HiGHS
 
 solution = parsimonious_flux_balance_analysis(
     model;
-    optimizer = GLPK.Optimizer, # GLPK is good for LP but cannot do QP
-    settings = [silence],
-    parsimonious_optimizer = Clarabel.Optimizer, # Clarabel is not very precise but can solve QP
+    optimizer = HiGHS.Optimizer, # HiGHS is used only for LP here
+    parsimonious_optimizer = Clarabel.Optimizer, # Clarabel is great for solving QPs
 )
 
 @test isapprox(solution.objective, 0.873922; atol = TEST_TOLERANCE) #src
@@ -91,7 +86,7 @@ solution = parsimonious_flux_balance_analysis(
 # at all:
 
 linear_solution =
-    linear_parsimonious_flux_balance_analysis(model; optimizer = GLPK.Optimizer)
+    linear_parsimonious_flux_balance_analysis(model; optimizer = HiGHS.Optimizer)
 
 #
 

docs/src/examples/03c-envelopes.jl

@@ -32,7 +32,7 @@ download_model(
 )
 
 import JSONFBCModels
-import GLPK
+import HiGHS
 
 model = load_model("e_coli_core.json")
 
@@ -46,7 +46,7 @@ envelope = objective_production_envelope(
     model,
     ["EX_o2_e", "EX_co2_e"];
     breaks = 5,
-    optimizer = GLPK.Optimizer,
+    optimizer = HiGHS.Optimizer,
 )
 
 @test count(isnothing, envelope.objective_values) == 18 #src

docs/src/examples/04-community-models.jl

@@ -31,7 +31,7 @@ download_model(
 # and a few supporting packages.
 
 import JSONFBCModels
-import GLPK
+import HiGHS
 import AbstractFBCModels.CanonicalModel as CM
 import ConstraintTrees as C
 
@@ -62,7 +62,7 @@ my_community = Dict("bug1" => (ecoli1, 0.2), "bug2" => (ecoli2, 0.8))
 solution = community_flux_balance_analysis(
     my_community,
     ["EX_glc__D_e" => (-10.0, 0.0)],
-    optimizer = GLPK.Optimizer,
+    optimizer = HiGHS.Optimizer,
 )
 
 @test isapprox(solution.community_biomass, 0.5237157737585179, atol = TEST_TOLERANCE) #src
@@ -92,7 +92,7 @@ screen(0.0:0.1:1.0) do ratio2
         [("bug1" => (ecoli1, ratio1)), ("bug2" => (ecoli2, ratio2))],
         ["EX_glc__D_e" => (-10.0, 0.0)],
         interface = :sbo, # usually more reproducible
-        optimizer = GLPK.Optimizer,
+        optimizer = HiGHS.Optimizer,
     )
     (ratio1, ratio2) => (isnothing(res) ? nothing : res.community_biomass)
 end
@@ -181,7 +181,7 @@ custom_solution = optimized_values(
     custom_community,
     objective = custom_community.community_biomass.value,
     output = custom_community.community_biomass,
-    optimizer = GLPK.Optimizer,
+    optimizer = HiGHS.Optimizer,
 )
 
 @test isapprox(custom_solution, solution.community_biomass, atol = TEST_TOLERANCE) #src

docs/src/examples/05a-minimization-of-metabolic-adjustment.jl

@@ -36,7 +36,7 @@ download_model(
 # reference solution" typically refers to Euclidean ("L2") distance which
 # requires a QP solver, but Manhattan ("L1") distance is also demonstrated
 # below.
-import Clarabel, GLPK
+import Clarabel, HiGHS
 
 # Because we will have to perform some perturbation, we import the model in
 # canonical Julia structures:
@@ -60,11 +60,7 @@ limited_ecoli.reactions["CYTBD"].upper_bound = 10.0
 # just right. For later comparison, we first get the optimal parsimonious flux
 # distribution in the perturbed model:
 
-solution = parsimonious_flux_balance_analysis(
-    limited_ecoli,
-    optimizer = Clarabel.Optimizer,
-    settings = [silence],
-)
+solution = parsimonious_flux_balance_analysis(limited_ecoli, optimizer = Clarabel.Optimizer)
 
 # Now, how much is the flux going to differ if we assume the bacterium did only
 # minimal adjustment from the previous state with unlimited CYTBD?
@@ -73,7 +69,6 @@ moma_solution = metabolic_adjustment_minimization_analysis(
     limited_ecoli, # the model to be examined
     ecoli; # the model that gives the reference flux
     optimizer = Clarabel.Optimizer,
-    settings = [silence],
 )
 
 @test isapprox(moma_solution.objective, 0.241496699165187, atol = TEST_TOLERANCE) #src
@@ -104,17 +99,12 @@ sort(collect(difference.fluxes), by = last)
 # from a known reaction flux. We can supply it manually as the second argument
 # (instead of the reference model).
 
-ref = parsimonious_flux_balance_analysis(
-    ecoli,
-    optimizer = Clarabel.Optimizer,
-    settings = [silence],
-)
+ref = parsimonious_flux_balance_analysis(ecoli, optimizer = Clarabel.Optimizer)
 
 ref_closest_solution = metabolic_adjustment_minimization_analysis(
     limited_ecoli,
     ref.fluxes;
     optimizer = Clarabel.Optimizer,
-    settings = [silence],
 )
 
 @test isapprox( #src
@@ -137,7 +127,6 @@ solution_close_to_measurement = metabolic_adjustment_minimization_analysis(
     limited_ecoli,
     measured_fluxes;
     optimizer = Clarabel.Optimizer,
-    settings = [silence],
 )
 
 @test isapprox(solution_close_to_measurement.objective, 0.272247, atol = TEST_TOLERANCE) #src
@@ -151,8 +140,7 @@ solution_close_to_measurement = metabolic_adjustment_minimization_analysis(
 linear_moma_solution = linear_metabolic_adjustment_minimization_analysis(
     limited_ecoli,
     ecoli;
-    optimizer = GLPK.Optimizer,
-    settings = [silence],
+    optimizer = HiGHS.Optimizer,
 )
 
 sort(collect(linear_moma_solution.fluxes), by = last)