diff --git a/HISTORY.md b/HISTORY.md index 0e818a04cf..96f46a5603 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -2,6 +2,86 @@ ## Unreleased (on master) +#### New: Alias elimination for `ReactionSystem`s + +- **New `aliases` field on `ReactionSystem`** for declaring symbol equivalences + (`lhs ~ rhs` means lhs is eliminated in favor of rhs). Aliases are automatically + eliminated during model conversion (`ode_model`, `sde_model`, `jump_model`, + `hybrid_model`, `ss_ode_model`) and problem construction (`ODEProblem`, `SDEProblem`, + `JumpProblem`, `HybridProblem`, `NonlinearProblem`, `SteadyStateProblem`). + + Supported alias types (v1): + - Ordinary species to ordinary species + - BC-species to BC-species + - Constant species to constant species + - Non-species unknowns (variables) to non-species unknowns + - Parameters to parameters (with strict metadata matching) + + Not supported in v1: Brownian, Poissonian, compound species, binding key, and + observable aliasing. Cross-category aliasing (e.g. species to parameter) is rejected. + +- **New `@aliases` DSL option** for declaring aliases in `@reaction_network` and + `@network_component`: + ```julia + rn = @reaction_network begin + @aliases begin + A ~ B # Species alias: A eliminated, B survives + k1 ~ k2 # Parameter alias: k1 eliminated, k2 survives + end + k1, A + C --> D + k2, B --> E + end + ``` + +- **New `eliminate_aliases` keyword argument** (`true` by default) on all model + conversion functions and problem constructors. Set to `false` to keep aliases as + algebraic constraints (ODE/SDE paths only; jump systems require elimination). + +- **New `remap_alias_inputs` function** for remapping user-provided `u0`/`p` maps + when using the `*_model` → `*Problem(::System, ...)` workflow. Works on any + `AbstractSystem` with alias maps in metadata. + +- **Aliases are handled through composition**: `compose(...; aliases=...)` accepts + cross-subsystem aliases, `extend` unions aliases from both systems, and `flatten` + collects aliases from subsystems with namespacing. + +- **`isequivalent` compares aliases** (set equality). + +- **Serialization of aliases is not yet supported** and will error if attempted. + +- **Non-symbolic (function-based) jumps** (`ConstantRateJump`/`VariableRateJump` with + plain Julia functions for rate/affect) are now rejected at `ReactionSystem` + construction time with an informative error. These were never supported by MTKBase. + + **Known limitations (v1):** + + The following alias types are rejected with an informative error: + - Brownian ↔ Brownian aliasing + - Poissonian ↔ Poissonian aliasing + - Compound species aliasing (either side) + - Aliases involving binding keys (either side) + - Aliases involving observable symbols directly + - Cross-category aliasing (e.g. species ↔ parameter, BC ↔ ordinary species) + - `MassActionJump` in systems with aliases (stoichiometry remapping deferred) + + The following metadata mismatches are rejected: + - Parameter type mismatch (e.g. `Int64` ↔ `Float64`) + - Parameter or species/variable unit mismatch + - Time-dependent ↔ non-time-dependent parameter mismatch + + Other known limitations: + - Serialization of systems with aliases errors (deferred alongside tstops/brownians/poissonians) + - `eliminate_aliases=false` is not supported for jump systems (algebraic constraints + cannot be carried by jump systems) + - When using the `*_model` → `*Problem(::System, ...)` workflow, users must call + `remap_alias_inputs` manually on `u0`/`p` maps. The automatic remapping only + applies to the `*Problem(::ReactionSystem, ...)` path. + - Substitution is unconditional per-expression (no expression-level early-exit guard); + mitigated by category-level fast paths when alias submaps are empty + - `system_to_reactionsystem` reverse conversion produces an empty `aliases` field + (the original alias declarations are cleared after elimination, though the + substitution maps survive in metadata for `remap_alias_inputs`) + ## Catalyst 16.1 - Added `use_jump_ratelaws` keyword argument to `ode_model`, `sde_model`, `hybrid_model`, diff --git a/docs/src/devdocs/dev_guide.md b/docs/src/devdocs/dev_guide.md index e57937f8ba..858e2be1d7 100644 --- a/docs/src/devdocs/dev_guide.md +++ b/docs/src/devdocs/dev_guide.md @@ -53,3 +53,57 @@ To build the Catalyst documentation locally: ### [Spellchecking in your code](@id devdocs_advice_codespellchecker) Especially when writing documentation, but also when writing normal code, it can be useful to have a spellchecker running through your texts. While code can be copied into a spellchecker and checked there (which is still useful to check grammar), it can also be very useful to (for users of VSCode) run the [Code Spell Checker](https://marketplace.visualstudio.com/items?itemName=streetsidesoftware.code-spell-checker) extension. This will automatically provide simple spell-checking for code and documentation as you write it. + +## [Adding a new field to `ReactionSystem`](@id devdocs_new_field) + +When adding a new field to the `ReactionSystem` struct, the following locations +must all be updated. Missing any of these can cause subtle bugs. + +### Core struct and constructors +- **`reactionsystem_fields` constant** (`src/reactionsystem.jl`) — add the field + name in the correct position. +- **Struct definition** (`src/reactionsystem.jl`) — add the field with docstring. +- **Inner constructor** (`src/reactionsystem.jl`) — add as positional argument and + in the `new{...}(...)` call. +- **5-argument constructor** (`src/reactionsystem.jl`) — add as keyword argument + with default, pass to inner constructor. +- **`make_ReactionSystem_internal`** (`src/reactionsystem.jl`) — add as keyword + argument, pass through to 5-arg constructor. + +### Composition and flattening +- **`flatten`** (`src/reactionsystem.jl`) — collect field from subsystems + (recursively if needed) and pass to constructor. +- **`compose`** (`src/reactionsystem.jl`) — add keyword argument if the field + should be settable at compose time; merge via `@set!`. +- **`extend`** (`src/reactionsystem.jl`) — union/merge field from both systems. + +### Accessors +- **`get_*` accessor** (`src/reactionsystem.jl`) — top-level only (uses `getfield`). +- **Recursive accessor** (`src/reactionsystem.jl`) — collects from subsystems with + namespacing (if applicable). +- **Exports** (`src/Catalyst.jl`) — export public accessors. + +### Equality and serialization +- **`isequivalent`** (`src/reactionsystem.jl`) — add comparison for the field. +- **`save_reactionsystem`** (`src/reactionsystem_serialisation/serialise_reactionsystem.jl`) + — add serialization support or an error guard if not yet supported. +- **`reactionsystem_uptodate_check`** — automatically covered by `reactionsystem_fields`. + +### Conversion pipeline +- **`hybrid_model`** (`src/reactionsystem_conversions.jl`) — pass field through to + `MT.System` constructor if applicable. +- **`ss_ode_model`** (`src/reactionsystem_conversions.jl`) — same (separate code path). +- **`sde_model` legacy path** (`src/reactionsystem_conversions.jl`) — same. +- **`eliminate_aliases`** (`src/alias_elimination.jl`) — substitute through the + field if it contains symbolic expressions. + +### DSL +- **`option_keys`** (`src/dsl.jl`) — add if the field has a DSL option. +- **`read_*_option`** (`src/dsl.jl`) — implement the option reader. +- **`make_reaction_system`** (`src/dsl.jl`) — wire into the generated code. + +### Other +- **`system_to_reactionsystem`** (`src/reactionsystem_conversions.jl`) — handle in + reverse conversion (or document as lost). +- **Tests** — add to appropriate test file. +- **HISTORY.md** — document for users. diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 2339490f58..0ec55a8cb4 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -27,7 +27,7 @@ import Symbolics: SymbolicT using Symbolics: iscall, sorted_arguments, value using ModelingToolkitBase: get_unknowns, get_ps, get_iv, get_systems, get_eqs, toparam, get_var_to_name, get_observed, - getvar, has_iv, JumpType + getvar, has_iv, JumpType, get_noise_eqs import ModelingToolkitBase: get_variables, namespace_expr, namespace_equation, modified_unknowns!, namespace_variables, @@ -108,6 +108,7 @@ const CatalystEqType = Union{Reaction, Equation} include("reactionsystem.jl") export ReactionSystem, isspatial export species, nonspecies, reactions, nonreactions, speciesmap, paramsmap +export aliases, has_aliases, aliases_present export numspecies, numreactions, numparams export make_empty_network export dependants, dependents, substoichmat, prodstoichmat, netstoichmat @@ -126,6 +127,10 @@ include("reactionsystem_metadata.jl") @public has_u0_map, get_u0_map, set_u0_map @public has_parameter_map, get_parameter_map, set_parameter_map +# Alias classification, validation, resolution, and elimination. +include("alias_elimination.jl") +export AliasClass + # Conversions of the `ReactionSystem` structure. include("reactionsystem_conversions.jl") export ODEProblem, SDEProblem, JumpProblem, NonlinearProblem, diff --git a/src/alias_elimination.jl b/src/alias_elimination.jl new file mode 100644 index 0000000000..68b24d1893 --- /dev/null +++ b/src/alias_elimination.jl @@ -0,0 +1,838 @@ +### Alias Elimination ### + +# This file contains the alias classification, validation, resolution, and elimination +# logic for ReactionSystems. Aliases declare that two symbols represent the same entity +# (e.g. species A is the same as species B). Alias elimination substitutes eliminated +# symbols throughout the system, producing a simplified system. + +### AliasClass EnumX ### + +""" + @enumx AliasClass + +Classification of a symbolic variable for alias validation purposes. +Used to enforce same-category aliasing and reject unsupported combinations. +""" +@enumx AliasClass begin + Species # ordinary species (non-BC, non-constant, non-compound) + BCSpecies # boundary condition species (in species/unknowns, has BC metadata) + ConstantSpecies # constant species (in parameters, has ParameterConstantSpecies metadata) + CompoundSpecies # compound species (components reference other species) + Unknown # non-species unknowns (variables) + Parameter # ordinary parameters (excludes constant species) + Brownian + Poissonian + Bound # binding key + Observable # LHS of an observed equation + Unsupported # anything else (e.g. independent variable) +end + +### Classification Context ### + +# Precomputed sets for efficient classification. Built once per validation/elimination call. +struct AliasClassificationContext + species_set::Set{SymbolicT} + unknown_set::Set{SymbolicT} + param_set::Set{SymbolicT} + brownian_set::Set{SymbolicT} + poissonian_set::Set{SymbolicT} + binding_keys::Set{SymbolicT} + obs_lhs_set::Set{SymbolicT} +end + +function _make_classification_context(rs::ReactionSystem) + AliasClassificationContext( + Set(get_species(rs)), + Set(get_unknowns(rs)), + Set(get_ps(rs)), + Set(MT.get_brownians(rs)), + Set(MT.get_poissonians(rs)), + Set(keys(MT.get_bindings(rs))), + Set(eq.lhs for eq in MT.get_observed(rs))) +end + +### Classification Function (flattened systems only) ### + +""" + alias_class(x, ctx::AliasClassificationContext) -> AliasClass.T + +Classify a symbolic variable for alias validation. Requires a flattened system +(uses non-recursive `get_*` accessors via the precomputed context). +Order matters — more specific checks come first. +""" +function alias_class(x, ctx::AliasClassificationContext) + x in ctx.binding_keys && return AliasClass.Bound + x in ctx.obs_lhs_set && return AliasClass.Observable + x in ctx.brownian_set && return AliasClass.Brownian + x in ctx.poissonian_set && return AliasClass.Poissonian + if x in ctx.species_set + isbc(x) && return AliasClass.BCSpecies + hasmetadata(x, CompoundSpecies) && return AliasClass.CompoundSpecies + return AliasClass.Species + end + x in ctx.unknown_set && return AliasClass.Unknown + if x in ctx.param_set + isconstant(x) && return AliasClass.ConstantSpecies + return AliasClass.Parameter + end + return AliasClass.Unsupported +end + +### Compatibility Predicates ### + +# V1 allowed pairings — must be identical class on both sides. +const _ALIAS_COMPATIBLE_CLASSES = Set([ + AliasClass.Species, + AliasClass.BCSpecies, + AliasClass.ConstantSpecies, + AliasClass.Unknown, + AliasClass.Parameter, +]) + +""" + alias_compatible(lhs_class::AliasClass.T, rhs_class::AliasClass.T) -> Bool + +Check whether two alias classes are compatible for aliasing in v1. +Both sides must be the same supported class. +""" +function alias_compatible(lhs_class::AliasClass.T, rhs_class::AliasClass.T) + return lhs_class == rhs_class && lhs_class in _ALIAS_COMPATIBLE_CLASSES +end + +""" + alias_compatibility_error(lhs, rhs, lhs_class, rhs_class) -> ErrorException + +Construct an informative error message for incompatible alias classes. +""" +function alias_compatibility_error(lhs, rhs, lhs_class::AliasClass.T, rhs_class::AliasClass.T) + # Same class but unsupported + if lhs_class == rhs_class + return ErrorException( + "Cannot alias $(lhs_class) variables: $lhs ~ $rhs. " * + "$(lhs_class) aliasing is not yet supported in Catalyst.") + end + # Different classes + return ErrorException( + "Cannot alias $lhs ($(lhs_class)) to $rhs ($(rhs_class)). " * + "Only same-type aliasing is allowed: Species↔Species, BCSpecies↔BCSpecies, " * + "ConstantSpecies↔ConstantSpecies, Unknown↔Unknown, Parameter↔Parameter.") +end + +### Metadata Compatibility ### + +""" + alias_metadata_compatible(lhs, rhs, alias_type::AliasClass.T) + +Check metadata compatibility for an alias pair. Throws informative error on mismatch. +V1 uses strict same-category matching — no metadata is ever propagated or mutated. +""" +function alias_metadata_compatible(lhs, rhs, alias_type::AliasClass.T) + # Unit check — applies to all supported types + _check_alias_units(lhs, rhs) + + if alias_type == AliasClass.Parameter || alias_type == AliasClass.ConstantSpecies + _check_alias_parameter_metadata(lhs, rhs) + end + nothing +end + +# Check unit metadata compatibility. +function _check_alias_units(lhs, rhs) + lhs_unit = MT.getmetadata(lhs, MT.VariableUnit, nothing) + rhs_unit = MT.getmetadata(rhs, MT.VariableUnit, nothing) + if lhs_unit !== nothing && rhs_unit !== nothing + lhs_unit == rhs_unit || error( + "Unit mismatch in alias $lhs ~ $rhs: " * + "$lhs has unit $lhs_unit but $rhs has unit $rhs_unit.") + elseif lhs_unit !== nothing || rhs_unit !== nothing + error("Unit mismatch in alias $lhs ~ $rhs: " * + "one side has units and the other does not. " * + "Both sides must have matching units or neither should have units.") + end +end + +# Check parameter-specific metadata: time-dependent, discrete, declared type. +function _check_alias_parameter_metadata(lhs, rhs) + # Time-dependent (called parameter): both must match + lhs_called = MT.iscalledparameter(lhs) + rhs_called = MT.iscalledparameter(rhs) + if lhs_called != rhs_called + error("Time-dependency mismatch in alias $lhs ~ $rhs: " * + "$(lhs_called ? "$lhs is time-dependent" : "$lhs is not time-dependent") but " * + "$(rhs_called ? "$rhs is time-dependent" : "$rhs is not time-dependent"). " * + "Both must match.") + end + + # Declared type: if both have type metadata, must match + lhs_type = MT.symtype(unwrap(lhs)) + rhs_type = MT.symtype(unwrap(rhs)) + if lhs_type != rhs_type + error("Type mismatch in alias $lhs ~ $rhs: " * + "$lhs has type $lhs_type but $rhs has type $rhs_type. Both must match.") + end +end + +### Constructor-Time Checks (Cheap, Structural Only) ### + +""" + check_aliases(aliases::Vector{Equation}) + +Cheap structural checks for aliases, called during constructor when `checks=true`. +Does NOT require a flattened system or symbol lookups. Full validation +(type classification, metadata, cycles) is deferred to `eliminate_aliases`. +""" +function check_aliases(aliases::Vector{Equation}) + seen_lhs = Set{Any}() + for eq in aliases + lhs, rhs = eq.lhs, eq.rhs + # Self-alias check + isequal(lhs, rhs) && error("Self-alias not allowed: $lhs ~ $rhs") + # Duplicate LHS check + lhs in seen_lhs && error("Duplicate alias: $lhs appears as LHS in multiple aliases") + push!(seen_lhs, lhs) + end +end + +### Alias Resolution ### + +# Visitor states for cycle detection. +@enum VisitState UNSEEN VISITING DONE + +# Result of alias resolution — substitution maps for unknowns and parameters. +struct AliasResolutionResult + unknown_submap::Dict{SymbolicT, SymbolicT} + param_submap::Dict{SymbolicT, SymbolicT} +end + +""" + build_alias_graphs(aliases, rs) -> (unknown_graph, param_graph) + +Build directed alias graphs from alias equations. Performs full validation: +type classification, compatibility, and metadata checks. Requires flattened system. +""" +function build_alias_graphs(aliases::Vector{Equation}, rs::ReactionSystem) + ctx = _make_classification_context(rs) + + unknown_graph = Dict{SymbolicT, SymbolicT}() + param_graph = Dict{SymbolicT, SymbolicT}() + seen_lhs = Set{SymbolicT}() + + for eq in aliases + lhs, rhs = eq.lhs, eq.rhs + + # Self-alias check + isequal(lhs, rhs) && error("Self-alias not allowed: $lhs ~ $rhs") + + # Duplicate LHS check + lhs in seen_lhs && error("Duplicate alias: $lhs appears as LHS in multiple aliases") + push!(seen_lhs, lhs) + + # Classify both sides + lhs_class = alias_class(lhs, ctx) + rhs_class = alias_class(rhs, ctx) + + # Compatibility check + alias_compatible(lhs_class, rhs_class) || + throw(alias_compatibility_error(lhs, rhs, lhs_class, rhs_class)) + + # Metadata compatibility + alias_metadata_compatible(lhs, rhs, lhs_class) + + # Insert into appropriate graph + if lhs_class in (AliasClass.Parameter, AliasClass.ConstantSpecies) + param_graph[lhs] = rhs + else + unknown_graph[lhs] = rhs + end + end + + return unknown_graph, param_graph +end + +""" + resolve_canonical!(graph, canonical, visit_state) + +Resolve alias chains to canonical representatives with cycle detection. +Uses iterative DFS with path compression. O(n) complexity. +""" +function resolve_canonical!( + graph::Dict{SymbolicT, SymbolicT}, + canonical::Dict{SymbolicT, SymbolicT}, + visit_state::Dict{SymbolicT, VisitState} +) + path = SymbolicT[] # reused across iterations + for start in keys(graph) + haskey(canonical, start) && continue + + empty!(path) + current = start + local canon + + while true + state = get(visit_state, current, UNSEEN) + + if state == DONE + canon = get(canonical, current, current) + break + elseif state == VISITING + cycle_start = findfirst(isequal(current), path) + cycle_syms = path[cycle_start:end] + error("Alias cycle detected involving: $cycle_syms") + else # UNSEEN + visit_state[current] = VISITING + push!(path, current) + + if haskey(graph, current) + current = graph[current] + else + canon = current + break + end + end + end + + # Path compression — only store eliminated nodes (not the canonical root itself) + for node in path + visit_state[node] = DONE + isequal(node, canon) || (canonical[node] = canon) + end + end +end + +""" + validate_and_resolve_aliases(aliases, rs) -> AliasResolutionResult + +Full validation and resolution of alias equations on a flattened system. +Returns substitution maps for unknowns and parameters, or `nothing` if empty. +""" +function validate_and_resolve_aliases(aliases::Vector{Equation}, rs::ReactionSystem) + isempty(aliases) && return nothing + + unknown_graph, param_graph = build_alias_graphs(aliases, rs) + + unknown_canonical = Dict{SymbolicT, SymbolicT}() + param_canonical = Dict{SymbolicT, SymbolicT}() + visit_state = Dict{SymbolicT, VisitState}() + + resolve_canonical!(unknown_graph, unknown_canonical, visit_state) + resolve_canonical!(param_graph, param_canonical, visit_state) + + return AliasResolutionResult(unknown_canonical, param_canonical) +end + +### Substitution Helpers ### + +""" + substitute_reaction(rx, unknown_submap, combined, has_unknown_aliases, + has_param_aliases, eliminated_unknown_set) + +Substitute aliases through a reaction. Returns `nothing` if the reaction becomes a no-op. +Uses cheap structural guards to skip work when possible. +""" +function substitute_reaction(rx::Reaction, unknown_submap, combined, + has_unknown_aliases, has_param_aliases, eliminated_unknown_set) + # Rate substitution is unconditional (expression-level guards deferred; see plan). + new_rate = Symbolics.substitute(rx.rate, combined) + + # Skip species list merging if no unknown aliases affect this reaction's species + if has_unknown_aliases && + ((!isnothing(rx.substrates) && any(in(eliminated_unknown_set), rx.substrates)) || + (!isnothing(rx.products) && any(in(eliminated_unknown_set), rx.products))) + new_subs, new_substoich = merge_species_list(rx.substrates, rx.substoich, unknown_submap) + new_prods, new_prodstoich = merge_species_list(rx.products, rx.prodstoich, unknown_submap) + else + new_subs, new_substoich = rx.substrates, rx.substoich + new_prods, new_prodstoich = rx.products, rx.prodstoich + end + + # No-op detection + if is_noop_reaction(new_subs, new_substoich, new_prods, new_prodstoich) + return nothing + end + + # Substitute through symbolic reaction metadata (e.g. noise_scaling) + new_metadata = isempty(combined) ? rx.metadata : + substitute_reaction_metadata(rx.metadata, combined) + + Reaction(new_rate, new_subs, new_prods, new_substoich, new_prodstoich; + only_use_rate = rx.only_use_rate, metadata = new_metadata) +end + +""" + merge_species_list(species_list, stoich_list, submap) + +Substitute and merge duplicate species in a substrate/product list. +Preserves first-seen order deterministically. +""" +function merge_species_list(species_list, stoich_list, submap) + isnothing(species_list) && return (nothing, nothing) + + merged_species = SymbolicT[] + merged_stoich = eltype(stoich_list)[] + seen_index = Dict{SymbolicT, Int}() + + for (sp, st) in zip(species_list, stoich_list) + canonical = get(submap, sp, sp) + idx = get(seen_index, canonical, 0) + if idx == 0 + push!(merged_species, canonical) + push!(merged_stoich, st) + seen_index[canonical] = length(merged_species) + else + merged_stoich[idx] += st + end + end + + return (merged_species, merged_stoich) +end + +""" + is_noop_reaction(subs, substoich, prods, prodstoich) + +Return `true` if a reaction has zero net stoichiometry after merging. +""" +function is_noop_reaction(subs, substoich, prods, prodstoich) + (isnothing(subs) || isempty(subs)) && (isnothing(prods) || isempty(prods)) && return true + + net = Dict{SymbolicT, Any}() + if !isnothing(subs) + for (sp, st) in zip(subs, substoich) + net[sp] = get(net, sp, 0) - st + end + end + if !isnothing(prods) + for (sp, st) in zip(prods, prodstoich) + net[sp] = get(net, sp, 0) + st + end + end + + return all(iszero, values(net)) +end + +""" + substitute_reaction_metadata(metadata, combined_submap) + +Substitute aliases through symbolic reaction metadata values (e.g. noise_scaling). +Non-symbolic values are passed through unchanged. +""" +function substitute_reaction_metadata(metadata, combined_submap) + return [key => (Symbolics.issym(val) || Symbolics.iscall(val) ? + Symbolics.substitute(val, combined_submap) : val) + for (key, val) in metadata] +end + +""" + substitute_equation(eq, combined_submap) + +Substitute aliases through both sides of an equation. +""" +function substitute_equation(eq::Equation, combined_submap) + new_lhs = Symbolics.substitute(eq.lhs, combined_submap) + new_rhs = Symbolics.substitute(eq.rhs, combined_submap) + return new_lhs ~ new_rhs +end + +# Substitute aliases through a SymbolicContinuousCallback. +function substitute_continuous_event(evt::MT.SymbolicContinuousCallback, combined_submap) + new_conds = [Symbolics.substitute(c, combined_submap) for c in evt.conditions] + new_affect = _substitute_affect(evt.affect, combined_submap) + new_affect_neg = _substitute_affect(evt.affect_neg, combined_submap) + MT.SymbolicContinuousCallback(new_conds, new_affect; affect_neg = new_affect_neg) +end + +# Substitute aliases through a SymbolicDiscreteCallback. +function substitute_discrete_event(evt::MT.SymbolicDiscreteCallback, combined_submap) + new_conds = Symbolics.substitute(evt.conditions, combined_submap) + new_affect = _substitute_affect(evt.affect, combined_submap) + MT.SymbolicDiscreteCallback(new_conds, new_affect) +end + +# Substitute through a SymbolicAffect (or nothing). +function _substitute_affect(aff, combined_submap) + isnothing(aff) && return nothing + new_eqs = [substitute_equation(eq, combined_submap) for eq in aff.affect] + # Remap discrete_parameters through alias submap (aliased discrete params stay stale otherwise). + new_disc_ps = [get(combined_submap, dp, dp) for dp in aff.discrete_parameters] + MT.SymbolicAffect(new_eqs; discrete_parameters = new_disc_ps) +end + +# Substitute through a jump. Per-type handling: +# - ConstantRateJump: substitute rate and affect! +# - VariableRateJump: use @set! to preserve all extra fields (save_positions, bounds, etc.) +# - MassActionJump: error guard (stoichiometry remapping deferred to future version) +function substitute_jump(j::ConstantRateJump, combined_submap) + new_rate = Symbolics.substitute(j.rate, combined_submap) + new_affect = Symbolics.substitute(j.affect!, combined_submap) + ConstantRateJump(new_rate, new_affect) +end + +function substitute_jump(j::VariableRateJump, combined_submap) + j_new = j + @set! j_new.rate = Symbolics.substitute(j.rate, combined_submap) + @set! j_new.affect! = Symbolics.substitute(j.affect!, combined_submap) + return j_new +end + +function substitute_jump(j::MassActionJump, combined_submap) + error("Alias elimination through MassActionJump is not yet supported. " * + "Please eliminate aliases before adding MassActionJumps, or use " * + "ConstantRateJump/VariableRateJump instead.") +end + +""" + substitute_variable_defaults(syms, combined_submap) + +Substitute aliases through VariableDefaultValue metadata on surviving symbols. +The MTKBase System constructor re-extracts VariableDefaultValue via process_variables!, +so stale references must be fixed. Only creates new symbol objects when the default +actually references an eliminated symbol (isequal short-circuit). +""" +function substitute_variable_defaults(syms::Vector{SymbolicT}, combined_submap) + return map(syms) do sym + MT.hasdefault(sym) || return sym + old_def = MT.getdefault(sym) + new_def = Symbolics.substitute(old_def, combined_submap) + isequal(old_def, new_def) && return sym + return MT.setdefault(sym, new_def) + end +end + +### Initial Conditions Transfer ### + +""" + transfer_initial_conditions(ics, unknown_submap, param_submap, combined) + +Transfer initial conditions from eliminated symbols to canonical ones. +Errors on conflicting values. Also substitutes through values. +""" +function transfer_initial_conditions(ics, unknown_submap, param_submap, combined) + new_ics = copy(ics) + + # Phase 1: Remap keys (eliminated → canonical) with conflict detection + for submap in (unknown_submap, param_submap) + for (eliminated, canonical) in submap + elim_val = get(new_ics, eliminated, nothing) + canon_val = get(new_ics, canonical, nothing) + + if elim_val !== nothing && canon_val !== nothing + if !isequal(elim_val, canon_val) + error("Initial condition conflict for alias $eliminated ~ $canonical: " * + "$eliminated has value $elim_val but $canonical has value $canon_val. " * + "These must agree or only one should be specified.") + end + elseif elim_val !== nothing + new_ics[canonical] = elim_val + end + + delete!(new_ics, eliminated) + end + end + + # Phase 2: Substitute through values (they can be symbolic expressions + # referencing aliased parameters, e.g. [B => 2*k1] with k1 aliased) + for (key, val) in new_ics + new_ics[key] = Symbolics.substitute(val, combined) + end + + return new_ics +end + +### Alias Map Metadata Storage ### + +""" + AliasSubstitutionMaps + +Metadata key type for storing alias substitution maps on a system. +Following the pattern from `reactionsystem_metadata.jl` (U0Map, ParameterMap). +""" +struct AliasSubstitutionMaps end + +function store_alias_maps_in_metadata(metadata, unknown_submap, param_submap) + # MetadataT is ImmutableDict{DataType, Any} — extend by constructing a new one. + return Base.ImmutableDict(metadata, + AliasSubstitutionMaps => (; + unknown_submap = Dict(unknown_submap), + param_submap = Dict(param_submap))) +end + +""" + get_alias_maps(sys) + +Retrieve alias substitution maps from system metadata, or `nothing` if not present. +""" +function get_alias_maps(sys) + MT.getmetadata(sys, AliasSubstitutionMaps, nothing) +end + +### Main Elimination Function ### + +""" + eliminate_aliases(rs::ReactionSystem; name = nameof(rs)) + +Return a new `ReactionSystem` with all alias relationships resolved by substitution. +Eliminated unknowns become observables. Eliminated parameters are removed. +Requires a flattened system (no subsystems). + +Supported alias types (v1): ordinary species, BC-species, constant species, +non-species unknowns (variables), ordinary parameters. + +Not supported: Brownians, Poissonians, compounds, bindings, observables. +""" +function eliminate_aliases(rs::ReactionSystem; name = nameof(rs)) + # Preconditions + !isempty(get_systems(rs)) && + error("System must be flattened before alias elimination. Call `flatten` first.") + alias_eqs = get_aliases(rs) + isempty(alias_eqs) && return rs + + # Validate and resolve + result = validate_and_resolve_aliases(alias_eqs, rs) + unknown_submap = result.unknown_submap + param_submap = result.param_submap + combined = merge(unknown_submap, param_submap) + + # Pre-build sets for cheap membership checks (avoids per-expression allocation). + has_unknown_aliases = !isempty(unknown_submap) + has_param_aliases = !isempty(param_submap) + eliminated_unknown_set = has_unknown_aliases ? Set(keys(unknown_submap)) : Set{SymbolicT}() + + # Substitute through reactions + new_rxs = Reaction[] + for rx in get_rxs(rs) + new_rx = substitute_reaction(rx, unknown_submap, combined, + has_unknown_aliases, has_param_aliases, eliminated_unknown_set) + new_rx === nothing && continue + push!(new_rxs, new_rx) + end + + # Substitute through non-reaction equations (skip if no aliases touch expressions) + new_noneq_eqs = Equation[] + for eq in get_eqs(rs) + eq isa Reaction && continue + push!(new_noneq_eqs, substitute_equation(eq, combined)) + end + new_eqs = CatalystEqType[new_rxs; new_noneq_eqs] + + # Rewrite existing observables (RHS only) + new_obs = Equation[] + for obs_eq in MT.get_observed(rs) + new_rhs = Symbolics.substitute(obs_eq.rhs, combined) + push!(new_obs, obs_eq.lhs ~ new_rhs) + end + + # Add eliminated unknowns as new observables (pointing to canonical) + for (eliminated, canonical) in unknown_submap + push!(new_obs, eliminated ~ canonical) + end + + # Substitute through events (skip if empty) + cevents = MT.get_continuous_events(rs) + new_cevents = isempty(cevents) ? cevents : + [substitute_continuous_event(evt, combined) for evt in cevents] + devents = MT.get_discrete_events(rs) + new_devents = isempty(devents) ? devents : + [substitute_discrete_event(evt, combined) for evt in devents] + + # Substitute through jumps (skip if empty) + old_jumps = MT.get_jumps(rs) + new_jumps = isempty(old_jumps) ? old_jumps : + [substitute_jump(j, combined) for j in old_jumps] + + # Substitute through tstops (skip if empty or no param aliases) + old_tstops = MT.get_tstops(rs) + new_tstops = (isempty(old_tstops) || !has_param_aliases) ? collect(old_tstops) : + [Symbolics.substitute(ts, combined) for ts in old_tstops] + + # Handle initial_conditions + new_ics = transfer_initial_conditions( + MT.initial_conditions(rs), unknown_submap, param_submap, combined) + + # Rebuild bindings with substituted values (skip if empty) + old_bindings = MT.get_bindings(rs) + if isempty(old_bindings) + new_bindings_dict = SymmapT() + else + new_bindings_dict = SymmapT() + for (key, val) in old_bindings + new_bindings_dict[key] = Symbolics.substitute(val, combined) + end + end + + # Remove eliminated symbols + new_unknowns = has_unknown_aliases ? + filter(u -> u ∉ eliminated_unknown_set, get_unknowns(rs)) : + collect(get_unknowns(rs)) + eliminated_params = has_param_aliases ? Set(keys(param_submap)) : Set{SymbolicT}() + new_params = has_param_aliases ? + filter(p -> p ∉ eliminated_params, get_ps(rs)) : + collect(get_ps(rs)) + + # Substitute through VariableDefaultValue on surviving symbols + new_unknowns = substitute_variable_defaults(new_unknowns, combined) + new_params = substitute_variable_defaults(new_params, combined) + + # Store alias maps in metadata for downstream input remapping + new_metadata = store_alias_maps_in_metadata( + MT.get_metadata(rs), unknown_submap, param_submap) + + # Construct new system with cleared aliases + new_rs = ReactionSystem(new_eqs, get_iv(rs), new_unknowns, new_params, MT.get_brownians(rs); + poissonians = MT.get_poissonians(rs), + jumps = new_jumps, + observed = new_obs, + name, + bindings = new_bindings_dict, + initial_conditions = new_ics, + checks = false, + combinatoric_ratelaws = get_combinatoric_ratelaws(rs), + balanced_bc_check = false, + spatial_ivs = get_sivs(rs), + continuous_events = new_cevents, + discrete_events = new_devents, + tstops = new_tstops, + metadata = new_metadata, + aliases = Equation[]) + + # Preserve completeness status from the input system. + if MT.iscomplete(rs) + new_rs = complete(new_rs) + end + return new_rs +end + +### Input Remapping ### + +""" + remap_alias_inputs(u0_or_p, sys) + +Remap eliminated symbol keys in user-provided `u0` or `p` maps to their canonical +equivalents. Works on any system with alias maps in metadata. +Handles both `Dict` and `Vector{Pair}` inputs (preserving format and ordering). +Errors if both eliminated and canonical keys are provided with different values. +""" +function remap_alias_inputs(u0_or_p, sys) + maps = get_alias_maps(sys) + isnothing(maps) && return u0_or_p + + combined = merge(maps.unknown_submap, maps.param_submap) + isempty(combined) && return u0_or_p + + return _remap_inputs(u0_or_p, combined) +end + +# No-op for NullParameters (common default for p in problem constructors). +remap_alias_inputs(p::DiffEqBase.NullParameters, sys) = p + +# Dict path: uses get/setindex!/delete! directly. +function _remap_inputs(d::AbstractDict, combined) + result = copy(d) + for (key, val) in d + canonical = get(combined, key, nothing) + isnothing(canonical) && continue + _check_remap_conflict(result, key, canonical, val) + result[canonical] = val + delete!(result, key) + end + return result +end + +# Vector{Pair} path: rebuild pairs list preserving stable ordering. +# Two-pass: first collect all values by canonical key, then emit deduplicated pairs. +function _remap_inputs(pairs::AbstractVector{<:Pair}, combined) + # First pass: collect values keyed by canonical target, detect conflicts. + canonical_vals = Dict{Any, Any}() + for (key, val) in pairs + canon = get(combined, key, nothing) + target = isnothing(canon) ? key : canon + prev = get(canonical_vals, target, nothing) + if prev !== nothing && !isequal(prev, val) + error("Conflict in input map: $key has value $val " * + "but $target already has value $prev. " * + "These must agree or only one should be specified.") + end + canonical_vals[target] = val + end + + # Second pass: emit pairs in first-seen order, skipping duplicates. + result = similar(pairs, 0) + sizehint!(result, length(canonical_vals)) + emitted = Set{Any}() + for (key, _) in pairs + canon = get(combined, key, nothing) + target = isnothing(canon) ? key : canon + if target ∉ emitted + push!(result, target => canonical_vals[target]) + push!(emitted, target) + end + end + return result +end + +# Helper: check for conflict when remapping a Dict entry. +function _check_remap_conflict(result::AbstractDict, key, canonical, val) + existing = get(result, canonical, nothing) + if existing !== nothing && !isequal(existing, val) + error("Conflict in input map: eliminated symbol $key has value $val " * + "but canonical symbol $canonical has value $existing. " * + "These must agree or only one should be specified.") + end +end + +### Conversion Pipeline Helper ### + +""" + prepare_aliases_for_conversion(flatrs; eliminate_aliases=true, allow_constraints=false) + +Process aliases on a flattened `ReactionSystem` for conversion. Hot path: when no aliases +are present, returns immediately with zero overhead (single `isempty` check). + +When `eliminate_aliases=true` (default): calls `eliminate_aliases(flatrs)`. +When `eliminate_aliases=false` and `allow_constraints=true`: materializes unknown aliases +as algebraic constraints (`lhs ~ rhs`) via `@set!`. Parameter aliases always error. +When `eliminate_aliases=false` and `allow_constraints=false` (jump path): errors. +""" +function prepare_aliases_for_conversion(flatrs::ReactionSystem; + eliminate_aliases::Bool = true, + allow_constraints::Bool = false) + !aliases_present(flatrs) && return flatrs + + if eliminate_aliases + return Catalyst.eliminate_aliases(flatrs) + end + + # Not eliminating — check constraints + if !isempty(get_parameter_aliases(flatrs)) + error("Parameter aliases must be eliminated before conversion. " * + "Use `eliminate_aliases=true` (default).") + end + + if !allow_constraints + error("This conversion path does not support algebraic constraints from " * + "uneliminated aliases. Use `eliminate_aliases=true` (default).") + end + + # Materialize unknown aliases as algebraic constraints (lhs ~ rhs kept as-is). + # Only two fields change, so use @set! instead of full reconstruction. + result = flatrs + @set! result.eqs = CatalystEqType[get_eqs(flatrs); get_unknown_aliases(flatrs)] + @set! result.aliases = Equation[] + return result +end + +# Helper: filter aliases to just unknown (species/variable) aliases. +function get_unknown_aliases(rs::ReactionSystem) + ctx = _make_classification_context(rs) + filter(get_aliases(rs)) do eq + cls = alias_class(eq.lhs, ctx) + cls in (AliasClass.Species, AliasClass.BCSpecies, AliasClass.Unknown) + end +end + +# Helper: filter aliases to just parameter aliases. +function get_parameter_aliases(rs::ReactionSystem) + ctx = _make_classification_context(rs) + filter(get_aliases(rs)) do eq + cls = alias_class(eq.lhs, ctx) + cls in (AliasClass.Parameter, AliasClass.ConstantSpecies) + end +end diff --git a/src/dsl.jl b/src/dsl.jl index 88a2c8b6c3..203909a8e7 100644 --- a/src/dsl.jl +++ b/src/dsl.jl @@ -10,7 +10,8 @@ const pure_rate_arrows = Set{Symbol}([:(=>), :(<=), :⇐, :⟽, :⇒, :⟾, :⇔ # Declares the keys used for various options. const option_keys = (:species, :parameters, :variables, :discretes, :ivs, :compounds, :observables, :default_noise_scaling, :differentials, :equations, :continuous_events, :discrete_events, - :tstops, :brownians, :poissonians, :combinatoric_ratelaws, :require_declaration, :unit_checks) + :tstops, :brownians, :poissonians, :combinatoric_ratelaws, :require_declaration, :unit_checks, + :aliases) ### `@species` Macro ### @@ -291,6 +292,7 @@ function make_reaction_system(ex::Expr, name) default_reaction_metadata = read_default_noise_scaling_option(options) combinatoric_ratelaws = read_combinatoric_ratelaws_option(options) unit_checks = read_unit_checks_option(options) + aliases_expr = read_aliases_option(options) # Creates expressions corresponding to actual code from the internal DSL representation. psexpr_init = get_psexpr(ps_inferred, stoich_ps, options) @@ -334,13 +336,14 @@ function make_reaction_system(ex::Expr, name) _tstops = $tstops_expr _combinatoric_ratelaws = $combinatoric_ratelaws _default_reaction_metadata = $default_reaction_metadata + _aliases = $aliases_expr remake_ReactionSystem_internal( make_ReactionSystem_internal(rx_eq_vec, $tiv, us, ps, $brownsvar; poissonians = $poissvar, name, spatial_ivs, observed = _observed, continuous_events = _continuous_events, discrete_events = _discrete_events, tstops = _tstops, combinatoric_ratelaws = _combinatoric_ratelaws, - unit_checks = _unit_checks); + unit_checks = _unit_checks, aliases = _aliases); default_reaction_metadata = _default_reaction_metadata) end)) end @@ -1066,6 +1069,26 @@ function read_unit_checks_option(options) end end +# Reads the `@aliases` option. Parses `A ~ B` equations and generates an expression +# that creates a vector of Equation objects at runtime. +function read_aliases_option(options) + !haskey(options, :aliases) && return :(Equation[]) + block = option_block_form(get_block_option(options[:aliases])) + alias_exprs = Expr[] + for line in block.args + line isa LineNumberNode && continue + # Each alias must be a binary ~ expression + if line isa Expr && line.head == :call && length(line.args) == 3 && line.args[1] == :(~) + lhs = line.args[2] + rhs = line.args[3] + push!(alias_exprs, :($lhs ~ $rhs)) + else + error("Invalid alias syntax: `$line`. Expected form: `A ~ B`.") + end + end + return :(Equation[$(alias_exprs...)]) +end + ### `@reaction` Macro & its Internals ### """ diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 8db6feb94c..bb34ff7216 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -253,7 +253,7 @@ all such code and updating it appropriately (e.g. serialization). Please use a s # structure have been updated (in the `reactionsystem_uptodate_check` function). const reactionsystem_fields = ( :eqs, :rxs, :iv, :sivs, :unknowns, :species, :ps, :var_to_name, - :observed, :name, :systems, :bindings, :initial_conditions, + :observed, :aliases, :name, :systems, :bindings, :initial_conditions, :networkproperties, :combinatoric_ratelaws, :continuous_events, :discrete_events, :tstops, :brownians, :poissonians, :jumps, :metadata, :complete, :parent) @@ -321,6 +321,8 @@ struct ReactionSystem{V <: NetworkProperties} <: MT.AbstractSystem var_to_name::Dict{Symbol, SymbolicT} """Equations for observed variables.""" observed::Vector{Equation} + """Alias equations declaring symbol equivalences. `lhs ~ rhs` means lhs is eliminated in favor of rhs.""" + aliases::Vector{Equation} """The name of the system""" name::Symbol """Internal sub-systems""" @@ -380,9 +382,9 @@ struct ReactionSystem{V <: NetworkProperties} <: MT.AbstractSystem # inner constructor is considered private and may change between non-breaking releases. function ReactionSystem(eqs, rxs, iv, sivs, unknowns, spcs, ps, var_to_name, observed, - name, systems, bindings, initial_conditions, nps, cls, cevs, devs, tstops, - brownians, poissonians, jumps, metadata, complete = false, parent = nothing; - checks::Bool = true, unit_checks::Bool = false) + aliases, name, systems, bindings, initial_conditions, nps, cls, cevs, devs, + tstops, brownians, poissonians, jumps, metadata, complete = false, + parent = nothing; checks::Bool = true, unit_checks::Bool = false) # Structural checks (fast, always on by default). if checks && isempty(sivs) @@ -400,9 +402,25 @@ struct ReactionSystem{V <: NetworkProperties} <: MT.AbstractSystem (hasnode(is_species_diff, eq.lhs) || hasnode(is_species_diff, eq.rhs)) && error("An equation ($eq) contains a differential with respect to a species. This is currently not supported. If this is a functionality you require, please raise an issue on the Catalyst GitHub page and we can consider the best way to implement it.") end + # Cheap structural alias checks (no symbol lookups — full validation deferred to elimination). + if checks && !isempty(aliases) + check_aliases(aliases) + end + + # Reject non-symbolic (imperative/function-based) jumps — not supported by MTKBase. + if checks && !isempty(jumps) + for j in jumps + if (j isa ConstantRateJump || j isa VariableRateJump) && + (j.rate isa Function || j.affect! isa Function) + error("Non-symbolic (function-based) jumps are not supported in ReactionSystem. " * + "Use symbolic rate and affect expressions for ConstantRateJump/VariableRateJump.") + end + end + end + rs = new{typeof(nps)}( eqs, rxs, iv, sivs, unknowns, spcs, ps, var_to_name, observed, - name, systems, bindings, initial_conditions, nps, cls, cevs, + aliases, name, systems, bindings, initial_conditions, nps, cls, cevs, devs, tstops, brownians, poissonians, jumps, metadata, complete, parent) unit_checks && assert_valid_units(rs; info = string("ReactionSystem constructor for ", name)) rs @@ -426,6 +444,7 @@ function ReactionSystem(eqs, iv, unknowns, ps, brownians = SymbolicT[]; poissonians = SymbolicT[], jumps = JumpType[], observed = Equation[], + aliases = Equation[], systems = [], name = nothing, bindings = SymmapT(), @@ -540,8 +559,8 @@ function ReactionSystem(eqs, iv, unknowns, ps, brownians = SymbolicT[]; tstops′ = Any[Symbolics.value(ts) for ts in tstops] ReactionSystem( - eqs′, rxs, iv′, sivs′, unknowns′, spcs, ps′, var_to_name, observed, name, - systems, bindings, initial_conditions, nps, combinatoric_ratelaws, + eqs′, rxs, iv′, sivs′, unknowns′, spcs, ps′, var_to_name, observed, aliases, + name, systems, bindings, initial_conditions, nps, combinatoric_ratelaws, continuous_events, discrete_events, tstops′, brownians′, poissonians′, jumps′, metadata; checks, unit_checks) end @@ -591,7 +610,8 @@ end # but carried out at a later stage. function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in, brownians = SymbolicT[]; poissonians = SymbolicT[], spatial_ivs = nothing, continuous_events = [], - discrete_events = [], tstops = [], observed = [], jumps = JumpType[], kwargs...) + discrete_events = [], tstops = [], observed = [], jumps = JumpType[], + aliases = Equation[], kwargs...) # Error if any observables have been declared a species or variable obs_vars = Set(obs_eq.lhs for obs_eq in observed) @@ -657,6 +677,11 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in, bro MT.collect_vars!(us, ps, eq.rhs, t) end + # Discover parameters/unknowns from alias equations. + for eq in aliases + MT.collect_vars!(us, ps, eq, t) + end + # Converts the found unknowns and parameters to vectors. usv = collect(us) @@ -679,7 +704,7 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in, bro # Passes the processed input into the next `ReactionSystem` call. # Note: brownians are passed as the 5th positional argument. ReactionSystem(fulleqs, t, usv, psv, brownians; poissonians, spatial_ivs, - continuous_events, discrete_events, tstops, observed, jumps, kwargs...) + continuous_events, discrete_events, tstops, observed, jumps, aliases, kwargs...) end ### Base Function Dispatches ### @@ -780,6 +805,8 @@ function isequivalent(rn1::ReactionSystem, rn2::ReactionSystem; ignorenames = tr debug_comparer( issetequal, MT.get_observed(rn1), MT.get_observed(rn2), "observed"; debug) || return false + debug_comparer(issetequal, get_aliases(rn1), get_aliases(rn2), "aliases"; debug) || + return false debug_comparer(issetequal, get_eqs(rn1), get_eqs(rn2), "eqs"; debug) || return false # Use custom event comparison functions to work around MTK issue #3907 debug_comparer(continuous_events_equal, MT.get_continuous_events(rn1), @@ -853,6 +880,44 @@ Return the system's `Reaction` vector (toplevel system only). get_rxs(sys::ReactionSystem) = getfield(sys, :rxs) has_rxs(sys::ReactionSystem) = isdefined(sys, :rxs) +""" + get_aliases(sys::ReactionSystem) + +Return the alias equations for the toplevel system only. +""" +get_aliases(sys::ReactionSystem) = getfield(sys, :aliases) + +""" + aliases(sys::ReactionSystem) + +Return all alias equations, recursively collecting from subsystems with namespacing. +""" +function aliases(sys::ReactionSystem) + systems = get_systems(sys) + isempty(systems) && return get_aliases(sys) + alias_eqs = copy(get_aliases(sys)) + for subsys in systems + for eq in aliases(subsys) + push!(alias_eqs, MT.renamespace(nameof(subsys), eq.lhs) ~ MT.renamespace(nameof(subsys), eq.rhs)) + end + end + alias_eqs +end + +""" + has_aliases(sys::ReactionSystem) + +Return `true` if the system type supports alias equations. Always `true` for `ReactionSystem`. +""" +has_aliases(sys::ReactionSystem) = true + +""" + aliases_present(sys::ReactionSystem) + +Return `true` if the system or any of its subsystems contains alias equations. +""" +aliases_present(sys::ReactionSystem) = !isempty(get_aliases(sys)) || any(aliases_present, get_systems(sys)) + """ get_sivs(sys::ReactionSystem) @@ -1535,6 +1600,7 @@ function MT.flatten(rs::ReactionSystem; name = nameof(rs)) poissonians = MT.poissonians(rs), jumps = MT.jumps(rs), observed = MT.observed(rs), + aliases = aliases(rs), name, initial_conditions = MT.initial_conditions(rs), checks = false, @@ -1565,7 +1631,8 @@ Notes: - By default, the new `ReactionSystem` will have the same name as `sys`. - Brownians and jumps from subsystems are collected at flatten time via recursive accessors. """ -function MT.compose(sys::ReactionSystem, systems::AbstractArray; name = nameof(sys)) +function MT.compose(sys::ReactionSystem, systems::AbstractArray; name = nameof(sys), + aliases = Equation[]) complete_check(sys, "MT.compose") foreach(s -> complete_check(s, "MT.compose"), systems) @@ -1595,6 +1662,11 @@ function MT.compose(sys::ReactionSystem, systems::AbstractArray; name = nameof(s @set! sys.ps = union(get_ps(sys), newparams) end + # Merge aliases (provided at composed-system namespace level, no renamespacing). + if !isempty(aliases) + @set! sys.aliases = [get_aliases(sys); aliases] + end + return sys end @@ -1633,15 +1705,17 @@ function MT.extend(sys::ReactionSystem, rs::ReactionSystem; Catalyst.get_combinatoric_ratelaws(rs) sivs = union(get_sivs(sys), get_sivs(rs)) - # Union brownians, poissonians, and jumps from both systems + # Union brownians, poissonians, jumps, and aliases from both systems new_brownians = union(MT.get_brownians(rs), MT.get_brownians(sys)) new_poissonians = union(MT.get_poissonians(rs), MT.get_poissonians(sys)) new_jumps = union(MT.get_jumps(rs), MT.get_jumps(sys)) + new_aliases = union(get_aliases(rs), get_aliases(sys)) ReactionSystem(eqs, t, sts, ps, collect(new_brownians); poissonians = collect(new_poissonians), jumps = collect(new_jumps), observed = obs, + aliases = collect(new_aliases), systems = syss, name, initial_conditions = defs, diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl index cba5a24b5b..b59846732e 100644 --- a/src/reactionsystem_conversions.jl +++ b/src/reactionsystem_conversions.jl @@ -631,6 +631,7 @@ function hybrid_model(rs::ReactionSystem; checks = false, initial_conditions = Dict(), use_jump_ratelaws = false, + eliminate_aliases = true, kwargs...) # Error checks. @@ -639,6 +640,17 @@ function hybrid_model(rs::ReactionSystem; flatrs = Catalyst.flatten(rs) + # Alias handling: eliminate or materialize as constraints. + # Hot path: returns immediately when no aliases present. + flatrs = prepare_aliases_for_conversion(flatrs; eliminate_aliases, + allow_constraints = true) + + # Conservation law interaction: must eliminate aliases first. + if remove_conserved && aliases_present(flatrs) + error("Cannot remove conserved species on a system with uneliminated aliases. " * + "Call `eliminate_aliases` first or use `eliminate_aliases=true`.") + end + # Resolve scales: _override_all_scales (internal) takes precedence over everything. if _override_all_scales !== nothing scales = fill(_override_all_scales, length(reactions(flatrs))) @@ -740,7 +752,7 @@ function hybrid_model(rs::ReactionSystem; continuous_events = MT.get_continuous_events(flatrs), discrete_events = MT.get_discrete_events(flatrs), tstops = MT.get_tstops(flatrs), - metadata = MT.get_metadata(rs), + metadata = MT.get_metadata(flatrs), kwargs...) end @@ -898,7 +910,7 @@ function ss_ode_model(rs::ReactionSystem; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), remove_conserved = false, checks = false, initial_conditions = Dict(), all_differentials_permitted = false, expand_catalyst_funs = true, - include_cl_as_eqs = false, kwargs...) + include_cl_as_eqs = false, eliminate_aliases = true, kwargs...) # Error checks. iscomplete(rs) || error(COMPLETENESS_ERROR) spatial_convert_err(rs::ReactionSystem, System) @@ -906,6 +918,14 @@ function ss_ode_model(rs::ReactionSystem; name = nameof(rs), # Generates system equations. fullrs = Catalyst.flatten(rs) + + # Alias handling. + fullrs = prepare_aliases_for_conversion(fullrs; eliminate_aliases, + allow_constraints = true) + if remove_conserved && aliases_present(fullrs) + error("Cannot remove conserved species on a system with uneliminated aliases. " * + "Call `eliminate_aliases` first or use `eliminate_aliases=true`.") + end remove_conserved && conservationlaws(fullrs) ists, ispcs = get_indep_sts(fullrs, (remove_conserved && !include_cl_as_eqs)) eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved, @@ -913,13 +933,13 @@ function ss_ode_model(rs::ReactionSystem; name = nameof(rs), eqs, us, ps, obs, ics, initeqs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved, compute_cl_initeqs = !include_cl_as_eqs, include_cl_as_eqs) - # Comoutes the correct initial conditions and bindings. - initial_conditions, bindings = MT.convert_bindings_for_time_independent_system(rs) + # Computes the correct initial conditions and bindings. + initial_conditions, bindings = MT.convert_bindings_for_time_independent_system(fullrs) initial_conditions = merge(initial_conditions, ics) # Throws a warning if there are differential equations in non-standard format. # Next, sets all differential terms to `0`. - all_differentials_permitted || nonlinear_convert_differentials_check(rs) + all_differentials_permitted || nonlinear_convert_differentials_check(fullrs) eqs = Equation[remove_diffs(eq.lhs) ~ remove_diffs(eq.rhs) for eq in eqs] System(eqs, us, ps; name, @@ -927,7 +947,7 @@ function ss_ode_model(rs::ReactionSystem; name = nameof(rs), bindings, initial_conditions, checks, - metadata = MT.get_metadata(rs), + metadata = MT.get_metadata(fullrs), kwargs...) end @@ -992,11 +1012,16 @@ function sde_model(rs::ReactionSystem; initial_conditions = Dict(), expand_catalyst_funs = true, use_legacy_noise = true, use_jump_ratelaws = false, + eliminate_aliases = true, kwargs...) # Flatten once upfront and check for constraints. flatrs = Catalyst.flatten(rs) + # Alias handling. + flatrs = prepare_aliases_for_conversion(flatrs; eliminate_aliases, + allow_constraints = true) + # Error if ReactionSystem has coupled jumps or poissonians (SDE + jumps hybrid not supported yet). if !isempty(MT.jumps(flatrs)) error("""Cannot convert ReactionSystem with coupled jumps to a pure SDE system. @@ -1016,8 +1041,8 @@ function sde_model(rs::ReactionSystem; # use legacy noise_eqs matrix approach (avoids mtkcompile overhead). # If user brownians are present, use the new Brownian-based path. if use_legacy_noise && !has_constraints && !has_user_brownians - iscomplete(rs) || error(COMPLETENESS_ERROR) - spatial_convert_err(rs, MT.System) + iscomplete(flatrs) || error(COMPLETENESS_ERROR) + spatial_convert_err(flatrs, MT.System) remove_conserved && conservationlaws(flatrs) ists, ispcs = get_indep_sts(flatrs, remove_conserved) @@ -1037,7 +1062,7 @@ function sde_model(rs::ReactionSystem; continuous_events = MT.get_continuous_events(flatrs), discrete_events = MT.get_discrete_events(flatrs), tstops = MT.get_tstops(flatrs), - metadata = MT.get_metadata(rs), + metadata = MT.get_metadata(flatrs), kwargs...) else # New path: Brownians via hybrid_model (requires mtkcompile for SDEProblem). @@ -1127,7 +1152,7 @@ function jump_model(rs::ReactionSystem; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), remove_conserved = nothing, checks = false, initial_conditions = Dict(), expand_catalyst_funs = true, save_positions = (true, true), - physical_scales = nothing, kwargs...) + physical_scales = nothing, eliminate_aliases = true, kwargs...) (remove_conserved !== nothing) && throw(ArgumentError("Catalyst does not support removing conserved species when converting to jump Systems.")) @@ -1141,6 +1166,11 @@ function jump_model(rs::ReactionSystem; name = nameof(rs), Use `HybridProblem` instead, which will process poissonians via mtkcompile.""") end + # Handle aliases before checking equations or computing scales. + # Jump systems cannot carry algebraic constraints, so allow_constraints=false. + flatrs = prepare_aliases_for_conversion(flatrs; eliminate_aliases, + allow_constraints = false) + # Error on non-reaction ODE/algebraic/SDE equations (pure Jump only supports reactions). # This also catches brownians since they appear in SDE equations. non_rxn_eqs = filter(eq -> !(eq isa Reaction), equations(flatrs)) @@ -1178,19 +1208,23 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan, combinatoric_ratelaws = get_combinatoric_ratelaws(rs), include_zero_odes = true, remove_conserved = false, checks = false, expand_catalyst_funs = true, mtkcompile = false, - use_jump_ratelaws = false, kwargs...) + use_jump_ratelaws = false, eliminate_aliases = true, kwargs...) osys = ode_model(rs; name, combinatoric_ratelaws, include_zero_odes, checks, - remove_conserved, expand_catalyst_funs, use_jump_ratelaws) + remove_conserved, expand_catalyst_funs, use_jump_ratelaws, eliminate_aliases) # Handles potential differential algebraic equations (which requires `mtkcompile`). + # Check the processed osys (not the original rs) so alias-generated constraints are detected. if mtkcompile osys = MT.mtkcompile(osys) - elseif has_alg_equations(rs) - error("The input ReactionSystem has algebraic equations. This requires setting `mtkcompile = true` within `ODEProblem` call.") + elseif has_alg_equations(osys) + error("The system has algebraic equations (possibly from uneliminated aliases). " * + "This requires setting `mtkcompile = true` within the `ODEProblem` call.") else osys = complete(osys) end + u0 = remap_alias_inputs(u0, osys) + p = remap_alias_inputs(p, osys) prob_cond = (p isa DiffEqBase.NullParameters) ? u0 : merge(Dict{Any,Any}(u0), Dict{Any,Any}(p)) return ODEProblem(osys, prob_cond, tspan, args...; check_length, kwargs...) end @@ -1228,10 +1262,13 @@ function DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0, p = DiffEqBase.NullParameters(), args...; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), remove_conserved = false, checks = false, check_length = false, expand_catalyst_funs = true, - mtkcompile = false, all_differentials_permitted = false, kwargs...) + mtkcompile = false, all_differentials_permitted = false, + eliminate_aliases = true, kwargs...) nlsys = ss_ode_model(rs; name, combinatoric_ratelaws, checks, all_differentials_permitted, - remove_conserved, expand_catalyst_funs) + remove_conserved, expand_catalyst_funs, eliminate_aliases) nlsys = mtkcompile ? MT.mtkcompile(nlsys) : complete(nlsys) + u0 = remap_alias_inputs(u0, nlsys) + p = remap_alias_inputs(p, nlsys) prob_cond = (p isa DiffEqBase.NullParameters) ? u0 : merge(Dict{Any,Any}(u0), Dict{Any,Any}(p)) return NonlinearProblem(nlsys, prob_cond, args...; check_length, kwargs...) @@ -1244,30 +1281,19 @@ function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan, include_zero_odes = true, checks = false, check_length = false, remove_conserved = false, mtkcompile = false, expand_catalyst_funs = true, use_legacy_noise = true, - use_jump_ratelaws = false, kwargs...) + use_jump_ratelaws = false, eliminate_aliases = true, kwargs...) - # Flatten once upfront and pass to sde_model. - flatrs = Catalyst.flatten(rs) - sde_sys = sde_model(flatrs; name, combinatoric_ratelaws, expand_catalyst_funs, - include_zero_odes, checks, remove_conserved, use_legacy_noise, use_jump_ratelaws) - - # Determine if we need mtkcompile: - # - If user brownians are present, sde_model routes through hybrid_model which requires mtkcompile - # - If using Brownian-based approach (not legacy), mtkcompile extracts the noise matrix - # - If there are algebraic equations, mtkcompile handles structural simplification - # - If mtkcompile is requested explicitly - has_constraints = has_alg_equations(flatrs) || any(isbc, get_unknowns(flatrs)) - has_user_brownians = !isempty(MT.brownians(flatrs)) - needs_mtkcompile = mtkcompile || - has_alg_equations(flatrs) || - !use_legacy_noise || - has_user_brownians || - has_constraints + sde_sys = sde_model(rs; name, combinatoric_ratelaws, expand_catalyst_funs, + include_zero_odes, checks, remove_conserved, use_legacy_noise, use_jump_ratelaws, + eliminate_aliases) + u0 = remap_alias_inputs(u0, sde_sys) + p = remap_alias_inputs(p, sde_sys) prob_cond = (p isa DiffEqBase.NullParameters) ? u0 : merge(Dict{Any,Any}(u0), Dict{Any,Any}(p)) + needs_mtkcompile = mtkcompile || (get_noise_eqs(sde_sys) === nothing) if needs_mtkcompile - if !mtkcompile && has_alg_equations(flatrs) + if !mtkcompile && has_alg_equations(sde_sys) error("The input ReactionSystem has algebraic equations. This requires setting `mtkcompile = true` within `SDEProblem` call.") end sde_sys = MT.mtkcompile(sde_sys) @@ -1275,7 +1301,7 @@ function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan, else # Legacy path: complete + noise_rate_prototype sde_sys = complete(sde_sys) - p_matrix = zeros(length(get_unknowns(sde_sys)), numreactions(flatrs)) + p_matrix = zeros(size(get_noise_eqs(sde_sys))...) return SDEProblem(sde_sys, prob_cond, tspan, args...; check_length, noise_rate_prototype = p_matrix, kwargs...) end @@ -1321,10 +1347,14 @@ function JumpProcesses.JumpProblem(rs::ReactionSystem, u0, tspan, expand_catalyst_funs = true, save_positions = (true, true), checks = false, + eliminate_aliases = true, kwargs...) # Pure jump system - use HybridProblem for hybrid ODE+SDE+Jump systems. jsys = complete(jump_model(rs; name, combinatoric_ratelaws, checks, - expand_catalyst_funs, save_positions)) + expand_catalyst_funs, save_positions, eliminate_aliases)) + # Remap eliminated symbols in user inputs. + u0 = remap_alias_inputs(u0, jsys) + p = remap_alias_inputs(p, jsys) # Use Dict{Any,Any} to prevent type promotion during merge (MTK converts to this anyway). op = (p isa DiffEqBase.NullParameters) ? u0 : merge(Dict{Any,Any}(u0), Dict{Any,Any}(p)) @@ -1416,9 +1446,13 @@ function HybridProblem(rs::ReactionSystem, u0, tspan, checks = false, mtkcompile = false, use_jump_ratelaws = false, + eliminate_aliases = true, kwargs...) # Determine which scale types are present. flatrs = Catalyst.flatten(rs) + # Handle aliases before scale resolution (alias elimination may change reaction count). + flatrs = prepare_aliases_for_conversion(flatrs; eliminate_aliases, + allow_constraints = true) resolved_scales = merge_physical_scales(reactions(flatrs), physical_scales, default_scale) has_ode, has_sde, has_jump = detect_scale_types(resolved_scales) @@ -1434,8 +1468,11 @@ function HybridProblem(rs::ReactionSystem, u0, tspan, # Build the unified System from the flattened ReactionSystem. sys = hybrid_model(flatrs; name, physical_scales, default_scale, combinatoric_ratelaws, expand_catalyst_funs, save_positions, checks, - use_jump_ratelaws) + use_jump_ratelaws, eliminate_aliases) + # Remap eliminated symbols in user inputs. + u0 = remap_alias_inputs(u0, sys) + p = remap_alias_inputs(p, sys) # Build problem conditions (u0 + p merged). # Use Dict{Any,Any} to prevent type promotion during merge (MTK converts to this anyway). prob_cond = (p isa DiffEqBase.NullParameters) ? u0 : merge(Dict{Any,Any}(u0), Dict{Any,Any}(p)) @@ -1482,19 +1519,24 @@ function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0, check_length = false, name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), remove_conserved = false, include_zero_odes = true, checks = false, - expand_catalyst_funs = true, mtkcompile = false, kwargs...) + expand_catalyst_funs = true, mtkcompile = false, + eliminate_aliases = true, kwargs...) osys = ode_model(rs; name, combinatoric_ratelaws, include_zero_odes, checks, - remove_conserved, expand_catalyst_funs) + remove_conserved, expand_catalyst_funs, eliminate_aliases) # Handles potential differential algebraic equations (which requires `mtkcompile`). + # Check the processed osys (not the original rs) so alias-generated constraints are detected. if mtkcompile (osys = MT.mtkcompile(osys)) - elseif has_alg_equations(rs) - error("The input ReactionSystem has algebraic equations. This requires setting `mtkcompile = true` within `ODEProblem` call.") + elseif has_alg_equations(osys) + error("The system has algebraic equations (possibly from uneliminated aliases). " * + "This requires setting `mtkcompile = true` within the `SteadyStateProblem` call.") else osys = complete(osys) end + u0 = remap_alias_inputs(u0, osys) + p = remap_alias_inputs(p, osys) prob_cond = (p isa DiffEqBase.NullParameters) ? u0 : merge(Dict{Any,Any}(u0), Dict{Any,Any}(p)) return SteadyStateProblem(osys, prob_cond, args...; check_length, kwargs...) end diff --git a/src/reactionsystem_serialisation/serialise_reactionsystem.jl b/src/reactionsystem_serialisation/serialise_reactionsystem.jl index 4d8a9e7ac2..c5f7e14219 100644 --- a/src/reactionsystem_serialisation/serialise_reactionsystem.jl +++ b/src/reactionsystem_serialisation/serialise_reactionsystem.jl @@ -34,6 +34,11 @@ function save_reactionsystem(filename::String, rn::ReactionSystem; annotate = true, safety_check = true) # Error and warning checks. reactionsystem_uptodate_check() + if aliases_present(rn) + error("Serialization of ReactionSystems with aliases is not yet supported. " * + "Please eliminate aliases before saving, or track " * + "https://github.com/SciML/Catalyst.jl/issues for updates.") + end if !isempty(MT.get_tstops(rn)) || (!isempty(get_systems(rn)) && !isempty(MT.symbolic_tstops(rn))) error("Serialization of ReactionSystems with tstops is not yet supported. " * "Please remove tstops before saving, or track " * diff --git a/test/miscellaneous_tests/explicit_imports.jl b/test/miscellaneous_tests/explicit_imports.jl index b630a9dbba..910299b65b 100644 --- a/test/miscellaneous_tests/explicit_imports.jl +++ b/test/miscellaneous_tests/explicit_imports.jl @@ -4,13 +4,13 @@ using Catalyst, ExplicitImports, Test @testset "Explicit Imports" begin # Test that there are no implicit imports - # allow_unanalyzable is needed because PhysicalScale enum causes parsing issues + # allow_unanalyzable is needed because EnumX-generated submodules cause parsing issues @test check_no_implicit_imports(Catalyst; - allow_unanalyzable=(Catalyst, Catalyst.PhysicalScale)) === nothing + allow_unanalyzable=(Catalyst, Catalyst.PhysicalScale, Catalyst.AliasClass)) === nothing # Test that there are no stale explicit imports @test check_no_stale_explicit_imports(Catalyst; - allow_unanalyzable=(Catalyst, Catalyst.PhysicalScale)) === nothing + allow_unanalyzable=(Catalyst, Catalyst.PhysicalScale, Catalyst.AliasClass)) === nothing # Test that all explicit imports are from owner modules @test check_all_explicit_imports_via_owners(Catalyst) === nothing diff --git a/test/reactionsystem_core/alias_elimination.jl b/test/reactionsystem_core/alias_elimination.jl new file mode 100644 index 0000000000..6b1cf75df4 --- /dev/null +++ b/test/reactionsystem_core/alias_elimination.jl @@ -0,0 +1,507 @@ +# Tests for alias elimination functionality. +using Catalyst, Test +import ModelingToolkitBase as MT +import DiffEqBase +using JumpProcesses + +const t = Catalyst.default_t() + +### AliasClass and Validation Tests ### + +@testset "AliasClass enum" begin + @test AliasClass.Species isa AliasClass.T + @test AliasClass.BCSpecies isa AliasClass.T + @test AliasClass.ConstantSpecies isa AliasClass.T + @test AliasClass.CompoundSpecies isa AliasClass.T + @test AliasClass.Unknown isa AliasClass.T + @test AliasClass.Parameter isa AliasClass.T + @test AliasClass.Brownian isa AliasClass.T + @test AliasClass.Poissonian isa AliasClass.T + @test AliasClass.Bound isa AliasClass.T + @test AliasClass.Observable isa AliasClass.T + @test AliasClass.Unsupported isa AliasClass.T +end + +@testset "alias_compatible" begin + @test Catalyst.alias_compatible(AliasClass.Species, AliasClass.Species) + @test Catalyst.alias_compatible(AliasClass.BCSpecies, AliasClass.BCSpecies) + @test Catalyst.alias_compatible(AliasClass.ConstantSpecies, AliasClass.ConstantSpecies) + @test Catalyst.alias_compatible(AliasClass.Unknown, AliasClass.Unknown) + @test Catalyst.alias_compatible(AliasClass.Parameter, AliasClass.Parameter) + @test !Catalyst.alias_compatible(AliasClass.Species, AliasClass.Parameter) + @test !Catalyst.alias_compatible(AliasClass.Brownian, AliasClass.Brownian) + @test !Catalyst.alias_compatible(AliasClass.CompoundSpecies, AliasClass.CompoundSpecies) +end + +@testset "check_aliases (constructor checks)" begin + @species A(t) B(t) + @parameters k1 k2 + Catalyst.check_aliases(Equation[]) + Catalyst.check_aliases([A ~ B]) + Catalyst.check_aliases([A ~ B, k1 ~ k2]) + @test_throws ErrorException Catalyst.check_aliases([A ~ A]) + @test_throws ErrorException Catalyst.check_aliases([A ~ B, A ~ B]) +end + +### Alias Resolution Tests ### + +@testset "Alias resolution" begin + @species A(t) B(t) C(t) D(t) + @parameters k1 k2 + + # Basic species alias + @named rn = ReactionSystem([Reaction(k1, [A], [B])], t, [A, B], [k1, k2]; + aliases = [A ~ B]) + rn_c = complete(rn) + flat = Catalyst.flatten(rn_c) + result = Catalyst.validate_and_resolve_aliases(Catalyst.get_aliases(flat), flat) + @test isequal(result.unknown_submap[A], B) + + # Chain aliasing + @named rn_chain = ReactionSystem( + [Reaction(k1, [A], [C])], t, [A, B, C], [k1]; + aliases = [A ~ B, B ~ C]) + flat_chain = Catalyst.flatten(complete(rn_chain)) + result_chain = Catalyst.validate_and_resolve_aliases(Catalyst.get_aliases(flat_chain), flat_chain) + @test isequal(result_chain.unknown_submap[A], C) + @test isequal(result_chain.unknown_submap[B], C) + + # Fan-in + chain + @named rn_fan = ReactionSystem( + [Reaction(k1, [A], [C])], t, [A, B, C, D], [k1]; + aliases = [A ~ C, B ~ C, D ~ B]) + flat_fan = Catalyst.flatten(complete(rn_fan)) + result_fan = Catalyst.validate_and_resolve_aliases(Catalyst.get_aliases(flat_fan), flat_fan) + @test isequal(result_fan.unknown_submap[A], C) + @test isequal(result_fan.unknown_submap[B], C) + @test isequal(result_fan.unknown_submap[D], C) + + # Cycle detection + @named rn_cycle = ReactionSystem( + [Reaction(k1, [A], [C])], t, [A, B, C], [k1]; + aliases = [A ~ B, B ~ A]) + flat_cycle = Catalyst.flatten(complete(rn_cycle)) + @test_throws ErrorException Catalyst.validate_and_resolve_aliases( + Catalyst.get_aliases(flat_cycle), flat_cycle) + + # Parameter alias + @named rn_param = ReactionSystem( + [Reaction(k1, [A], [B])], t, [A, B], [k1, k2]; + aliases = [k1 ~ k2]) + flat_param = Catalyst.flatten(complete(rn_param)) + result_param = Catalyst.validate_and_resolve_aliases( + Catalyst.get_aliases(flat_param), flat_param) + @test isequal(result_param.param_submap[k1], k2) +end + +### Alias Elimination Tests ### + +@testset "Alias elimination" begin + @species A(t) B(t) C(t) + @parameters k1 k2 + + # Basic species elimination + @named rn = ReactionSystem( + [Reaction(k1, [A], [B]), Reaction(k2, [B], [C])], t, [A, B, C], [k1, k2]; + aliases = [A ~ B]) + flat = Catalyst.flatten(complete(rn)) + elim = Catalyst.eliminate_aliases(flat) + @test !any(isequal(A), Catalyst.get_unknowns(elim)) + @test any(eq -> isequal(eq.lhs, A), Catalyst.get_observed(elim)) + @test isempty(Catalyst.get_aliases(elim)) + + # Parameter elimination + @named rn_p = ReactionSystem( + [Reaction(k1, [A], [B])], t, [A, B], [k1, k2]; + aliases = [k1 ~ k2]) + flat_p = Catalyst.flatten(complete(rn_p)) + elim_p = Catalyst.eliminate_aliases(flat_p) + @test !any(isequal(k1), Catalyst.get_ps(elim_p)) + @test any(isequal(k2), Catalyst.get_ps(elim_p)) + + # No-op reaction removal (A→B with A~B becomes B→B, dropped) + @named rn_noop = ReactionSystem( + [Reaction(k1, [A], [B]), Reaction(k2, [B], [C])], t, [A, B, C], [k1, k2]; + aliases = [A ~ B]) + elim_noop = Catalyst.eliminate_aliases(Catalyst.flatten(complete(rn_noop))) + @test length(Catalyst.get_rxs(elim_noop)) == 1 + + # Stoichiometry merging: 2A + B → C with A~B becomes 3B → C + @named rn_merge = ReactionSystem( + [Reaction(k1, [A, B], [C], [2, 1], [1])], t, [A, B, C], [k1]; + aliases = [A ~ B]) + elim_merge = Catalyst.eliminate_aliases(Catalyst.flatten(complete(rn_merge))) + rx = first(Catalyst.get_rxs(elim_merge)) + @test length(rx.substrates) == 1 + @test isequal(rx.substrates[1], B) + @test rx.substoich[1] == 3 + + # Initial conditions transfer + @named rn_ic = ReactionSystem( + [Reaction(k1, [A], [B])], t, [A, B], [k1]; + aliases = [A ~ B], initial_conditions = Dict(A => 5.0)) + elim_ic = Catalyst.eliminate_aliases(Catalyst.flatten(complete(rn_ic))) + ics = MT.initial_conditions(elim_ic) + @test !haskey(ics, A) + @test haskey(ics, B) + + # Initial conditions conflict + @named rn_conflict = ReactionSystem( + [Reaction(k1, [A], [B])], t, [A, B], [k1]; + aliases = [A ~ B], initial_conditions = Dict(A => 5.0, B => 10.0)) + @test_throws ErrorException Catalyst.eliminate_aliases( + Catalyst.flatten(complete(rn_conflict))) + + # Alias maps in metadata + @named rn_meta = ReactionSystem( + [Reaction(k1, [A], [B])], t, [A, B], [k1, k2]; + aliases = [A ~ B, k1 ~ k2]) + elim_meta = Catalyst.eliminate_aliases(Catalyst.flatten(complete(rn_meta))) + maps = Catalyst.get_alias_maps(elim_meta) + @test maps !== nothing + @test isequal(maps.unknown_submap[A], B) + + # Input remapping + u0 = Dict(A => 1.0) + remapped = Catalyst.remap_alias_inputs(u0, elim_meta) + @test !haskey(remapped, A) + @test haskey(remapped, B) +end + +### DSL Tests ### + +@testset "DSL @aliases" begin + rn = @reaction_network begin + @aliases A ~ B + @species A(t) B(t) C(t) + @parameters k + k, A --> C + k, B --> C + end + @test aliases_present(rn) + @test length(Catalyst.get_aliases(rn)) == 1 + + rn2 = @reaction_network begin + @aliases begin + A ~ B + k1 ~ k2 + end + @species A(t) B(t) C(t) + @parameters k1 k2 + k1, A --> C + k2, B --> C + end + @test length(Catalyst.get_aliases(rn2)) == 2 +end + +### Conversion Pipeline Tests ### + +@testset "Conversion pipeline with aliases" begin + @species A(t) B(t) C(t) + @parameters k1 k2 + + @named rn = ReactionSystem( + [Reaction(k1, [A], [B]), Reaction(k2, [B], [C])], t, [A, B, C], [k1, k2]; + aliases = [A ~ B]) + rn_c = complete(rn) + + # Model conversions + @test ode_model(rn_c) isa MT.System + @test sde_model(rn_c) isa MT.System + @test jump_model(rn_c) isa MT.System + @test hybrid_model(rn_c; default_scale = PhysicalScale.ODE) isa MT.System + + # Problem construction + prob = ODEProblem(rn_c, [B => 1.0, C => 0.0], (0.0, 1.0), [k1 => 0.5, k2 => 0.3]) + @test prob isa ODEProblem +end + +### Compose / Extend / Flatten Tests ### + +@testset "Compose/extend/flatten with aliases" begin + @species A(t) B(t) C(t) D(t) + @parameters k1 k2 k3 + + # Flatten collects aliases + @named sub1 = ReactionSystem([Reaction(k1, [A], [B])], t, [A, B], [k1]; + aliases = [A ~ B]) + @named sub2 = ReactionSystem([Reaction(k2, [C], [D])], t, [C, D], [k2]) + @named parent = ReactionSystem(Reaction[], t; systems = [sub1, sub2]) + flat = Catalyst.flatten(parent) + @test length(Catalyst.get_aliases(flat)) == 1 + + # Compose with aliases keyword + @named sys1 = ReactionSystem([Reaction(k1, [A], [B])], t, [A, B], [k1]) + @named sys2 = ReactionSystem([Reaction(k2, [C], [D])], t, [C, D], [k2]) + composed = compose(sys1, [sys2]; aliases = [A ~ C]) + @test length(Catalyst.get_aliases(composed)) == 1 + + # Extend unions aliases + @named ext1 = ReactionSystem([Reaction(k1, [A], [B])], t, [A, B], [k1]; + aliases = [A ~ B]) + @named ext2 = ReactionSystem([Reaction(k2, [C], [D])], t, [C, D], [k2, k3]; + aliases = [k2 ~ k3]) + extended = extend(ext1, ext2) + @test length(Catalyst.get_aliases(extended)) == 2 +end + +### isequivalent Tests ### + +@testset "isequivalent with aliases" begin + @species A(t) B(t) + @parameters k1 k2 + + @named rn1 = ReactionSystem([Reaction(k1, [A], [B])], t, [A, B], [k1, k2]; + aliases = [A ~ B]) + @named rn2 = ReactionSystem([Reaction(k1, [A], [B])], t, [A, B], [k1, k2]; + aliases = [A ~ B]) + @test Catalyst.isequivalent(rn1, rn2) + + @named rn3 = ReactionSystem([Reaction(k1, [A], [B])], t, [A, B], [k1, k2]; + aliases = [k1 ~ k2]) + @test !Catalyst.isequivalent(rn1, rn3) + + @named rn4 = ReactionSystem([Reaction(k1, [A], [B])], t, [A, B], [k1, k2]) + @test !Catalyst.isequivalent(rn1, rn4) +end + +### Serialization Guard Test ### + +@testset "Serialization guard" begin + @species A(t) B(t) + @parameters k1 + @named rn = ReactionSystem([Reaction(k1, [A], [B])], t, [A, B], [k1]; + aliases = [A ~ B]) + @test_throws ErrorException Catalyst.save_reactionsystem(tempname() * ".jl", rn) +end + +### Input Remapping Tests ### + +@testset "remap_alias_inputs formats" begin + @species A(t) B(t) C(t) + @parameters k1 k2 + + @named rn = ReactionSystem( + [Reaction(k1, [A], [B])], t, [A, B], [k1, k2]; + aliases = [A ~ B, k1 ~ k2]) + elim = Catalyst.eliminate_aliases(Catalyst.flatten(complete(rn))) + + # Vector{Pair} input + u0_vec = [A => 1.0, C => 2.0] + remapped_vec = Catalyst.remap_alias_inputs(u0_vec, elim) + @test remapped_vec isa Vector + @test any(p -> isequal(p[1], B), remapped_vec) + @test any(p -> isequal(p[1], C), remapped_vec) + + # Dict input + u0_dict = Dict(A => 1.0, C => 2.0) + remapped_dict = Catalyst.remap_alias_inputs(u0_dict, elim) + @test remapped_dict isa Dict + @test haskey(remapped_dict, B) + @test !haskey(remapped_dict, A) + + # NullParameters no-op + @test Catalyst.remap_alias_inputs( + DiffEqBase.NullParameters(), elim) isa DiffEqBase.NullParameters + + # Conflict: Dict with both eliminated and canonical keys (different values) + @test_throws ErrorException Catalyst.remap_alias_inputs(Dict(A => 1.0, B => 5.0), elim) + + # Conflict: Vector{Pair} with both eliminated and canonical keys (different values) + @test_throws ErrorException Catalyst.remap_alias_inputs([A => 1.0, B => 5.0], elim) +end + +### Constructor-Level Remapping Tests ### + +@testset "Constructor remapping with eliminated symbols" begin + @species A(t) B(t) C(t) + @parameters k1 k2 + + @named rn = ReactionSystem( + [Reaction(k1, [A], [B]), Reaction(k2, [B], [C])], t, [A, B, C], [k1, k2]; + aliases = [A ~ B]) + rn_c = complete(rn) + + # ODEProblem with eliminated symbol A in u0 + prob = ODEProblem(rn_c, [A => 1.0, C => 0.0], (0.0, 1.0), [k1 => 0.5, k2 => 0.3]) + @test prob isa ODEProblem + + # ODEProblem with NullParameters (system with no free parameters) + @named rn_np = ReactionSystem( + [Reaction(1.0, [A], [B])], t, [A, B], []; + aliases = [A ~ B]) + prob_np = ODEProblem(complete(rn_np), [A => 1.0], (0.0, 1.0)) + @test prob_np isa ODEProblem +end + +### Non-Symbolic Jump Guard ### + +@testset "Non-symbolic jump constructor guard" begin + @species A(t) B(t) + @parameters k1 + # Imperative (function-based) CRJ should be rejected at construction + crj = ConstantRateJump((u, p, t) -> p[1], integrator -> (integrator.u[1] -= 1)) + @test_throws ErrorException ReactionSystem( + Reaction[], Catalyst.default_t(), [A, B], [k1]; jumps = [crj], name = :test) +end + +### MassActionJump + Aliases Guard ### + +@testset "MassActionJump + aliases error guard" begin + @species A(t) B(t) + @parameters k1 + maj = MassActionJump(k1, [A => 1], [B => 1]) + @test_throws ErrorException Catalyst.substitute_jump(maj, Dict(A => B)) +end + +### Additional Constructor-Level Remapping Tests ### + +@testset "SDEProblem with eliminated-symbol p" begin + @species A(t) B(t) + @parameters k1 k2 + # Use parameter alias (no species alias avoids no-op reaction / noise matrix size issue) + @named rn = ReactionSystem( + [Reaction(k1, [A], [B]), Reaction(k2, [B], [A])], t, [A, B], [k1, k2]; + aliases = [k1 ~ k2]) + prob = SDEProblem(complete(rn), [A => 1.0, B => 0.0], (0.0, 1.0), [k1 => 0.5]) + @test prob isa SDEProblem +end + +@testset "JumpProblem with eliminated-symbol u0" begin + @species A(t) B(t) C(t) + @parameters k1 k2 + @named rn = ReactionSystem( + [Reaction(k1, [A], [B]), Reaction(k2, [B], [C])], t, [A, B, C], [k1, k2]; + aliases = [A ~ B]) + prob = JumpProblem(complete(rn), [A => 10, C => 0], (0.0, 1.0), [k1 => 0.5, k2 => 0.3]) + @test prob isa JumpProblem +end + +@testset "HybridProblem with eliminated-symbol u0" begin + @species A(t) B(t) C(t) + @parameters k1 k2 + @named rn = ReactionSystem( + [Reaction(k1, [A], [B]), Reaction(k2, [B], [C])], t, [A, B, C], [k1, k2]; + aliases = [A ~ B]) + prob = HybridProblem(complete(rn), [A => 1.0, C => 0.0], (0.0, 1.0), + [k1 => 0.5, k2 => 0.3]; default_scale = PhysicalScale.ODE) + @test prob isa ODEProblem +end + +@testset "SteadyStateProblem with eliminated-symbol u0" begin + @species A(t) B(t) C(t) + @parameters k1 k2 + @named rn = ReactionSystem( + [Reaction(k1, [A], [B]), Reaction(k2, [B], [C])], t, [A, B, C], [k1, k2]; + aliases = [A ~ B]) + prob = SteadyStateProblem(complete(rn), [A => 1.0, C => 0.0], [k1 => 0.5, k2 => 0.3]) + @test prob isa SteadyStateProblem +end + +### eliminate_aliases=false Tests ### + +@testset "eliminate_aliases=false" begin + @species A(t) B(t) C(t) + @parameters k1 k2 + @named rn = ReactionSystem( + [Reaction(k1, [A], [B]), Reaction(k2, [B], [C])], t, [A, B, C], [k1, k2]; + aliases = [A ~ B]) + rn_c = complete(rn) + + # Without mtkcompile should error (alias-generated constraints need DAE handling) + @test_throws ErrorException ODEProblem(rn_c, [B => 1.0, C => 0.0], (0.0, 1.0), + [k1 => 0.5, k2 => 0.3]; eliminate_aliases = false) +end + +### Binding Value Substitution Test ### + +@testset "Binding value substitution through aliases" begin + @species A(t) B(t) + @parameters k1 k2 + + # Test initial_conditions value substitution + @named rn_ic = ReactionSystem( + [Reaction(k1, [A], [B])], t, [A, B], [k1, k2]; + aliases = [k1 ~ k2], + initial_conditions = Dict(A => 2 * k1)) + elim_ic = Catalyst.eliminate_aliases(Catalyst.flatten(complete(rn_ic))) + ics = MT.initial_conditions(elim_ic) + @test haskey(ics, A) + ic_vars = Symbolics.get_variables(ics[A]) + @test !any(isequal(k1), ic_vars) + @test any(isequal(k2), ic_vars) + + # Test actual bindings-map rebuilding: create a species with a symbolic default + # (which MTKBase extracts into bindings via process_variables!). + @species X(t) = 3 * k1 + @named rn_bind = ReactionSystem( + [Reaction(k1, [X], [B])], t, [X, B], [k1, k2]; + aliases = [k1 ~ k2]) + flat_bind = Catalyst.flatten(complete(rn_bind)) + elim_bind = Catalyst.eliminate_aliases(flat_bind) + bindings = MT.get_bindings(elim_bind) + # X's binding should exist and reference k2 (canonical), not k1 (eliminated) + @test haskey(bindings, X) + bind_vars = Symbolics.get_variables(bindings[X]) + @test !any(isequal(k1), bind_vars) + @test any(isequal(k2), bind_vars) +end + +### Discrete Event Alias Substitution Test ### + +@testset "Discrete event with aliased parameter" begin + @species A(t) B(t) + @parameters k1 k2 + + # Create a system with a discrete event that references an aliased parameter. + # The event affect should have k1 substituted to k2 after elimination. + devt = MT.SymbolicDiscreteCallback( + 1.0, [A ~ A + k1]) + @named rn = ReactionSystem( + [Reaction(k1, [A], [B])], t, [A, B], [k1, k2]; + aliases = [k1 ~ k2], + discrete_events = [devt]) + flat = Catalyst.flatten(complete(rn)) + elim = Catalyst.eliminate_aliases(flat) + + # After elimination, discrete events should reference k2, not k1 + devents = MT.get_discrete_events(elim) + @test length(devents) == 1 + # Check the affect equation: should be A ~ A + k2 (not k1) + affect_eqs = devents[1].affect.affect + affect_vars = union([Symbolics.get_variables(eq.rhs) for eq in affect_eqs]...) + @test !any(isequal(k1), affect_vars) + @test any(isequal(k2), affect_vars) +end + +### Alias Autodiscovery Tests ### + +@testset "Autodiscovery of species from aliases" begin + # B only appears in the alias, not in any reaction or the explicit species list. + # It should be autodiscovered as a species from the alias equation. + @species A(t) B(t) C(t) + @parameters k + @named rn = ReactionSystem( + [Reaction(k, [A], [C])]; + aliases = [A ~ B]) + rn = complete(Catalyst.flatten(rn)) + @test any(isequal(B), Catalyst.get_species(rn)) + elim = Catalyst.eliminate_aliases(rn) + @test !any(isequal(A), Catalyst.get_species(elim)) + @test any(isequal(B), Catalyst.get_species(elim)) +end + +@testset "Autodiscovery of parameters from aliases" begin + # k2 only appears in the alias, not in any reaction or the explicit parameter list. + # It should be autodiscovered as a parameter from the alias equation. + @species A(t) B(t) + @parameters k1 k2 + @named rn = ReactionSystem( + [Reaction(k1, [A], [B])]; + aliases = [k1 ~ k2]) + rn = complete(Catalyst.flatten(rn)) + @test any(isequal(k2), parameters(rn)) + elim = Catalyst.eliminate_aliases(rn) + @test !any(isequal(k1), parameters(elim)) + @test any(isequal(k2), parameters(elim)) +end diff --git a/test/runtests.jl b/test/runtests.jl index 6a02c8a964..5d0bfafd5b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -34,6 +34,7 @@ end # Tests compositional and hierarchical modelling. @time @safetestset "ReactionSystem Components Based Creation" begin include("compositional_modelling/component_based_model_creation.jl") end # hierarchical modelling broken due to https://github.com/SciML/ModelingToolkit.jl/pull/4101 @time @safetestset "Brownians, Jumps, and Poissonians Composition" begin include("compositional_modelling/brownians_and_jumps_composition.jl") end + @time @safetestset "Alias Elimination" begin include("reactionsystem_core/alias_elimination.jl") end # Tests various miscellaneous features. @time @safetestset "API" begin include("miscellaneous_tests/api.jl") end