diff --git a/.github/workflows/CI-test.yml b/.github/workflows/CI-test.yml index f0480a4..6d6c010 100644 --- a/.github/workflows/CI-test.yml +++ b/.github/workflows/CI-test.yml @@ -69,14 +69,6 @@ jobs: echo 'export PATH="${{ env.pythonLocation }}/bin:$PATH"' >> ~/.bash_profile - name: Check Python Version run: python --version - - name: Java Setup - uses: actions/setup-java@v5 - with: - distribution: 'temurin' # See 'Supported distributions' for available options - java-version: '21' - - name: Set JAVA_HOME for mac OS - if: runner.os == 'macOS' - run: echo "JAVA_HOME=$(/usr/libexec/java_home)" >> "$GITHUB_ENV" - name: Build for PyPI (if pip) if: matrix.package-manager == 'pip' run: | @@ -97,10 +89,6 @@ jobs: if: matrix.package-manager == 'pip' run: | pip install pyscipopt - - name: Install jpype1 (pip) - if: matrix.package-manager == 'pip' - run: | - pip install jpype1 - name: Install pytest (pip) if: matrix.package-manager == 'pip' run: | @@ -139,10 +127,6 @@ jobs: if: matrix.package-manager == 'conda' run: | micromamba install -c conda-forge pyscipopt scip - - name: Install jpype1 (conda) - if: matrix.package-manager == 'conda' - run: | - pip install jpype1 - name: Install pytest (conda) if: matrix.package-manager == 'conda' run: | diff --git a/README.md b/README.md index a2a0d1e..b8ba489 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,7 @@ To get started, check out the [StrainDesign documentation](https://straindesign. Flux-space animation

