diff --git a/README.md b/README.md index f225cf1..a7c6d14 100644 --- a/README.md +++ b/README.md @@ -19,6 +19,9 @@ Deployed at https://meadows.apoena.dev substrate; the domain Model is the source of truth. - [ADR-0003](./docs/adr/0003-flow-as-node-materialised-clouds.md) — a Flow is a node; Source/Sink clouds are materialised nodes. + - [ADR-0004](./docs/adr/0004-rate-rules-not-formulas.md) — simulation behaviour + comes from a small fixed vocabulary of rules over Information Links, not + free-form formulas. ## Stack diff --git a/docs/adr/0004-rate-rules-not-formulas.md b/docs/adr/0004-rate-rules-not-formulas.md new file mode 100644 index 0000000..7a9b34c --- /dev/null +++ b/docs/adr/0004-rate-rules-not-formulas.md @@ -0,0 +1,74 @@ +# Behaviour comes from rules over Information Links, not free-form formulas + +_Part of [meadows](../../README.md) · see [DESIGN.md](../../DESIGN.md)._ + +Numeric simulation (phase 2) makes a **Model** _alive_: Stocks accumulate over +time, Flows and Converters recompute each instant. Two coupled decisions shape +how a Model carries the numbers, and both trade expressive power for **valid by +construction** — the right trade for a tool that exists to _popularise_ systems +thinking, not to compete with Vensim. + +## Decision + +**1. The Information Link _is_ the declared dependency.** A Flow's or Converter's +inputs are exactly the elements that link into it. There is no separate "equation +references" namespace to keep in sync — the wiring you draw _is_ the wiring the +simulator reads. The same signed graph the loop detector walks (ADR-0001) is the +graph the simulator integrates, so the loops you _see_ classified R/B are the +loops you _run_. They can never disagree. + +**2. A Flow/Converter computes from a small fixed vocabulary of `Rule`s, not a +typed-in formula.** Each instantaneous element picks one rule and a plain number +or two — never an expression: + +| Rule | Value | Reads (via Information Links) | Emergent behaviour | +| ---------------- | --------------------------- | ----------------------------------------- | -------------------------- | +| **Constant** | a fixed number | nothing | linear Stock change | +| **Proportional** | `factor × (its `+` inputs)` | the `+`-polarity inputs | exponential growth / decay | +| **Gap** | `factor × (level − target)` | the `+` input is _level_, `−` is _target_ | goal-seeking / asymptotic | + +The famous curves are _compositions_ of these over the structure — a logistic +S-curve is Proportional growth meeting a Gap-driven ceiling (limits-to-growth); +goal-seeking decay is a lone Gap (coffee cooling). The user sets up a local rule; +the global shape **emerges**. That emergence _is_ the lesson. + +**Polarity does double duty.** The `+`/`−` already captured for loop +classification (ADR-0001) also selects each operand's role: Proportional reads +its `+` inputs; Gap reads its `+` input as the level and its `−` input as the +target. One gesture, two payoffs — no new per-link data. + +## Considered Options + +- **Free-form expression strings** (`birth rate = Population × fertility`) — + maximally expressive, and what `equation?: string` originally anticipated. + Rejected: needs a parser + a sandbox (never `eval`/`new Function`), invites + broken-formula and name-resolution errors (auto-names contain spaces), and lets + a learner _paint_ a curve instead of discovering it from structure. +- **Pick the output curve** (label a Stock "exponential" / "logarithmic") — + rejected: it is the answer, not the cause, and it breaks the moment feedback + decides the shape. "Logarithmic" in particular has no honest local rule; what + people mean by it is asymptotic approach — which _is_ the Gap rule. +- **Rules over Information Links (chosen)** — no parser, valid by construction, + and it teaches structure → behaviour. + +## Consequences + +- The domain types gain a `Rule` union on Flow/Converter (replacing the unused + `equation?: string`), an optional `initialValue` on Stock, and an optional + `SimSpec` (`start` / `stop` / `dt`) on the Model. All optional and additive, so + existing saved Models still load (F8); they are simply not _simulatable_ until + equipped. +- **Algebraic loops are an error.** A cycle in the wiring is legitimate feedback + **iff it passes through a Stock** — the Stock supplies last-step state and so + breaks the within-step dependency. A cycle among only Flows/Converters has no + Stock to break it: the simulator cannot order it and rejects it. The simulator + reuses the cycle machinery to detect this; it is a new _sim-readiness_ check, + distinct from structural validity (validation.ts). +- A reader seeing `{ kind: "gap", factor }` on a Flow and wondering where its + operands come from should look here: they are the Flow's inbound Information + Links, picked by Polarity. +- The vocabulary starts deliberately small (Constant / Proportional / Gap — + enough for linear, exponential, and goal-seeking, and for the coffee and + savings samples). Growing it is additive: a new `kind` in the union plus a case + in the evaluator. Multi-input products (e.g. `Population × fertility`) are a + later increment, not a phase-2 blocker. diff --git a/src/model/io.ts b/src/model/io.ts index 7f74793..9210e79 100644 --- a/src/model/io.ts +++ b/src/model/io.ts @@ -21,6 +21,8 @@ import { type NodeKind, type Polarity, type Position, + type Rule, + type SimSpec, } from "./types" /** Pretty JSON so an exported Model is human-readable and diff-friendly (F8). */ @@ -49,6 +51,24 @@ function isPolarity(value: unknown): value is Polarity { return value === "+" || value === "-" } +/** A Rule must be one of the fixed kinds with its numeric parameter (ADR-0004). */ +function isRule(value: unknown): value is Rule { + if (!isObject(value)) return false + if (value.kind === "constant") return isFiniteNumber(value.value) + if (value.kind === "proportional" || value.kind === "gap") return isFiniteNumber(value.factor) + return false +} + +function isSimSpec(value: unknown): value is SimSpec { + return ( + isObject(value) && + isFiniteNumber(value.start) && + isFiniteNumber(value.stop) && + isFiniteNumber(value.dt) && + value.dt > 0 + ) +} + /** Validate one node by kind, returning an error string or null when valid. */ function nodeError(value: unknown, index: number): string | null { if (!isObject(value)) return `nodes[${index}] is not an object` @@ -65,6 +85,17 @@ function nodeError(value: unknown, index: number): string | null { if (typeof value.source !== "string") return `${at}.source must be a string` if (typeof value.target !== "string") return `${at}.target must be a string` } + // Simulation fields (ADR-0004) are optional, but if present must be well-formed. + if (kind === "stock" && value.initialValue !== undefined && !isFiniteNumber(value.initialValue)) { + return `${at}.initialValue must be a finite number` + } + if ( + (kind === "flow" || kind === "converter") && + value.rule !== undefined && + !isRule(value.rule) + ) { + return `${at}.rule must be a valid rule (constant/proportional/gap)` + } return null } @@ -104,6 +135,9 @@ export function parseModel(text: string): ParseResult { if (typeof data.name !== "string") return { ok: false, error: "name must be a string" } if (!Array.isArray(data.nodes)) return { ok: false, error: "nodes must be an array" } if (!Array.isArray(data.infoLinks)) return { ok: false, error: "infoLinks must be an array" } + if (data.sim !== undefined && !isSimSpec(data.sim)) { + return { ok: false, error: "sim must be {start, stop, dt} numbers with dt > 0" } + } for (let i = 0; i < data.nodes.length; i++) { const error = nodeError(data.nodes[i], i) @@ -134,6 +168,13 @@ export function parseModel(text: string): ParseResult { return { ok: true, - model: { version: MODEL_VERSION, id: data.id, name: data.name, nodes, infoLinks: links }, + model: { + version: MODEL_VERSION, + id: data.id, + name: data.name, + nodes, + infoLinks: links, + ...(data.sim !== undefined && { sim: data.sim as SimSpec }), + }, } } diff --git a/src/model/simulation.ts b/src/model/simulation.ts new file mode 100644 index 0000000..e54948b --- /dev/null +++ b/src/model/simulation.ts @@ -0,0 +1,228 @@ +/** + * Simulation engine (ADR-0004) — what makes a Model *alive*. Pure data in, time + * series out; it knows nothing of Vue or the canvas, so it is trivially testable + * and runs off any reactive frame. + * + * The system-dynamics loop, one step of `dt`: + * 1. Evaluate the instantaneous network — every Converter and Flow — in + * dependency order, reading the *current* Stock values. A Converter/Flow + * depends on the elements that link into it (ADR-0004: the Information Link + * *is* the dependency); Stocks are state, available without ordering. + * 2. Integrate every Stock *simultaneously* (forward Euler): + * `stock += dt × (Σ inflow rates − Σ outflow rates)`. All net rates are read + * from the same pre-update state, so updating one Stock never feeds another + * within the same step. + * + * A cycle in the wiring is legitimate feedback only if it passes through a Stock + * (the Stock supplies last-step state and breaks the within-step dependency). A + * cycle among only Converters/Flows is an **algebraic loop** — unorderable, and + * rejected by `evaluationOrder` / surfaced by `checkSimReady`. + */ +import { + type ConverterNode, + DEFAULT_SIM_SPEC, + type FlowNode, + type InformationLink, + type Model, + type Rule, + type SimSpec, +} from "./types" + +/** Thrown when a Model cannot be simulated as wired (e.g. an algebraic loop). */ +export class SimulationError extends Error { + constructor(message: string) { + super(message) + this.name = "SimulationError" + } +} + +/** A completed run: aligned `times` and per-element value tracks. */ +export interface Run { + /** The time at each recorded sample, `start … stop` in steps of `dt`. */ + times: number[] + /** nodeId → its value at each time index (Stocks, Flows, Converters; not Clouds). */ + series: Map + /** True if the run ran past representable numbers and stopped early (see below). */ + diverged: boolean +} + +/** A Converter or Flow — the stateless, instantaneous elements an order applies to. */ +type Instant = ConverterNode | FlowNode + +/** + * Order the instantaneous network so each element is computed after everything + * it reads. Throws `SimulationError` on an algebraic loop (a Converter/Flow cycle + * with no Stock to break it). Stocks are excluded — they are state, not computed. + */ +export function evaluationOrder(model: Model): Instant[] { + const instant = model.nodes.filter( + (node): node is Instant => node.kind === "flow" || node.kind === "converter", + ) + const ids = new Set(instant.map((node) => node.id)) + + // deps[x] = the instantaneous elements x reads (links from another Flow/Converter). + const deps = new Map>(instant.map((node) => [node.id, new Set()])) + for (const link of model.infoLinks) { + if (ids.has(link.source) && ids.has(link.target)) deps.get(link.target)?.add(link.source) + } + + const order: Instant[] = [] + const resolved = new Set() + while (order.length < instant.length) { + const next = instant.find( + (node) => + !resolved.has(node.id) && [...(deps.get(node.id) ?? [])].every((d) => resolved.has(d)), + ) + if (!next) { + const names = instant + .filter((node) => !resolved.has(node.id)) + .map((node) => node.name) + .join(", ") + throw new SimulationError( + `Algebraic loop: ${names} depend on each other with no Stock to break the cycle.`, + ) + } + resolved.add(next.id) + order.push(next) + } + return order +} + +/** Evaluate one Rule given a way to read the values feeding in via `links`. */ +function evalRule(rule: Rule, links: InformationLink[], valueOf: (id: string) => number): number { + switch (rule.kind) { + case "constant": + return rule.value + case "proportional": { + // factor × the product of every `+`-polarity input (one input → factor × it). + let value = rule.factor + for (const link of links) if (link.polarity === "+") value *= valueOf(link.source) + return value + } + case "gap": { + // factor × (level − target): the `+` input is the level, the `−` the target. + const level = links.find((link) => link.polarity === "+") + const target = links.find((link) => link.polarity === "-") + return ( + rule.factor * ((level ? valueOf(level.source) : 0) - (target ? valueOf(target.source) : 0)) + ) + } + } +} + +/** + * Run the Model and return aligned time series. Assumes a sim-ready Model + * (see `checkSimReady`); a missing rule evaluates to 0 and a missing initial + * value to 0 rather than throwing, so a half-built Model still produces a plot. + * Throws `SimulationError` only on an algebraic loop, which has no defined order. + */ +export function simulate(model: Model, spec: SimSpec = model.sim ?? DEFAULT_SIM_SPEC): Run { + const order = evaluationOrder(model) + const nodeById = new Map(model.nodes.map((node) => [node.id, node])) + + const inbound = new Map() + for (const link of model.infoLinks) { + const list = inbound.get(link.target) + if (list) list.push(link) + else inbound.set(link.target, [link]) + } + + const stocks = model.nodes.filter((node) => node.kind === "stock") + const flows = model.nodes.filter((node): node is FlowNode => node.kind === "flow") + + const stockValues = new Map(stocks.map((s) => [s.id, s.initialValue ?? 0])) + + const times: number[] = [] + const series = new Map() + for (const node of model.nodes) if (node.kind !== "cloud") series.set(node.id, []) + let diverged = false + + // dt ≤ 0 would never advance (and 0 diverges); a non-positive step means no run. + const steps = spec.dt > 0 ? Math.max(0, Math.floor((spec.stop - spec.start) / spec.dt)) : 0 + for (let i = 0; i <= steps; i++) { + // A run-away Reinforcing loop over a long horizon can exceed what a float + // holds. Stop at the last valid sample and flag it rather than plotting NaN. + if (!stocks.every((s) => Number.isFinite(stockValues.get(s.id) ?? 0))) { + diverged = true + break + } + + // 1. Evaluate the instantaneous network from the current Stock values. + const computed = new Map() + const valueOf = (id: string): number => { + const node = nodeById.get(id) + if (!node) return 0 + if (node.kind === "stock") return stockValues.get(id) ?? 0 + if (node.kind === "cloud") return 0 + return computed.get(id) ?? 0 + } + for (const node of order) { + computed.set( + node.id, + node.rule ? evalRule(node.rule, inbound.get(node.id) ?? [], valueOf) : 0, + ) + } + + // Record this sample (Stocks at their current value, Flows/Converters as just computed). + times.push(spec.start + i * spec.dt) + for (const s of stocks) series.get(s.id)?.push(stockValues.get(s.id) ?? 0) + for (const node of order) series.get(node.id)?.push(computed.get(node.id) ?? 0) + + if (i >= steps) break + + // 2. Integrate every Stock simultaneously (forward Euler) — but with + // non-negative stocks: an outflow can't drain more than its source holds this + // step. Scale a stock's competing outflows together if they would overdraw, + // and apply the scaled rate to both ends so quantity is conserved. You can't + // infect more people than are susceptible — and that floor is exactly what + // stops a bilinear model (S × I) from flipping sign and diverging. + const rate = new Map(flows.map((f) => [f.id, computed.get(f.id) ?? 0])) + for (const s of stocks) { + const available = stockValues.get(s.id) ?? 0 + const drains = flows.filter((f) => f.source === s.id && (rate.get(f.id) ?? 0) > 0) + const totalOut = drains.reduce((sum, f) => sum + (rate.get(f.id) ?? 0), 0) + if (totalOut * spec.dt > available) { + const scale = available / (totalOut * spec.dt) + for (const f of drains) rate.set(f.id, (rate.get(f.id) ?? 0) * scale) + } + } + + const next = new Map(stockValues) + for (const s of stocks) { + let net = 0 + for (const flow of flows) { + const r = rate.get(flow.id) ?? 0 + if (flow.target === s.id) net += r // an inflow fills it + if (flow.source === s.id) net -= r // an outflow drains it + } + next.set(s.id, (stockValues.get(s.id) ?? 0) + spec.dt * net) + } + for (const [id, value] of next) stockValues.set(id, value) + } + + return { times, series, diverged } +} + +/** + * What stands between this Model and a run, as human-readable lines (empty array + * = ready). Distinct from structural validity (validation.ts): a Model can be a + * perfectly valid diagram yet not carry the numbers a simulation needs. + */ +export function checkSimReady(model: Model): string[] { + const problems: string[] = [] + for (const node of model.nodes) { + if (node.kind === "stock" && node.initialValue === undefined) { + problems.push(`${node.name} has no initial value.`) + } + if ((node.kind === "flow" || node.kind === "converter") && !node.rule) { + problems.push(`${node.name} has no rule yet.`) + } + } + try { + evaluationOrder(model) + } catch (error) { + if (error instanceof SimulationError) problems.push(error.message) + else throw error + } + return problems +} diff --git a/src/model/types.ts b/src/model/types.ts index 89db7ff..51fa0ad 100644 --- a/src/model/types.ts +++ b/src/model/types.ts @@ -31,6 +31,37 @@ export interface Position { y: number } +/** + * How a Flow's rate or a Converter's value is computed — a small fixed + * vocabulary, *not* a free-form formula, so a Model stays valid by construction + * and teaches structure → behaviour (ADR-0004). Operands are read from the + * element's inbound Information Links; Polarity selects each operand's role. + * + * - `constant` — a fixed number; reads nothing. (→ linear Stock change) + * - `proportional` — `factor ×` the product of its `+`-polarity inputs. + * (→ exponential growth/decay) + * - `gap` — `factor × (level − target)`, where the `+` input is the + * level and the `−` input the target. (→ goal-seeking) + */ +export type Rule = + | { kind: "constant"; value: number } + | { kind: "proportional"; factor: number } + | { kind: "gap"; factor: number } + +/** + * The run parameters for a simulation: integrate from `start` to `stop` in steps + * of `dt` (forward Euler). Optional on a Model; absent means the diagram phase or + * a not-yet-simulated Model — the engine falls back to `DEFAULT_SIM_SPEC`. + */ +export interface SimSpec { + start: number + stop: number + dt: number +} + +/** Sensible default run when a Model carries no `sim` of its own. */ +export const DEFAULT_SIM_SPEC: SimSpec = { start: 0, stop: 100, dt: 1 } + interface BaseNode { id: string position: Position @@ -60,8 +91,8 @@ export interface FlowNode extends BaseNode { source: string /** Node id of the Stock or Cloud the Flow feeds into. */ target: string - /** Rate expression, recomputed each instant. Optional in the diagram phase. */ - equation?: string + /** How its rate is computed each instant (ADR-0004). Optional in the diagram phase. */ + rule?: Rule } /** @@ -71,8 +102,8 @@ export interface FlowNode extends BaseNode { export interface ConverterNode extends BaseNode { kind: "converter" name: string - /** Expression or constant. Optional in the diagram phase. */ - equation?: string + /** How its value is computed each instant (ADR-0004). Optional in the diagram phase. */ + rule?: Rule } /** @@ -107,4 +138,6 @@ export interface Model { name: string nodes: ModelNode[] infoLinks: InformationLink[] + /** Run parameters for simulation (ADR-0004). Optional; absent → not yet simulated. */ + sim?: SimSpec }