diff --git a/README.md b/README.md index a355d23f..897d214c 100644 --- a/README.md +++ b/README.md @@ -81,6 +81,29 @@ Further examples can be found directly within this repository in the [Examples P * How to run a [stochastic simulation](https://github.com/draeger-lab/SBSCL/blob/master/src/main/java/fern/Start.java) +## Using OptSolvX (LP) in SBSCL + +### Run the OptSolvX demo + +* In your IDE, run: `org.simulator.optsolvx.OptSolvXDemo`. +* The demo builds a tiny LP and solves it via OptSolvX (CommonsMath backend for now). +* Make sure the OptSolvX **jdk8** artifact is on the classpath (declared in SBSCL’s `pom.xml`). + +### Enable debug logs + +Create the adapter with `debug = true`: + +```java +LPSolverAdapter solver = + new OptSolvXSolverAdapter(new CommonsMathSolver(), true); +``` + +### (Preview) SBML/FBC → LP + +An experimental entry point is available at `org.simulator.optsolvx.FbaToOptSolvX`. + + + ### Comparison to Similar Libraries To compare SBSCL to other simulation engines and to benchmark its predictions and results, a separate project, [SBSCL simulator comparison](https://github.com/matthiaskoenig/sbscl-simulator-comparison), is available. diff --git a/pom.xml b/pom.xml index e162e1e8..de26c1f9 100644 --- a/pom.xml +++ b/pom.xml @@ -471,6 +471,14 @@ + + + + org.optsolvx + optsolvx + 0.1.0-SNAPSHOT + jdk8 + diff --git a/src/main/java/org/simulator/optsolvx/FbaToOptSolvX.java b/src/main/java/org/simulator/optsolvx/FbaToOptSolvX.java new file mode 100644 index 00000000..300ef4bb --- /dev/null +++ b/src/main/java/org/simulator/optsolvx/FbaToOptSolvX.java @@ -0,0 +1,194 @@ +package org.simulator.optsolvx; + +import org.optsolvx.model.AbstractLPModel; +import org.optsolvx.model.Constraint; +import org.optsolvx.model.OptimizationDirection; +import org.sbml.jsbml.*; +import org.sbml.jsbml.ext.fbc.*; + +import java.text.MessageFormat; +import java.util.LinkedHashMap; +import java.util.Map; +import java.util.logging.Logger; + +/** + * SBML/FBC -> OptSolvX bridge (minimal functional version). + * Maps reactions to variables (with bounds), builds S·v=0 constraints, and sets the active objective. + * NOTE: This intentionally ignores InitialAssignments and StoichiometryMath for now (TODO). + */ +public final class FbaToOptSolvX { + private static final Logger LOGGER = Logger.getLogger(FbaToOptSolvX.class.getName()); + + private FbaToOptSolvX() { /* no instances */ } + + /** Build an OptSolvX LP model from SBML/FBC (v1 or v2). */ + public static AbstractLPModel fromSBML(SBMLDocument doc) { + if (doc == null || !doc.isSetModel()) { + throw new IllegalArgumentException("SBMLDocument must contain a Model."); + } + final Model m = doc.getModel(); + + // Prepare LP model + final AbstractLPModel lp = new AbstractLPModel(); + + // Resolve FBC namespaces for level/version + final int level = doc.getLevel(); + final int version = doc.getVersion(); + final String fbcNSv1 = FBCConstants.getNamespaceURI(level, version, 1); + final String fbcNSv2 = FBCConstants.getNamespaceURI(level, version, 2); + + // --- 1) Variables (reactions) with bounds --------------------------------------------- + // Default bounds if FBC has no values: 0 .. +inf (conventional for FBA) + final Map lb = new LinkedHashMap<>(); + final Map ub = new LinkedHashMap<>(); + + for (Reaction r : m.getListOfReactions()) { + String rid = r.getId(); + double lower = 0.0d; + double upper = Double.POSITIVE_INFINITY; + + // FBC v2: bounds via FBCReactionPlugin lower/upper "Parameter" instances + if (r.isSetPlugin(fbcNSv2)) { + FBCReactionPlugin rp = (FBCReactionPlugin) r.getPlugin(fbcNSv2); + Parameter lpi = rp.getLowerFluxBoundInstance(); + Parameter upi = rp.getUpperFluxBoundInstance(); + if (lpi != null) lower = valueOf(lpi); + if (upi != null) upper = valueOf(upi); + } + // Store preliminary bounds; FBC v1 may override below via FluxBounds list + lb.put(rid, lower); + ub.put(rid, upper); + } + + // FBC v1: FluxBounds at model level (override) + FBCModelPlugin mpV1 = (FBCModelPlugin) m.getPlugin(fbcNSv1); + if (mpV1 != null && mpV1.isSetListOfFluxBounds()) { + for (FluxBound fb : mpV1.getListOfFluxBounds()) { + String rid = fb.getReaction(); + if (rid == null) continue; + switch (fb.getOperation()) { // Java 8 compatible switch + case GREATER_EQUAL: + lb.put(rid, fb.getValue()); + break; + case LESS_EQUAL: + ub.put(rid, fb.getValue()); + break; + case EQUAL: + lb.put(rid, fb.getValue()); + ub.put(rid, fb.getValue()); + break; + default: + LOGGER.warning("FBC v1: Unsupported FluxBound operation on " + rid); + } + } + } + + // Add variables to LP + for (Reaction r : m.getListOfReactions()) { + String rid = r.getId(); + double lower = nvl(lb.get(rid), 0.0d); + double upper = nvl(ub.get(rid), Double.POSITIVE_INFINITY); + lp.addVariable(rid, lower, upper); + } + + // --- 2) Mass-balance constraints S·v = 0 (ignore boundary species) --------------------- + for (Species s : m.getListOfSpecies()) { + boolean isBoundary = s.isSetBoundaryCondition() && s.getBoundaryCondition(); + if (isBoundary) continue; + + Map coeffs = new LinkedHashMap<>(); + // For each reaction, accumulate stoichiometry of species s + for (Reaction r : m.getListOfReactions()) { + double sum = 0.0d; + // Reactants contribute negative stoichiometry + for (SpeciesReference sr : r.getListOfReactants()) { + if (s.getId().equals(sr.getSpecies())) { + sum += -stoich(sr); + } + } + // Products contribute positive stoichiometry + for (SpeciesReference sr : r.getListOfProducts()) { + if (s.getId().equals(sr.getSpecies())) { + sum += +stoich(sr); + } + } + if (sum != 0.0d) { + coeffs.put(r.getId(), sum); + } + } + if (!coeffs.isEmpty()) { + lp.addConstraint("mb_" + s.getId(), coeffs, Constraint.Relation.EQ, 0.0d); + } else { + // Optional: many species simply don't appear; keep quiet to avoid noisy logs + } + } + + // --- 3) Objective (active FBC objective) ----------------------------------------------- + Map obj = new LinkedHashMap<>(); + OptimizationDirection dir = OptimizationDirection.MAXIMIZE; // sensible default + + // Prefer FBC v2 at model level; fall back to v1 + FBCModelPlugin mpV2 = (FBCModelPlugin) m.getPlugin(fbcNSv2); + FBCModelPlugin mp = (mpV2 != null) ? mpV2 : mpV1; + + if (mp != null && mp.isSetActiveObjective()) { + Objective o = mp.getActiveObjectiveInstance(); + if (o != null) { + // Direction + Objective.Type t = o.getType(); + if (t == Objective.Type.MINIMIZE) { + dir = OptimizationDirection.MINIMIZE; + } else if (t == Objective.Type.MAXIMIZE) { + dir = OptimizationDirection.MAXIMIZE; + } + // Coefficients + for (FluxObjective fo : o.getListOfFluxObjectives()) { + String rid = fo.getReaction(); + if (rid != null) obj.put(rid, fo.getCoefficient()); + } + } else { + LOGGER.warning("No active FBC Objective instance set; using empty objective."); + } + } else { + LOGGER.warning("FBC Objective not set; using empty objective (objective=0)."); + } + + if (!obj.isEmpty()) { + lp.setObjective(obj, dir); + } else { + // Keep objective empty → objective value will be 0.0; direction still recorded if needed + lp.setObjective(new LinkedHashMap(), dir); + } + + // --- 4) Finalize & Log ----------------------------------------------------------------- + lp.build(); + LOGGER.info(MessageFormat.format( + "FbaToOptSolvX: built LP (vars={0}, cons={1}, objectiveVars={2}, dir={3})", + lp.getVariables().size(), + lp.getConstraints().size(), + obj.size(), + dir + )); + return lp; + } + + // --- Helpers ----------------------------------------------------------------------------- + + /** Get numeric value from a Parameter (defaults to 0.0 if unset). */ + private static double valueOf(Parameter p) { + if (p == null) return 0.0d; + // Prefer explicit value; ignoring units & possible initial assignments for now (TODO) + return p.isSetValue() ? p.getValue() : 0.0d; + } + + /** Null-safe value with default. */ + private static double nvl(Double v, double def) { + return (v != null) ? v.doubleValue() : def; + } + + /** Basic stoichiometry extraction (ignores StoichiometryMath & InitialAssignments for now). */ + private static double stoich(SpeciesReference sr) { + if (sr == null) return 0.0d; + return sr.isSetStoichiometry() ? sr.getStoichiometry() : 1.0d; + } +} diff --git a/src/main/java/org/simulator/optsolvx/OptSolvXDemo.java b/src/main/java/org/simulator/optsolvx/OptSolvXDemo.java new file mode 100644 index 00000000..f5ed8676 --- /dev/null +++ b/src/main/java/org/simulator/optsolvx/OptSolvXDemo.java @@ -0,0 +1,50 @@ +package org.simulator.optsolvx; + +import org.optsolvx.model.AbstractLPModel; +import org.optsolvx.model.Constraint; +import org.optsolvx.model.OptimizationDirection; +import org.optsolvx.solver.CommonsMathSolver; +import org.optsolvx.solver.LPSolution; +import org.optsolvx.solver.LPSolverAdapter; + +import java.util.HashMap; +import java.util.Map; + +public class OptSolvXDemo { + public static void main(String[] args) { + // Build a tiny LP: maximize x + y + AbstractLPModel model = new AbstractLPModel(); + model.addVariable("x", 0.0d, Double.POSITIVE_INFINITY); + model.addVariable("y", 0.0d, Double.POSITIVE_INFINITY); + + // Objective: max x + y (Java 8 style map) + Map objective = new HashMap<>(); + objective.put("x", 1.0d); + objective.put("y", 1.0d); + model.setObjective(objective, OptimizationDirection.MAXIMIZE); + + // Constraints (Java 8 style maps) + Map c1 = new HashMap<>(); + c1.put("x", 1.0d); + c1.put("y", 2.0d); + model.addConstraint("c1", c1, Constraint.Relation.LEQ, 4.0d); + + Map c2 = new HashMap<>(); + c2.put("x", 1.0d); + model.addConstraint("c2", c2, Constraint.Relation.LEQ, 3.0d); + + // Finalize model + model.build(); + + // Solve via adapter, using CommonsMath as backend + LPSolverAdapter backend = new CommonsMathSolver(); + LPSolverAdapter solver = new OptSolvXSolverAdapter(backend, /*debug=*/true); + + LPSolution sol = solver.solve(model); + + // Print result (expected optimum: x=3.0, y=0.5, objective=3.5) + System.out.println("Variables: " + sol.getVariableValues()); + System.out.println("Objective: " + sol.getObjectiveValue()); + System.out.println("Feasible: " + sol.isFeasible()); + } +} diff --git a/src/main/java/org/simulator/optsolvx/OptSolvXSolverAdapter.java b/src/main/java/org/simulator/optsolvx/OptSolvXSolverAdapter.java new file mode 100644 index 00000000..b423f51d --- /dev/null +++ b/src/main/java/org/simulator/optsolvx/OptSolvXSolverAdapter.java @@ -0,0 +1,66 @@ +package org.simulator.optsolvx; + +import org.optsolvx.model.AbstractLPModel; +import org.optsolvx.solver.LPSolution; +import org.optsolvx.solver.LPSolverAdapter; + +import java.util.Objects; +import java.text.MessageFormat; +import java.util.logging.Logger; + +public final class OptSolvXSolverAdapter implements LPSolverAdapter { + private static final Logger LOGGER = + Logger.getLogger(OptSolvXSolverAdapter.class.getName()); + + private final LPSolverAdapter backend; + private final boolean debug; + + /** + * @param backend concrete OptSolvX solver (e.g. CommonsMathSolver) + * @param debug enable verbose logging + */ + public OptSolvXSolverAdapter(LPSolverAdapter backend, boolean debug) { + this.backend = Objects.requireNonNull(backend, "backend must not be null"); + this.debug = debug; + } + + /** Convenience: debug disabled by default. */ + public OptSolvXSolverAdapter(LPSolverAdapter backend) { + this(backend, false); + } + + @Override + public LPSolution solve(AbstractLPModel model) { + if (model == null) { + throw new IllegalArgumentException("model must not be null"); + } + if (!model.isBuilt()) { + throw new IllegalStateException("Model must be built() before solving."); + } + + if (debug) { + LOGGER.info(MessageFormat.format( + "{0}: solving with {1} (vars={2}, cons={3})", + getClass().getSimpleName(), + backend.getClass().getSimpleName(), + model.getVariables().size(), + model.getConstraints().size() + )); + } + + long t0 = System.nanoTime(); + LPSolution sol = backend.solve(model); + long dtMs = (System.nanoTime() - t0) / 1_000_000L; + + if (debug) { + LOGGER.info(MessageFormat.format( + "{0}: result feasible={1}, objective={2}, time={3} ms", + getClass().getSimpleName(), + sol.isFeasible(), + sol.getObjectiveValue(), + dtMs + )); + } + return sol; + } +} diff --git a/src/test/java/org/simulator/optsolvx/BridgeObjectiveTest.java b/src/test/java/org/simulator/optsolvx/BridgeObjectiveTest.java new file mode 100644 index 00000000..a52d6194 --- /dev/null +++ b/src/test/java/org/simulator/optsolvx/BridgeObjectiveTest.java @@ -0,0 +1,79 @@ +package org.simulator.optsolvx; + +import org.junit.Test; +import static org.junit.Assert.*; + +import org.optsolvx.model.AbstractLPModel; +import org.optsolvx.solver.LPSolution; +import org.optsolvx.solver.CommonsMathSolver; +import org.sbml.jsbml.*; +import org.sbml.jsbml.ext.fbc.*; + +public class BridgeObjectiveTest { + + @Test + public void maps_fbc_v2_bounds_and_objective() { + SBMLDocument doc = new SBMLDocument(3, 2); + Model m = doc.createModel("toy_fbc2"); + Compartment c = m.createCompartment("c"); c.setConstant(true); c.setSize(1.0); + + Species x = m.createSpecies("X"); + x.setCompartment("c"); + x.setBoundaryCondition(false); + + Reaction rin = m.createReaction("R_in"); + rin.setReversible(false); + rin.createProduct().setSpecies("X"); + + Reaction rout = m.createReaction("R_out"); + rout.setReversible(false); + rout.createReactant().setSpecies("X"); + + // --- Attach FBC v2 plugin to model and reactions + final String fbcNS = FBCConstants.getNamespaceURI(doc.getLevel(), doc.getVersion(), 2); + FBCModelPlugin fbcModel = (FBCModelPlugin) m.getPlugin(fbcNS); + if (fbcModel == null) { + fbcModel = new FBCModelPlugin(m); + m.addExtension(FBCConstants.shortLabel, fbcModel); + } + + FBCReactionPlugin fbcIn = (FBCReactionPlugin) rin.getPlugin(fbcNS); + if (fbcIn == null) { + fbcIn = new FBCReactionPlugin(rin); + rin.addExtension(FBCConstants.shortLabel, fbcIn); + } + FBCReactionPlugin fbcOut = (FBCReactionPlugin) rout.getPlugin(fbcNS); + if (fbcOut == null) { + fbcOut = new FBCReactionPlugin(rout); + rout.addExtension(FBCConstants.shortLabel, fbcOut); + } + + // Bounds via Parameters referenced by FBC (v2 style) + Parameter lbIn = m.createParameter("LB_IN"); lbIn.setConstant(true); lbIn.setValue(0.0); + Parameter ubIn = m.createParameter("UB_IN"); ubIn.setConstant(true); ubIn.setValue(10.0); + Parameter lbOut= m.createParameter("LB_OUT");lbOut.setConstant(true);lbOut.setValue(0.0); + Parameter ubOut= m.createParameter("UB_OUT");ubOut.setConstant(true);ubOut.setValue(10.0); + + fbcIn.setLowerFluxBound(lbIn.getId()); + fbcIn.setUpperFluxBound(ubIn.getId()); + fbcOut.setLowerFluxBound(lbOut.getId()); + fbcOut.setUpperFluxBound(ubOut.getId()); + + // Objective: maximize v_out + Objective obj = fbcModel.createObjective("obj"); + obj.setType(Objective.Type.MAXIMIZE); + FluxObjective fo = obj.createFluxObjective(); + fo.setReaction(rout.getId()); + fo.setCoefficient(1.0); + fbcModel.setActiveObjective(obj.getId()); + + // --- Bridge & solve + AbstractLPModel lp = FbaToOptSolvX.fromSBML(doc); + LPSolution sol = new OptSolvXSolverAdapter(new CommonsMathSolver()).solve(lp); + + assertTrue(sol.isFeasible()); + assertEquals(10.0, sol.getObjectiveValue(), 1e-6); // max v_out = 10 + assertEquals(10.0, sol.getVariableValues().get("R_in"), 1e-6); + assertEquals(10.0, sol.getVariableValues().get("R_out"), 1e-6); + } +} diff --git a/src/test/java/org/simulator/optsolvx/BridgeSmokeTest.java b/src/test/java/org/simulator/optsolvx/BridgeSmokeTest.java new file mode 100644 index 00000000..970df053 --- /dev/null +++ b/src/test/java/org/simulator/optsolvx/BridgeSmokeTest.java @@ -0,0 +1,79 @@ +package org.simulator.optsolvx; + +import org.junit.Test; +import static org.junit.Assert.*; + +import org.optsolvx.model.AbstractLPModel; +import org.optsolvx.solver.LPSolution; +import org.optsolvx.solver.CommonsMathSolver; +import org.sbml.jsbml.*; +import org.sbml.jsbml.ext.fbc.*; + +public class BridgeSmokeTest { + + @Test + public void maps_fbc_v2_bounds_and_objective() { + SBMLDocument doc = new SBMLDocument(3, 2); + Model m = doc.createModel("toy_fbc2"); + Compartment c = m.createCompartment("c"); c.setConstant(true); c.setSize(1.0); + + Species x = m.createSpecies("X"); + x.setCompartment("c"); + x.setBoundaryCondition(false); + + Reaction rin = m.createReaction("R_in"); + rin.setReversible(false); + rin.createProduct().setSpecies("X"); + + Reaction rout = m.createReaction("R_out"); + rout.setReversible(false); + rout.createReactant().setSpecies("X"); + + // --- Attach FBC v2 plugin to model and reactions + final String fbcNS = FBCConstants.getNamespaceURI(doc.getLevel(), doc.getVersion(), 2); + FBCModelPlugin fbcModel = (FBCModelPlugin) m.getPlugin(fbcNS); + if (fbcModel == null) { + fbcModel = new FBCModelPlugin(m); + m.addExtension(FBCConstants.shortLabel, fbcModel); + } + + FBCReactionPlugin fbcIn = (FBCReactionPlugin) rin.getPlugin(fbcNS); + if (fbcIn == null) { + fbcIn = new FBCReactionPlugin(rin); + rin.addExtension(FBCConstants.shortLabel, fbcIn); + } + FBCReactionPlugin fbcOut = (FBCReactionPlugin) rout.getPlugin(fbcNS); + if (fbcOut == null) { + fbcOut = new FBCReactionPlugin(rout); + rout.addExtension(FBCConstants.shortLabel, fbcOut); + } + + // Bounds via Parameters referenced by FBC (v2 style) + Parameter lbIn = m.createParameter("LB_IN"); lbIn.setConstant(true); lbIn.setValue(0.0); + Parameter ubIn = m.createParameter("UB_IN"); ubIn.setConstant(true); ubIn.setValue(10.0); + Parameter lbOut= m.createParameter("LB_OUT");lbOut.setConstant(true);lbOut.setValue(0.0); + Parameter ubOut= m.createParameter("UB_OUT");ubOut.setConstant(true);ubOut.setValue(10.0); + + fbcIn.setLowerFluxBound(lbIn.getId()); + fbcIn.setUpperFluxBound(ubIn.getId()); + fbcOut.setLowerFluxBound(lbOut.getId()); + fbcOut.setUpperFluxBound(ubOut.getId()); + + // Objective: maximize v_out + Objective obj = fbcModel.createObjective("obj"); + obj.setType(Objective.Type.MAXIMIZE); + FluxObjective fo = obj.createFluxObjective(); + fo.setReaction(rout.getId()); + fo.setCoefficient(1.0); + fbcModel.setActiveObjective(obj.getId()); + + // --- Bridge & solve + AbstractLPModel lp = FbaToOptSolvX.fromSBML(doc); + LPSolution sol = new OptSolvXSolverAdapter(new CommonsMathSolver()).solve(lp); + + assertTrue(sol.isFeasible()); + assertEquals(10.0, sol.getObjectiveValue(), 1e-6); // max v_out = 10 + assertEquals(10.0, sol.getVariableValues().get("R_in"), 1e-6); + assertEquals(10.0, sol.getVariableValues().get("R_out"), 1e-6); + } +} diff --git a/src/test/java/org/simulator/optsolvx/OptSolvXSolverAdapterTest.java b/src/test/java/org/simulator/optsolvx/OptSolvXSolverAdapterTest.java new file mode 100644 index 00000000..c36b9716 --- /dev/null +++ b/src/test/java/org/simulator/optsolvx/OptSolvXSolverAdapterTest.java @@ -0,0 +1,55 @@ +package org.simulator.optsolvx; + +import org.junit.Test; +import static org.junit.Assert.*; +import org.optsolvx.model.*; +import org.optsolvx.solver.*; + +import java.util.HashMap; +import java.util.Map; + +public class OptSolvXSolverAdapterTest { + + @Test(expected = IllegalArgumentException.class) + public void solve_requires_non_null_model() { + LPSolverAdapter s = new OptSolvXSolverAdapter(new CommonsMathSolver(), false); + s.solve(null); // must throw + } + + @Test(expected = IllegalStateException.class) + public void solve_requires_build() { + AbstractLPModel m = new AbstractLPModel(); + m.addVariable("x", 0.0d, 10.0d); + LPSolverAdapter s = new OptSolvXSolverAdapter(new CommonsMathSolver(), false); + s.solve(m); // not built -> must throw + } + + @Test + public void smoke_minimize_with_eq_and_bounds() { + // Minimize 2x + y, s.t. x + y = 5, 0 <= x,y <= 5 + AbstractLPModel m = new AbstractLPModel(); + m.addVariable("x", 0.0d, 5.0d); + m.addVariable("y", 0.0d, 5.0d); + + Map obj = new HashMap<>(); + obj.put("x", 2.0d); + obj.put("y", 1.0d); + m.setObjective(obj, OptimizationDirection.MINIMIZE); + + Map eq = new HashMap<>(); + eq.put("x", 1.0d); + eq.put("y", 1.0d); + m.addConstraint("sum", eq, Constraint.Relation.EQ, 5.0d); + + m.build(); + + LPSolverAdapter s = new OptSolvXSolverAdapter(new CommonsMathSolver(), false); + LPSolution sol = s.solve(m); + + assertTrue(sol.isFeasible()); + // Optimum at x=0, y=5 -> objective = 5 + assertEquals(5.0d, sol.getObjectiveValue(), 1e-6d); + assertEquals(0.0d, sol.getVariableValues().get("x"), 1e-6d); + assertEquals(5.0d, sol.getVariableValues().get("y"), 1e-6d); + } +}