-Parts of the compression routine are done by efmtool's compression function ([csb.ethz.ch/tools/software/efmtool.html](https://csb.ethz.ch/tools/software/efmtool.html)[6]). Therefore some source code from the [efmtool_link](https://github.com/cnapy-org/efmtool_link) package was adopted. +The compression routine follows the approach of efmtool's compression function ([csb.ethz.ch/tools/software/efmtool.html](https://csb.ethz.ch/tools/software/efmtool.html)[6]), reimplemented in pure Python with exact rational arithmetic. ## Installation @@ -51,10 +51,6 @@ pip install -e . in the main folder. Through the installation with `-e`, updates from a `git pull` are at once available in your Python environment without the need for a reinstallation. -### Legacy Java backend (optional) - -Java is not required for the default compression (`compression_backend='sparse_rref'`). A legacy Java-based EFMTool backend (`compression_backend='efmtool_rref'`) is optionally available via `pip install straindesign[java]`. For setup help see the [Legacy Methods](https://straindesign.readthedocs.io/en/latest/legacy_methods.html) documentation page. - ## Install additional solvers The cobra package is shipped with the GLPK solver. The more powerful commercial solvers IBM CPLEX and Gurobi may be used by cobra and the straindesign package. This makes sense in particular when using strain design algorithms like MCS, OptKnock etc. As another alternative solver, SCIP may be used. In the following, you will find installation instructions for the individual solvers. diff --git a/docs/source/developers_guide.md b/docs/source/developers_guide.md index 662622d..33de71a 100644 --- a/docs/source/developers_guide.md +++ b/docs/source/developers_guide.md @@ -29,7 +29,7 @@ which drifts with every edit. Grep for the symbol. 1. [**Orientation & the strain-design problem**](#ch1) — the MCS problem, SUPPRESS/PROTECT/bilevel semantics, interventions & cost, the binary `z` vector, invocation, and the master notation table. 2. [**The constraint-based foundation**](#ch2) — `Sv=0`, the flux polytope/cone, FBA & FVA as LPs, the internal standard form, and the convex geometry needed for duality. -3. [**Network compression**](#ch3) — exact rational nullspace compression; parallel, coupled, conservation and blocked reductions; lump scaling; GPR propagation and simplification; compression maps; and the legacy efmtool backend. +3. [**Network compression**](#ch3) — exact rational nullspace compression; parallel, coupled, conservation and blocked reductions; lump scaling; GPR propagation and simplification; compression maps;. 4. [**GPR integration**](#ch4) — GPR reduction and Boolean simplification, `extend_model_gpr`, reversible splitting, module remapping, and the two-compression-pass boundary. 5. [**FVA in preprocessing**](#ch5) — pre-compression sign classification, desired-region essentiality, final bound/module FVA, the single-classical-module fold, and size-1 MCS extraction. 6. [**Dualization (the mathematical core)**](#ch6) — LP duality, Farkas certificates, and the strong-duality encodings shared by the supported module types. @@ -67,9 +67,7 @@ straindesign/ ├── gurobi_interface.py # Gurobi backend (Gurobi_MILP_LP) ├── scip_interface.py # SCIP backend (SCIP_MILP_LP) ├── glpk_interface.py # GLPK backend (GLPK_MILP_LP) -├── efmtool_cmp_interface.py # EFMtool JAR interface (legacy compression backend) ├── pool.py # SDPool: cross-platform multiprocessing pool -└── efmtool.jar # Bundled EFMtool binary ``` Which chapter covers which module: compression → [Ch 3](#ch3); GPR / networktools → [Ch 4](#ch4), [Ch 12](#ch12); FVA / lptools / @@ -1082,7 +1080,7 @@ nonzero flux in any steady state — a *contradicting* group. Then the master *a and `contradicting_removed` is set, which is the flag that triggers a re-iteration of the whole pass (→): removing a contradicting group changes the flux space and may make previously-uncoupled reactions coupled. A consistent (nonempty) group removes only the slaves -. This bound-intersection logic replaced a Java-era behaviour that could drop +. This bound-intersection logic replaced an earlier behaviour that could drop reactions incorrectly; getting the translate-and-intersect direction right (especially the `λ<0` flip and the `±inf` handling) is exactly the subject of the closed issue #44 cautionary tale in [Ch 10](#ch10). @@ -1107,9 +1105,7 @@ Two design points. First, this is a **row-rank reduction**, complementary to the §3.4 — together they push `S` toward full rank (the §3.1 hypothesis). Second, the *ordering* matters: conservation removal runs *before* the expensive coupled step in each cycle (`compress_model`, –). Fewer metabolite rows means the nullspace RREF that drives coupling detection operates -on a smaller matrix, so removing dependent rows first makes the costliest stage cheaper. (There is a -legacy Java oracle, `_remove_conservation_relations_java` at, selectable via the -`efmtool_rref` backend; the default `sparse_rref` path uses the pure-Python exact RREF above.) +on a smaller matrix, so removing dependent rows first makes the costliest stage cheaper. ### 3.6 Blocked and zero-flux removal @@ -1266,202 +1262,6 @@ re-injection, and gene translation — are owned by [Ch 9](#ch9); this section o map that [Ch 9](#ch9) consumes. -### 3.11 The legacy efmtool (Java) backend - -Everything in §3.2–§3.10 describes the **default** compression engine: the pure-Python, exact -integer/rational `sparse_rref` backend. That engine is a *reimplementation*. The original backend — -and the one every pre-1.15 release actually ran — was **efmtool**, Marco Terzer's Java tool for -elementary-flux-mode enumeration and network compression (the compression stage of efmtool is exactly -the coupled/zero/contradicting reduction that §3.4/§3.6 now do in Python). It is still shipped and -still reachable, selected with `compression_backend='efmtool_rref'`, and this section documents how the -bridge works and *why* it has been demoted to legacy. Reading it also explains the vocabulary the -Python code inherited: the Python `CompressionMethod` enum (`compression.py`), the Python class -name `StoichMatrixCompressor` (`compression.py`), and the `CoupledZero`/`CoupledCombine`/ -`CoupledContradicting` method names are all deliberate echoes of the efmtool Java API they replaced. - -#### 3.11.1 What efmtool is and how straindesign reaches it - -efmtool is a Java library (namespace `ch.javasoft.*`, packaged as `efmtool.jar` alongside the Python -sources at `straindesign/efmtool.jar`). straindesign uses only its *compression* half — not its EFM -enumeration — through the classes loaded in `efmtool_cmp_interface.py`–: -`ch.javasoft.smx.impl.DefaultBigIntegerRationalMatrix` (an arbitrary-precision rational matrix), -`ch.javasoft.smx.ops.Gauss` (rational Gaussian elimination), `ch.javasoft.metabolic.compress. -StoichMatrixCompressor` and `CompressionMethod`, and `ch.javasoft.math.BigFraction` / -`java.math.BigInteger`. The bridge is **JPype**: `_start_jvm` (`efmtool_cmp_interface.py`) starts an -in-process JVM, adds `efmtool.jar` to the classpath, and imports the Java classes via -`jpype.imports` so they become callable Python objects. - -The routing has three layers. - -1. **Import time.** `__init__.py`– calls `_start_jvm` *eagerly* at `import straindesign`. - This is a no-op when jpype1 or a JVM is absent (neither is a package dependency), so a normal install - never touches Java. When Java *is* present the JVM must be started here — before NumPy/OpenBLAS spins - up worker threads — or JNI calls later crash with SIGBUS/SIGSEGV (`__init__.py`–; the code is - littered with such mitigations, see §3.11.4). -2. **Backend selection.** `compute_strain_designs` reads the kwarg - `compression_backend = kwargs.get('compression_backend', 'sparse_rref')` - (`compute_strain_designs.py`) and threads it into both `compress_model` calls - . `compress_model` sets `use_java = (compression_backend == 'efmtool_rref')` - (`compression.py`). -3. **Dispatch inside the fixpoint.** Crucially, `efmtool_rref` does **not** replace the whole - compression pipeline — only two of its three reducers. Inside the alternating fixpoint (§3.7, - `compression.py`–): - - **Parallel merge** (step 1, §3.8) is **always** the Python hash-based `compress_model_parallel` — - efmtool has no equivalent and it is never routed to Java. - - **Conservation removal** (step 2, §3.5) forks on `use_java` : Java goes through - `_remove_conservation_relations_java`, Python through `remove_conservation_relations`. - - **Coupled merge** (step 3, §3.4) forks inside `compress_model_coupled`: Java calls - `compress_model_java` (`efmtool_cmp_interface.py`), Python calls `compress_cobra_model`. - - So `efmtool_rref` is really a **hybrid**: Python parallel-merge + Java conservation-removal + Java - coupled-merge, iterated by the same Python fixpoint driver. The two backends differ only in the - *nullspace/rank algorithm* used for steps 2 and 3. - -#### 3.11.2 Data marshalling: cobra model → Java → cobra model - -The coupled step, `compress_model_java` (`efmtool_cmp_interface.py`), is where the interesting -marshalling lives. It mutates the cobra model in place and returns the same -`{compressed_id: {orig_id: factor}}` reaction map that the Python backend produces, so the rest of the -pipeline (module remapping, cost compression, decompression in [Ch 9](#ch9)) is backend-agnostic. - -**Into Java.** -- `stoichmat_coeff_to_fraction(model)` first converts every stoichiometric coefficient to an - exact `Fraction`/sympy-`Rational` — the same exactness discipline as §3.2.1, done *before* any Java - call. -- All gene rules are cleared, `r.gene_reaction_rule = ''`, matching the Python coupled path - (§3.9); GPR is re-attached afterward (below). -- A `DefaultBigIntegerRationalMatrix(num_met, num_active)` is allocated and filled column by - column. Reactions whose upper bound is `≤ 0` are **flipped** to the forward direction - (`model.reactions[mi] *= -1`,–) and their index recorded in `flipped`; efmtool's - compressor assumes a canonical orientation. Each coefficient `v` is converted by - `sympyRat2jBigIntegerPair` into a Java `BigInteger` numerator/denominator pair — using - `BigInteger.valueOf` for values that fit in 63 bits and `BigInteger(str(...))` otherwise — and set as - a `BigFraction(n, d)`. This path is **exact**: efmtool's `DefaultBigIntegerRational - Matrix` is arbitrary-precision, so the Java core does *not* overflow. -- A `StoichMatrixCompressor(subset_compression)` is built, where `subset_compression = - [CoupledZero, CoupledCombine, CoupledContradicting]` : remove structurally - zero-flux reactions, combine coupled groups, and drop contradicting groups — the Java analogues of - §3.3's three removal kinds. `smc.compress(stoich_mat, reversible, …, reacNames, None)` - returns a `comprec` whose `post` matrix is the reaction transformation (the Java counterpart of the - Python `post` in §3.3, `v_original = post · v_compressed`). - -**Back to Python.** Here is the seam that matters for correctness: - -```python -subset_matrix = jpypeArrayOfArrays2numpy_mat(comprec.post.getDoubleRows()) # :424 — DOUBLES -``` - -The *structure* of the compression (which original reaction maps into which compressed column, and the -zero pattern) is read back as a **double-precision** numpy matrix via `getDoubleRows`. The -per-reaction merge then: -- flags a reaction zero-flux iff its `subset_matrix` row is all-zero; -- for each compressed column `j`, gathers members from `subset_matrix[:,j].nonzero`, scales - each member's stoichiometry by the **exact** factor `jBigFraction2sympyRat(comprec.post. - getBigFractionValueAt(ai, j))` (–, exact `BigFraction → sympy.Rational`), and **rescales - its bounds by `/= abs(subset_matrix[ai, j])`** (–, i.e. by a **double**); -- merges member reactions into the group representative, concatenating ids with `*` and truncating past - ~220 chars to `...` — the same naming convention as the parallel backend (§3.8); -- records `subset_rxns`/`subset_stoich` per representative (negating the stoich for `flipped` - reactions,–) and finally assembles `rational_map` from them. - -So the *factors* are exact rationals, but the *pattern detection and the bound rescaling* pass through -double precision. The `suppressed_reactions` argument — reaction ids that must survive -because a strain-design module references them — are excluded from the active set entirely and re-added -as standalone identity entries, a workaround for efmtool's `CoupledContradicting` step, -which will otherwise delete reactions it deems inconsistent (contrast the Python backend, which keeps -them via the exact bounds-intersection of §3.4.4). Back in `compress_model_coupled` the Java branch -then sweeps up any leftover `(0,0)` reactions (`compression.py`–) and — identically to the -Python branch — re-attaches the **AND-combined GPR** from the pre-merge snapshot -(`compression.py`–). GPR propagation is therefore the *same* for both backends on the -coupled step. - -**The conservation path.** `_remove_conservation_relations_java` (`compression.py`) builds `S` as -a LIL matrix, **densifies its transpose** (`stoich_mat.transpose.toarray`), and hands it -to `basic_columns_rat_java` (`efmtool_cmp_interface.py`). That function wraps the dense array into a -`DefaultBigIntegerRationalMatrix` via `numpy_mat2jpypeArrayOfArrays` — which builds a **`JDouble[rows, -cols]`** — then runs `Gauss.getRationalInstance.rowEchelon(...)` and returns the -pivot columns, i.e. the independent metabolite rows; the non-pivot metabolites are dependent -(conservation relations) and removed (`compression.py`–). This is the exact-RREF -independence oracle of §3.5, but computed in Java — and note it marshals the stoichiometry through a -**dense double** array, both memory-heavy on genome-scale models and lossy for large coefficients. - -#### 3.11.3 Why it is legacy - -The pure-Python `sparse_rref` engine (§3.2) was written to replace efmtool for four concrete reasons, -each a decisive advantage on a genome-scale correctness/performance workload: - -1. **No JVM / JPype dependency.** efmtool needs a JVM, the `efmtool.jar`, `jpype1`, and `sympy` all - present and version-compatible (`_init_java`, `efmtool_cmp_interface.py`, raises `ImportError` - for any missing piece). The Python backend needs only NumPy/SciPy, which straindesign already - depends on. A default that requires a working Java toolchain is a default that fails on many - installs. -2. **Native-crash fragility.** The bridge is defensive to a degree that itself signals the risk: - eager JVM startup ordered before OpenBLAS threads (§3.11.1); `gc.disable` wrapped around *every* - JNI block (`efmtool_cmp_interface.py`–,–) because Python's garbage collector - finalizing a JPype proxy mid-call causes Bus error / SIGSEGV; an `atexit` JVM-shutdown hook to dodge - a JPype teardown race. None of this can occur in a pure-Python engine. -3. **Big-integer safety at the interface.** efmtool's Java core is arbitrary-precision (`DefaultBig - IntegerRationalMatrix`), so the *internal* arithmetic does not overflow. The hazard is at the - **marshalling boundary**: the compression structure and bound rescaling are read back through - `getDoubleRows` and `abs(subset_matrix[...])` in double precision (§3.11.2), and conservation - removal pushes `S` through a dense `JDouble` array. On models whose exact subdeterminants are huge — - the verified extreme is **yeast-GEM, needing ~263-bit coefficients** (§3.2.5) — a double cannot - represent those magnitudes, so bound rescaling and pattern detection silently lose precision. The - Python engine keeps *everything* in Python big integers / `Fraction` end to end and switches to a - dict-of-`Fraction` store above int64 (§3.2.5), so it is exact even on yeast-GEM. This is the single - most important reason the Python path is the default. -4. **It is the default and the tested path.** The measured pipeline numbers (§3.1, and the iML1515 - timings in CONTEXT) are all on `sparse_rref`; that is the code that receives ongoing correctness - work (e.g. the bounds-intersection fix of §3.4.4 / issue #44). - -**The trade-off, honestly stated.** efmtool is not bad code — it is a mature, well-tested Java library -whose fraction-free rational Gauss elimination is fast compiled code, and for a decade it *was* the -compression engine for this and related tools. If you have a JVM handy and a model whose coefficients -stay comfortably inside double range, `efmtool_rref` will produce a correct compression at competitive -speed. Its costs are the heavy dependency stack, the native-crash surface, and the double-precision -marshalling seam. Given a pure-Python alternative that is exact to arbitrary precision, needs no JVM, -and is the maintained default, the Java backend earns its "legacy" label: **there is essentially no -production reason to select it.** The realistic remaining uses are (a) cross-validation — regression- -testing the Python engine's output against the historical efmtool result on a model both can handle — -and (b) a fallback if a bug were ever found in the Python RREF. For everyday strain design, leave -`compression_backend` at its default. - -#### 3.11.4 Behavioral differences to be aware of - -The two backends are *intended* to produce the same lossless flux-space compression, but they are not -byte-identical and a few divergences are worth knowing: - -- **GPR propagation is identical on the coupled step.** Both backends clear gene rules before merging - and re-attach the AND-combined GPR from the saved AST snapshot in `compress_model_coupled` - (`compression.py`–), and the parallel OR-combine is always the Python - `compress_model_parallel` (§3.9). So GPR handling does *not* diverge between backends. -- **Protected reactions are honored only by the Python backend.** `compress_model` passes gene- - controlled reactions as `protected_reactions` (`no_coupled_compress_reacs`, `compression.py`– - ) so they survive COMPRESS #1 un-merged and gene multiplicity is preserved for GPR - integration (§3.4.2, [Ch 4](#ch4)). `compress_model_java` **ignores `protected_reactions`** — it reads only - `suppressed_reactions`, which `compress_model` never populates on this path. On the Java backend those - reactions can therefore be lumped in COMPRESS #1, a genuine semantic divergence in the gene-KO - pipeline. -- **Contradicting groups are handled differently.** efmtool's `CoupledContradicting` deletes groups it - finds inconsistent (the reason `suppressed_reactions` exists as a shield). The Python backend instead - computes the exact **bounds intersection** of the coupled group and removes only genuinely - empty/zero groups (§3.4.4). This is precisely the logic whose Java-era version "could drop reactions - incorrectly" — the cautionary tale of closed issue #44 ([Ch 10](#ch10)). The two backends can thus disagree on - which reactions a contradicting group costs you. -- **Direction bookkeeping differs.** The Java path physically flips `ub ≤ 0` reactions (`*= -1`) and - negates their recorded stoich (`efmtool_cmp_interface.py`–,–); the Python - coupled backend carries sign inside the exact `ratios` (§3.4.3). Same flux space, different maps — - which is fine because decompression ([Ch 9](#ch9)) consumes whichever map its backend produced. -- **Bound rescaling precision.** Java rescales merged-reaction bounds by a **double** - (`efmtool_cmp_interface.py`–); the Python backend intersects bounds using exact rationals - (§3.4.4). On well-scaled models this is invisible; on large-coefficient models it is another place the - Java path can drift. - -The safe reading: `efmtool_rref` is preserved for provenance and cross-checking, exercises the same -fixpoint and produces the same *kind* of map, but the exact-arithmetic Python backend is the one whose -compression you should trust for correctness-sensitive strain design. - - (ch4)= ## 4. GPR integration @@ -5926,9 +5726,7 @@ pytest tests -v --log-cli-level=INFO --junit-xml=test-results.xml ``` **CI matrix** (`.github/workflows/CI-test.yml`): OS `ubuntu-latest` / `windows-latest`; Python -`3.10`–`3.13`; both `pip` and `conda`. CPLEX is excluded for Python 3.13 (max supported: 3.12). A -JPype/JVM-shutdown segfault on Ubuntu is tolerated via a JUnit-XML exit-code check rather than the raw -process exit code. +`3.10`–`3.13`; both `pip` and `conda`. CPLEX is excluded for Python 3.13 (max supported: 3.12). **Correctness gates.** The canonical known-answer tests are the ones to keep green after any change to the pipeline: gene-level MCS on `e_coli_core` = **455** solutions, and on `iML1515` = **393**. These diff --git a/docs/source/examples/JN_08_compression.ipynb b/docs/source/examples/JN_08_compression.ipynb index 7ff6762..1800c97 100644 --- a/docs/source/examples/JN_08_compression.ipynb +++ b/docs/source/examples/JN_08_compression.ipynb @@ -3,7 +3,7 @@ { "cell_type": "markdown", "metadata": {}, - "source": "# Standalone network compression\n\nAn effective network compression is essential to any strain design computation. Since it may also be of interest outside the context of strain design, this example may help you using the network compression routine independently. Likewise, StrainDesign also offers the integration of GPR rules into the metabolic networks as a separate function.\n\nThe network compression routine removes blocked reactions, removes conservation relations and then alternately lumps **coupled** reactions (compress_model_coupled, using a sparse integer RREF nullspace algorithm by default) and **parallel** reactions (compress_model_parallel). The compression returns a compressed network and a list of so-called \"compression maps\". Each map consists of a dictionary that contains complete information for reversing the compression steps successively and expand information obtained from the compressed model to the full model. Each entry of each map contains the id of a compressed reaction, associated with the original reaction names and their factor (provided as a rational number) with which they were lumped.\n\nThe default is compression_backend='sparse_rref' (pure Python, no extra dependencies). A legacy Java-based backend, compression_backend='efmtool_rref', is also available via pip install straindesign[java].\n\nFurthermore, the user can select reactions that should be exempt from the parallel compression. In the following, we provide the code snippet that can be used to call the compression." + "source": "# Standalone network compression\n\nAn effective network compression is essential to any strain design computation. Since it may also be of interest outside the context of strain design, this example may help you using the network compression routine independently. Likewise, StrainDesign also offers the integration of GPR rules into the metabolic networks as a separate function.\n\nThe network compression routine removes blocked reactions, removes conservation relations and then alternately lumps **coupled** reactions (compress_model_coupled, using a sparse integer RREF nullspace algorithm) and **parallel** reactions (compress_model_parallel). The compression returns a compressed network and a list of so-called \"compression maps\". Each map consists of a dictionary that contains complete information for reversing the compression steps successively and expand information obtained from the compressed model to the full model. Each entry of each map contains the id of a compressed reaction, associated with the original reaction names and their factor (provided as a rational number) with which they were lumped.\n\nFurthermore, the user can select reactions that should be exempt from the parallel compression. In the following, we provide the code snippet that can be used to call the compression." }, { "cell_type": "code", @@ -257,4 +257,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} \ No newline at end of file +} diff --git a/docs/source/index.rst b/docs/source/index.rst index ce5e24c..a053158 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -75,7 +75,7 @@ The comprehensive StrainDesign package for MILP-based strain design computation :width: 40% :alt: Plot animation -The default compression uses a pure Python sparse RREF implementation. A legacy Java-based compression via EFMTool :html:`[6]` is optionally available (see :doc:`legacy_methods`). Note that the Java backend (via JPype) is known to conflict with CPLEX's native library when both are loaded in the same Python session. If you use CPLEX, we recommend the default Python compression backend. +Network compression uses a pure Python sparse RREF implementation with exact rational arithmetic, following the approach of EFMTool :html:`[6]`. It has no dependencies beyond NumPy/SciPy. :html:``\ Installation: ================================================ @@ -136,7 +136,6 @@ How to cite: examples/JN_08_compression.ipynb 9_cnapy_integration api_reference - legacy_methods developers_guide .. diff --git a/docs/source/legacy_methods.rst b/docs/source/legacy_methods.rst deleted file mode 100644 index 7a2cfd6..0000000 --- a/docs/source/legacy_methods.rst +++ /dev/null @@ -1,42 +0,0 @@ -Legacy Methods -============== - -This page documents optional legacy functionality that requires additional -dependencies beyond the core StrainDesign installation. - -Java-based EFMTool compression (``compression_backend='efmtool_rref'``) ------------------------------------------------------------------------ - -The default compression backend is ``compression_backend='sparse_rref'``, a pure Python -implementation with no extra dependencies. A legacy Java-based backend is -available for comparison or reproducibility purposes. - -To use it, install the optional Java dependency:: - - pip install straindesign[java] - -or:: - - pip install jpype1 - -Then pass ``compression_backend='efmtool_rref'`` to :func:`~straindesign.compress_model` -or to ``compute_strain_designs`` via the ``compression_backend`` keyword argument. - -JAVA_HOME path --------------- - -In some cases, using the ``efmtool_rref`` backend may fail with: - -``JVMNotFoundException: No JVM shared library file (libjli.dylib) found. Try setting up the JAVA_HOME environment variable.`` - -In this case, make sure Java is installed correctly and the JAVA_HOME variable -is set. See `JAVA_HOME environment variable `_ -for platform-specific instructions. - -If you're on OS X and get the error - -``OSError: [Errno 0] JVM DLL not found`` - -check that your `Java and the JPype library is set up correctly `_. -The easiest way to avoid this error is to use conda to install StrainDesign and -Java together. diff --git a/pyproject.toml b/pyproject.toml index de5717b..d4265de 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,9 +31,6 @@ dependencies = [ "psutil", ] -[project.optional-dependencies] -java = ["jpype1"] - [project.urls] Homepage = "https://github.com/klamt-lab/straindesign" Documentation = "https://straindesign.readthedocs.io/en/latest/index.html" @@ -43,5 +40,3 @@ Source = "https://github.com/klamt-lab/straindesign/" [tool.setuptools.packages.find] include = ["straindesign*"] -[tool.setuptools.package-data] -straindesign = ["efmtool.jar"] diff --git a/straindesign/__init__.py b/straindesign/__init__.py index 34e0f36..ac79c2f 100644 --- a/straindesign/__init__.py +++ b/straindesign/__init__.py @@ -43,15 +43,6 @@ def __exit__(self, exit_type, exit_value, exit_traceback): if module_exists("pyscipopt"): avail_solvers.add(SCIP) -# Conditional eager JVM startup — required for stable JPype operation. -# The JVM must start before NumPy/OpenBLAS spawns worker threads, otherwise -# JNI calls crash with SIGBUS/SIGSEGV (jpype#808, jpype#934). -# No-op when jpype1 or Java is not installed (neither is a dependency). -# See developers_guide.md "efmtool_cmp_interface.py — JPype/JVM Initialization". -from .efmtool_cmp_interface import _start_jvm as _start_jvm -_start_jvm() -del _start_jvm - from .solver_interface import * from .indicatorConstraints import * from .pool import * @@ -64,4 +55,4 @@ def __exit__(self, exit_type, exit_value, exit_traceback): from .strainDesignProblem import * from .strainDesignMILP import * from .compute_strain_designs import * -from .compression import sparse_nullspace, sparse_nullspace as nullspace, RationalMatrix, ExactCOO +from .compression import sparse_nullspace, RationalMatrix, ExactCOO diff --git a/straindesign/compression.py b/straindesign/compression.py index 262a308..8d450c9 100644 --- a/straindesign/compression.py +++ b/straindesign/compression.py @@ -2092,7 +2092,7 @@ def simplify_model_gprs(model, budget=50000): logging.info(' GPR rule simplification: %d rules, %d rewritten.' % (n, nchg)) -def compress_model(model, no_par_compress_reacs=set(), compression_backend='sparse_rref', propagate_gpr=False, +def compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False, no_coupled_compress_reacs=set()): """Compress a metabolic model using multiple techniques. @@ -2108,11 +2108,6 @@ def compress_model(model, no_par_compress_reacs=set(), compression_backend='spar un-merged through COMPRESS#1 so that gene multiplicity is preserved exactly once GPR rules are integrated (correct gene-regulatory semantics under compression). To also exempt them from parallel merging, include them in no_par_compress_reacs. - compression_backend: Compression backend to use: - - 'sparse_rref' (default): Pure Python sparse integer RREF. - No external dependencies beyond NumPy/SciPy. - - 'efmtool_rref' (legacy): Java-based EFMTool via JPype. - Requires a JVM and the jpype1 package. propagate_gpr: If True, propagate and simplify GPR rules through compression (AND for coupled, OR for parallel merges). Empty GPR rules are correctly handled: skipped in AND (always @@ -2126,16 +2121,6 @@ def compress_model(model, no_par_compress_reacs=set(), compression_backend='spar no_coupled_compress_reacs = set(no_coupled_compress_reacs) with suppress_lp_context(model): cmp_mapReac = [] - use_java = (compression_backend == 'efmtool_rref') - if use_java: - # The Python compressor re-expresses each lump in one member's units (see - # StoichMatrixCompressor._restore_group_scale); the legacy Java backend does not, so a - # lump can come out at an extreme scale. The returned map carries the factor, so - # expanding a design stays exact -- but a bound stated on a lumped reaction is read in - # the lump's units, which is how 'biomass >= 0.001' can end up below feasibility tolerance. - LOG.warning(' Compression backend "efmtool_rref" does not normalize lumped-reaction ' - 'scales; bounds and constraints on lumped reactions are expressed in the ' - 'lump\'s units. Use "sparse_rref" if you constrain lumped reactions.') LOG.info(' Removing blocked reactions.') remove_blocked_reactions(model) LOG.info(' Converting coefficients to rationals.') @@ -2155,10 +2140,7 @@ def compress_model(model, no_par_compress_reacs=set(), compression_backend='spar cmp_mapReac.append({"reac_map_exp": reac_map_exp, "parallel": True}) # 2. Conservation relation removal (reduces S rows for RREF) - if use_java: - _remove_conservation_relations_java(model) - else: - remove_conservation_relations(model) + remove_conservation_relations(model) # 3. Exit if either parallel or coupled found nothing (after # at least one full cycle). If one step found nothing, @@ -2171,7 +2153,7 @@ def compress_model(model, no_par_compress_reacs=set(), compression_backend='spar # 4. Coupled (expensive — nullspace/RREF) numr_pre = len(model.reactions) LOG.info(f' Compression {run}: Lumping coupled reactions.') - reac_map_exp = compress_model_coupled(model, compression_backend, + reac_map_exp = compress_model_coupled(model, propagate_gpr=propagate_gpr, protected_reactions=no_coupled_compress_reacs) for new_reac, old_reac_val in reac_map_exp.items(): @@ -2198,38 +2180,20 @@ def compress_model(model, no_par_compress_reacs=set(), compression_backend='spar return cmp_mapReac -def _remove_conservation_relations_java(model) -> None: - """Remove conservation relations using Java efmtool.""" - from . import efmtool_cmp_interface as efm - stoich_mat = create_stoichiometric_matrix(model, array_type='lil') - basic_mets = efm.basic_columns_rat_java(stoich_mat.transpose().toarray(), tolerance=0) - dependent = [model.metabolites[i] for i in set(range(len(model.metabolites))) - set(basic_mets)] - if dependent: - model.remove_metabolites(dependent) - - -def compress_model_coupled(model, compression_backend='sparse_rref', propagate_gpr=False, - suppressed_reactions=set(), protected_reactions=set()): +def compress_model_coupled(model, propagate_gpr=False, protected_reactions=set()): """Compress by lumping stoichiometrically coupled (dependent) reactions. Identifies groups of reactions whose flux vectors are proportional in every steady state (i.e. they share a common nullspace direction) and merges each - group into a single lumped reaction. Both the pure-Python and legacy Java - backends perform this operation; the compression_backend controls the nullspace algorithm. + group into a single lumped reaction, via the sparse integer RREF nullspace. Args: model: COBRA model to compress in-place - compression_backend: 'sparse_rref' (default, Python) or 'efmtool_rref' (Java legacy) propagate_gpr: If True, AND-combine GPR rules of merged reactions (with sympy simplification). Empty GPRs are skipped. Default False. - suppressed_reactions: Set of reaction IDs to exclude from compression - (Java backend only). Used to protect reactions referenced in strain - design constraints from being deleted by the Java compressor's - CoupledContradicting logic. Ignored for the Python backend (which - handles contradicting groups correctly via bounds intersection). protected_reactions: Set of reaction IDs to exempt from coupled merging (kept as their own reactions; the rest of their coupled group still - merges). Python (sparse_rref) backend only. Used to keep gene-controlled + merges). Used to keep gene-controlled reactions intact through compression before GPR integration so that the gene multiplicity is preserved (correct gene-regulatory semantics). @@ -2239,27 +2203,18 @@ def compress_model_coupled(model, compression_backend='sparse_rref', propagate_g # Compression is pure linear algebra; keep it off the optlang solver. from straindesign.networktools import suppress_lp_context with suppress_lp_context(model): - # Save GPR AST bodies before either backend clears them + # Save GPR AST bodies before compression clears them if propagate_gpr: saved_gpr_bodies = {r.id: r.gpr.body for r in model.reactions} - if compression_backend == 'efmtool_rref': - from .efmtool_cmp_interface import compress_model_java - reaction_map = compress_model_java(model, suppressed_reactions=suppressed_reactions) - # Clean up any remaining zero-flux reactions that the Java compressor created. - zero_flux = {r for r in model.reactions if r.lower_bound == 0 and r.upper_bound == 0} - for r in zero_flux: - reaction_map.pop(r.id, None) - if zero_flux: - model.remove_reactions(list(zero_flux), remove_orphans=True) - else: - # Clear gene rules to match Java behavior - for r in model.reactions: - r.gene_reaction_rule = '' + # Gene rules are cleared here and re-derived below from the saved ASTs, so a lumped + # reaction's rule is the AND-combination of its members rather than one member's. + for r in model.reactions: + r.gene_reaction_rule = '' - result = compress_cobra_model(model, methods=CompressionMethod.standard(), in_place=True, - protected_reactions=protected_reactions) - reaction_map = result.reaction_map + result = compress_cobra_model(model, methods=CompressionMethod.standard(), in_place=True, + protected_reactions=protected_reactions) + reaction_map = result.reaction_map # Propagate GPR rules: AND-combine contributing reactions' GPR ASTs if propagate_gpr: @@ -2274,11 +2229,6 @@ def compress_model_coupled(model, compression_backend='sparse_rref', propagate_g return reaction_map -# Backward-compatibility alias (old name referenced efmtool, but the function -# is backend-agnostic — the new name compress_model_coupled is preferred). -compress_model_efmtool = compress_model_coupled - - def compress_model_parallel(model, protected_rxns=set(), propagate_gpr=False): """Compress by lumping parallel reactions. @@ -2413,7 +2363,6 @@ def _parallel_key(i): # High-level API 'compress_model', 'compress_model_coupled', - 'compress_model_efmtool', # backward-compat alias 'compress_model_parallel', # GPR propagation helpers '_gpr_ast_to_expr', diff --git a/straindesign/compute_strain_designs.py b/straindesign/compute_strain_designs.py index 7e53bc8..5cfd258 100644 --- a/straindesign/compute_strain_designs.py +++ b/straindesign/compute_strain_designs.py @@ -386,7 +386,7 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: """ allowed_keys = { MODULES, SETUP, SOLVER, MAX_COST, MAX_SOLUTIONS, 'M', 'compress', 'gene_kos', KOCOST, KICOST, GKOCOST, GKICOST, REGCOST, - SOLUTION_APPROACH, 'advanced', 'use_scenario', T_LIMIT, SEED, MILP_THREADS, 'compression_backend', 'dump_preprocessed' + SOLUTION_APPROACH, 'advanced', 'use_scenario', T_LIMIT, SEED, MILP_THREADS, 'dump_preprocessed' } logging.info('Preparing strain design computation.') if SETUP in kwargs: @@ -564,7 +564,6 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: # also exempt them from parallel merging so their names stay stable across the # compression passes (keeps the coupled-exemption matching them by name) no_par_compress_reacs.update(no_coupled_compress_reacs) - compression_backend = kwargs.get('compression_backend', 'sparse_rref') # --- Reversibility pre-tightening (BEFORE compress #1) --- # Sign-only FVA (cheaper than full FVA): fix lb/ub to 0 for directions carrying no flux in the # base polytope. Design-neutral (a base-infeasible direction stays infeasible under any module @@ -586,7 +585,6 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: logging.info('Compressing Network (' + str(len(cmp_model.reactions)) + ' reactions).') t0 = time.time() cmp_mapReac_1 = compress_model(cmp_model, no_par_compress_reacs, - compression_backend=compression_backend, propagate_gpr=True, no_coupled_compress_reacs=no_coupled_compress_reacs) sd_modules = compress_modules(sd_modules, cmp_mapReac_1) @@ -671,7 +669,7 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: t0 = time.time() no_par_compress_reacs = _collect_no_par_compress_reacs(sd_modules) cmp_mapReac_2 = compress_model(cmp_model, no_par_compress_reacs, - compression_backend=compression_backend) +) sd_modules = compress_modules(sd_modules, cmp_mapReac_2) cmp_ko_cost, cmp_ki_cost, cmp_mapReac_2 = compress_ki_ko_cost( cmp_ko_cost, cmp_ki_cost, cmp_mapReac_2) diff --git a/straindesign/efmtool.jar b/straindesign/efmtool.jar deleted file mode 100644 index c189479..0000000 Binary files a/straindesign/efmtool.jar and /dev/null differ diff --git a/straindesign/efmtool_cmp_interface.py b/straindesign/efmtool_cmp_interface.py deleted file mode 100644 index 190df4d..0000000 --- a/straindesign/efmtool_cmp_interface.py +++ /dev/null @@ -1,529 +0,0 @@ -#!/usr/bin/env python3 -# -# Copyright 2022 Max Planck Insitute Magdeburg -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -"""EFMtool compression interface for straindesign. - -This module provides compression utilities for metabolic networks. -By default, the pure Python 'sparse_rref' compression backend is used. -The Java EFMTool backend is available via compression_backend='efmtool_rref' (requires jpype1). - -For the documentation of the compression API provided by StrainDesign, -refer to straindesign.compression.compress_model. -""" - -import logging -import numpy as np -import os -import sys - -# ============================================================================= -# Pure Python Implementation (Default) -# ============================================================================= - - -def basic_columns_rat(mx, tolerance=0): - """Find basic columns using exact rational arithmetic (FLINT or sympy).""" - from .compression import basic_columns_from_numpy - return basic_columns_from_numpy(mx) - - -# ============================================================================= -# Lazy Java Initialization -# ============================================================================= - -_JAVA_INITIALIZED = False -_JPYPE_AVAILABLE = None - -# Java classes (populated by _init_java) -DefaultBigIntegerRationalMatrix = None -Gauss = None -CompressionMethod = None -StoichMatrixCompressor = None -BigFraction = None -BigInteger = None -subset_compression = None -jTrue = None -jSystem = None - - -def _check_jpype_available(): - """Check if jpype is available without importing it.""" - global _JPYPE_AVAILABLE - if _JPYPE_AVAILABLE is None: - import importlib.util - _JPYPE_AVAILABLE = importlib.util.find_spec("jpype") is not None - return _JPYPE_AVAILABLE - - -def _check_sympy_available(): - """Check if sympy is available without importing it.""" - import importlib.util - return importlib.util.find_spec("sympy") is not None - - -def _search_for_jvm(): - """Search for JVM in common locations.""" - common_java_paths = [ - "C:\\Program Files\\Java", # Windows - "/usr/lib/jvm", # Linux - "/Library/Java/JavaVirtualMachines", # macOS - os.path.dirname(sys.executable) - ] - for base in common_java_paths: - if os.path.exists(base): - for root, _dirs, files in os.walk(base): - if any(lib in files for lib in ["jvm.dll", "libjvm.so", "libjvm.dylib"]): - return root - return None - - -def _start_jvm(): - """Start the JVM and load Java classes, matching the v1.14 eager init pattern. - - Called eagerly at package import time (from __init__.py) when jpype and Java - are available. This replicates the v1.14 behaviour where ``import straindesign`` - immediately started the JVM and loaded all efmtool Java classes via - ``import jpype.imports``. That approach is the only one known to be stable - on Linux/macOS CI runners; deferred JVM startup or JClass-based loading - causes SIGBUS/SIGSEGV on larger matrices (iMLcore+). - - No-op when jpype or Java is not installed. - """ - global _JAVA_INITIALIZED - global DefaultBigIntegerRationalMatrix, Gauss, CompressionMethod - global StoichMatrixCompressor, BigFraction, BigInteger - global subset_compression, jTrue, jSystem - - if _JAVA_INITIALIZED: - return - - if not _check_jpype_available(): - return # jpype not installed — nothing to do - - import jpype - import io - from contextlib import redirect_stdout, redirect_stderr - - # Add efmtool.jar to classpath - efmtool_jar = os.path.join(os.path.dirname(__file__), 'efmtool.jar') - if os.path.exists(efmtool_jar): - jpype.addClassPath(efmtool_jar) - - if not jpype.isJVMStarted(): - # Look up JVM at different locations - if not os.environ.get("JAVA_HOME"): - candidate = _search_for_jvm() - if candidate: - os.environ["JAVA_HOME"] = candidate - - try: - # Suppress faulthandler during JVM startup to prevent ugly - # "Windows fatal exception: access violation" messages. - import faulthandler as _fh - _fh_was_enabled = _fh.is_enabled() - if _fh_was_enabled: - _fh.disable() - try: - with redirect_stdout(io.StringIO()), redirect_stderr(io.StringIO()): - # --enable-native-access: allow JPype's System.load() calls - # without warnings on Java 17+ and prevent blocking on 24+. - jpype.startJVM("--enable-native-access=ALL-UNNAMED") - finally: - if _fh_was_enabled: - _fh.enable() - # Register explicit JVM shutdown to avoid SIGSEGV during Python - # exit (known JPype race condition in _JTerminate, see - # https://github.com/jpype-project/jpype/issues/842). - import atexit - def _shutdown_jvm(): - try: - import jpype as _jp - if _jp.isJVMStarted(): - _jp.shutdownJVM() - except Exception: - pass - atexit.register(_shutdown_jvm) - except Exception: - return # JVM startup failed — will raise a clear error in _init_java() - - # Load Java classes via import (v1.14 pattern) — this activates JPype's - # import hook which sets up JNI references differently from JClass(). - try: - import jpype.imports # noqa: F401 — activates the import hook - - import ch.javasoft.smx.impl.DefaultBigIntegerRationalMatrix as _DBIRM - import ch.javasoft.smx.ops.Gauss as _Gauss - import ch.javasoft.metabolic.compress.CompressionMethod as _CM - import ch.javasoft.metabolic.compress.StoichMatrixCompressor as _SMC - import ch.javasoft.math.BigFraction as _BF - import java.math.BigInteger as _BI - - DefaultBigIntegerRationalMatrix = _DBIRM - Gauss = _Gauss - CompressionMethod = _CM - StoichMatrixCompressor = _SMC - BigFraction = _BF - BigInteger = _BI - - subset_compression = CompressionMethod[:]( - [CompressionMethod.CoupledZero, CompressionMethod.CoupledCombine, - CompressionMethod.CoupledContradicting]) - jTrue = jpype.JBoolean(True) - jSystem = jpype.JClass("java.lang.System") - - _JAVA_INITIALIZED = True - except Exception: - pass # class loading failed — _init_java() will retry or raise - - -def _init_java(): - """Ensure Java classes are loaded. Usually a no-op (classes loaded eagerly - by _start_jvm). Falls back to JClass loading if eager init was skipped.""" - global _JAVA_INITIALIZED - global DefaultBigIntegerRationalMatrix, Gauss, CompressionMethod - global StoichMatrixCompressor, BigFraction, BigInteger - global subset_compression, jTrue, jSystem - - if _JAVA_INITIALIZED: - return - - if not _check_jpype_available(): - raise ImportError("jpype1 is not installed. Legacy Java compression requires jpype1.\n" - "Install with: pip install jpype1\n" - "Or use the default Python compression (compression_backend='sparse_rref').") - - if not _check_sympy_available(): - raise ImportError("sympy is not installed. Legacy Java compression requires sympy.\n" - "Install with: pip install sympy\n" - "Or use the default Python compression (compression_backend='sparse_rref').") - - import jpype - - # Try eager init first (may have been skipped if _start_jvm wasn't called) - _start_jvm() - if _JAVA_INITIALIZED: - return - - # Fallback: start JVM and load classes via JClass - if not jpype.isJVMStarted(): - extra_info = "" - if not os.environ.get("JAVA_HOME"): - extra_info = " JAVA_HOME is not defined." - raise RuntimeError( - "Failed to start JVM. Please ensure that Java (OpenJDK) is installed." + extra_info + - " If using conda, install openjdk from conda-forge and set JAVA_HOME to the OpenJDK installation path.") - - efmtool_jar = os.path.join(os.path.dirname(__file__), 'efmtool.jar') - if not os.path.exists(efmtool_jar): - raise FileNotFoundError(f"efmtool.jar not found at {efmtool_jar}. " - "Legacy Java compression requires the efmtool.jar file.") - jpype.addClassPath(efmtool_jar) - - try: - DefaultBigIntegerRationalMatrix = jpype.JClass('ch.javasoft.smx.impl.DefaultBigIntegerRationalMatrix') - Gauss = jpype.JClass('ch.javasoft.smx.ops.Gauss') - CompressionMethod = jpype.JClass('ch.javasoft.metabolic.compress.CompressionMethod') - StoichMatrixCompressor = jpype.JClass('ch.javasoft.metabolic.compress.StoichMatrixCompressor') - BigFraction = jpype.JClass('ch.javasoft.math.BigFraction') - BigInteger = jpype.JClass('java.math.BigInteger') - except Exception as e: - raise RuntimeError( - "Failed to load EFMTool Java classes. The JVM started but the efmtool.jar " - "classes could not be loaded. Use compression_backend='sparse_rref' instead.") from e - - subset_compression = CompressionMethod[:]( - [CompressionMethod.CoupledZero, CompressionMethod.CoupledCombine, CompressionMethod.CoupledContradicting]) - jTrue = jpype.JBoolean(True) - jSystem = jpype.JClass("java.lang.System") - - _JAVA_INITIALIZED = True - - -# ============================================================================= -# Java Conversion Utilities -# ============================================================================= - - -def numpy_mat2jpypeArrayOfArrays(npmat): - """Convert numpy matrix to jpype array of arrays (requires Java init).""" - _init_java() - import jpype - - rows = npmat.shape[0] - cols = npmat.shape[1] - jmat = jpype.JDouble[rows, cols] - for r in range(rows): - for c in range(cols): - jmat[r][c] = npmat[r, c] - return jmat - - -def jpypeArrayOfArrays2numpy_mat(jmat): - """Convert jpype array of arrays to numpy matrix.""" - rows = len(jmat) - cols = len(jmat[0]) - npmat = np.zeros((rows, cols)) - for r in range(rows): - for c in range(cols): - npmat[r, c] = jmat[r][c] - return npmat - - -def sympyRat2jBigIntegerPair(val): - """Convert Fraction or sympy Rational to Java BigInteger pair (requires Java init).""" - _init_java() - - # Support both fractions.Fraction (.numerator/.denominator) and sympy.Rational (.p/.q) - numer = val.numerator if hasattr(val, 'numerator') else val.p - if numer.bit_length() <= 63: - numer = BigInteger.valueOf(numer) - else: - numer = BigInteger(str(numer)) - - denom = val.denominator if hasattr(val, 'denominator') else val.q - if denom.bit_length() <= 63: - denom = BigInteger.valueOf(denom) - else: - denom = BigInteger(str(denom)) - - return (numer, denom) - - -def jBigFraction2sympyRat(val): - """Convert Java BigFraction to sympy Rational (requires Java init).""" - return jBigIntegerPair2sympyRat(val.getNumerator(), val.getDenominator()) - - -def jBigFraction2fraction(val): - """Convert Java BigFraction to fractions.Fraction.""" - from fractions import Fraction - r = jBigIntegerPair2sympyRat(val.getNumerator(), val.getDenominator()) - return Fraction(int(r.p), int(r.q)) - - -def jBigIntegerPair2sympyRat(numer, denom): - """Convert Java BigInteger pair to sympy Rational (requires sympy).""" - import sympy - - if numer.bitLength() <= 63: - numer = numer.longValue() - else: - numer = str(numer.toString()) - - if denom.bitLength() <= 63: - denom = denom.longValue() - else: - denom = str(denom.toString()) - - return sympy.Rational(numer, denom) - - -# ============================================================================= -# Legacy Java Compression Functions -# ============================================================================= - - -def basic_columns_rat_java(mx, tolerance=0): - """ - Find basic columns using Java Gaussian elimination. - - Legacy implementation using jpype and Java efmtool. - Requires jpype1 and sympy to be installed. - - Args: - mx: Matrix (numpy array or Java matrix) - tolerance: Tolerance (unused in exact arithmetic) - - Returns: - Array of indices of basic columns - """ - _init_java() - import gc - import jpype - - if isinstance(mx, np.ndarray): - mx = DefaultBigIntegerRationalMatrix(numpy_mat2jpypeArrayOfArrays(mx), jTrue, jTrue) - - row_map = jpype.JInt[mx.getRowCount()] - col_map = jpype.JInt[:](range(mx.getColumnCount())) - # Disable GC during the Java call — Python's garbage collector can - # attempt to finalize JPype proxy objects mid-computation, causing - # Bus error / SIGSEGV on macOS and Linux. - gc.disable() - try: - rank = Gauss.getRationalInstance().rowEchelon(mx, False, row_map, col_map) - finally: - gc.enable() - - return col_map[0:rank] - - -def compress_model_java(model, suppressed_reactions=set()): - """Legacy Java compression using jpype (requires jpype and sympy). - - Args: - model: COBRA model (will be modified in place) - suppressed_reactions: Set of reaction IDs to exclude from compression. - These reactions are kept as standalone entries with identity mapping. - Used to protect reactions referenced in strain design constraints - from being deleted by the Java compressor's CoupledContradicting logic. - - Returns: - dict: Reaction map from compressed to original reactions with scaling factors - """ - import jpype - from .networktools import stoichmat_coeff_to_fraction - - # Initialize Java if not already done - _init_java() - - # Convert to rational coefficients for Java - stoichmat_coeff_to_fraction(model) - - for r in model.reactions: - r.gene_reaction_rule = '' - - suppressed_set = set(suppressed_reactions) if suppressed_reactions else set() - num_met = len(model.metabolites) - num_reac = len(model.reactions) - old_reac_ids = [r.id for r in model.reactions] - - # Build mapping between active (non-suppressed) indices and model indices - active_to_model = [i for i in range(num_reac) if old_reac_ids[i] not in suppressed_set] - num_active = len(active_to_model) - - # Disable GC for the entire Java interaction block — Python's garbage - # collector can finalize JPype proxy objects mid-JNI call, causing - # Bus error / SIGSEGV (non-deterministic, see jpype-project/jpype#934). - import gc - gc.disable() - try: - stoich_mat = DefaultBigIntegerRationalMatrix(num_met, num_active) - reversible = jpype.JBoolean[:]([model.reactions[active_to_model[ai]].reversibility for ai in range(num_active)]) - flipped = set() - for ai in range(num_active): - mi = active_to_model[ai] - if model.reactions[mi].upper_bound <= 0: - model.reactions[mi] *= -1 - flipped.add(ai) - logging.debug("Flipped " + model.reactions[mi].id) - for k, v in model.reactions[mi]._metabolites.items(): - n, d = sympyRat2jBigIntegerPair(v) - stoich_mat.setValueAt(model.metabolites.index(k.id), ai, BigFraction(n, d)) - - # Compress active reactions only - smc = StoichMatrixCompressor(subset_compression) - reacNames = jpype.JString[:]([old_reac_ids[active_to_model[ai]] for ai in range(num_active)]) - comprec = smc.compress(stoich_mat, reversible, jpype.JString[num_met], reacNames, None) - subset_matrix = jpypeArrayOfArrays2numpy_mat(comprec.post.getDoubleRows()) - finally: - gc.enable() - - # subset_matrix shape: (num_active, num_compressed) - del_model = np.zeros(num_reac, dtype=bool) - - # Mark zero-flux active reactions for deletion - for ai in range(num_active): - if not np.any(subset_matrix[ai, :]): - del_model[active_to_model[ai]] = True - - for j in range(subset_matrix.shape[1]): - rxn_ai = subset_matrix[:, j].nonzero()[0] - if len(rxn_ai) == 0: - continue - r0_mi = active_to_model[rxn_ai[0]] - model.reactions[r0_mi].subset_rxns = [] - model.reactions[r0_mi].subset_stoich = [] - for ai in rxn_ai: - mi = active_to_model[ai] - factor = jBigFraction2fraction(comprec.post.getBigFractionValueAt(ai, j)) - model.reactions[mi] *= factor - if model.reactions[mi].lower_bound not in (0, -float('inf')): - model.reactions[mi].lower_bound /= abs(subset_matrix[ai, j]) - if model.reactions[mi].upper_bound not in (0, float('inf')): - model.reactions[mi].upper_bound /= abs(subset_matrix[ai, j]) - model.reactions[r0_mi].subset_rxns.append(mi) - if ai in flipped: - model.reactions[r0_mi].subset_stoich.append(-factor) - else: - model.reactions[r0_mi].subset_stoich.append(factor) - for ai in rxn_ai[1:]: - mi = active_to_model[ai] - if len(model.reactions[r0_mi].id) + len(model.reactions[mi].id) < 220 and model.reactions[r0_mi].id[-3:] != '...': - model.reactions[r0_mi].id += '*' + model.reactions[mi].id - elif not model.reactions[r0_mi].id[-3:] == '...': - model.reactions[r0_mi].id += '...' - model.reactions[r0_mi] += model.reactions[mi] - if model.reactions[mi].lower_bound > model.reactions[r0_mi].lower_bound: - model.reactions[r0_mi].lower_bound = model.reactions[mi].lower_bound - if model.reactions[mi].upper_bound < model.reactions[r0_mi].upper_bound: - model.reactions[r0_mi].upper_bound = model.reactions[mi].upper_bound - del_model[mi] = True - # Update stored objective through compression factors - if len(rxn_ai) > 1: - _obj = getattr(model, '_suppressed_obj', None) - if _obj is not None: - merged_obj = 0.0 - for ai in rxn_ai: - mi = active_to_model[ai] - factor = jBigFraction2fraction(comprec.post.getBigFractionValueAt(ai, j)) - merged_obj += _obj.pop(old_reac_ids[mi], 0.0) * float(factor) - if merged_obj != 0: - _obj[model.reactions[r0_mi].id] = merged_obj - - # Add suppressed reactions as standalone entries - from fractions import Fraction - for mi in range(num_reac): - if old_reac_ids[mi] in suppressed_set: - model.reactions[mi].subset_rxns = [mi] - model.reactions[mi].subset_stoich = [Fraction(1)] - - # Delete reactions (reverse order to preserve indices) - del_indices = np.where(del_model)[0] - for i in range(len(del_indices) - 1, -1, -1): - model.reactions[del_indices[i]].remove_from_model(remove_orphans=True) - - # Build rational_map - rational_map = {} - for j in range(len(model.reactions)): - rational_map[model.reactions[j].id] = { - old_reac_ids[mi]: v - for mi, v in zip(model.reactions[j].subset_rxns, model.reactions[j].subset_stoich) - } - return rational_map - - -# ============================================================================= -# Exports -# ============================================================================= - -__all__ = [ - # Pure Python - 'basic_columns_rat', - # Java initialization - '_start_jvm', - '_init_java', - '_check_jpype_available', - # Java compression - 'basic_columns_rat_java', - 'compress_model_java', - # Java conversion utilities - 'numpy_mat2jpypeArrayOfArrays', - 'jpypeArrayOfArrays2numpy_mat', - 'sympyRat2jBigIntegerPair', - 'jBigFraction2sympyRat', - 'jBigIntegerPair2sympyRat', -] diff --git a/straindesign/lptools.py b/straindesign/lptools.py index d70f79d..663869f 100644 --- a/straindesign/lptools.py +++ b/straindesign/lptools.py @@ -389,52 +389,6 @@ def fva_legacy(model, **kwargs) -> DataFrame: ) -def remove_redundant_bounds(model, **kwargs) -> DataFrame: - """Remove non-binding bounds from a model using FVA. - - Runs FVA and relaxes bounds that never bind at steady state: - - If fva_min > lb + tol: set lb = -inf (lower bound is not binding) - - If fva_max < ub - tol: set ub = +inf (upper bound is not binding) - - Modifies the model IN-PLACE. Returns the FVA DataFrame. - - Args: - model (cobra.Model): - A metabolic model. Modified in-place. - - solver (optional (str)): - Solver for FVA. - - constraints (optional): - Constraints passed through to fva(). - - compress (optional (bool)): - Compress before FVA (passed through). - - threads (optional (int)): - Parallel threads for FVA (passed through). - - tol (optional (float)): (Default: 1e-6) - Tolerance for considering a bound as binding. - - Returns: - (pandas.DataFrame): - FVA results with 'minimum' and 'maximum' columns. - """ - tol = kwargs.pop('tol', 1e-6) - fva_result = fva(model, **kwargs) - - for rxn in model.reactions: - fva_min = fva_result.loc[rxn.id, 'minimum'] - fva_max = fva_result.loc[rxn.id, 'maximum'] - if fva_min > rxn.lower_bound + tol: - rxn.lower_bound = -float('inf') - if fva_max < rxn.upper_bound - tol: - rxn.upper_bound = float('inf') - - return fva_result - - def fba(model, **kwargs) -> Solution: """Flux Balance Analysis (FBA), parsimonius Flux Balance Analysis (pFBA), diff --git a/straindesign/networktools.py b/straindesign/networktools.py index 120cdf7..f7b0356 100644 --- a/straindesign/networktools.py +++ b/straindesign/networktools.py @@ -484,7 +484,6 @@ def _silent_io(): from straindesign.compression import ( compress_model, compress_model_coupled, - compress_model_efmtool, # backward-compat alias compress_model_parallel, remove_blocked_reactions, remove_ext_mets, @@ -1384,21 +1383,6 @@ def modules_coeff_to_fraction(sd_modules): return sd_modules -def modules_coeff2float(sd_modules): - """Convert coefficients occurring in SDModule objects to floats""" - for i, module in enumerate(sd_modules): - for param in [CONSTRAINTS, INNER_OBJECTIVE, OUTER_OBJECTIVE, PROD_ID]: - if param in module and module[param] is not None: - if param == CONSTRAINTS: - for constr in module[CONSTRAINTS]: - for reac in constr[0].keys(): - constr[0][reac] = float(constr[0][reac]) - if param in [INNER_OBJECTIVE, OUTER_OBJECTIVE, PROD_ID]: - for reac in module[param].keys(): - module[param][reac] = float(module[param][reac]) - return sd_modules - - def bound_blocked_or_irrevers_fva(model, **kwargs): """Use FVA to determine flux ranges, update model bounds, and return the ranges. diff --git a/straindesign/parse_constr.py b/straindesign/parse_constr.py index f1db034..cb7b89e 100644 --- a/straindesign/parse_constr.py +++ b/straindesign/parse_constr.py @@ -86,58 +86,6 @@ def parse_linexpr(expr, reaction_ids) -> List: return [linexpr2dict(e, reaction_ids) if type(e) is str else e for e in expr] -def lineq2mat(equations, reaction_ids) -> Tuple[sparse.csr_matrix, Tuple, sparse.csr_matrix, Tuple]: - """Translates *linear* (in)equalities to matrices - - Input inequalities in the form of strings is translated into matrices and vectors. The reaction - list defines the order of variables and thus the columns of the resulting matrices, the order - of (in)equalities will be preserved in the output matrices. As an example, take the input: - - equations = ["2 c - b +3 a <= 2","c - b = 0","2 b -a >=-2"], reaction_ids = ["a","b","c"] - - This will be translated to the form A_ineq * x <= b_ineq, A_eq * x = b_eq and hence to - - A_ineq = sparse.csr_matrix([[3,-1,2],[1,-2,0]]), b_ineq = [2,2], - A_eq = sparse.csr_matrix([[1,-2,0]]), b_eq = [0] - - Args: - equations (list of str): - (List of) (in)equalities in string form equations=["r1 + 3 r2 = 0.3", "-5 r3 -r4 <= -0.5"] - - reaction_ids (list of str): - List of reaction identifiers or variable names that are used to recognize variables in - the provided (in)equalities - - Returns: - (Tuple): - A_ineq, b_ineq, A_eq, b_eq. Coefficient matrices and right hand sides that represent the input - (in)equalities as matrix-vector multiplications - """ - numr = len(reaction_ids) - A_ineq = sparse.csr_matrix((0, numr)) - b_ineq = [] - A_eq = sparse.csr_matrix((0, numr)) - b_eq = [] - for equation in equations: - try: - lhs, rhs = re.split(r"<=|=|>=", equation) - eq_sign = re.search(r"<=|>=|=", equation)[0] - rhs = float(rhs) - except: - raise Exception("Equations must contain exactly one (in)equality sign: <=,=,>=. Right hand side must be a float number.") - A = linexpr2mat(lhs, reaction_ids) - if eq_sign == "=": - A_eq = sparse.vstack((A_eq, A)) - b_eq += [rhs] - elif eq_sign == "<=": - A_ineq = sparse.vstack((A_ineq, A)) - b_ineq += [rhs] - elif eq_sign == ">=": - A_ineq = sparse.vstack((A_ineq, -A)) - b_ineq += [-rhs] - return A_ineq, b_ineq, A_eq, b_eq - - def lineq2list(equations, reaction_ids) -> List: """Translates *linear* (in)equalities to list format: [lhs,sign,rhs] @@ -177,28 +125,6 @@ def lineq2list(equations, reaction_ids) -> List: return D -def lineqlist2str(D): - """Translates a *linear* (in)equality from the list format [lhs,sign,rhs] to a string - - E.g. input: D=[{"a":3.0,"b":-1.0,"c":2.0},"<=",2.0]] is translated to: out="3.0 a - 1.0 b + 2.0 c <= 2" - - Args: - D (list): - (In)equality in list form, e.g.: D=[{"a":3.0,"b":-1.0,"c":2.0},"<=",2.0]] - - Returns: - (str): - A list of (in)equalities in string form - - """ - if D[0]: - return linexprdict2str(D[0]) + " " + D[1] + " " + str(D[2]) - elif D[1] and D[2]: - return D[1] + " " + str(D[2]) - else: - return "" - - def lineqlist2mat(D, reaction_ids) -> Tuple[sparse.csr_matrix, Tuple, sparse.csr_matrix, Tuple]: """Translates *linear* (in)equalities presented in the list of lists format to matrices @@ -402,30 +328,3 @@ def linexprdict2str(D): else: return "" - -def get_rids(expr, reaction_ids): - """Get reaction identifiers that are present in string - - E.g.: input: D={"R1":-1.0, "R3": 2.0}, translates to the string: "- 1.0 R1 + 2.0 R3" - - Args: - expr (str): - A character string - - reaction_ids (list of str): - List of reaction identifiers or variable names - - Returns: - (list of str): - A list of strings containing the reaction/variable strings present in the input string - """ - expr_parts = [re.sub(r"^(\s|-|\+|\()*|(\s|-|\+|\|<|\=|>)*$", "", part) for part in expr.split()] - reacIDs = [] - for part in expr_parts: - if part in reaction_ids: - reacIDs += [part] - continue - if re.match(r"^\d*\.{0,1}\d*$", part) is not None: - continue - raise Exception("Expression invalid. Unknown identifier " + part + ".") - return reacIDs diff --git a/tests/conftest.py b/tests/conftest.py index f1dbaa7..f4556fc 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,4 +1,3 @@ -import platform import pytest from importlib.util import find_spec from cobra import Configuration @@ -13,7 +12,6 @@ def pytest_addoption(parser): for name, help_text in [ ("--medium", "Run iMLcore genome-scale benchmarks (~4 min total)."), ("--large", "Run iML1515 large-model benchmarks (several min/solver)."), - ("--java", "Run JPype/JVM tests on Linux/macOS too (flaky, see jpype#934)."), ]: try: parser.addoption(name, action="store_true", default=False, help=help_text) @@ -25,7 +23,6 @@ def pytest_configure(config): for marker, desc in [ ("medium", "genome-scale benchmark; enable with --medium"), ("large", "large-model benchmark; enable with --large"), - ("java", "requires JPype/JVM; skipped on non-Windows CI (jpype#934)"), ]: config.addinivalue_line("markers", f"{marker}: {desc}") # Suppress known third-party warnings @@ -43,23 +40,6 @@ def pytest_collection_modifyitems(config, items): if marker in item.keywords: item.add_marker(skip) - # --------------------------------------------------------------------------- - # Platform-based skips (centralized here for visibility) - # --------------------------------------------------------------------------- - # JPype's JNI bridge crashes non-deterministically (~1-in-20) on Linux/macOS - # CI runners due to a GC finalization race (jpype#934). Windows is unaffected. - # Tested jpype1==1.5.0 pinning — no improvement (still segfaults, plus no - # Python 3.13 wheel causing build failures on macOS ARM64). - # --java forces them on anyway, which is how a Java-backend change gets verified without - # round-tripping through the Windows CI leg. - if platform.system() != 'Windows' and not config.getoption("--java", default=False): - skip_java = pytest.mark.skip( - reason="JPype JNI crashes non-deterministically on Linux/macOS (jpype#934); " - "pass --java to run anyway") - for item in items: - if "java" in item.keywords: - item.add_marker(skip_java) - cobra_conf = Configuration() bound_thres = max((abs(cobra_conf.lower_bound), abs(cobra_conf.upper_bound))) diff --git a/tests/test_04_preprocessing.py b/tests/test_04_preprocessing.py index 32893d3..a0f60d9 100644 --- a/tests/test_04_preprocessing.py +++ b/tests/test_04_preprocessing.py @@ -268,76 +268,3 @@ def test_coupled_group_r3_rpex(self, gpr_model): assert simplify_logic(result_sympy ^ expected) == False, \ f"GPR mismatch. Got: {result_sympy}, expected: {expected}" - -@pytest.mark.java -class TestEfmtoolBackendGpr: - """Test GPR propagation with efmtool_rref backend (if Java available).""" - - @pytest.fixture - def java_available(self): - try: - from straindesign.efmtool_cmp_interface import _check_jpype_available - if not _check_jpype_available(): - pytest.skip("jpype not installed") - from straindesign.efmtool_cmp_interface import _init_java - _init_java() - except Exception as e: - pytest.skip(f"Java/jpype not available: {e}") - - def test_efmtool_coupled_gpr_matches_sparse_rref(self, gpr_model, java_available): - """Both backends should produce semantically equivalent GPR rules.""" - model_java = copy.deepcopy(gpr_model) - - # Sparse RREF path - remove_blocked_reactions(gpr_model) - stoichmat_coeff_to_fraction(gpr_model) - remove_conservation_relations(gpr_model) - rref_map = compress_model_coupled(gpr_model, compression_backend='sparse_rref', propagate_gpr=True) - - # Efmtool path - remove_blocked_reactions(model_java) - stoichmat_coeff_to_fraction(model_java) - remove_conservation_relations(model_java) - java_map = compress_model_coupled(model_java, compression_backend='efmtool_rref', propagate_gpr=True) - - def gpr_by_group(model, reac_map): - result = {} - for rid, orig_map in reac_map.items(): - key = frozenset(orig_map.keys()) - rxn = model.reactions.get_by_id(rid) - result[key] = rxn.gene_reaction_rule - return result - - rref_gprs = gpr_by_group(gpr_model, rref_map) - java_gprs = gpr_by_group(model_java, java_map) - - assert rref_gprs.keys() == java_gprs.keys(), \ - f"Different groups: rref={rref_gprs.keys()}, java={java_gprs.keys()}" - - for group_key in rref_gprs: - gpr_rref = rref_gprs[group_key] - gpr_java = java_gprs[group_key] - - if not gpr_rref and not gpr_java: - continue - - from cobra.core.gene import GPR - sym_rref = GPR.from_string(gpr_rref).as_symbolic() - sym_java = GPR.from_string(gpr_java).as_symbolic() - assert simplify_logic(sym_rref ^ sym_java) == False, \ - f"GPR mismatch for group {sorted(group_key)}: rref='{gpr_rref}', java='{gpr_java}'" - - def test_efmtool_full_compression_gpr(self, gpr_model, java_available): - """Full compression with efmtool backend should preserve gene info.""" - orig_genes = set() - for r in gpr_model.reactions: - orig_genes.update(g.id for g in r.genes) - - compress_model(gpr_model, compression_backend='efmtool_rref', propagate_gpr=True) - - compressed_genes = set() - for r in gpr_model.reactions: - compressed_genes.update(g.id for g in r.genes) - - assert compressed_genes <= orig_genes - assert len(compressed_genes) > 0, "All gene information was lost" diff --git a/tests/test_07_compression.py b/tests/test_07_compression.py index 2397847..a60961b 100644 --- a/tests/test_07_compression.py +++ b/tests/test_07_compression.py @@ -1,4 +1,4 @@ -"""Compression tests: unit tests, compression_backend parity, FVA equivalence, and MCS validation.""" +"""Compression tests: unit tests, map correctness, FVA equivalence, and MCS validation.""" import sys import pytest import numpy as np @@ -33,31 +33,10 @@ def model_small_example(): # ============================================================================= -# Unit tests (Python-only, no Java required) +# Unit tests # ============================================================================= -def test_no_jpype_loaded(): - """Verify that jpype is not loaded when importing straindesign.""" - jpype_before = [m for m in sys.modules if 'jpype' in m.lower()] - # Save and remove straindesign modules to test a fresh import - saved_modules = {m: sys.modules[m] for m in list(sys.modules) if m.startswith('straindesign')} - for m in saved_modules: - del sys.modules[m] - try: - import straindesign as sd_fresh - jpype_after = [m for m in sys.modules if 'jpype' in m.lower()] - new_jpype = set(jpype_after) - set(jpype_before) - assert len(new_jpype) == 0, f"straindesign loaded jpype modules: {new_jpype}" - finally: - # Restore original modules so function identity is preserved for - # multiprocessing pickle (SDPool serialises fva_worker_init by reference). - for m in list(sys.modules): - if m.startswith('straindesign'): - del sys.modules[m] - sys.modules.update(saved_modules) - - def test_python_compression_basic(model_gpr): """Compression reduces reaction count and returns a non-empty map.""" sd.extend_model_gpr(model_gpr, use_names=False) @@ -68,10 +47,10 @@ def test_python_compression_basic(model_gpr): def test_python_compression_coupled_function(model_small_example): - """compress_model_coupled with compression_backend='sparse_rref' returns a dict.""" + """compress_model_coupled returns a dict.""" nt.stoichmat_coeff_to_fraction(model_small_example) nt.remove_conservation_relations(model_small_example) - reac_map = nt.compress_model_coupled(model_small_example, compression_backend='sparse_rref') + reac_map = nt.compress_model_coupled(model_small_example) assert isinstance(reac_map, dict) @@ -79,7 +58,7 @@ def test_compression_coefficient_type(model_small_example): """Compression coefficients are exact rational number types.""" nt.stoichmat_coeff_to_fraction(model_small_example) nt.remove_conservation_relations(model_small_example) - reac_map = nt.compress_model_coupled(model_small_example, compression_backend='sparse_rref') + reac_map = nt.compress_model_coupled(model_small_example) for new_reac, old_reacs in reac_map.items(): for old_reac, coeff in old_reacs.items(): assert is_rational_type(coeff), (f"Coefficient for {old_reac} in {new_reac}: expected rational, got {type(coeff)}") @@ -93,11 +72,11 @@ def test_stoichmat_coeff_to_fraction_uses_rational_type(model_small_example): assert is_rational_type(coeff), (f"Coefficient for {metabolite.id} in {reaction.id}: expected rational, got {type(coeff)}") -def test_basic_columns_rat_python(): - """basic_columns_rat returns correct pivot count for a rank-2 matrix.""" - import straindesign.efmtool_cmp_interface as efm +def test_basic_columns_from_numpy(): + """basic_columns_from_numpy returns correct pivot count for a rank-2 matrix.""" + from straindesign.compression import basic_columns_from_numpy mx = np.array([[1.0, 0.0, 1.0], [0.0, 1.0, 1.0], [1.0, 1.0, 2.0]]) - basic_cols = efm.basic_columns_rat(mx) + basic_cols = basic_columns_from_numpy(mx) assert len(basic_cols) == 2, f"Expected 2 basic columns, got {len(basic_cols)}" @@ -119,48 +98,26 @@ def test_compression_preserves_flux_space(model_small_example): @pytest.mark.timeout(30) -def test_full_strain_design_without_java(model_gpr): - """A full strain design computation completes using only the sparse backend.""" +def test_full_strain_design_compressed(model_gpr): + """A full strain design computation completes with compression enabled.""" from straindesign.names import SUPPRESS, ANY sd.extend_model_gpr(model_gpr, use_names=False) module = sd.SDModule(model_gpr, module_type=SUPPRESS, constraints='r_bm >= 0.1') - try: - sd.compute_strain_designs( - model_gpr, - sd_modules=[module], - max_solutions=1, - max_cost=2, - compress=True, - solution_approach=ANY, - ) - except ImportError as e: - if 'jpype' in str(e).lower(): - pytest.fail(f"Java/jpype was required but should not be: {e}") - raise + sd.compute_strain_designs( + model_gpr, + sd_modules=[module], + max_solutions=1, + max_cost=2, + compress=True, + solution_approach=ANY, + ) # ============================================================================= -# Backend parity tests (Java required; skipped if jpype unavailable) +# Compression correctness through the reaction map # ============================================================================= -@pytest.fixture -def jpype_available(): - jpype = pytest.importorskip("jpype", reason="jpype not installed; skipping Java parity tests") - return jpype - - -@pytest.mark.java -def test_compression_parity_reaction_count(jpype_available): - """Both compression backends compress e_coli_core to the same number of reactions.""" - model_py = load_model("e_coli_core") - nt.compress_model(model_py, compression_backend='sparse_rref') - model_java = load_model("e_coli_core") - nt.compress_model(model_java, compression_backend='efmtool_rref') - assert len(model_py.reactions) == len( - model_java.reactions), (f"Reaction count mismatch: sparse_rref={len(model_py.reactions)}, efmtool_rref={len(model_java.reactions)}") - - def _trace_lump(cmp_maps, orig_id): """Follow an original reaction through the compression rounds. @@ -176,36 +133,28 @@ def _trace_lump(cmp_maps, orig_id): return cur, factor -@pytest.mark.java -def test_fba_equivalence(jpype_available): - """Both backends preserve the uncompressed optimum once the lump factor is applied. +def test_fba_optimum_recovered_through_map(): + """Compression preserves the uncompressed optimum once the lump factor is applied. A lump's overall scale is free: only its ratios are fixed, so the raw objective value of a - lumped reaction is backend-specific and not a meaningful thing to compare. sparse_rref - re-expresses each lump in one member's units, efmtool_rref does not, so their biomass columns - differ by a constant factor. What must agree -- and what a caller actually relies on -- is the - flux recovered through the compression map. + lumped reaction is not meaningful on its own. What a caller relies on is the flux recovered + through the compression map, which must reproduce the uncompressed optimum exactly. This also + exercises the map itself -- it would catch factors drifting out of step with the column + scaling applied when a lump is re-expressed in one member's units. """ base = load_model("e_coli_core") biomass = next((r.id for r in base.reactions if 'biomass' in r.id.lower()), None) assert biomass, "Could not find biomass reaction" ref = sd.fba(base, obj={biomass: 1}, obj_sense='maximize').objective_value - recovered = {} - for backend in ('sparse_rref', 'efmtool_rref'): - model = load_model("e_coli_core") - cmp_maps = nt.compress_model(model, compression_backend=backend) - cmp_id, factor = _trace_lump(cmp_maps, biomass) - assert cmp_id in [r.id for r in model.reactions], ( - f"{backend}: compression map names {cmp_id}, which is not in the compressed model") - val = sd.fba(model, obj={cmp_id: 1}, obj_sense='maximize').objective_value - recovered[backend] = factor * val - - for backend, val in recovered.items(): - assert abs(val - ref) < 1e-6, ( - f"{backend}: recovered optimum {val} != uncompressed {ref}") - assert abs(recovered['sparse_rref'] - recovered['efmtool_rref']) < 1e-6, ( - f"Backend mismatch after mapping back: {recovered}") + model = load_model("e_coli_core") + cmp_maps = nt.compress_model(model) + cmp_id, factor = _trace_lump(cmp_maps, biomass) + assert cmp_id in [r.id for r in model.reactions], ( + f"compression map names {cmp_id}, which is not in the compressed model") + val = sd.fba(model, obj={cmp_id: 1}, obj_sense='maximize').objective_value + assert abs(factor * val - ref) < 1e-6, ( + f"recovered optimum {factor * val} != uncompressed {ref}") def test_cobra_optimize_after_compression(): @@ -223,7 +172,7 @@ def test_cobra_optimize_after_compression(): val_orig = model_orig.optimize().objective_value model_cmp = load_model("e_coli_core") - cmp_map = nt.compress_model(model_cmp, compression_backend='sparse_rref') + cmp_map = nt.compress_model(model_cmp) # Find biomass in compressed model via compression map biomass_cmp_id = None @@ -243,34 +192,6 @@ def test_cobra_optimize_after_compression(): f"Expanded objective mismatch: original={val_orig}, expanded={val_expanded}") -@pytest.mark.java -def test_fva_equivalence(jpype_available): - """Both compression backends produce flux spaces with no true FVA mismatches. - - Sign-convention differences (Python = -Java for some lumped reactions) are - mathematically equivalent and are not counted as mismatches. - """ - model_py = load_model("e_coli_core") - nt.compress_model(model_py, compression_backend='sparse_rref') - model_java = load_model("e_coli_core") - nt.compress_model(model_java, compression_backend='efmtool_rref') - - fva_py = flux_variability_analysis(model_py, fraction_of_optimum=0.0, processes=1) - fva_java = flux_variability_analysis(model_java, fraction_of_optimum=0.0, processes=1) - - common = set(fva_py.index) & set(fva_java.index) - true_mismatches = [] - for r_id in common: - py_min, py_max = fva_py.loc[r_id, 'minimum'], fva_py.loc[r_id, 'maximum'] - java_min, java_max = fva_java.loc[r_id, 'minimum'], fva_java.loc[r_id, 'maximum'] - direct = abs(py_min - java_min) < 1e-6 and abs(py_max - java_max) < 1e-6 - flipped = abs(py_min - (-java_max)) < 1e-6 and abs(py_max - (-java_min)) < 1e-6 - if not direct and not flipped: - true_mismatches.append(r_id) - - assert len(true_mismatches) == 0, (f"True FVA mismatches between sparse_rref and efmtool_rref backends: {true_mismatches}") - - # ============================================================================= # FVA back-mapping test (sparse only) # ============================================================================= @@ -283,7 +204,7 @@ def test_fva_expansion(): fva_orig = flux_variability_analysis(model_orig, fraction_of_optimum=0.0, processes=1) model_cmp = load_model("e_coli_core") - cmp_map = nt.compress_model(model_cmp, compression_backend='sparse_rref') + cmp_map = nt.compress_model(model_cmp) fva_cmp = flux_variability_analysis(model_cmp, fraction_of_optimum=0.0, processes=1) # Build inverse map: orig_id -> (compressed_id, coefficient) @@ -317,21 +238,12 @@ def test_fva_expansion(): # ============================================================================= -@pytest.mark.parametrize("compression_backend", [ - "sparse_rref", - pytest.param("efmtool_rref", marks=pytest.mark.java), -]) -def test_mcs_e_coli_core(compression_backend): +def test_mcs_e_coli_core(): """MCS computation on e_coli_core returns the expected 455 solutions. - Parametrized over both compression backends so regressions in either - are caught. The efmtool_rref variant is skipped when jpype is not installed. - Requires a strong MILP solver (Gurobi, CPLEX, or SCIP). GLPK cannot reliably enumerate all solutions via POPULATE and is excluded. """ - if compression_backend == "efmtool_rref": - pytest.importorskip("jpype", reason="jpype not installed; skipping efmtool backend") from straindesign.names import SUPPRESS, POPULATE, GLPK, SCIP, GUROBI, CPLEX # Solver priority: SCIP (no size limit) > CPLEX > GUROBI (both have free-tier limits) strong_solvers = sd.avail_solvers - {GLPK} @@ -346,51 +258,7 @@ def test_mcs_e_coli_core(compression_backend): max_cost=3, gene_kos=True, solver=solver, - compression_backend=compression_backend) +) assert len(sols.reaction_sd) == 455, ( - f"Expected 455 MCS for e_coli_core (compression_backend={compression_backend}), got {len(sols.reaction_sd)}") - + f"Expected 455 MCS for e_coli_core, got {len(sols.reaction_sd)}") -@pytest.mark.timeout(300) -@pytest.mark.large -def test_imlcore_compression_parity(jpype_available): - """Both compression backends compress iMLcore to the same number of reactions. - - Marked --large: JPype's JNI bridge crashes (SIGBUS/SIGSEGV) on GitHub Actions - runners when processing iMLcore-sized matrices through the Java RREF. - The e_coli_core parity tests above cover the same code path on a smaller matrix. - Run locally with: pytest --large -k test_imlcore_compression_parity - """ - model_py = read_sbml_model(dirname(abspath(__file__)) + r"/iMLcore.xml") - nt.compress_model(model_py, compression_backend='sparse_rref') - model_java = read_sbml_model(dirname(abspath(__file__)) + r"/iMLcore.xml") - nt.compress_model(model_java, compression_backend='efmtool_rref') - assert len(model_py.reactions) == len( - model_java.reactions), (f"iMLcore reaction count mismatch: sparse_rref={len(model_py.reactions)}, efmtool_rref={len(model_java.reactions)}") - - -@pytest.mark.timeout(600) -@pytest.mark.large -def test_mcs_imlcore_parity(jpype_available): - """MCS on iMLcore returns the same solutions with both compression backends. - - Marked --large: JPype's JNI bridge crashes on CI runners with iMLcore-sized - matrices (see test_imlcore_compression_parity docstring). - Run locally with: pytest --large -k test_mcs_imlcore_parity - """ - from straindesign.names import SUPPRESS, POPULATE, GLPK - strong_solvers = sd.avail_solvers - {GLPK} - if not strong_solvers: - pytest.skip("iMLcore MCS parity test requires Gurobi, CPLEX, or SCIP") - solver = next(iter(strong_solvers)) - results = {} - for backend in ['sparse_rref', 'efmtool_rref']: - model = read_sbml_model(dirname(abspath(__file__)) + r"/iMLcore.xml") - modules = [sd.SDModule(model, SUPPRESS, - constraints='BIOMASS_Ec_iML1515_core_75p37M >= 0.001')] - sols = sd.compute_strain_designs(model, sd_modules=modules, solution_approach=POPULATE, - max_cost=3, solver=solver, - compression_backend=backend) - results[backend] = len(sols.reaction_sd) - assert results['sparse_rref'] == results['efmtool_rref'], ( - f"iMLcore MCS count mismatch: sparse_rref={results['sparse_rref']}, efmtool_rref={results['efmtool_rref']}")