From 84be8f3473d95a544e01ba754d956cefe5bc49b8 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 08:40:07 +0100 Subject: [PATCH 1/6] fix: read_yaml_model reads empty enzymes omaps and doubled-quote scalars An ec reaction without enzymes is written as '- enzymes: !!omap' with no content; the reader now gives it no enzymes, as RAVEN's readYAMLmodel does. A scalar written as ""text"" is repaired line by line only after the strict parse fails, and reads as "text". --- CHANGELOG.md | 2 + src/raven_toolbox/io/yaml.py | 68 ++++++++++++++++++++++++++--- tests/test_io_yaml_ec_data.py | 80 +++++++++++++++++++++++++++++++++++ 3 files changed, 143 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b629dd22..ff60a30d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,8 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased +* `io.read_yaml_model` reads an ec reaction written as `- enzymes: !!omap` with nothing under the tag, the form GECKO and RAVEN write for a reaction without enzymes, as a reaction with no enzymes, as `readYAMLmodel` does. It raised `AttributeError: 'ScalarNode' object has no attribute 'flow_style'`. +* `io.read_yaml_model` reads a scalar written as `""text""` (invalid YAML, written by older RAVEN and GECKO for a note that carried its own quotes) as `"text"`, as `readYAMLmodel` does. Files that parse as YAML are read unchanged. * `init.ftinit` and `init.run_ftinit` give each scored reaction its own big-M, its flux capacity `max(min(ub, big_m), min(-lb, big_m), force_on)` with an infinite bound counted as `big_m`, as RAVEN's `ftINITInternalAlg` does. `big_m` is the cap on that capacity; its diff --git a/src/raven_toolbox/io/yaml.py b/src/raven_toolbox/io/yaml.py index bf5b52bb..011d4180 100644 --- a/src/raven_toolbox/io/yaml.py +++ b/src/raven_toolbox/io/yaml.py @@ -49,8 +49,10 @@ import cobra from cobra.io.dict import model_from_dict, model_to_dict -from cobra.io.yaml import yaml as _cobra_yaml # ruamel round-trip YAML (handles !!omap) from ruamel.yaml import YAML +from ruamel.yaml.constructor import RoundTripConstructor +from ruamel.yaml.error import YAMLError +from ruamel.yaml.nodes import ScalarNode from ruamel.yaml.scalarstring import DoubleQuotedScalarString from raven_toolbox.io.ec_data import ( @@ -59,11 +61,10 @@ ec_data_to_yaml_sections, ) -# A dedicated instance for *writing* only (reading still goes through -# _cobra_yaml, whose parsing is unaffected by any of this). Mutating -# _cobra_yaml's own width would also change cobra.io.save_yaml_model for -# the rest of the process, so this module gets its own ruamel round-trip -# instance instead. `width` is set far past any real line length so +# A dedicated instance for *writing* only (reading uses _read_yaml below). +# Mutating cobra's own YAML instance would change its width for +# cobra.io.save_yaml_model in the rest of the process, so this module gets its own +# ruamel round-trip instance instead. `width` is set far past any real line length so # ruamel never folds a long scalar onto a continuation line, matching # writeYAMLmodel.m. `indent` is pinned explicitly (to the values ruamel # already used implicitly, verified byte-identical) rather than left for @@ -83,6 +84,59 @@ _write_yaml.indent(mapping=2, sequence=2, offset=0) +class _OmapConstructor(RoundTripConstructor): + """Round-trip constructor that reads a content-less ``!!omap`` as an empty map. + + GECKO and RAVEN write an ec reaction without enzymes as ``- enzymes: !!omap`` + with nothing under the tag. YAML parses that as a null scalar carrying the omap + tag, which ruamel's own constructor rejects; RAVEN's ``readYAMLmodel.m`` reads + it as an empty enzyme list. + """ + + def construct_yaml_omap(self, node): + if isinstance(node, ScalarNode) and not node.value: + yield OrderedDict() + return + yield from super().construct_yaml_omap(node) + + +_OmapConstructor.add_constructor("tag:yaml.org,2002:omap", _OmapConstructor.construct_yaml_omap) + +# A dedicated instance for *reading*, so the constructor above does not change +# cobra.io.load_yaml_model for the rest of the process. +_read_yaml = YAML(typ="rt") +_read_yaml.Constructor = _OmapConstructor + +# A scalar written as `""text""` (two quotes on each side), the shape older RAVEN +# and GECKO writers gave to a note that already carried its own quotes. It is not +# valid YAML. `readYAMLmodel.m` strips one pair of quotes per line and unescapes +# `\\` and `\"` in what is left, so the value is `"text"`. +_DOUBLED_QUOTES = re.compile( + r"""^(?P\s*-\s+|\s*(?:-\s+)?[^\s:\-][^:]*:\s+)"(?P".*")"\s*$""" +) + + +def _fix_doubled_quotes(text: str) -> str: + """Rewrite every ``""text""`` scalar as the single-quoted scalar ``'"text"'``.""" + + def fix(match: re.Match) -> str: + inner = re.sub(r'\\(["\\])', r"\1", match["body"]) + return f"{match['head']}'{inner.replace(chr(39), chr(39) * 2)}'" + + return "\n".join(_DOUBLED_QUOTES.sub(fix, line) for line in text.split("\n")) + + +def _load_yaml_text(text: str): + """Parse ``text`` strictly; when that fails, retry with doubled quotes repaired.""" + try: + return _read_yaml.load(text) + except YAMLError: + repaired = _fix_doubled_quotes(text) + if repaired == text: + raise + return _read_yaml.load(repaired) + + def _open_text(path: str | Path, mode: str): """Open ``path`` as a text handle, transparently gzipping when it ends ``.gz``.""" if str(path).endswith(".gz"): @@ -166,7 +220,7 @@ def read_yaml_model(path: str | Path) -> cobra.Model: latter is merged into one mapping before parsing. """ with _open_text(path, "r") as handle: - raw = _to_plain(_cobra_yaml.load(handle)) + raw = _to_plain(_load_yaml_text(handle.read())) if isinstance(raw, list): # Very old RAVEN files: a sequence of one-key mappings instead of diff --git a/tests/test_io_yaml_ec_data.py b/tests/test_io_yaml_ec_data.py index f236b851..4295927e 100644 --- a/tests/test_io_yaml_ec_data.py +++ b/tests/test_io_yaml_ec_data.py @@ -432,3 +432,83 @@ def test_validate_catches_coupling_matrix_shape_drift(): def test_empty_gecko_light_flag_propagates(): assert EcData.empty(1, 1, gecko_light=True).gecko_light is True assert EcData.empty(1, 1).gecko_light is False + + +# --------------------------------------------------------------------------- # +# Files as GECKO and RAVEN write them +# --------------------------------------------------------------------------- # + +_GECKO_FILE = """\ +--- +!!omap +- metabolites: + - !!omap + - id: "a" + - name: "A" + - compartment: "c" +- reactions: + - !!omap + - id: "R1" + - lower_bound: 0 + - upper_bound: 1000 + - metabolites: !!omap + - a: -1 +- genes: [] +- compartments: !!omap + - c: "cytoplasm" +- ec-rxns: + - !!omap + - id: "std" + - kcat: 10 + - source: "custom" + - notes: {note} + - enzymes: !!omap{empty} + - !!omap + - id: "enz" + - kcat: 5 + - enzymes: !!omap + - P1: 1 +- ec-enzymes: + - !!omap + - genes: "G1" + - enzymes: "P1" + - mw: 30000 + - sequence: "MAGIC" +""" + + +@pytest.mark.parametrize("empty", ["", " []"]) +def test_load_reaction_with_empty_enzymes_omap(tmp_path, empty): + # GECKO and RAVEN write `- enzymes: !!omap` with nothing under it for an ec + # reaction without enzymes; raven-toolbox writes `!!omap []`. Both read as no enzymes. + path = tmp_path / "m.yml" + path.write_text(_GECKO_FILE.format(note='"ok"', empty=empty), encoding="utf-8") + model = read_yaml_model(path) + assert model.ec.rxns == ["std", "enz"] + assert model.ec.enzymes == ["P1"] + assert model.ec.rxn_enz_mat.toarray().tolist() == [[0.0], [1.0]] + + +def test_load_notes_in_doubled_quotes_keep_one_pair(tmp_path): + # An older writer put quotes around a note that carried its own, giving `""text""`, + # which is not valid YAML. readYAMLmodel.m strips one pair and keeps the other. + path = tmp_path / "m.yml" + path.write_text(_GECKO_FILE.format(note='""Putative: lowest kcat.""', empty=""), encoding="utf-8") + model = read_yaml_model(path) + assert model.ec.notes[0] == '"Putative: lowest kcat."' + assert model.ec.rxns == ["std", "enz"] + + +def test_doubled_quotes_unescape_like_matlab_and_keep_apostrophes(tmp_path): + path = tmp_path / "m.yml" + note = '""it\'s a \\\\ and \\"inner\\" one""' + path.write_text(_GECKO_FILE.format(note=note, empty=""), encoding="utf-8") + model = read_yaml_model(path) + assert model.ec.notes[0] == '"it\'s a \\ and "inner" one"' + + +def test_invalid_yaml_other_than_doubled_quotes_still_raises(tmp_path): + path = tmp_path / "m.yml" + path.write_text(_GECKO_FILE.format(note='"ok"', empty="").replace("kcat: 10", "kcat: [10"), encoding="utf-8") + with pytest.raises(Exception): + read_yaml_model(path) From 9ae970da873fceced166b05995d734f8dd156785 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 13:55:17 +0100 Subject: [PATCH 2/6] fix: read_yaml_model keeps an unescaped quote inside a double-quoted value The fallback after a failed strict parse now repairs any double-quoted value with an unescaped interior quote, not only one that starts and ends with a doubled quote. The outer pair is dropped and \\ and \" are unescaped, as readYAMLmodel does. --- CHANGELOG.md | 1 + src/raven_toolbox/io/yaml.py | 43 ++++++++++++++++++++++++++--------- tests/test_io_yaml_ec_data.py | 30 ++++++++++++++++++++++++ 3 files changed, 63 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ff60a30d..ed08ec15 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,7 @@ Milestones in the raven-toolbox port. For function-level status see * `io.read_yaml_model` reads an ec reaction written as `- enzymes: !!omap` with nothing under the tag, the form GECKO and RAVEN write for a reaction without enzymes, as a reaction with no enzymes, as `readYAMLmodel` does. It raised `AttributeError: 'ScalarNode' object has no attribute 'flow_style'`. * `io.read_yaml_model` reads a scalar written as `""text""` (invalid YAML, written by older RAVEN and GECKO for a note that carried its own quotes) as `"text"`, as `readYAMLmodel` does. Files that parse as YAML are read unchanged. +* `io.read_yaml_model` reads a double-quoted scalar with an unescaped quote inside, such as `- notes: "a "b" c"`, as `a "b" c`: it drops the outer pair of quotes and unescapes `\\` and `\"`, as `readYAMLmodel` does. This repair runs only when the strict parse fails. * `init.ftinit` and `init.run_ftinit` give each scored reaction its own big-M, its flux capacity `max(min(ub, big_m), min(-lb, big_m), force_on)` with an infinite bound counted as `big_m`, as RAVEN's `ftINITInternalAlg` does. `big_m` is the cap on that capacity; its diff --git a/src/raven_toolbox/io/yaml.py b/src/raven_toolbox/io/yaml.py index 011d4180..72627e54 100644 --- a/src/raven_toolbox/io/yaml.py +++ b/src/raven_toolbox/io/yaml.py @@ -107,27 +107,48 @@ def construct_yaml_omap(self, node): _read_yaml = YAML(typ="rt") _read_yaml.Constructor = _OmapConstructor -# A scalar written as `""text""` (two quotes on each side), the shape older RAVEN -# and GECKO writers gave to a note that already carried its own quotes. It is not -# valid YAML. `readYAMLmodel.m` strips one pair of quotes per line and unescapes -# `\\` and `\"` in what is left, so the value is `"text"`. -_DOUBLED_QUOTES = re.compile( - r"""^(?P\s*-\s+|\s*(?:-\s+)?[^\s:\-][^:]*:\s+)"(?P".*")"\s*$""" +# A double-quoted scalar that holds a double quote no backslash escapes, such as +# `""text""` (the shape older RAVEN and GECKO writers gave to a note that already +# carried its own quotes) or `"a "b" c"`. It is not valid YAML. `readYAMLmodel.m` strips +# the outer pair of quotes of the line and unescapes `\\` and `\"` in what is left, so the +# values are `"text"` and `a "b" c`. +_QUOTED_VALUE = re.compile( + r"""^(?P\s*-\s+|\s*(?:-\s+)?[^\s:\-][^:]*:\s+)"(?P.*)"\s*$""" ) +def _has_unescaped_quote(body: str) -> bool: + """Whether ``body``, the text between the outer quotes, has a ``"`` that no backslash escapes.""" + i = 0 + while i < len(body): + if body[i] == "\\": + i += 2 + elif body[i] == '"': + return True + else: + i += 1 + return False + + def _fix_doubled_quotes(text: str) -> str: - """Rewrite every ``""text""`` scalar as the single-quoted scalar ``'"text"'``.""" + """Rewrite each double-quoted scalar with an unescaped interior quote as a single-quoted scalar. + + The value is the text between the outer quotes with escaped backslashes and quotes + unescaped, as ``readYAMLmodel.m`` reads it; quoted scalars without such a quote are left alone. + """ def fix(match: re.Match) -> str: - inner = re.sub(r'\\(["\\])', r"\1", match["body"]) - return f"{match['head']}'{inner.replace(chr(39), chr(39) * 2)}'" + body = match["body"] + if not _has_unescaped_quote(body): + return match[0] + value = re.sub(r'\\(["\\])', r"\1", body) + return f"{match['head']}'{value.replace(chr(39), chr(39) * 2)}'" - return "\n".join(_DOUBLED_QUOTES.sub(fix, line) for line in text.split("\n")) + return "\n".join(_QUOTED_VALUE.sub(fix, line) for line in text.split("\n")) def _load_yaml_text(text: str): - """Parse ``text`` strictly; when that fails, retry with doubled quotes repaired.""" + """Parse ``text`` strictly; when that fails, retry with quotes inside double-quoted scalars repaired.""" try: return _read_yaml.load(text) except YAMLError: diff --git a/tests/test_io_yaml_ec_data.py b/tests/test_io_yaml_ec_data.py index 4295927e..bf494e03 100644 --- a/tests/test_io_yaml_ec_data.py +++ b/tests/test_io_yaml_ec_data.py @@ -507,6 +507,36 @@ def test_doubled_quotes_unescape_like_matlab_and_keep_apostrophes(tmp_path): assert model.ec.notes[0] == '"it\'s a \\ and "inner" one"' +@pytest.mark.parametrize( + "note, expected", + [ + ('"a "b" c"', 'a "b" c'), + ('"say "hi""', 'say "hi"'), + ('"a \\"b\\" and "c" end"', 'a "b" and "c" end'), + ('"x \\\\ "y" it\'s"', 'x \\ "y" it\'s'), + ('"trailing "q""', 'trailing "q"'), + ], +) +def test_interior_unescaped_quote_is_kept_like_matlab(tmp_path, note, expected): + # readYAMLmodel.m strips the outer pair of quotes and keeps every quote inside. + path = tmp_path / "m.yml" + path.write_text(_GECKO_FILE.format(note=note, empty=""), encoding="utf-8") + model = read_yaml_model(path) + assert model.ec.notes[0] == expected + assert model.ec.rxns == ["std", "enz"] + + +def test_valid_quoted_scalars_are_read_unchanged_next_to_a_repaired_one(tmp_path): + path = tmp_path / "m.yml" + text = _GECKO_FILE.format(note='"a "b" c"', empty="").replace( + '- source: "custom"', '- source: "cu\\\\st\\"om"' + ).replace("- id: \"enz\"", "- id: \"enz\"\n - notes: \"fine\\n\"") + path.write_text(text, encoding="utf-8") + model = read_yaml_model(path) + assert model.ec.notes == ['a "b" c', "fine\n"] + assert model.ec.source[0] == 'cu\\st"om' + + def test_invalid_yaml_other_than_doubled_quotes_still_raises(tmp_path): path = tmp_path / "m.yml" path.write_text(_GECKO_FILE.format(note='"ok"', empty="").replace("kcat: 10", "kcat: [10"), encoding="utf-8") From a0d6e087389d0279773e348f893932b128a8ced4 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 14:25:54 +0100 Subject: [PATCH 3/6] fix: fseof holds the objective at 0.999999 of its optimum in the pFBA step RAVEN's solveLP minimises the flux sum with the objective fixed at 0.999999 of the optimum it just found. fseof demanded the full value, which the HiGHS solver of the scipy interface rejects at the top level of ecYeastGEM because its optimum is accurate only to the feasibility tolerance. All enforced levels are now returned on scipy, as on Gurobi and GLPK. --- CHANGELOG.md | 1 + src/raven_toolbox/analysis/fseof.py | 13 ++++++++++++- tests/test_analysis_fseof.py | 27 +++++++++++++++++++++++++++ 3 files changed, 40 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ed08ec15..637ed9b5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased +* `analysis.fseof` (and so `geckopy.ec_fseof`) solves each level with pFBA that holds the objective at 0.999999 of its optimum at that level, as RAVEN's `solveLP` does, where it held the full value. On the scipy solver (HiGHS) the optimum is only accurate to the solver's feasibility tolerance, so the second step was infeasible at the top level of ecYeastGEM and the scan returned 9 of 10 levels; it now returns all 10, as on Gurobi and GLPK. A level at which the biomass maximisation is infeasible still ends the scan. * `io.read_yaml_model` reads an ec reaction written as `- enzymes: !!omap` with nothing under the tag, the form GECKO and RAVEN write for a reaction without enzymes, as a reaction with no enzymes, as `readYAMLmodel` does. It raised `AttributeError: 'ScalarNode' object has no attribute 'flow_style'`. * `io.read_yaml_model` reads a scalar written as `""text""` (invalid YAML, written by older RAVEN and GECKO for a note that carried its own quotes) as `"text"`, as `readYAMLmodel` does. Files that parse as YAML are read unchanged. * `io.read_yaml_model` reads a double-quoted scalar with an unescaped quote inside, such as `- notes: "a "b" c"`, as `a "b" c`: it drops the outer pair of quotes and unescapes `\\` and `\"`, as `readYAMLmodel` does. This repair runs only when the strict parse fails. diff --git a/src/raven_toolbox/analysis/fseof.py b/src/raven_toolbox/analysis/fseof.py index b3e8f605..72305d92 100644 --- a/src/raven_toolbox/analysis/fseof.py +++ b/src/raven_toolbox/analysis/fseof.py @@ -27,6 +27,13 @@ from scipy.stats import linregress +# pFBA keeps the objective at this fraction of the optimum found at the same level, as +# RAVEN's solveLP does when it minimises the flux sum (``solution.f*0.999999``). A solver +# reports its optimum to within its feasibility tolerance, so a second solve that demands +# the full value can be infeasible on that solver (HiGHS at the top level of ecYeastGEM). +_PFBA_FRACTION = 0.999999 + + @dataclass class FSEOFResult: """FSEOF output. @@ -102,6 +109,10 @@ def fseof( it is ``knockout`` when the slope is negative and the last absolute flux is at most 1e-8, ``knockdown`` for a negative slope otherwise. + Each level is solved in two steps as in RAVEN: the objective is maximised, then the + flux sum is minimised with the objective held at 0.999999 of its optimum at that + level. A level at which the first step is infeasible ends the scan. + Differences from ``FSEOF.m``: the scan stops at the first infeasible enforced level, so the last absolute flux used for the knockout label is that of the last feasible level; ``targets`` is sorted by ``|slope| * |r|`` and not grouped by label; and @@ -133,7 +144,7 @@ def fseof( model.objective = biomass_rxn model.reactions.get_by_id(target_rxn).lower_bound = level try: - columns[level] = pfba(model).fluxes + columns[level] = pfba(model, fraction_of_optimum=_PFBA_FRACTION).fluxes except OptimizationError: break # enforced flux became infeasible: stop scanning enforced.append(level) diff --git a/tests/test_analysis_fseof.py b/tests/test_analysis_fseof.py index 0ee9b3de..20deadee 100644 --- a/tests/test_analysis_fseof.py +++ b/tests/test_analysis_fseof.py @@ -171,3 +171,30 @@ def test_knockout_threshold_is_1e_minus_8(): table = _classify(m, scan, enforced, corr_threshold=0.5).set_index("reaction") assert table.loc["r_down", "target_type"] == "knockdown" assert table.loc["r_out", "target_type"] == "knockout" + + +def test_pfba_holds_objective_at_matlab_fraction(model, monkeypatch): + # RAVEN's solveLP minimises the flux sum with the objective fixed at 0.999999 of its optimum, + # so a solver whose optimum overshoots by its tolerance still has a feasible second step. + import importlib + + module = importlib.import_module("raven_toolbox.analysis.fseof") + + fractions = [] + real = module.pfba + + def spy(m, **kwargs): + fractions.append(kwargs.get("fraction_of_optimum")) + return real(m, **kwargs) + + monkeypatch.setattr(module, "pfba", spy) + res = fseof(model, "EX_P", n_steps=4) + assert len(fractions) >= len(res.enforced) >= 2 + assert set(fractions) == {0.999999} + + +def test_level_with_infeasible_maximisation_ends_the_scan(model): + # The theoretical maximum of EX_P is 10; levels above it cannot be enforced and are not + # reported, however the second pFBA step is tolerated. + res = fseof(model, "EX_P", n_steps=4, max_fraction=2.0) + assert res.enforced == pytest.approx([5.0, 10.0]) From 5caa922511c3bdbe555e3395a6aa908e1b7431f4 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 18:10:21 +0100 Subject: [PATCH 4/6] fix: read_yaml_model keeps an apostrophe inside a single-quoted value The fallback after a failed strict parse now also repairs a single-quoted value with an apostrophe that no second apostrophe pairs with, such as 'it's'. The outer pair is dropped and each '' is read as ', as readYAMLmodel does. --- CHANGELOG.md | 1 + src/raven_toolbox/io/yaml.py | 42 ++++++++++++++++++++++++++++++----- tests/test_io_yaml_ec_data.py | 28 +++++++++++++++++++++++ 3 files changed, 66 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 637ed9b5..878495ba 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ Milestones in the raven-toolbox port. For function-level status see * `io.read_yaml_model` reads an ec reaction written as `- enzymes: !!omap` with nothing under the tag, the form GECKO and RAVEN write for a reaction without enzymes, as a reaction with no enzymes, as `readYAMLmodel` does. It raised `AttributeError: 'ScalarNode' object has no attribute 'flow_style'`. * `io.read_yaml_model` reads a scalar written as `""text""` (invalid YAML, written by older RAVEN and GECKO for a note that carried its own quotes) as `"text"`, as `readYAMLmodel` does. Files that parse as YAML are read unchanged. * `io.read_yaml_model` reads a double-quoted scalar with an unescaped quote inside, such as `- notes: "a "b" c"`, as `a "b" c`: it drops the outer pair of quotes and unescapes `\\` and `\"`, as `readYAMLmodel` does. This repair runs only when the strict parse fails. +* `io.read_yaml_model` reads a single-quoted scalar with an apostrophe inside, such as `- name: 'it's'`, as `it's`: it drops the outer pair of quotes and reads each `''` as `'`, as `readYAMLmodel` does. This repair runs only when the strict parse fails; a comment after the closing quote is not handled. * `init.ftinit` and `init.run_ftinit` give each scored reaction its own big-M, its flux capacity `max(min(ub, big_m), min(-lb, big_m), force_on)` with an infinite bound counted as `big_m`, as RAVEN's `ftINITInternalAlg` does. `big_m` is the cap on that capacity; its diff --git a/src/raven_toolbox/io/yaml.py b/src/raven_toolbox/io/yaml.py index 72627e54..9ccad62c 100644 --- a/src/raven_toolbox/io/yaml.py +++ b/src/raven_toolbox/io/yaml.py @@ -117,6 +117,27 @@ def construct_yaml_omap(self, node): ) +# A single-quoted scalar that holds an apostrophe no second apostrophe pairs with, such as +# `'it's'`. It is not valid YAML. `readYAMLmodel.m` strips the outer pair of quotes of the line +# and turns each `''` in what is left into `'`, so the value is `it's`. +_SINGLE_QUOTED_VALUE = re.compile( + r"""^(?P\s*-\s+|\s*(?:-\s+)?[^\s:\-][^:]*:\s+)'(?P.*)'\s*$""" +) + + +def _has_lone_apostrophe(body: str) -> bool: + """Whether ``body``, the text between the outer quotes, has a ``'`` that is not half of a ``''``.""" + i = 0 + while i < len(body): + if body[i] != "'": + i += 1 + elif body[i : i + 2] == "''": + i += 2 + else: + return True + return False + + def _has_unescaped_quote(body: str) -> bool: """Whether ``body``, the text between the outer quotes, has a ``"`` that no backslash escapes.""" i = 0 @@ -131,10 +152,13 @@ def _has_unescaped_quote(body: str) -> bool: def _fix_doubled_quotes(text: str) -> str: - """Rewrite each double-quoted scalar with an unescaped interior quote as a single-quoted scalar. + """Rewrite each quoted scalar with an unpaired interior quote as a valid single-quoted scalar. - The value is the text between the outer quotes with escaped backslashes and quotes - unescaped, as ``readYAMLmodel.m`` reads it; quoted scalars without such a quote are left alone. + A double-quoted scalar with an unescaped ``"`` takes the text between the outer quotes with + escaped backslashes and quotes unescaped as its value. A single-quoted scalar with an + apostrophe that no second apostrophe pairs with takes the text between the outer quotes + with each ``''`` read as ``'``. Both are the values ``readYAMLmodel.m`` reads; quoted + scalars without such a quote are left alone. """ def fix(match: re.Match) -> str: @@ -144,11 +168,19 @@ def fix(match: re.Match) -> str: value = re.sub(r'\\(["\\])', r"\1", body) return f"{match['head']}'{value.replace(chr(39), chr(39) * 2)}'" - return "\n".join(_QUOTED_VALUE.sub(fix, line) for line in text.split("\n")) + def fix_single(match: re.Match) -> str: + body = match["body"] + if not _has_lone_apostrophe(body): + return match[0] + return f"{match['head']}'{body.replace(chr(39) * 2, chr(39)).replace(chr(39), chr(39) * 2)}'" + + return "\n".join( + _SINGLE_QUOTED_VALUE.sub(fix_single, _QUOTED_VALUE.sub(fix, line)) for line in text.split("\n") + ) def _load_yaml_text(text: str): - """Parse ``text`` strictly; when that fails, retry with quotes inside double-quoted scalars repaired.""" + """Parse ``text`` strictly; when that fails, retry with unpaired quotes inside quoted scalars repaired.""" try: return _read_yaml.load(text) except YAMLError: diff --git a/tests/test_io_yaml_ec_data.py b/tests/test_io_yaml_ec_data.py index bf494e03..5ff1cdae 100644 --- a/tests/test_io_yaml_ec_data.py +++ b/tests/test_io_yaml_ec_data.py @@ -526,6 +526,34 @@ def test_interior_unescaped_quote_is_kept_like_matlab(tmp_path, note, expected): assert model.ec.rxns == ["std", "enz"] +@pytest.mark.parametrize( + "value, expected", + [ + ("'it's'", "it's"), + ("'a 'b' c'", "a 'b' c"), + ("'it's ''x'' z'", "it's 'x' z"), + ("'say \"hi\" it's'", 'say "hi" it\'s'), + ], +) +def test_interior_apostrophe_in_single_quoted_value_is_kept_like_matlab(tmp_path, value, expected): + # readYAMLmodel.m strips the outer pair of quotes and turns each '' into ', keeping a lone '. + path = tmp_path / "m.yml" + text = _GECKO_FILE.format(note='"ok"', empty="").replace('- name: "A"', f"- name: {value}") + path.write_text(text, encoding="utf-8") + model = read_yaml_model(path) + assert model.metabolites.get_by_id("a").name == expected + assert model.ec.rxns == ["std", "enz"] + + +def test_valid_single_quoted_values_are_read_unchanged(tmp_path): + path = tmp_path / "m.yml" + text = _GECKO_FILE.format(note="'a ''b'' c'", empty="").replace('- name: "A"', "- name: 'it''s'") + path.write_text(text, encoding="utf-8") + model = read_yaml_model(path) + assert model.metabolites.get_by_id("a").name == "it's" + assert model.ec.notes[0] == "a 'b' c" + + def test_valid_quoted_scalars_are_read_unchanged_next_to_a_repaired_one(tmp_path): path = tmp_path / "m.yml" text = _GECKO_FILE.format(note='"a "b" c"', empty="").replace( From a54cad9256d1b210256437996028ad9feaa17e13 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 18:10:22 +0100 Subject: [PATCH 5/6] fix: fseof takes its top level from the pFBA maximum, sorts stably and can fix the target The top of the scan is max_fraction of the target flux in the pFBA solution that maximises the target, as RAVEN's FSEOF and GECKO's ecFSEOF do. Targets of equal score keep model order. The new equality option fixes the target flux at each level, as ecFSEOF does. --- CHANGELOG.md | 3 ++ src/raven_toolbox/analysis/fseof.py | 28 ++++++++++--- tests/test_analysis_fseof.py | 64 +++++++++++++++++++++++++++++ 3 files changed, 90 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 878495ba..3b8810aa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,9 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased * `analysis.fseof` (and so `geckopy.ec_fseof`) solves each level with pFBA that holds the objective at 0.999999 of its optimum at that level, as RAVEN's `solveLP` does, where it held the full value. On the scipy solver (HiGHS) the optimum is only accurate to the solver's feasibility tolerance, so the second step was infeasible at the top level of ecYeastGEM and the scan returned 9 of 10 levels; it now returns all 10, as on Gurobi and GLPK. A level at which the biomass maximisation is infeasible still ends the scan. +* `analysis.fseof` takes the top of the scan as `max_fraction` of the target flux in the pFBA solution that maximises the target, as RAVEN's `FSEOF` and GECKO's `ecFSEOF` do, where it took the plain maximum (about 1e-6 relative difference). +* `analysis.fseof` lists targets of equal score in model order (a stable sort), as MATLAB's `sort` does, where their order was arbitrary. +* `analysis.fseof` gained `equality`: `True` fixes the target flux at each level as GECKO's `ecFSEOF` does, `False` (the default) keeps the lower bound of RAVEN's `FSEOF`. * `io.read_yaml_model` reads an ec reaction written as `- enzymes: !!omap` with nothing under the tag, the form GECKO and RAVEN write for a reaction without enzymes, as a reaction with no enzymes, as `readYAMLmodel` does. It raised `AttributeError: 'ScalarNode' object has no attribute 'flow_style'`. * `io.read_yaml_model` reads a scalar written as `""text""` (invalid YAML, written by older RAVEN and GECKO for a note that carried its own quotes) as `"text"`, as `readYAMLmodel` does. Files that parse as YAML are read unchanged. * `io.read_yaml_model` reads a double-quoted scalar with an unescaped quote inside, such as `- notes: "a "b" c"`, as `a "b" c`: it drops the outer pair of quotes and unescapes `\\` and `\"`, as `readYAMLmodel` does. This repair runs only when the strict parse fails. diff --git a/src/raven_toolbox/analysis/fseof.py b/src/raven_toolbox/analysis/fseof.py index 72305d92..d038d73c 100644 --- a/src/raven_toolbox/analysis/fseof.py +++ b/src/raven_toolbox/analysis/fseof.py @@ -88,11 +88,13 @@ def fseof( correlation_threshold: float = 0.9, flux_eps: float = 1e-6, min_target: float | None = None, + equality: bool = False, ) -> FSEOFResult: """Run FSEOF for over-production of ``target_rxn``'s product. Enforces target flux from ``max_fraction/n_steps`` up to ``max_fraction`` of the - theoretical maximum in ``n_steps`` steps, maximising growth (``biomass_rxn`` or the + theoretical maximum (the target flux of the pFBA solution that maximises it) in + ``n_steps`` steps, maximising growth (``biomass_rxn`` or the model's current objective) with pFBA at each step. Returns an :class:`FSEOFResult`. Parameters @@ -113,9 +115,13 @@ def fseof( flux sum is minimised with the objective held at 0.999999 of its optimum at that level. A level at which the first step is infeasible ends the scan. + The top of the scan is ``max_fraction`` of the target flux in the pFBA solution that + maximises the target, as in RAVEN's ``FSEOF`` and GECKO's ``ecFSEOF``. + Differences from ``FSEOF.m``: the scan stops at the first infeasible enforced level, so the last absolute flux used for the knockout label is that of the last feasible - level; ``targets`` is sorted by ``|slope| * |r|`` and not grouped by label; and + level; ``targets`` is sorted by ``|slope| * |r|`` (equal scores in model order) and is + neither grouped by label nor split into reaction and transport targets; and ``flux_eps`` only gates the check of the theoretical maximum. min_target: Floor of the enforced-flux range, evenly spaced up to the scan's ceiling @@ -123,6 +129,10 @@ def fseof( Defaults to ``None``, which floors the range at ``max_fraction of the theoretical maximum / n_steps`` instead. Pass e.g. the target reaction's flux at growth-optimum to scan from there. + equality: + How the enforced flux is imposed. ``False`` (the default) sets it as the lower bound + of the target reaction, as RAVEN's ``FSEOF`` does; ``True`` fixes the target flux at + the level (lower and upper bound), as GECKO's ``ecFSEOF`` does. """ with model: # find the theoretical maximum target flux model.objective = target_rxn @@ -130,7 +140,10 @@ def fseof( # slim_optimize returns NaN on an infeasible model; np.isfinite catches that too. if target_opt is None or not np.isfinite(target_opt) or target_opt <= flux_eps: raise ValueError(f"{target_rxn!r} cannot carry positive flux; nothing to scan.") - target_max = target_opt * max_fraction + with model: # the maximum is read from the pFBA solution, as RAVEN's solveLP(model, 1) gives it + model.objective = target_rxn + target_pfba = pfba(model, fraction_of_optimum=_PFBA_FRACTION).fluxes[target_rxn] + target_max = float(target_pfba) * max_fraction if min_target is None: levels = [target_max * (i + 1) / n_steps for i in range(n_steps)] else: @@ -142,7 +155,11 @@ def fseof( with model: if biomass_rxn is not None: model.objective = biomass_rxn - model.reactions.get_by_id(target_rxn).lower_bound = level + target = model.reactions.get_by_id(target_rxn) + if equality: + target.bounds = (level, level) + else: + target.lower_bound = level try: columns[level] = pfba(model, fraction_of_optimum=_PFBA_FRACTION).fluxes except OptimizationError: @@ -198,4 +215,5 @@ def _classify(model, scan, enforced, corr_threshold) -> pd.DataFrame: "reaction", "name", "subsystem", "gene_reaction_rule", "genes", "target_type", "slope", "correlation", "initial_flux", "final_flux", "score", ]) - return table.sort_values("score", ascending=False, ignore_index=True) + # A stable sort keeps reactions with equal scores in model order, as MATLAB's sort does. + return table.sort_values("score", ascending=False, ignore_index=True, kind="stable") diff --git a/tests/test_analysis_fseof.py b/tests/test_analysis_fseof.py index 20deadee..a48f5450 100644 --- a/tests/test_analysis_fseof.py +++ b/tests/test_analysis_fseof.py @@ -198,3 +198,67 @@ def test_level_with_infeasible_maximisation_ends_the_scan(model): # reported, however the second pFBA step is tolerated. res = fseof(model, "EX_P", n_steps=4, max_fraction=2.0) assert res.enforced == pytest.approx([5.0, 10.0]) + + +def test_equality_fixes_the_target_flux_at_each_level(model): + res = fseof(model, "EX_P", n_steps=4, equality=True) + assert res.scan.loc["EX_P"].to_numpy() == pytest.approx(res.enforced) + # the default lower bound lets a level exceed the enforced flux when that suits growth; + # here the product branch is the only sink of the shared intermediate at the top level. + lower = fseof(model, "EX_P", n_steps=4) + assert (lower.scan.loc["EX_P"].to_numpy() >= lower.enforced).all() + + +def test_equality_differs_from_lower_bound_when_the_target_is_coupled_to_growth(): + # S -> B + P is the only route to biomass, so growth alone exports 10 of P; a lower bound + # below that leaves the flux at 10, while an equality holds it at the level. + m = cobra.Model("m") + S, P, B = (cobra.Metabolite(x, compartment="c") for x in ("S", "P", "B")) + m.add_metabolites([S, P, B]) + sup = cobra.Reaction("sup", lower_bound=0, upper_bound=10) + sup.add_metabolites({S: 1}) + v1 = cobra.Reaction("v1", lower_bound=0, upper_bound=1000) + v1.add_metabolites({S: -1, B: 1, P: 1}) + ex_p = cobra.Reaction("EX_P", lower_bound=0, upper_bound=1000) + ex_p.add_metabolites({P: -1}) + ex_b = cobra.Reaction("EX_B", lower_bound=0, upper_bound=1000) + ex_b.add_metabolites({B: -1}) + m.add_reactions([sup, v1, ex_p, ex_b]) + m.objective = "EX_B" + lower = fseof(m, "EX_P", n_steps=4, min_target=2.0) + fixed = fseof(m, "EX_P", n_steps=4, min_target=2.0, equality=True) + assert lower.scan.loc["EX_P"].to_numpy() == pytest.approx(10.0) + assert fixed.scan.loc["EX_P"].to_numpy() == pytest.approx(fixed.enforced) + + +def test_top_level_is_a_fraction_of_the_pfba_maximum(model): + # pFBA holds the target at 0.999999 of its maximum while it minimises the flux sum, + # so the top level is 0.9 of that flux and not of the plain maximum. + from cobra.flux_analysis import pfba + + with model: + model.objective = "EX_P" + plain = model.slim_optimize() + in_pfba = pfba(model, fraction_of_optimum=0.999999).fluxes["EX_P"] + assert in_pfba < plain + res = fseof(model, "EX_P", n_steps=4) + assert res.enforced[-1] == pytest.approx(0.9 * in_pfba, rel=1e-8) + assert res.enforced[-1] != pytest.approx(0.9 * plain, rel=1e-8) + + +def test_equal_scores_keep_model_order(): + # Reactions of three score classes are interleaved in the model; within a class the + # scores tie exactly, and the table lists each class in model order, highest score first. + import numpy as np + import pandas as pd + + from raven_toolbox.analysis.fseof import _classify + + m = cobra.Model("m") + ids = [f"r{i:03d}" for i in range(600)] + m.add_reactions([cobra.Reaction(rid) for rid in ids]) + enforced = np.array([1.0, 2.0, 3.0, 4.0]) + scan = pd.DataFrame([enforced * (1 + i % 3) for i in range(600)], index=ids, columns=enforced) + table = _classify(m, scan, enforced, 0.9) + expected = [ids[i] for c in (2, 1, 0) for i in range(600) if i % 3 == c] + assert list(table["reaction"]) == expected From b0adf72d742eaef8b6190f07bee5e63571cd1231 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 30 Oct 2026 04:54:00 +0100 Subject: [PATCH 6/6] fix: sort the fseof imports and name the exception the invalid-YAML test expects --- src/raven_toolbox/analysis/fseof.py | 1 - tests/test_io_yaml_ec_data.py | 3 ++- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/raven_toolbox/analysis/fseof.py b/src/raven_toolbox/analysis/fseof.py index d038d73c..029bedbb 100644 --- a/src/raven_toolbox/analysis/fseof.py +++ b/src/raven_toolbox/analysis/fseof.py @@ -26,7 +26,6 @@ from cobra.flux_analysis import pfba from scipy.stats import linregress - # pFBA keeps the objective at this fraction of the optimum found at the same level, as # RAVEN's solveLP does when it minimises the flux sum (``solution.f*0.999999``). A solver # reports its optimum to within its feasibility tolerance, so a second solve that demands diff --git a/tests/test_io_yaml_ec_data.py b/tests/test_io_yaml_ec_data.py index 5ff1cdae..5c4feef0 100644 --- a/tests/test_io_yaml_ec_data.py +++ b/tests/test_io_yaml_ec_data.py @@ -18,6 +18,7 @@ import numpy as np import pytest from cobra.io.yaml import yaml as cobra_yaml +from ruamel.yaml import YAMLError from scipy import sparse from raven_toolbox.io import EcData, read_yaml_model, write_yaml_model @@ -568,5 +569,5 @@ def test_valid_quoted_scalars_are_read_unchanged_next_to_a_repaired_one(tmp_path def test_invalid_yaml_other_than_doubled_quotes_still_raises(tmp_path): path = tmp_path / "m.yml" path.write_text(_GECKO_FILE.format(note='"ok"', empty="").replace("kcat: 10", "kcat: [10"), encoding="utf-8") - with pytest.raises(Exception): + with pytest.raises(YAMLError): read_yaml_model(path)