"""Implements simulation of models."""
import os
import shutil
import tempfile
import random
import warnings
import csv
import subprocess
from collections import defaultdict
from dataclasses import dataclass
from functools import cached_property
from types import MappingProxyType
from typing import Optional, Iterable, Mapping, Self, TYPE_CHECKING
import cloudpickle
if TYPE_CHECKING:
from graphviz import Source
from pykappa.mixture import Mixture
from pykappa.rule import Rule
from pykappa.pattern import Component, Pattern, Site
from pykappa.analysis import Monitor
from pykappa.expression import Expression
from pykappa._utils import str_table
[docs]
@dataclass(frozen=True)
class RuleTally:
"""Counts outcomes of stochastic attempts to apply a rule."""
applied: int = 0
failed: int = 0
@property
def attempts(self) -> int:
"""The total number of attempted rule applications."""
return self.applied + self.failed
[docs]
class System:
"""A Kappa system containing agents, rules, observables, and variables for simulation."""
_mixture: Mixture
_rules: dict[str, Rule]
_observables: dict[str, Expression]
_variables: dict[str, Expression]
_site_defaults: dict[str, Mapping[str, str]]
_tokens: dict[str, float]
_monitor: Optional["Monitor"]
_time: float
_next_update_time: Optional[float]
_reactivity_cache: Optional[tuple[float, ...]]
_tallies: dict[str, RuleTally]
_rng: random.Random # Random number generator for reproducibility of updates
[docs]
@classmethod
def read_ka(cls, filepath: str, seed: Optional[int] = None) -> Self:
"""Read and parse a Kappa .ka file to create a System.
Args:
filepath: Path to the Kappa file.
seed: Random seed for reproducibility.
"""
with open(filepath) as f:
return cls.from_ka(f.read(), seed=seed)
[docs]
@classmethod
def from_ka(cls, ka_str: str, seed: Optional[int] = None) -> Self:
"""Create a System from a Kappa (.ka style) string.
Args:
ka_str: Kappa language string containing a system definition.
seed: Random seed for reproducibility.
"""
from pykappa._parsing import (
kappa_parser,
KappaTransformer,
ExpressionTransformer,
)
input_tree = kappa_parser.parse(ka_str)
assert input_tree.data == "kappa_input"
variables: dict[str, Expression] = {}
observables: dict[str, Expression] = {}
token_inits: list[tuple[str, Expression]] = []
rules: list[Rule] = []
inits: list[tuple[Expression, Pattern]] = []
for child in input_tree.children:
tag = child.data
if tag in ["f_rule", "fr_rule", "ambi_rule", "ambi_fr_rule"]:
new_rules = KappaTransformer().transform(child)
rules.extend(new_rules)
elif tag == "variable_declaration":
name_tree = child.children[0]
assert name_tree.data == "declared_variable_name"
name = name_tree.children[0].value.strip("'\"")
expr_tree = child.children[1]
assert expr_tree.data == "algebraic_expression"
value = ExpressionTransformer.from_tree(expr_tree)
variables[name] = value
elif tag == "plot_declaration":
raise NotImplementedError
elif tag == "observable_declaration":
label_tree = child.children[0]
assert isinstance(label_tree, str)
name = label_tree.strip("'\"")
expr_tree = child.children[1]
assert expr_tree.data == "algebraic_expression"
value = ExpressionTransformer.from_tree(expr_tree)
observables[name] = value
elif tag == "signature_declaration":
pass # ignore agent signatures
elif tag == "declared_token":
pass # token declarations are handled via init_declaration
elif tag == "init_declaration":
amount = ExpressionTransformer.from_tree(child.children[0])
target = child.children[1]
if target.data == "declared_token_name":
token_inits.append((str(target.children[0]), amount))
else:
pattern = KappaTransformer().transform(target)
inits.append((amount, pattern))
elif tag == "definition":
pass # %def: directives not used
elif tag == "pattern":
raise NotImplementedError
else:
raise TypeError(f"Unsupported input type: {tag}")
system = cls(None, rules, observables, variables, seed=seed)
for init in inits:
system.add(init[1], int(init[0].evaluate(system)))
for token_name, amount_expr in token_inits:
system._tokens[token_name] = float(amount_expr.evaluate(system))
return system
[docs]
@classmethod
def from_kappa(
cls,
mixture: Optional[dict[str, int]] = None,
rules: Optional[Iterable[str]] = None,
observables: Optional[list[str] | dict[str, str]] = None,
variables: Optional[dict[str, str]] = None,
*args,
**kwargs,
) -> Self:
"""Create a System from Kappa strings.
Args:
mixture: Dictionary mapping agent patterns to initial counts.
rules: Iterable of rule strings in Kappa format.
observables: List of observable expressions or dict mapping names to expressions.
variables: Dictionary mapping variable names to expressions.
*args: Additional arguments passed to System constructor.
**kwargs: Additional keyword arguments passed to System constructor.
"""
real_rules = []
if rules is not None:
for rule in rules:
real_rules.extend(Rule.list_from_kappa(rule))
if observables is None:
real_observables = {}
elif isinstance(observables, list):
real_observables = {
f"o{i}": Expression.from_kappa(obs) for i, obs in enumerate(observables)
}
else:
real_observables = {
name: Expression.from_kappa(obs) for name, obs in observables.items()
}
real_variables = (
{}
if variables is None
else {name: Expression.from_kappa(var) for name, var in variables.items()}
)
system = cls(
None, real_rules, real_observables, real_variables, *args, **kwargs
)
if mixture is not None:
for pattern_str, count in mixture.items():
system.add(pattern_str, count)
return system
def __init__(
self,
mixture: Optional[Mixture] = None,
rules: Optional[Iterable[Rule]] = None,
observables: Optional[dict[str, Expression]] = None,
variables: Optional[dict[str, Expression]] = None,
tokens: Optional[dict[str, float]] = None,
site_defaults: Optional[dict[str, dict[str, str]]] = None,
monitor: bool = True,
seed: Optional[int] = None,
):
"""
Args:
mixture: Initial mixture state.
rules: Collection of rules to apply.
observables: Dictionary of observable expressions.
variables: Dictionary of variable expressions.
tokens: Dictionary of token names to initial values.
site_defaults: Maps agent types to site default states.
monitor: Whether to enable monitoring of simulation history.
seed: Random seed for reproducibility.
"""
self._rng = random.Random() if seed is None else random.Random(seed)
self._rules = (
{} if rules is None else {f"r{i}": rule for i, rule in enumerate(rules)}
)
if mixture is None:
mixture = Mixture(
track_components=any(
rule.component_constraint != "any" for rule in self._rules.values()
)
)
self._observables = {} if observables is None else dict(observables)
self._variables = {} if variables is None else dict(variables)
self._site_defaults = {
agent_type: MappingProxyType(dict(defaults))
for agent_type, defaults in (site_defaults or {}).items()
}
self._set_mixture(mixture)
self._time = 0
self._next_update_time = None
self._reactivity_cache = None
self._tokens = {} if tokens is None else dict(tokens)
self._tallies = {}
self._monitor = Monitor(self) if monitor else None
def __str__(self):
return self.kappa_str
def __getitem__(self, name: str) -> int | float:
"""Get the value of an observable or variable.
Raises:
KeyError: If name doesn't correspond to any observable or variable.
"""
if name in self._observables:
return self._observables[name].evaluate(self)
elif name in self._variables:
return self._variables[name].evaluate(self)
else:
raise KeyError(
f"Name {name} doesn't correspond to a declared observable or variable"
)
def __setitem__(self, name: str, value: float) -> None:
"""Update an existing variable to a new numeric value.
Args:
name: Name of a declared variable.
value: New numeric value.
Raises:
KeyError: If the name is not a declared variable.
ValueError: If the declared variable is not a numeric literal.
"""
if name not in self._variables:
raise KeyError(f"'{name}' is not a declared variable")
if self._variables[name]._type != "literal":
raise ValueError(
f"'{name}' is not a numeric literal and cannot be reassigned"
)
self._variables[name] = Expression("literal", value=value)
self._invalidate_next_event()
[docs]
def set_token(self, name: str, value: float) -> None:
"""Set a token's value."""
self._tokens[name] = value
self._invalidate_next_event()
@property
def mixture(self) -> Mixture:
"""The current state of agents and their connections."""
return self._mixture
@property
def time(self) -> float:
"""The current simulation time."""
return self._time
@property
def next_update_time(self) -> Optional[float]:
"""The time of the next update, or ``None`` if the system is nonreactive."""
if self._next_update_time is None:
if (reactivity := self.reactivity) == 0:
return None
self._next_update_time = self._time + self._rng.expovariate(reactivity)
return self._next_update_time
@property
def monitor(self) -> Optional["Monitor"]:
"""The monitor tracking simulation history, if enabled."""
return self._monitor
@property
def rules(self) -> Mapping[str, Rule]:
"""Maps rule names to rules."""
return MappingProxyType(self._rules)
@cached_property
def _signatures(self) -> dict[str, frozenset[str]]:
sites_by_type: dict[str, set[str]] = defaultdict(set)
for rule in self._rules.values():
for pattern in (rule.left, rule.right):
for agent in pattern.agents:
if agent is not None:
sites_by_type[agent.type].update(site.label for site in agent)
return {
agent_type: frozenset(sites) for agent_type, sites in sites_by_type.items()
}
@property
def signatures(self) -> Mapping[str, frozenset[str]]:
"""The complete site interface for each agent type as inferred from the rule set."""
return MappingProxyType(self._signatures)
@property
def site_defaults(self) -> Mapping[str, Mapping[str, str]]:
"""Maps agent types to their default site states."""
return MappingProxyType(self._site_defaults)
@property
def observables(self) -> Mapping[str, Expression]:
"""Maps observable names to expressions."""
return MappingProxyType(self._observables)
@property
def variables(self) -> Mapping[str, Expression]:
"""Maps variable names to expressions."""
return MappingProxyType(self._variables)
@property
def tokens(self) -> Mapping[str, float]:
"""Maps token names to their current values."""
return MappingProxyType(self._tokens)
@property
def tallies(self) -> Mapping[str, RuleTally]:
"""Maps rule names to counts of stochastic application outcomes."""
return MappingProxyType(self._tallies)
@property
def tally_totals(self) -> RuleTally:
"""Counts of all stochastic application attempts, combined across rules."""
return RuleTally(
applied=sum(tally.applied for tally in self._tallies.values()),
failed=sum(tally.failed for tally in self._tallies.values()),
)
@property
def tallies_table(self) -> str:
"""A formatted summary of stochastic rule application outcomes."""
totals = self.tally_totals
return str_table(
[
[str(rule), tally.applied, tally.failed, tally.attempts]
for rule, tally in self._tallies.items()
]
+ [["Total", totals.applied, totals.failed, totals.attempts]],
header=["Rule", "Applied", "Failed", "Attempts"],
)
@cached_property
def _reversible_rules(self) -> list[tuple[str, str]]:
"""Find forward/reverse rule pairs by checking pattern symmetry."""
names = list(self._rules.keys())
pairs = []
used = set()
for i, name_a in enumerate(names):
if name_a in used:
continue
rule_a = self._rules[name_a]
for name_b in names[i + 1 :]:
if name_b in used:
continue
rule_b = self._rules[name_b]
if rule_a.left.n_isomorphisms(
rule_b.right
) and rule_a.right.n_isomorphisms(rule_b.left):
pairs.append((name_a, name_b))
used.add(name_a)
used.add(name_b)
break
return pairs
@property
def kappa_str(self) -> str:
"""The system representation in Kappa (.ka style) format."""
kappa_list = []
# Append the inferred agent signature at the top
for agent, sites in self._signatures.items():
sig = ", ".join(sites)
kappa_list.append(f"%agent: {agent}({sig})")
# Format reversible rules with <-> notation
pairs = self._reversible_rules
paired = {name for pair in pairs for name in pair}
for fwd_name, rev_name in pairs:
fwd = self._rules[fwd_name]
rev = self._rules[rev_name]
kappa_list.append(
f"{fwd.left.kappa_str} <-> {fwd.right.kappa_str} "
f"@ {fwd._rate_str}, {rev._rate_str}"
)
# Otherwise format with -> notation
for name, rule in self._rules.items():
if name not in paired:
kappa_list.append(rule.kappa_str)
for var_name, var in self._variables.items():
kappa_list.append(f"%var: '{var_name}' {var.kappa_str}")
for obs_name, obs in self._observables.items():
obs_str = (
f"|{obs.kappa_str}|" if isinstance(obs, Component) else obs.kappa_str
)
kappa_list.append(f"%obs: '{obs_name}' {obs_str}")
kappa_list.append(self._mixture.kappa_str)
return "\n".join(kappa_list)
[docs]
def write_ka(self, filepath: str) -> None:
"""Write system information to a Kappa file."""
with open(filepath, "w") as f:
f.write(self.kappa_str)
[docs]
def save(self, filepath: str) -> None:
"""Save a checkpoint that can be continued with :meth:`System.load`."""
with open(filepath, "wb") as f:
cloudpickle.dump(self, f)
[docs]
@classmethod
def load(cls, filepath: str) -> Self:
"""Load a trusted checkpoint created by :meth:`System.save`.
Note:
Checkpoints must only be loaded from a trusted source and are intended
for use with the same PyKappa and Python versions.
"""
with open(filepath, "rb") as f:
system = cloudpickle.load(f)
if not isinstance(system, cls):
raise TypeError(
f"Checkpoint contains {type(system).__name__}, not {cls.__name__}"
)
return system
def _set_mixture(self, mixture: Mixture) -> None:
"""Set the system's mixture and update tracking."""
self._mixture = mixture
for rule in self._rules.values():
for component in rule.left.components:
if component not in mixture._embeddings:
mixture._track_component(component)
for expr in [*self._observables.values(), *self._variables.values()]:
for component_expr in expr._filter("component_pattern"):
mixture._track_component(component_expr._attrs["value"])
def _invalidate_next_event(self) -> None:
"""Discard cached reactivities after a state or rate change."""
self._next_update_time = None
self._reactivity_cache = None
def _enforce_signature(self, agent: "Agent") -> None:
"""Validate agent type and sites against the inferred signature and fill missing sites.
Raises:
ValueError: If the agent type is unknown or has sites not in the signature.
"""
if not self._signatures:
return
known = self._signatures.get(agent.type)
if known is None:
raise ValueError(
f"Agent type '{agent.type}' is not declared by any rule. "
f"Known agent types: {set(self._signatures)}"
)
unknown = {s.label for s in agent} - known
if unknown:
raise ValueError(
f"Agent '{agent.type}' has unknown site(s) {unknown}. "
f"Known sites for this type: {known}"
)
for label in known - agent.interface.keys():
agent._add_site(
Site(
label, self._site_defaults.get(agent.type, {}).get(label, "?"), "."
)
)
[docs]
def add(self, pattern: Pattern | Component | str, n_copies: int = 1) -> None:
"""Add instances of a pattern or component to the mixture using inferred agent signatures."""
if isinstance(pattern, str):
pattern = Pattern.from_kappa(pattern)
components = [pattern] if isinstance(pattern, Component) else pattern.components
for component in components:
for _ in range(n_copies):
self._mixture._add_component(
component, prepare_agent=self._enforce_signature
)
self._invalidate_next_event()
[docs]
def remove(self, component: Component) -> None:
"""Remove a specific component from the current mixture."""
self._mixture._remove_component(component)
self._invalidate_next_event()
@property
def reactivity(self) -> float:
"""The total reactivity of the system."""
return sum(self._rule_reactivities())
def _rule_reactivities(self) -> tuple[float, ...]:
"""Return rule reactivities cached for the pending simulation event."""
if self._reactivity_cache is None:
self._reactivity_cache = tuple(
rule.reactivity(self) for rule in self._rules.values()
)
return self._reactivity_cache
[docs]
def advance_time_to(self, time: float) -> None:
"""Advance time without applying an update.
Raises:
ValueError: If ``time`` is outside the interval before the next update.
"""
if time < self._time:
raise ValueError("cannot run backwards in simulation time")
if self.next_update_time is not None and time >= self.next_update_time:
raise ValueError("cannot advance to or past the next update")
self._time = time
[docs]
def update(self) -> None:
"""Perform one simulation step."""
if self._monitor is not None and not self._monitor.history["time"]:
self._monitor.update()
if (next_update_time := self.next_update_time) is None:
warnings.warn("system has no reactivity", RuntimeWarning)
if self._monitor is not None:
self._monitor.update()
return
self._time = next_update_time
self._next_update_time = None
rule = self._rng.choices(
list(self._rules.values()),
weights=self._rule_reactivities(),
)[0]
# Apply the rule
update = rule._select(self._mixture, rng=self._rng)
name = str(rule)
tally = self._tallies.get(name, RuleTally())
if update is not None:
self._tallies[name] = RuleTally(
applied=tally.applied + 1, failed=tally.failed
)
for agent in update.agents_to_add:
self._enforce_signature(agent)
previous_components, current_components = self._mixture._apply_update(
update
)
for expr, name in rule.token_updates:
self._tokens[name] += expr.evaluate(self)
self._reactivity_cache = tuple(
(
candidate._reactivity_from_embeddings(
candidate.update_component_weights(
self._mixture, previous_components, current_components
),
self,
)
if candidate.component_constraint != "any"
else candidate.reactivity(self)
)
for candidate in self._rules.values()
)
else:
self._tallies[name] = RuleTally(
applied=tally.applied, failed=tally.failed + 1
)
if self._monitor is not None:
self._monitor.update()
if update is None:
self._reactivity_cache = None
[docs]
def apply(self, transformation: str, n: int = 1) -> None:
"""Apply a transformation immediately for a specified number of times.
Unlike `update`, this does not advance simulation time or use stochastic
selection — the rule fires exactly ``n`` times using randomly chosen
embeddings.
Args:
transformation: Kappa string representation of the rule.
n: Number of times to apply the rule.
"""
rule = Rule.from_kappa(transformation + " @ 0")
self._invalidate_next_event()
for _ in range(n):
update = rule._select(self._mixture, rng=self._rng)
if update is not None:
for agent in update.agents_to_add:
self._enforce_signature(agent)
self._mixture._apply_update(update)
[docs]
def update_via_kasim(self, time: float) -> None:
"""Simulate for a given amount of time using KaSim.
Note:
KaSim must be installed and in the PATH.
Some features are not compatible between PyKappa and KaSim.
"""
self._invalidate_next_event()
assert shutil.which("KaSim"), "KaSim not found in the PATH."
if any(rule.n_symmetries > 1 for rule in self._rules.values()):
warnings.warn(
"Some rules have multiple symmetries. "
"PyKappa normalizes reactivities accordingly: results may differ from KaSim."
)
history = None # the observable history
with tempfile.TemporaryDirectory() as tmpdirname:
snap_path = os.path.join(tmpdirname, "snap.ka")
out_path = os.path.join(tmpdirname, "out.ka")
in_path = os.path.join(tmpdirname, "in.ka")
output_lines = [
self.kappa_str,
f'%mod: alarm {time} do $SNAPSHOT "{snap_path}";',
]
if self._observables:
output_lines.append("%mod: [true] do $PLOTENTRY; repeat [true]")
# Run KaSim
with open(in_path, "w") as f:
f.write("\n".join(output_lines))
subprocess.run(
["KaSim", in_path, "-l", str(time), "-d", tmpdirname, "-o", out_path],
stdout=subprocess.DEVNULL,
stderr=subprocess.DEVNULL,
)
# Read KaSim output
with open(snap_path) as f:
content = f.read()
if self._observables:
with open(out_path) as f:
reader = csv.reader(f)
header = next(row for row in reader if row and row[0] == "[T]")
columns = ["time", *header[1:]]
history = {name: [] for name in columns}
for row in reader:
history["time"].append(self._time + float(row[0]))
for name, value in zip(columns[1:], row[1:]):
history[name].append(float(value))
content = content.replace(
",\n", ", "
) # KaSim splits long components across lines
output_kappa_str = "".join(
line.split("/")[0] + line.split("/")[-1]
for line in content.splitlines(keepends=True)
if line.startswith("%init")
)
# Apply the update
self._set_mixture(System.from_ka(output_kappa_str).mixture)
self._time += time
# Update the monitor
if self._monitor is not None and history is not None:
for name, values in history.items():
self._monitor.history[name].extend(values)
[docs]
def kd_table(self, volume: float = 1.0) -> str:
"""Summarize kinetic constants of two-component binding/unbinding rules
given volume in liters."""
from pykappa.analysis import _kd_table
return _kd_table(self, volume=volume)
[docs]
def rule_graph(self) -> "Source":
"""Visualize a ruleset as a site graph of local transformations.
Solid edges = bond formation; dashed edges = bond breaking. Sites that
change state show their transition as ``site {old→new}``. Creation and
degradation are shown as directed edges to/from a sink node.
Note:
This is a lossy projection that neglects conditions of transformations;
multiple rulesets can yield the same graph.
"""
from pykappa.analysis import _rule_graph
return _rule_graph(self)
