readYAMLmodel: fix position-dependent metCharges/metComps defaults - #709
Merged
Conversation
Both bugs come from the same shape: a per-metabolite field is filled
with a real default only for a *trailing* gap (the only case
emptyOrFill's numel-shortfall check reaches), while a gap earlier in
the list — followed by a metabolite that does have the field — was
already resolved to something else by an earlier conversion step, and
emptyOrFill never revisits it.
- metCharges: a mid-list gap already parses to NaN via str2double('')
and stays that way; only a trailing gap reached emptyOrFill's fill
value, which was 0 — so an unset charge was NaN or 0 depending
entirely on its position in the metabolite list. Moved into the same
NaN-fill group as metDeltaG so it is always NaN, matching what an
unset charge already meant everywhere except that one branch.
- metComps: a mid-list gap parses to '' and then ismember('', comps)
resolves it to 0 — not 1 ("assume first compartment"), and not a
valid MATLAB index, so anything indexing model.comps(0) for that
metabolite would error. Only a trailing gap (the array literally too
short for ismember to touch) reached emptyOrFill's pad-with-1
fallback. Replace every 0 left by the ismember lookup with 1
immediately after it, so the default applies regardless of position.
rxnComps/geneComps have the identical shape but are never populated by
writeYAMLmodel.m in the first place, so the bug path is not reachable
in practice; left alone.
Verified against tSyntax/tIO/tBiomass/tConditions/tQueries (54 passed,
0 failed) and two targeted probes confirming both fields are now
position-independent.
edkerk
added a commit
to SysBioChalmers/raven-toolbox
that referenced
this pull request
Aug 28, 2026
Matches readYAMLmodel.m's own "assume first compartment" convention (just fixed on the MATLAB side for a position-dependent bug in SysBioChalmers/RAVEN#709) for a metabolite with no explicit compartment: cobra's model_from_dict has no equivalent default and would otherwise leave met.compartment as None. Verified against the same probe file used to confirm the MATLAB fix: both readers now resolve it to the identical compartment assignment.
4 tasks
Function test results304 tests 277 ✅ 1m 11s ⏱️ Results for commit 1efaf16. |
edkerk
added a commit
to SysBioChalmers/raven-toolbox
that referenced
this pull request
Aug 28, 2026
#106) * io/yaml: match writeYAMLmodel.m's metaData order and MIRIAM collapsing Closes three concrete gaps found while byte-diffing MATLAB and Python output for the same model: - metaData now follows writeYAMLmodel.m's fixed field order (id, name, version, date, then the annotation-style fields), with the same defaults: id/name fall back to "blankID"/"blankName", date falls back to today. Previously the block just echoed whatever order model.notes['metaData'] happened to be in. - defaultLB/defaultUB are now recomputed from the model's current bounds and included, matching readYAMLmodel.m's own min(model.lb)/max(model.ub) derivation — previously never emitted at all unless already present in the source. - A singleton annotation value (e.g. one kegg.compound id) now collapses to a bare scalar, matching writeYAMLmodel.m's MIRIAM handling — except ec-code/smiles, which RAVEN's own writer always keeps as a list. The reader gained the matching read-side normalisation (wrap a bare scalar into a one-item list before model_from_dict), since cobra doesn't do that itself and a collapsed-then-reread value would otherwise come back as a bare string instead of list[str]. A fourth candidate (defaulting an unset metabolite charge to 0.0) was investigated and deliberately not replicated: readYAMLmodel.m's actual behaviour there is a position-dependent side effect of generic fill-missing-fields code (a mid-list gap parses as NaN and is omitted; only a gap at the very end of the metabolite list happens to hit a different tail-padding branch and become a real 0), not a "charge defaults to neutral" convention worth reproducing. Verified against real RAVEN-authored files: reading and rewriting tutorial/smallYeast.yml through this writer and through RAVEN MATLAB's writeYAMLmodel.m now produces byte-identical output (same SHA256). * io/yaml: default a metabolite's missing compartment to the first one Matches readYAMLmodel.m's own "assume first compartment" convention (just fixed on the MATLAB side for a position-dependent bug in SysBioChalmers/RAVEN#709) for a metabolite with no explicit compartment: cobra's model_from_dict has no equivalent default and would otherwise leave met.compartment as None. Verified against the same probe file used to confirm the MATLAB fix: both readers now resolve it to the identical compartment assignment.
edkerk
added a commit
to SysBioChalmers/raven-toolbox
that referenced
this pull request
Aug 28, 2026
* Add a release workflow (GitHub release + PyPI trusted publishing)
Builds sdist and wheel from the tag, verifies the built version matches it, takes the
release notes from the CHANGELOG section for that version, creates the GitHub release, and
publishes to PyPI via OIDC trusted publishing so no token is stored in the repo.
Triggers on a v* tag push, plus a workflow_dispatch that takes a tag input so a tag that was
already pushed can still be released -- push.tags cannot fire retroactively.
* CI: run the release artifact actions on Node 24
actions/upload-artifact@v4 and actions/download-artifact@v4 target Node.js 20, so every
release run is annotated with the Node 20 deprecation warning. The rest of the workflows moved
to Node 24 action majors in #35 (checkout@v5, setup-python@v6); this brings the release
workflow in line.
The two actions' majors are offset by one, so they cannot both take the same number. v5 (upload)
and v6 (download) only added preliminary Node 24 support and still default to Node 20, which
would leave the warning in place; upload v6 and download v7 are the first to set
'runs.using: node24'. This takes the latest matched pair -- upload v7 with download v8 -- since
upload v7's direct uploads are what download v8's direct downloads consume.
Behaviour notes:
* download-artifact v8 defaults 'digest-mismatch' to error rather than warn, so a corrupted
download now fails the run instead of publishing silently.
* Both require Actions Runner >= 2.327.1; the workflow runs on GitHub-hosted ubuntu-latest.
* The multi-path upload (dist/ plus notes.md) is unaffected: the paths still resolve through
their least common ancestor to the workspace root, so both downstream jobs see 'dist/*'
beside 'notes.md'. v7's single-file restriction applies only with 'archive: false'.
* docs: note diffModels grRule comparison to back-port into raven-toolbox
MATLAB RAVEN's diffModels compares grRules as logic (DNF + sort, so
"a and b" == "b and a"), while raven-toolbox's diff_models._normalise_gpr
only lowercases and collapses whitespace. Record the AST-based comparison as
a pending back-port into raven-toolbox.
* fix: make the placement master deterministic across runs
The score-maximising placement MILP iterated genes_in_scope, a set, so string
hash randomisation reordered its variables and constraints between processes and
a degenerate objective returned a different co-optimal placement on each run.
Iterate the genes sorted, and pin the Gurobi solve (single thread, fixed seed,
presolve level and zero MIP gap) so the same placement is returned every time.
A parameter sweep on the yeast master showed thread count, seed and presolve
level are the settings that change which co-optimal placement is chosen.
* KEGG FASTA/HMM path: align model output with MATLAB RAVEN
Bring the protein-FASTA reconstruction (get_kegg_model_from_sequences)
into byte-parity with RAVEN develop3 getKEGGModelForOrganism:
- Prune each gene-backed reaction's kegg.orthology annotation to the KOs
that actually matched a gene, via a new prune_orthology flag on
assemble_model_from_ko_genes (off by default, so the organism-annotation
path keeps the full KO list). The order-preserving intersection also
avoids the rxnMiriams index misalignment in the MATLAB pruning.
- Default model.id to the FASTA stem when no model_id is given, so the
draft never inherits the reference model's id.
- Use a neutral reaction note ("Included by KEGG HMM reconstruction").
Note: the hmmsearch -Z parity item is handled on the MATLAB side (Python's
-Z <nprofiles> matches the hmmscan scale the K15 cut-offs were calibrated
against).
* Compare grRules as logic in diff_models (order-insensitive)
Implements the back-port this branch previously only documented.
diff_models' GPR check DNF-expands each rule via the existing
manipulation.gpr_to_dnf, sorts genes within each isozyme clause and sorts
the clauses, so operand order no longer registers as a difference:
"a and b" == "b and a", "a or b" == "b or a". The previous _normalise_gpr
only lowercased and collapsed whitespace, flagging logically identical
rules that differed only in operand order.
Uses the reaction's already-parsed cobra .gpr (no re-parse) and reports
the original rule strings in the diff message while comparing canonically.
A rule cobra cannot parse falls back to the old string comparison, so
malformed rules are still compared rather than silently equated to empty.
Brings diff_models in line with MATLAB RAVEN's diffModels (RAVEN #686);
removes the now-resolved pending-back-port note.
* Fix load_delta_g_csv stamping the ΔG missing-value sentinel as a real value (#73)
* Fix load_delta_g_csv stamping the SEED "missing" sentinel as a real value
The ModelSEED-derived side-car tables encode "no valid dG" as the magic
value 10000000. load_delta_g_csv -- written for exactly these files, down
to its Var1/Var2 column defaults -- had no notion of it and stamped it
verbatim, so 777 of yeast-GEM's 4102 reactions (19.5%) carried a
physically impossible 1e7 kJ/mol presented as a measurement. Anything
reading notes["deltaG"] got garbage for a fifth of the model.
yeast-GEM's own checkrxnDirection.m gates on the same value:
if ~isequal(seed_rxnInfo{rxnIdx4(i),16},'10000000.0')
%check if database contains valid deltaG value
Treat the sentinel as missing, exactly as NaN already was, recognising it
whichever dtype the MATLAB/pandas round-trip produces (10000000,
10000000.0, "10000000.0"). The new keyword-only missing_value (default
SEED_DELTA_G_MISSING) tunes or disables it.
Real dG coverage of yeast-GEM is 78.2% (3207/4102), not the 97.1% the
loader previously implied.
* Close the confidence facet set and rewrite the remaining plan
The facet set is localization + equation + gene_association, and no
further facet is planned. Drop the planned fourth facet from the study
doc and from confidence.py's module docstring, which both still promised
one.
Replace the phasing section, whose P1/P2 entries only restated what had
already shipped, with what is actually left:
* Wire the facets together -- an annotate_confidence() umbrella, and let
curation_priority read the record so a mark_curated reaction stops
resurfacing in the review queue. That is the change that closes the
score -> review -> curate loop.
* Validate beyond one model. Every number in the doc is yeast-GEM's, and
three bands never fire there.
* Standards alignment for the paper (Thiele-Palsson / ECO); the SBO half
is already done.
* delta_g: drop unfounded ModelSEED provenance from the missing-ΔG sentinel
The 10000000 "no valid ΔG" sentinel handled by load_delta_g_csv was
described as ModelSEED's and the constant named SEED_DELTA_G_MISSING, but
the ΔG side-car tables (e.g. yeast-GEM's model_rxnDeltaG.csv) are not
ModelSEED-derived — the shared sentinel value does not establish that
provenance. Describe it neutrally as the tables' own missing-value marker
(still evidenced by yeast-GEM's checkrxnDirection.m gating on 10000000.0)
and rename the constant to DELTA_G_MISSING. No behaviour change.
* confidence: wire the facets together (#78)
* confidence: Thiele-Palsson score mapping + verified ECO ids (§10.3) (#80)
Add thiele_palsson_score(reaction): the Thiele & Palsson reconstruction
confidence score (0-4) derived from the gene_association facet's basis
(gpr+literature -> 3, gpr -> 2, no-gpr -> 1), the facet that captures
reaction-inclusion evidence. Grounded in the correspondences already
asserted in score_gene_association_confidence.
Respects the study doc's two cautions: the map runs forward only (basis ->
score), so the ambiguous Thiele-Palsson 2 is never reverse-mapped to an
evidence class, and curated is left unmapped (a curator's assertion does
not name its evidence class). The ECO evidence-class ids in §9 were each
verified against EBI's ECO (OLS) - ECO:0000044 sequence similarity,
ECO:0000015 mutant phenotype, ECO:0000002 direct assay - and the
assertion-method term ECO:0000305 is explicitly excluded as it is not an
evidence class.
Docs: §9 gains the concrete basis -> score -> ECO table; §10.3 marked
mapping-shipped. Tests for the mapping and the end-to-end basis path.
* confidence: cross-model validation of the facet bands (§10.2) (#81)
Run measure_confidence_facets.py over three further models beyond
yeast-GEM (Human-GEM, iYali, panAsp/pAo) and record the result:
- Two of the three bands §6 could not exercise DO fire on other real
models -- formula-generic on all three, formula-unparseable on both
drafts -- so they are calibrated, not yeast artefacts. Only
charge-unknown stays fixture-only across all four (a mass-balanced
reaction with a charge-less metabolite is genuinely rare).
- gpr+literature fires only where reactions carry a pubmed annotation.
- The gene-rubric-vs-recorded-Confidence-Level check needs a
per-reaction-curated model: iYali and pAo record a blanket value, so
they are not graded ground truth the way yeast-GEM is.
Also generalise the script to load SBML or cobra YAML/JSON by extension
(Human-GEM ships as YAML). Docs: §6 caveat resolved with a pointer to
§10.2; §10.2 written up and marked done. -W docs build clean.
* Deterministic score-aligned reaction placement in the assignment master
The placement master maximises a gene-localisation objective that never
mentions the per-reaction placement variable, so each reaction's compartment
was a free co-optimum: the pinned solver returned it reproducibly but
arbitrarily, giving 52.8% reaction agreement with curated yeast-GEM (an earlier
un-pinned build happened to land ~72%).
Add a lexicographic second pass: fix the gene layout to the primary optimum
(fixing the y solution, not the objective value, so there is no tolerance to
tune), then place each reaction in the compartment its own enzymes are
predicted to occupy -- the summed DeepLoc score of the reaction's genes, with a
small default_compartment prior so genes-free and score-tied reactions fall
there deterministically.
Yeast reaction agreement rises to 72.5% (1408/1943) and now rests on the
localisation evidence rather than a solver tie-break; gene agreement is
unchanged (88.7%, 716/807 -- the gene layout is untouched); coherent placement
adds fewer transports (1001 -> 967); growth and blocked fraction unchanged; the
warm-started second solve adds ~4s. Reproducible across independent runs.
Docs updated (yeast_validation, localization_redesign, multiorganism_validation
footnote); the multi-organism reaction-agreement rows await re-measurement.
* Refresh multi-organism validation on the deterministic placement (stacks on #83) (#84)
* Deterministic score-aligned reaction placement in the assignment master
The placement master maximises a gene-localisation objective that never
mentions the per-reaction placement variable, so each reaction's compartment
was a free co-optimum: the pinned solver returned it reproducibly but
arbitrarily, giving 52.8% reaction agreement with curated yeast-GEM (an earlier
un-pinned build happened to land ~72%).
Add a lexicographic second pass: fix the gene layout to the primary optimum
(fixing the y solution, not the objective value, so there is no tolerance to
tune), then place each reaction in the compartment its own enzymes are
predicted to occupy -- the summed DeepLoc score of the reaction's genes, with a
small default_compartment prior so genes-free and score-tied reactions fall
there deterministically.
Yeast reaction agreement rises to 72.5% (1408/1943) and now rests on the
localisation evidence rather than a solver tie-break; gene agreement is
unchanged (88.7%, 716/807 -- the gene layout is untouched); coherent placement
adds fewer transports (1001 -> 967); growth and blocked fraction unchanged; the
warm-started second solve adds ~4s. Reproducible across independent runs.
Docs updated (yeast_validation, localization_redesign, multiorganism_validation
footnote); the multi-organism reaction-agreement rows await re-measurement.
* Refresh multi-organism validation numbers on the deterministic placement
Re-measure Human-GEM, AraCore, and iCre1355 (scripts/benchmark_certified_multiorg.py)
on the deterministic score-aligned reaction placement, and refresh the yeast row
to match yeast_validation.md. The previous reaction-agreement figures came from
the old arbitrary co-optimum placement.
Reaction agreement rose or held (Human-GEM 52.8% -> 56.9%, the largest move;
AraCore 81.8% -> 82.1%; iCre1355 53.4% -> 54.1%; yeast 72.0% -> 72.5%), and
coherent placement needs fewer transports everywhere. Gene agreement is
unchanged (the tie-break leaves the gene layout untouched), confirmed per model.
One non-obvious shift: iCre1355's blocked fraction rose 19.0% -> 32.9% as tighter,
less transport-heavy placement leaves more reactions dead-ended (still certifies).
Removes the "await re-measurement" footnote; updates the prose reaction counts to
the current Human-GEM version (12854 -> 12877).
* assign_compartments gap-fill: reliable flux-based fill, not cobra's MILP (#82)
* assign_compartments gap-fill: reliable flux-based fill, not cobra's MILP
* gap-fill: run on a copy, batch-add candidates, warn on namespace mismatch
Harden the flux-based _gapfill (universal-DB fill in localization.certify):
operate on a model copy so the caller's model is never mutated; add candidates
in one batch (per-reaction adds are super-linear at scale); return the
flux-carrying set sorted, so the result does not depend on which co-optimal
vertex the solver picked; and warn when most universal candidates share no
metabolite id with the model, so a namespace mismatch is distinguishable from
a genuine empty fill instead of returning [] silently.
Knockout-recovery 60/60 (exact reaction each time) and 12/12 on realistic
incomplete drafts, vs cobra's 45% and 0/12.
* feat(ftinit): opt-in deterministic extraction (strict gap + canonical optimum) (#85)
Two keyword flags on ftinit() (and the underlying run_ftinit / fill_tasks), both default
off, so existing behaviour is unchanged.
* strict_gap — one solve per step proven to a fixed absolute gap (0.05, below the 0.1
reaction-score granularity), replacing the relative-gap escalation whose final
near-zero-objective run otherwise accepts an arbitrary within-gap incumbent.
* canonical — a lexicographic phase 2 over the removable reactions, applied to each
extraction step and to the task gap-fill: hold the primary objective at its optimum,
minimise the count of kept/added reactions, then their summed id rank.
* Both guarded by a solution-availability check: a solve that ends at the time limit with
no incumbent raises a clear error (or, for the canonical phase, falls back to the primary
optimum) instead of an opaque AttributeError on a null primal.
* The extraction seed becomes a module constant (_EXTRACT_SEED, still 1234) so the
determinism probe in the study can vary it without editing source.
These pin which of many equally data-consistent optima is returned; they are not a fix for
gene-essentiality determinism. On Human-GEM/DLD1 they cut the reaction-level seed-swing
125 -> 34 but worsened essential-gene flips 5 -> 19 at a 3-7x build-time cost. Measurements,
mechanism and limitations: docs/studies/ftinit_determinism.md.
* manipulation: export merge_compartments and copy_to_compartment (#86)
Both live in manipulation/compartments.py and have done since the package
was renamed, but manipulation/__init__.py never imported the module, so
neither name was public API. The package docstring already advertises
"compartment merge / copy", and docs/reference/migration.md documents both
as ported RAVEN functions (mergeCompartments, copyToComps) with source
links, so the omission contradicted both.
The tests reached into raven_toolbox.manipulation.compartments directly,
which is why nothing caught it. They now import through the package, plus
an explicit assertion that both names are in __all__.
No other module has the same problem: the remaining non-re-exported public
names are genuine internal helpers, and utils.parse_name_comp is documented
in IMPROVEMENTS.md as deliberately internal.
* convert_to_irreversible: carry annotations, subsystem and notes to _REV (#88)
MATLAB's convertToIrrev copies the per-reaction fields (eccodes,
rxnMiriams, subSystems, rxnNotes, ...) onto the reverse reaction when it
splits a reversible one; the Python port copied only the name and the
gene rule.
Downstream that silently cost the reverse reaction its EC code: in
geckopy, fill_eccodes_from_gem returned '' for every _REV reaction,
fuzzy BRENDA matching then found no kcat, and the reverse direction was
left unconstrained. GECKO's own unit-test suite expects kcat 10 on the
ecTestGEM reverse reactions where geckopy produced 0.
Annotations and notes are deep-copied, so editing one direction does not
write through to the other.
* Parity harness: make agreement with MATLAB RAVEN a test (#87)
* fix(io): read RAVEN YAML models that omit an empty section
RAVEN MATLAB writes no `genes:` block at all for a model without genes -- its
own tutorial/small.yml is such a model. cobra's model_from_dict indexes
obj["genes"] directly, so reading one raised KeyError: 'genes' and
raven-toolbox could not open a valid file written by the toolbox it is the
counterpart of.
Every YAML fixture in this repository happened to have genes, so no test
covered it. Found by reading RAVEN's own shipped models, which is now a
standing check (tests/parity/test_yaml_interop.py).
The three entity sections are defaulted to empty lists before handing the data
to cobra. A file that omits `metabolites:` while its reactions reference
metabolites is still an error, and still fails loudly -- that is a broken file,
not a RAVEN one.
* test(parity): make agreement with MATLAB RAVEN a test, not a claim
The Human-GEM, yeast and multi-organism validations live in study documents.
Nothing fails if the two implementations drift apart, so "validated against
MATLAB RAVEN" has been a statement about the past rather than a property of the
current code.
tests/parity/ starts the harness, run with `pytest -m parity` and in its own CI
job. Its README states the contract in three tiers, because "identical output"
is achievable for some functions and meaningless for others: exact for
deterministic transformations, set-level for MILP results that have many optima
of equal value, statistical for sampling. Putting a function in the wrong tier
gives either a flaky test or a vacuous one.
Reference values come from two places:
- RAVEN's own artefacts. Its repository ships models written by MATLAB RAVEN,
so reading them is a real cross-language check needing no MATLAB. RAVEN is
GPL and this package is MIT, so they are read from a checkout via RAVEN_ROOT
and never vendored here; the CI job checks RAVEN out separately. This is what
found the KeyError fixed in the previous commit -- with that fix reverted,
test_reads_every_raven_authored_model fails on small.yml.
- Recorded oracles, for behaviour no file can show: how RAVEN grades elemental
balance, what it normalises a gene rule to. scripts/parity/generate_oracles.m
runs RAVEN once over tests/data/parity/tiny.yml -- authored here, so no GPL
material is involved -- and writes JSON. Those tests skip with instructions
until someone generates them: a visible skip is honest, a fabricated oracle
is not. The comparison logic was verified end to end against throwaway
oracles, which were then deleted.
test_determinism.py is not cross-language but guards the property the tiers
rest on: comparing against MATLAB means nothing if raven-toolbox itself answers
differently each run. Placement and gap-filling were made deterministic in #76,
#83 and the assignment rework with no regression test. Both tests assert first
that the fixture still exercises a real choice, so they cannot quietly become
vacuous.
* docs: add the development roadmap, and rework the open-work backlog
todo.md was a 24-line pointer page. It is now the item-level backlog, and
roadmap.md is the plan it serves: four dependency-ordered phases with what each
delivers, what it depends on, and how you know it is finished.
The phases are driven by what stops a newcomer trusting the port rather than by
what is left to write: nothing tells a user how the two toolboxes differ (phase
0/1), nothing enforces the parity claims (phase 2, started in this branch), the
documentation site has no Python content at all (phase 3), and no one can tell
which functions are mature (phase 4).
Decisions taken and recorded: two new dual-language reconstruction tutorials
rather than porting the published hanpo-GEM protocol; SLIME lipid curation out
of scope while that method work proceeds; no external anchor for 1.0, so phases
3 and 4 may overlap. Open decisions are listed with what they block.
* Parity: add the set-level tier with a recorded baseline (#89)
* test(parity): add the set-level tier with a recorded baseline
PR #87 stated the tier-2 contract and left it without tests. Extraction is a
MILP, so there is no single right answer to assert -- but there is a wrong
outcome worth catching: the result quietly moving.
test_set_level.py runs run_init over RAVEN's smallYeast.yml and compares the
kept reactions against a baseline recorded by scripts/parity/record_baseline.py.
On failure it prints the Jaccard overlap and both directions of the diff, so a
one-reaction shift is immediately distinguishable from a rewrite, and says how
to re-record when the move is intended.
The assertion is exact set equality rather than an overlap band because that
was measured, not assumed: on this fixture GLPK and Gurobi return the same 13
of 53 reactions, and each is identical across three runs. A difference means
this package changed, not that the solver picked another optimum. On a fixture
where the solvers genuinely disagree the right form is a band with a measured
floor -- not a loosened threshold here.
The baseline is seeded from raven-toolbox, so this guards against regression
rather than proving agreement with MATLAB. The file's `source` field says so and
the failure message repeats it; when generate_oracles.m grows an extraction
oracle the same comparison becomes a real parity check.
A second test asserts the baseline still keeps some but not all reactions, so
the check cannot become vacuous if the fixture or scores change.
* test(parity): genome-scale checks and the nightly job that can run them
The Human-GEM, yeast and multi-organism validations are the strongest evidence
this package agrees with MATLAB RAVEN, and the least protected: measured once,
by hand, on models no free runner's solver licence will accept.
parity-nightly.yml runs them on a Gurobi Web License Service licence, supplied
as GUROBI_WLSACCESSID / GUROBI_WLSSECRET / GUROBI_LICENSEID. WLS is the licence
type meant for containers and CI; a named-user academic licence is tied to a
machine and does not belong on a hosted runner.
The job refuses to be quietly useless. If the secrets are missing it fails with
a message rather than skipping, and before any test runs, check_license.py
solves a 20,000-variable MIP -- far past the 2000-variable cap of the licence
bundled with gurobipy -- so a silent fallback to that licence cannot leave the
genome-scale checks skipped under a green tick.
test_genome_scale.py compares as a *band*, not exact equality, which is the
opposite choice from the small fixture and for a reason: at this size the MILP
has many optima of equal value and the solve is bounded (mip_gap, time_limit),
so a different optimum is legitimate. It also asserts a sanity floor that holds
whichever optimum turns up -- non-empty, not everything, no duplicates -- so
the run is checked even before a baseline exists. The band floor will be
recorded alongside the run that produces it rather than guessed now.
The pull-request job keeps running everything except these (`-m "parity and not
slow"`) and now installs gurobipy: its bundled licence is useless at genome
scale but ample for the 53-reaction fixture, which gives the cross-solver check
a second solver to compare against.
* Parity: keep the genome-scale numbers, budget the nightly, and check determinism across processes (#90)
* test(parity): keep the genome-scale numbers, and budget the job realistically
Ran the genome-scale check locally against the real Human-GEM (12,931
reactions) with a full Gurobi licence. It passed, and exposed two problems with
what #89 shipped.
**The numbers were thrown away.** pytest captures stdout for a passing test, so
the extraction's reaction count and timing -- the whole point of the run, and
the candidate for the next baseline -- were invisible in precisely the case
that matters. The fixture now writes a JSON report to
$RAVEN_PARITY_REPORT_DIR, appends a line to the GitHub job summary, and the
workflow uploads the report as an artefact even on failure, since a failure is
when the actual reaction set most needs inspecting.
**The job budget was wrong.** The run took 68 minutes on a fast workstation
against a 90-minute job timeout, which a two-core hosted runner would not have
met. Raised to 150, with the reasoning recorded: only ~2 of those 68 minutes
are parsing and copying (measured: 94s to read the YAML, 29s to copy the
model), and run_init's time_limit does not bound the call -- a 60-second limit
still ran for well over fifteen minutes. pytest-timeout is therefore the real
stop, and timeout-minutes only has to exceed it so an overrun is reported as a
test failure rather than a killed job.
That time_limit behaviour is a finding about run_init rather than about the
test, and is being measured separately.
* test(parity): check determinism across processes, not just within one
The existing determinism tests repeat a call in one process, which cannot see
the failure mode the placement fixes were actually about: Python randomises
string hashing per process, so a set iterated to build constraint rows or a
dict deciding a tie-break gives a stable answer within a run and a different
one in the next. Repeating in-process, however many times, is blind to it.
The new test runs the gap-fill in three subprocesses under PYTHONHASHSEED 0, 1
and 12345 and compares digests of the result -- which is what
scripts/determinism_probe.py was doing by hand while those fixes were being
made. All three agree today.
The worker sets cobra's process count to 1. The gap-filler runs FVA, which
otherwise starts a worker pool per invocation; with three workers doing that at
once the test took over ten minutes and spawned nineteen processes. Single
-process, on a four-reaction model, it takes 53 seconds.
Like the other determinism tests, it first asserts the fixture still requires a
fill, so it cannot pass by checking nothing.
* fix(homology): match RAVEN's BLAST e-value default (#91)
* fix(homology): match RAVEN's BLAST e-value default
run_blast defaulted to 1e-5 where MATLAB RAVEN's getBlast hardcodes
-evalue 10e-5 -- which is 1e-4, ten times looser. The same two proteomes
therefore gave the two toolboxes different hit sets, before any of the
homology logic ran.
Found while building the homology parity scenario in raven-gecko-parity: the
BLAST checkpoint had to pass RAVEN's value explicitly for the comparison to
mean anything. With the defaults aligned it does not.
RAVEN hardcodes the cutoff; here it stays an argument, so a stricter search is
still available to anyone who wants one.
run_diamond is left alone and documented instead: getDiamond passes no
--evalue at all, so DIAMOND's own default of 1e-3 applies on the MATLAB side --
a third value, looser than either. Picking one needs a decision about which is
right, not a quiet change.
* fix(homology): match RAVEN's DIAMOND e-value default too
run_diamond defaulted to 1e-5. getDiamond passes no --evalue at all, so
DIAMOND's own default of 1e-3 applies on the MATLAB side -- two hundred times
looser. Same reasoning as the BLAST default in the previous commit: MATLAB came
first, and nothing measured here says otherwise.
The two aligners still disagree with each other (1e-4 for BLAST, 1e-3 for
DIAMOND) because RAVEN's own defaults disagree; that is now stated in both
docstrings instead of being implicit in an unpassed flag.
Worth noting what this is *not*: the KO-assignment cutoffs in assign_kos do
diverge from RAVEN deliberately (1e-30 against RAVEN's 1e-50, gene ratio 0.9
against 0.8), and that divergence is backed by the measurements in
docs/studies/kegg_hmm_cutoff_calibration.md. No comparable study exists for the
homology e-value on either side, so there is nothing to weigh against following
MATLAB here.
* Homology cut-offs: measured, and min_align_len lowered to 100 (#92)
* docs: protocol for calibrating the homology cut-offs
RAVEN and raven-toolbox share three thresholds that decide which template
reactions transfer to a new organism -- max_evalue 1e-30, min_align_len 200,
min_identity 40 -- and neither records where they came from. This is the design
for measuring them, written before any numbers exist so the success criterion
cannot be chosen to fit the curves.
It retargets the question. The aligner's own cut-off is inert: run_blast and
getBlast collect at 1e-4, then get_model_from_homology filters again at 1e-30,
twenty-six orders stricter, so every hit that survives model building would have
survived collection. Tuning the collection threshold changes nothing; the levers
are the three above plus bidirectionality.
Ground truth comes from two sources whose weaknesses do not overlap. KEGG
cross-organism models share a reaction namespace, so the comparison needs no id
mapping and scales to a distance series -- but KEGG's own orthology is built on
all-against-all BLAST, so it measures agreement with KEGG's homology calls, not
correctness, and the protocol says so where the numbers will be read. Curated
GEM pairs (yeast-GEM against hanpo-GEM and rhto-GEM) supply the non-circular
check.
Measurement is at hit level first -- do the two genes share a KO -- because
reaction-level agreement is buffered: a reaction transfers if any one of its
template genes hits, so the sweep would look flat for reasons unrelated to the
parameter.
The sweep is cheap: align once per organism pair, cache the hit table, and every
threshold combination is post-processing of it.
The objective is fixed in the protocol, the very-distant pairs act as a
constraint rather than a target so recall cannot be maximised until everything
transfers, and "the current values are as good as anything measured" is named as
an acceptable outcome.
* docs: run the homology cut-off study -- defaults unchanged
The measurements are in. Under the criterion fixed before any data existed, no
change to max_evalue / min_align_len / min_identity is justified, and they stay
as they are. Three findings do not depend on that criterion.
max_evalue does nothing. From 1e-4 to 1e-50 the extracted model is identical --
on both ground truths, on every organism. Only 1e-100 changes anything, because
identity and alignment length have already excluded whatever a looser e-value
would admit. It is not a tuning knob and calibration effort spent on it is
wasted.
min_identity is the only real lever, and its best value moves with phylogenetic
distance: ~30 for a close relative, ~25 for medium and distant. The default of
40 costs 1 point of F1 on Kluyveromyces and 10 on Aspergillus. No single global
default is right for every reconstruction, which the documentation should say.
min_align_len is nearly inert below 150; the default 200 costs about 2 points.
The curated-GEM arm had to be retired, and the study can show why rather than
just suspect it. hanpo-GEM's draft was built by getModelFromHomology at
(1e-30, 150, 35), and the sweep's optimum lands exactly there. The
discriminating test: walk the thresholds loose one step at a time and ask what
fraction of newly admitted reactions are in the curated model. The rate holds
between 0.62 and 0.85 up to the build settings and collapses to 0.06 one step
past them. That cliff is the reference remembering its own construction, so any
optimisation against it returns the build settings whether or not they are good.
Most curated non-model fungal GEMs are RAVEN drafts from yeast-GEM, so this
applies to the class, not just this model.
The KEGG arm ran because UniProt carries KEGG cross-references, which solved the
gene-id problem that the missing (subscription-only) KEGG proteomes had created.
What is left open is a value judgement rather than a measurement: the criterion
blocked a change that improves F1 on all four organisms, including the
constraint organism, because it counted calls rather than correctness. Rewriting
it after seeing the curves is what fixing it in advance was meant to prevent, so
it stands, and the real question is stated instead -- what a wrongly transferred
reaction costs relative to a missed one. Gap-filling can recover the second; the
first pollutes the model silently. Answering that gives a loss function, and the
cached hit tables can be re-scored in minutes.
Both drivers are committed; the KEGG one reproduces every published figure from
the cached tables.
* docs: rescore the homology study at beta = 0.5 -- defaults confirmed
The loss function was the one thing the study deliberately left open, since it
is a value judgement rather than a measurement. It is now decided on the
asymmetry it turns on: a wrongly transferred reaction is worse than a missed
one, because gap-filling can recover the second while the first pollutes the
model and its gene associations silently. So precision is weighted above recall,
beta = 0.5.
Rescored, the conclusion reverses in the most useful direction: min_identity 40
is optimal on the close and very-distant organisms outright and within 0.01 of
the best on the other two. The 10-point deficit at distance that the F1 scoring
showed was entirely an artefact of counting a missed reaction as costly as a
wrong one. The default was right.
One candidate change misses narrowly and is recorded rather than adopted:
min_align_len 50 improves all four organisms (+0.009 to +0.020) at essentially
unchanged precision, but eco calls rise 137 -> 152, violating the constraint,
and the improvement of 0.011 is smaller than the 0.056 within-band spread the
margin rule compares against. Both gates are doing debatable work at this scale
-- comparing an improvement against the spread between organisms of different
difficulty is not a noise estimate -- so the evidence is stated and the decision
left to a maintainer rather than taken here.
The driver takes --beta, so the same cached hit tables can be re-scored under a
different weighting in a couple of minutes without realigning anything. That the
same measurements recommend an identity of 25 or 40 depending on beta alone is
the study's most transferable finding: a threshold recommendation without a
stated weighting is not a recommendation.
* feat(homology): lower min_align_len to 100, measured against KEGG orthology
The three filters that decide which template reactions transfer had never been
measured in either toolbox. Scored against KEGG's own orthology across a
distance series -- close yeast, medium, distant fungus, bacterium -- with
precision weighted above recall, because gap-filling can recover a missing
reaction while a wrongly transferred one is hard to find and harder to remove.
min_align_len 200 -> 100. Everything at or below 150 measures the same; the loss
appears between 150 and 200. Dropping to 100 recovers 3-4 points of recall on
every organism tested while precision moves by at most 0.006 -- real orthologs
that 200 was discarding for nothing. 50 and 100 are identical to three decimals,
so 100 is taken as the less permissive of two equivalent values.
min_identity 40 is confirmed: it wins outright on the closest and most distant
organisms and trails by under 0.01 on the two between. It is also the only filter
that materially decides anything.
max_evalue is inert. Any value from 1e-4 to 1e-50 yields the identical model,
because identity and length have already excluded whatever it would exclude. It
stays at 1e-30 for continuity with RAVEN, and the docstring now says it is not
worth tuning.
Worth recording: at beta = 1 the same measurements recommend identity 25 and
call the default a 10-point deficit at distance. At beta = 0.5 they confirm 40.
A cut-off recommendation without a stated weighting is not a recommendation.
min_align_len now diverges from RAVEN's 200, so a back-port is worth proposing
rather than leaving the toolboxes quietly different.
* feat(homology): report the near misses instead of discarding them
The filters are strict on purpose: a wrongly transferred reaction is hard to
spot and harder to remove, while a missing one can be gap-filled. That argument
justifies rejecting borderline hits. It does not justify throwing away the
evidence for them, which is what happened -- a reaction that missed the identity
threshold by half a percent left no trace.
get_model_from_homology now takes review_identity (e.g. 25). Reactions that
would have transferred at that looser identity are returned in
HomologyResult.candidates with the hit supporting each, sorted strongest first,
and are *not* added to the model. On a real reconstruction (yeast-GEM to
H. polymorpha) that is 303 candidates beside a 2,078-reaction draft: too many to
accept blindly, few enough to read.
The reported identity is the *limiting* one. Matching is bidirectional, so a
pair has to clear the threshold both ways; the first version of this reported
the better direction and the top candidates came out at 44 % against a threshold
of 40, which would have left a curator wondering why a comfortable match was
turned away. Now they come out at 39.7 %, 39.6 %, 39.5 % -- the near misses they
actually are. Covered by a regression test.
* docs: add the OMA and DIAMOND results, and trim the study to its findings
OMA infers orthologs without using BLAST, so it answers the obvious objection to
the KEGG numbers: that they partly measure agreement with the method being
tested. Both references put the best identity between 35 and 45 and neither
supports anything looser, so 40 stands on two independent sources rather than
one.
DIAMOND reaches the same optimum as BLAST on every organism, within 0.006. It
finds about half as many matches, but nearly all the missing ones are weak ones
these settings discard anyway -- after filtering the two agree on 87-92% of what
survives, while DIAMOND runs 10-20x faster. That matters because both toolboxes
apply the same three settings to whichever aligner ran, which until now was an
assumption.
The driver grew --aligner and --reference so all three combinations reproduce
from the cached alignments; verified against every published figure.
The write-up is now about what was found rather than how it was set up, and is
roughly half the length.
* Fix the nightly's Gurobi licence lookup
parity-nightly.yml had failed on its licence check every night since it was
added: it read GUROBI_WLSACCESSID / GUROBI_WLSSECRET / GUROBI_LICENSEID, but
only the first exists. The licence the organisation shares is GUROBI_EDUARD,
holding a licence file verbatim.
- Accept either shape, and report which kind was found without printing it.
A named-user licence cannot authenticate on a hosted runner.
- Still fails loudly when neither is configured.
- Add a test that parses every workflow run: block with bash -n; nothing else
checks a scheduled workflows shell until the night it runs.
* docs: give evidence_aware_transport_cost a Returns section (#93)
The closing paragraph of the docstring ("Returns every metabolite base
-> cost ...") sat directly against the Parameters list at the same
indentation as a parameter name, so the numpydoc parser read it as two
further parameters, named "Returns" and "is". Both then rendered as
parameter rows on the API page, and the generator warned four times:
No types or annotations for parameters ['Returns']
Parameter 'Returns' does not appear in the function signature
Move that paragraph into a proper Returns section, and document `model`,
`annotation` and `base_metabolite`, which the Parameters list had left
out -- six of the nine arguments were described.
The `name:` parameter style is unchanged: it is what the rest of the
package uses, and it resolves types from the signature without
complaint.
Documentation only; no functional change.
* Make every executable in a downloaded bundle executable
zipfile does not restore Unix permission bits, so extracted files arrive
non-executable on Linux and macOS. ensure_binary marked only the executable it
was asked for, and returned a cached bundle without marking anything at all.
A bundle can provide several tools: BLAST ships blastp and makeblastdb. Fetching
blastp then calling makeblastdb raised PermissionError, which is what broke the
parity nightly. Affects any Linux or macOS user, not only CI.
* io.yaml: adopt a shared RAVEN/raven-toolbox YAML layout (#97)
* io.yaml: adopt a shared RAVEN/raven-toolbox YAML layout
The two writers previously disagreed on line folding, number
formatting and how an empty subsystem was represented, which made
yeast-GEM / Human-GEM churn on every write regardless of which
implementation touched them last.
- write_yaml_model now dumps through its own ruamel instance with a
very large width instead of cobra's shared module-level one, so long
scalars are never folded onto a continuation line. Mutating the
shared instance's width would also change cobra.io.save_yaml_model
for the rest of the process.
- Every numeric leaf is coerced to float before dumping, so a whole
number always gets an explicit ".0" (matching RAVEN's writer, which
cannot distinguish an int-typed field from a float-typed one since
its model struct stores everything as double).
- A reaction's subsystem is omitted entirely when empty, like every
other optional field, instead of round-tripping cobra's own
round-trip preservation of older RAVEN files' "blank list entry"
convention. A non-empty subsystem is now always written as a list,
even a single one, matching writeYAMLmodel.m and the reader's
existing list-form support (a reaction can carry more than one).
Indentation and quoting were already ruamel's own defaults and needed
no change. Two tests asserted the previous round-trip shape for a
plain-string subsystem input; updated to expect the new single-item
list, which is what the format now guarantees regardless of how the
in-memory Reaction object represented it.
Verified byte-identical against RAVEN's writeYAMLmodel.m on both a
plain model and one exercising every optional field. Full test suite:
874 passed, 0 failed.
* io.yaml: trim comments to statements, not arguments
Comments were re-litigating why the layout diverges from cobrapy's
rather than just describing what it is; cut down to the latter.
* Default a new metabolite's name to its id in add_reactions_from_equations (#98)
The mets_by="id" path created a metabolite via Metabolite(token,
compartment=compartment), leaving name empty, where addRxns on the
MATLAB side names it after its id. A nameless metabolite is invisible
to anything that keys on name -- merge_compartments (whose own
docstring tells callers to pass base_metabolite=lambda m: m.name for
name-keyed models), add_reactions_from_model, and mergeCompartments on
the MATLAB side.
* Teach apply_condition's reset_exchanges the direction RAVEN reads it as (#99)
cobrapy has no concept of RAVEN's in/out exchange directions, so any
truthy reset_exchanges value reset every exchange reaction regardless
of what was actually named. applyCondition forwards the value straight
through to getExchangeRxns as a direction filter; a condition naming
"out" specifically -- exactly the value both docstrings use as their
own worked example -- silently reopened uptake reactions it never
asked to touch.
_exchange_direction replicates getExchangeRxns' own in/out rule
(nothing produced at all -> boundary met is implicitly the product,
"out"; nothing consumed -> "in"), not restricted to single-metabolite
reactions the way cobra's own Reaction.boundary is. "all"/"both", or
any other truthy value that isn't a recognised direction keyword,
falls back to resetting every exchange in either direction --
preserving prior behaviour for a bare `true`.
Verified against the exact scenario reproducer: reset_exchanges: "out"
on smallYeast now leaves glcIN/o2IN/ethIN untouched and reports
growth=0.0, matching RAVEN exactly (apply_condition_smallyeast MATCHes
end to end).
* Default check_tasks/find_task_essential_reactions to RAVEN's additive boundary semantics (#100)
close_boundaries defaulted to True, closing every exchange/sink/demand
reaction before applying a task's constraints -- the opposite of
RAVEN's own checkTasks, which relaxes metabolite balance for a task's
declared inputs/outputs and leaves the model's own reactions, including
open exchanges, untouched. A task the model can already satisfy through
its own boundary without the task naming it as an input/output passed
in RAVEN and failed here.
Confirmed on smallYeast: GLC_AER_TIGHT (declares only CO2 as output;
the model ships ethOUT/glyOUT/acOUT open) now agrees with RAVEN's
feasible verdict instead of reporting infeasible.
close_boundaries=True remains available for the stricter reading,
where a task's declared inputs/outputs are the complete boundary of
the system for that check.
* Hide RXNS bound cells matching the model's declared default in export_to_excel (#101)
export_to_excel always wrote LOWER BOUND / UPPER BOUND literally;
exportToExcelFormat leaves the cell blank when the bound equals the
model's own declared default (an irreversible reaction's lower bound
is hidden separately whenever it is exactly 0, regardless of the
declared default). Not a corner case: on smallYeast, whose bounds are
drawn entirely from {-1000, 0, 1000} (its own declared defaults plus
zero), nearly every reaction's bound cells were affected.
_rxn_bound_cells reads the same defaultLB/defaultUB metadata.notes
entries the MODEL sheet already reads, and replicates
exportToExcelFormat's hiding rule exactly.
Verified against the exact scenario reproducer: PGI/FBA/HXK now come
back fully blank and glcIN/o2IN blank-lb/literal-ub, matching RAVEN's
output cell for cell (export_to_excel_smallyeast MATCHes end to end).
* Make remove_duplicate_reactions' round trip with expand_model lossless (#103)
Matches RAVEN's contractModel, which this was a lossy counterpart to in
three ways: the survivor was the last-encountered reaction rather than
the first; its id kept an "_EXP_N" suffix rather than having it
stripped; and only the survivor's own GPR was kept rather than the
union of every duplicate's.
The isozyme relationship contractModel restores by design -- three
reactions that expand_model split out of one, each carrying one of the
original's isozyme genes -- was previously lost on the Python side:
two of the three genes came out unassociated with the reaction.
makeEcModel runs this exact expand-then-contract round trip, so the
gap reached the gecko pair too.
_top_level_or_clauses walks cobra's own GPR AST rather than
reimplementing RAVEN's bracket-matching string search; the union is
deduplicated and AND-clauses are parenthesised for readability.
Verified against the exact reproducer: expanding smallYeast's HXK (an
isozyme of three genes) and contracting back now returns HXK (not
HXK_EXP_3) with all three genes restored, matching RAVEN's own
documented output for this case exactly.
* Default load_delta_g_csv to stamping every value literally, matching deltaGCSV (#102)
* Default load_delta_g_csv to stamping every value literally, matching deltaGCSV
missing_value defaulted to DELTA_G_MISSING (1e7), so a matched CSV
value equal to yeast-GEM's own "no measurement" sentinel was silently
left unstamped instead of recorded. RAVEN's deltaGCSV has no sentinel
concept at all and stores whatever the CSV says, verbatim -- applying
the identical CSV to the identical model left the two implementations
disagreeing about whether an entity had a value on record at all.
Both sides are now symmetric: literal by default, with an explicit,
opt-in way to treat a chosen value as missing instead (deltaGCSV's new
missingValue argument; here, passing DELTA_G_MISSING explicitly rather
than relying on the default).
* test_annotation: opt these two into the sentinel default explicitly
load_delta_g_csv's missing_value default flipped to None (stamp
verbatim); these two tests were exercising the old skip-the-sentinel
default and need to request it explicitly now.
* fix: merge_compartments kept neither exchanges nor the objective (#96)
* fix: merge_compartments kept neither exchanges nor the objective
Two defects, both silent, both turning a working model into one that
cannot grow.
1. Every reaction left with a single metabolite after merging was
deleted, including the ones that had a single metabolite to begin
with -- the exchanges. On the small yeast model used by the
documentation site that removed all eight boundary reactions plus
biomassOUT: 11 reactions dropped where only 3 were genuine
transports. Only reactions that *become* trivial through the merge
are dropped now, which is what RAVEN does (it records the
single-metabolite reactions before merging and excludes them from
deletion).
2. The merged model is rebuilt from scratch and the objective was not
carried over, so the result optimised to 0.0 with nothing to
indicate why. It is copied now, with the direction, and a warning is
raised in the case where every objective-carrying reaction really
was removed.
Together these took that model from growth 0.1222 to 0.0000 with no
error and no warning; it now flattens to 0.1268, slightly higher, as
expected once compartmentalisation no longer constrains it.
The existing tests used a toy model with no exchange reactions and no
objective, so neither path was covered. Adds a test for each.
* compartments: fix import order (ruff I001)
* io/yaml: pin the writer's block-sequence indent explicitly (#105)
Everything about the writer's list layout already matched RAVEN's
writeYAMLmodel.m (each list entry two spaces past its parent key), but
only because ruamel's own indent defaults happened to produce that.
Set them explicitly so a future ruamel version can't silently change
the on-disk layout out from under either writer. Verified byte-identical
against the previous (implicit) output.
* manipulation: add set_exchange_bounds, RAVEN's setExchangeBounds port (#104)
cobra's model.medium covers the common single-direction, close-the-
rest case, but not everything setExchangeBounds does: independent lb
and ub per metabolite (medium sets only whichever bound is "active"
for a reaction, from one scalar), a media-only compartment filter, and
a check that every exchange reaction agrees on which sign of flux is
import before touching any of them (closing "import" on the rest
needs a consistent answer to that, or it closes the wrong bound).
Matches metabolites against the model's exchanged set by id or name
(case-insensitive), reports ones that weren't found, and warns (same
as RAVEN) when a metabolite is exchanged in more than one reaction --
RAVEN's own version of that specific warning has a latent bug (an
index-range mismatch between two differently-sized arrays means it
never actually fires); this one does.
Cross-validated directly against setExchangeBounds.m on smallYeast:
explicit bounds plus close_others under RAVEN's own mixed-direction
warning path, unused-metabolite reporting, the no-metabolites-given
default, and the media_only error path (smallYeast has no
extracellular compartment) -- all four cases match bound-for-bound.
* io/yaml: match writeYAMLmodel.m's metaData order and MIRIAM collapsing (#106)
* io/yaml: match writeYAMLmodel.m's metaData order and MIRIAM collapsing
Closes three concrete gaps found while byte-diffing MATLAB and Python
output for the same model:
- metaData now follows writeYAMLmodel.m's fixed field order (id, name,
version, date, then the annotation-style fields), with the same
defaults: id/name fall back to "blankID"/"blankName", date falls
back to today. Previously the block just echoed whatever order
model.notes['metaData'] happened to be in.
- defaultLB/defaultUB are now recomputed from the model's current
bounds and included, matching readYAMLmodel.m's own
min(model.lb)/max(model.ub) derivation — previously never emitted at
all unless already present in the source.
- A singleton annotation value (e.g. one kegg.compound id) now
collapses to a bare scalar, matching writeYAMLmodel.m's MIRIAM
handling — except ec-code/smiles, which RAVEN's own writer always
keeps as a list. The reader gained the matching read-side
normalisation (wrap a bare scalar into a one-item list before
model_from_dict), since cobra doesn't do that itself and a
collapsed-then-reread value would otherwise come back as a bare
string instead of list[str].
A fourth candidate (defaulting an unset metabolite charge to 0.0) was
investigated and deliberately not replicated: readYAMLmodel.m's actual
behaviour there is a position-dependent side effect of generic
fill-missing-fields code (a mid-list gap parses as NaN and is omitted;
only a gap at the very end of the metabolite list happens to hit a
different tail-padding branch and become a real 0), not a "charge
defaults to neutral" convention worth reproducing.
Verified against real RAVEN-authored files: reading and rewriting
tutorial/smallYeast.yml through this writer and through RAVEN
MATLAB's writeYAMLmodel.m now produces byte-identical output (same
SHA256).
* io/yaml: default a metabolite's missing compartment to the first one
Matches readYAMLmodel.m's own "assume first compartment" convention
(just fixed on the MATLAB side for a position-dependent bug in
SysBioChalmers/RAVEN#709) for a metabolite with no explicit
compartment: cobra's model_from_dict has no equivalent default and
would otherwise leave met.compartment as None. Verified against the
same probe file used to confirm the MATLAB fix: both readers now
resolve it to the identical compartment assignment.
* io/yaml: quote with double quotes, matching Prettier's default (#107)
write_yaml_model relied on ruamel's own default (single-quote a scalar
when quoting is required), which was never a deliberate choice, just
what the writer happened to inherit. Switched to double quotes to
match Prettier's YAML default instead, since byte-parity with a
particular tool's default isn't the actual constraint here, only
staying importable across RAVEN and raven-toolbox is.
ruamel's round-trip dumper has no global "prefer this quote character"
switch, only a per-value style override, so this needed its own
needs_quote check (ported from writeYAMLmodel.m's, so both writers
quote exactly the same set of values) and a recursive pass wrapping
exactly those strings in DoubleQuotedScalarString before dumping --
wrapping unconditionally would have quoted every scalar regardless of
whether YAML requires it.
No reader change needed: read_yaml_model goes through ruamel's real
parser, which already handles either quote style and its escaping
correctly.
* Prepare 0.4.0 release (#109)
Bump version 0.3.0 -> 0.4.0 and complete the CHANGELOG 0.4.0 section: six MATLAB RAVEN
parity fixes found by the new raven-gecko-parity harness, the YAML writer brought to
byte-parity with writeYAMLmodel.m, BLAST/homology defaults matched to RAVEN, deterministic
compartment placement, confidence-score interoperability (Thiele-Palsson, ECO), a new
set_exchange_bounds port, and merge_compartments made public and fixed.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Found while byte-comparing RAVEN's and raven-toolbox's YAML output for the same model. Both bugs share the same shape: a per-metabolite field is filled with a real default only for a trailing gap (the only case
emptyOrFill's numel-shortfall check reaches), while a gap earlier in the list — followed by a metabolite that does have the field — was already resolved to something else by an earlier conversion step thatemptyOrFillnever revisits.metCharges: a mid-list gap already parses toNaNviastr2double('')and stays that way; only a trailing gap reachedemptyOrFill's fill value, which was0— so an unset charge wasNaNor0depending entirely on its position in the metabolite list. Moved into the same NaN-fill group asmetDeltaGso it's alwaysNaN, matching what an unset charge already meant everywhere except that one branch.metComps: a mid-list gap parses to''and thenismember('', comps)resolves it to0— not1("assume first compartment"), and not a valid MATLAB index, so anything indexingmodel.comps(0)for that metabolite would error. Only a trailing gap (the array literally too short forismemberto touch) reachedemptyOrFill's pad-with-1fallback. Now every0left by theismemberlookup is replaced with1immediately after it, so the default applies regardless of position.rxnComps/geneCompshave the identical shape but are never populated bywriteYAMLmodel.min the first place, so the bug path isn't reachable in practice — left alone.Test plan
tSyntax,tIO,tBiomass,tConditions,tQueries— 54 passed, 0 failed, 2 filtered (unrelated missing dependencies)c, c, e)