diff --git a/experiments/dump_leaf_quality.pl b/experiments/dump_leaf_quality.pl
new file mode 100644
index 0000000..500214e
--- /dev/null
+++ b/experiments/dump_leaf_quality.pl
@@ -0,0 +1,99 @@
+#!/usr/bin/perl
+# Parity oracle for homemaker-py-gnw: per-leaf quality factors and per-storey
+# cost/value from Urb's programme-driven fitness.
+#
+# Stock urb-fitness.pl cannot emit the per-leaf DEBUG lines: Leaf.pm/Storey.pm
+# copy $Urb::Dom::Fitness::DEBUG into a package var at *compile* time, before
+# Fitness.pm's own `our $DEBUG = $ENV{DEBUG}` line has run, so the flag is
+# permanently false.  This wrapper sets those package vars after loading and
+# also wraps process_storey to print each storey's cost/value subtotal.
+#
+# Usage (from the corpus directory, URB_NO_OCCLUSION required):
+#   cd examples/programme-house
+#   URB_NO_OCCLUSION=1 perl -I$URB/lib dump_leaf_quality.pl a.dom b.dom ...
+#
+# Output (stdout):
+#   === FILE <name>
+#   ...Leaf.pm debug lines ('  quality perpendicular: <f>' x7, then
+#      'leaf level: L id: I type: T rate: R area: A width: W proportion: P')
+#   STOREY <level_id> cost: <c> value: <v>
+#   INITIAL_COST <c>
+#   FAIL <message>            (one per failure)
+#   SCORE <fitness>
+
+use strict;
+use warnings;
+use 5.010;
+use Urb::Dom;
+use Urb::Dom::Fitness;
+use Urb::Dom::Fitness::ProgrammeDriven;
+use YAML;
+use Carp;
+
+die "URB_NO_OCCLUSION=1 required: native fitness ports simple crinkliness only\n"
+    unless $ENV{URB_NO_OCCLUSION};
+
+# Re-enable the per-leaf/storey debug gates (see header) and route everything
+# to STDOUT so output interleaves deterministically.
+$Urb::Dom::Fitness::Leaf::DEBUG = 1;
+$Urb::Dom::Fitness::Storey::DEBUG = 1;
+$Urb::Dom::Fitness::DEBUG = 1;
+{
+    no warnings 'redefine';
+    *Urb::Dom::Fitness::debug = sub { print join(' ', @_), "\n"; return 1 };
+    *Urb::Dom::Fitness::Base::debug = sub {
+        my @text = @_;
+        shift @text if ref $text[0];    # method form: drop $self
+        print join(' ', @text), "\n";
+        return 1;
+    };
+}
+
+my $orig_process_storey = \&Urb::Dom::Fitness::Storey::process_storey;
+{
+    no warnings 'redefine';
+    *Urb::Dom::Fitness::Storey::process_storey = sub {
+        my @args = @_;
+        my $result = $orig_process_storey->(@args);
+        printf "STOREY %s cost: %.17g value: %.17g\n",
+            $args[5], $result->{cost}, $result->{value};
+        return $result;
+    };
+}
+
+# Config loading exactly as urb-fitness.pl (project-level then local).
+my $config = undef;
+for my $path ('../patterns.config', 'patterns.config')
+{
+    next unless -e $path;
+    my $config_temp = YAML::LoadFile ($path);
+    $config->{$_} = $config_temp->{$_} for keys %{$config_temp};
+}
+my $costs = undef;
+for my $path ('../costs.config', 'costs.config')
+{
+    next unless -e $path;
+    my $costs_temp = YAML::LoadFile ($path);
+    $costs->{$_} = $costs_temp->{$_} for keys %{$costs_temp};
+}
+
+die "no spaces config: wrapper only supports programme-driven mode\n"
+    unless $config && exists $config->{spaces};
+my $assessor = Urb::Dom::Fitness::ProgrammeDriven->new ($config, $costs);
+
+for my $path_yaml (@ARGV)
+{
+    my @items = YAML::LoadFile ($path_yaml);
+    next unless scalar @items;
+    my $dom = Urb::Dom->new;
+    $dom->Deserialise (@items);
+
+    say "=== FILE $path_yaml";
+    printf "INITIAL_COST %.17g\n", $assessor->Cost ('plot') * $dom->Area;
+    my $score = $assessor->_apply ($dom);
+    say "FAIL $_" for ($dom->Failures);
+    printf "SCORE %.17g\n", $score;
+    $dom->DESTROY;
+}
+
+0;
diff --git a/experiments/leaf_parity.py b/experiments/leaf_parity.py
new file mode 100644
index 0000000..1cc485e
--- /dev/null
+++ b/experiments/leaf_parity.py
@@ -0,0 +1,183 @@
+#!/usr/bin/env python3
+"""homemaker-py-gnw acceptance: per-leaf quality factors and per-storey
+cost/value parity between homemaker.fitness and Urb's programme-driven fitness.
+
+Runs experiments/dump_leaf_quality.pl (URB_NO_OCCLUSION oracle with per-leaf
+DEBUG re-enabled) over the corpus, parses its output, and diffs against the
+native evaluation: 7 quality factors + rate + area per leaf, cost/value per
+storey, initial (plot) cost, and the leaf-scope failure set.
+
+Usage: python3 experiments/leaf_parity.py [corpus_dir]
+"""
+
+from __future__ import annotations
+
+import math
+import re
+import subprocess
+import sys
+from collections import defaultdict
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
+from homemaker import dom, fitness, graph  # noqa: E402
+
+URB = Path("/home/bruno/src/urb")
+CORPUS = Path(sys.argv[1]) if len(sys.argv) > 1 else URB / "examples" / "programme-house"
+WRAPPER = Path(__file__).resolve().parent / "dump_leaf_quality.pl"
+
+REL_TOL = 1e-9
+
+# Failures emitted by the gnw scope (leaf quality + cost functions + the two
+# per-leaf structural checks in process_storey). All carry a "level/" prefix —
+# requiring it excludes building-level fails like 'no outside public access'
+# (homemaker-py-hgg scope), which would otherwise match the ' access' suffix.
+_LEAF_FAIL_RE = re.compile(
+    r"^\d+/.* (perpendicular|proportion|size|width|crinkliness|daylight|access"
+    r"|edge too long|unsupported covered outside|covered outside above ground)$"
+)
+
+_FACTOR_RE = re.compile(r"^  quality (\w+): (\S+)$")
+_LEAF_RE = re.compile(
+    r"^leaf level: (\d+) id: ([lr]*) type: (\S+) rate: (\S+) area: (\S+) "
+    r"width: (\S+) proportion: (\S+)$"
+)
+_STOREY_RE = re.compile(r"^STOREY (\d+) cost: (\S+) value: (\S+)$")
+
+
+def run_oracle(files: list[str]) -> dict[str, dict]:
+    """Run the wrapper over all files; return per-file parsed records."""
+    proc = subprocess.run(
+        ["perl", f"-I{URB}/lib", str(WRAPPER)] + files,
+        cwd=CORPUS,
+        env={"URB_NO_OCCLUSION": "1", "PATH": "/usr/bin:/bin"},
+        capture_output=True,
+        text=True,
+    )
+    if proc.returncode != 0:
+        sys.exit(f"wrapper failed:\n{proc.stderr}")
+
+    out: dict[str, dict] = {}
+    rec = None
+    pending: dict[str, float] = {}
+    for line in proc.stdout.splitlines():
+        if line.startswith("=== FILE "):
+            rec = out[line[len("=== FILE "):]] = {
+                "leaves": {}, "storeys": {}, "fails": [], "initial_cost": None,
+                "score": None,
+            }
+            pending = {}
+            continue
+        if rec is None:
+            continue
+        m = _FACTOR_RE.match(line)
+        if m:
+            pending[m.group(1)] = float(m.group(2))
+            continue
+        m = _LEAF_RE.match(line)
+        if m:
+            level, lid = int(m.group(1)), m.group(2)
+            rec["leaves"][(level, lid)] = {
+                "type": m.group(3),
+                "rate": float(m.group(4)),
+                "area": float(m.group(5)),
+                "factors": pending,
+            }
+            pending = {}
+            continue
+        m = _STOREY_RE.match(line)
+        if m:
+            rec["storeys"][int(m.group(1))] = (float(m.group(2)), float(m.group(3)))
+            continue
+        if line.startswith("INITIAL_COST "):
+            rec["initial_cost"] = float(line.split()[1])
+        elif line.startswith("FAIL "):
+            rec["fails"].append(line[len("FAIL "):])
+        elif line.startswith("SCORE "):
+            rec["score"] = float(line.split()[1])
+    return out
+
+
+def native_eval(path: Path, fit: fitness.Fitness) -> dict:
+    root = dom.load(str(path))
+    fails: list[str] = []
+    fit.preprocess_building(root)
+    dom.merge_divided(root)
+    graphs = graph.build_graphs(root)
+
+    rec: dict = {"leaves": {}, "storeys": {}, "fails": fails,
+                 "initial_cost": fit.plot_cost(root)}
+    for li, lvl in enumerate(dom.levels(root)):
+        se = fit.process_storey(lvl, graphs[li], li, fails.append)
+        rec["storeys"][li] = (se.cost, se.value)
+        for le in se.leaves:
+            rec["leaves"][(le.level, le.id)] = {
+                "type": le.type, "rate": le.rate, "area": le.area,
+                "factors": le.factors,
+            }
+    return rec
+
+
+def close(a: float, b: float) -> bool:
+    return math.isclose(a, b, rel_tol=REL_TOL, abs_tol=1e-12)
+
+
+def main() -> int:
+    files = sorted(p.name for p in CORPUS.glob("*.dom"))
+    conf, cost = fitness.load_config(CORPUS)
+    fit = fitness.Fitness(conf, cost)
+    oracle = run_oracle(files)
+
+    n_leaves = n_factors = 0
+    bad: dict[str, list[str]] = defaultdict(list)
+
+    for name in files:
+        o = oracle[name]
+        n = native_eval(CORPUS / name, fit)
+
+        if not close(o["initial_cost"], n["initial_cost"]):
+            bad[name].append(
+                f"initial cost {o['initial_cost']} != {n['initial_cost']}")
+
+        if set(o["leaves"]) != set(n["leaves"]):
+            bad[name].append(
+                f"leaf sets differ: oracle-only {sorted(set(o['leaves']) - set(n['leaves']))}, "
+                f"native-only {sorted(set(n['leaves']) - set(o['leaves']))}")
+        for key in sorted(set(o["leaves"]) & set(n["leaves"])):
+            ol, nl = o["leaves"][key], n["leaves"][key]
+            n_leaves += 1
+            for fld in ("rate", "area"):
+                if not close(ol[fld], nl[fld]):
+                    bad[name].append(f"leaf {key} {fld}: {ol[fld]} != {nl[fld]}")
+            if ol["type"] != nl["type"]:
+                bad[name].append(f"leaf {key} type: {ol['type']} != {nl['type']}")
+            for fac, ov in ol["factors"].items():
+                n_factors += 1
+                nv = nl["factors"].get(fac)
+                if nv is None or not close(ov, nv):
+                    bad[name].append(f"leaf {key} {fac}: {ov} != {nv}")
+
+        for li, (oc, ov) in o["storeys"].items():
+            nc, nv = n["storeys"].get(li, (None, None))
+            if nc is None or not close(oc, nc):
+                bad[name].append(f"storey {li} cost: {oc} != {nc}")
+            if nv is None or not close(ov, nv):
+                bad[name].append(f"storey {li} value: {ov} != {nv}")
+
+        o_fails = sorted(f for f in o["fails"] if _LEAF_FAIL_RE.search(f))
+        n_fails = sorted(n["fails"])
+        if o_fails != n_fails:
+            bad[name].append(
+                f"leaf fails differ:\n    oracle: {o_fails}\n    native: {n_fails}")
+
+    for name in sorted(bad):
+        print(f"MISMATCH {name}")
+        for msg in bad[name]:
+            print(f"  {msg}")
+    print(f"\n{len(files)} files, {n_leaves} leaves, {n_factors} factors compared; "
+          f"{len(bad)} files with mismatches")
+    return 1 if bad else 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())
diff --git a/src/homemaker/dom.py b/src/homemaker/dom.py
index 7a9f968..516fd98 100644
--- a/src/homemaker/dom.py
+++ b/src/homemaker/dom.py
@@ -235,6 +235,48 @@ def _below_leaves(n: Node) -> "list[Node]":
     return bm.leaves() if bm is not None else []
 
 
+def _above_node(n: Node) -> "Node | None":
+    """Matching node on the level above; mirrors ``Urb::Quad::Above``."""
+    lr = _level_root(n)
+    if lr.above is None:
+        return None
+    return lr.above.by_id(n.id)
+
+
+def _above_more(n: Node) -> "Node | None":
+    """Matching node on level above, walking up to parent if path absent;
+    mirrors ``Urb::Quad::Above_More``."""
+    if _level_root(n).above is None:
+        return None
+    a = _above_node(n)
+    if a is not None:
+        return a
+    if n.parent is None:
+        return None
+    return _above_more(n.parent)
+
+
+def _above_leaves(n: Node) -> "list[Node]":
+    am = _above_more(n)
+    return am.leaves() if am is not None else []
+
+
+def level_of(n: Node) -> int:
+    """Storey index of n (0 = ground); mirrors ``Urb::Quad::Level``
+    (= number of levels below n's level root)."""
+    lr = _level_root(n)
+    i = 0
+    while lr.below is not None:
+        lr = lr.below
+        i += 1
+    return i
+
+
+def is_covered(n: Node) -> bool:
+    """Any leaf above n is indoor; mirrors ``Urb::Dom::Is_Covered``."""
+    return any(not is_outside(lf) for lf in _above_leaves(n))
+
+
 def is_supported(n: Node) -> bool:
     """All leaves below n are indoor; mirrors ``Urb::Dom::Is_Supported``."""
     bl = _below_leaves(n)
@@ -247,6 +289,23 @@ def is_unsupported(n: Node) -> bool:
     return bool(bl) and all(is_outside(lf) for lf in bl)
 
 
+def is_usable(n: Node) -> bool:
+    """Indoors, or on the ground floor, or supported by building below;
+    mirrors ``Urb::Dom::Is_Usable``."""
+    if not is_outside(n):
+        return True
+    if level_of(n) == 0:
+        return True
+    return is_supported(n)
+
+
+def is_circulation(n: Node) -> bool:
+    """Usable and type 'c'/'s'; mirrors ``Urb::Dom::Is_Circulation``."""
+    if not is_usable(n):
+        return False
+    return n.type is not None and n.type[0].lower() in ("c", "s")
+
+
 # --------------------------------------------------------------------------- #
 # Merge_Divided (Urb::Dom::Merge_Divided)
 # --------------------------------------------------------------------------- #
diff --git a/src/homemaker/fitness.py b/src/homemaker/fitness.py
new file mode 100644
index 0000000..a982721
--- /dev/null
+++ b/src/homemaker/fitness.py
@@ -0,0 +1,488 @@
+"""Native port of Urb's programme-driven fitness: leaf quality terms + cost model.
+
+Scope (homemaker-py-gnw): per-leaf quality factors (perpendicular, proportion,
+size, width, crinkliness, daylight, access), the programme-driven parameter
+lookup chain (``get_space_params``), value rates, and the cost denominator
+(per-leaf area costs, interior/exterior wall edge costs, boundary costs).
+Storey/building checks, staircases, failure stacking and final assembly are
+homemaker-py-hgg; corpus-parity validation is homemaker-py-uxz.
+
+Source of truth: ``Urb::Dom::Fitness::{Base,Leaf,Storey,ProgrammeDriven}``.
+
+DESCOPE (DESIGN.md §6, decision 2026-06-12): this ports *simple* crinkliness —
+the CIEsky illumination factor is pinned to 1, exactly what Urb computes under
+``URB_NO_OCCLUSION=1``.  ``quality_daylight`` is likewise pinned to 1.  Parity
+targets the *flagged* oracle, never stock Urb.
+
+Call ``dom.merge_divided(root)`` and rebuild graphs before ``process_storey``
+— storey processing runs on the MERGED tree (two-phase pattern, see graph.py).
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass, field
+from pathlib import Path
+
+import networkx as nx
+import yaml
+
+from . import dom as dom_mod
+from . import geometry
+from .dom import Node
+
+FAIL_THRESHOLD = 0.1  # Urb::Dom::Fitness::Base
+
+# Urb::Dom::Fitness::Base $CONF — keep values byte-identical to the Perl
+# expressions (5.0/6 etc. evaluate to the same IEEE doubles in both languages).
+CONF_DEFAULTS: dict = {
+    "value_inside": 300.0,
+    "value_circulation": 50.0,
+    "value_outside": 100.0,
+    "value_supported": 300.0,
+    "storey_limit": 4,
+    "storey_minimum": 2,
+    "latitude": 53.3814,
+    "door_width": 1.2,
+    "plot_ratio": [2.00, 0.50],
+    "ratio_outside": [0.33, 0.15],
+    "ratio_circulation": [0.00, 0.20],
+    "uncrinkliness": [5.0 / 6, 1.1 / 3],
+    "uncrinkliness_circulation": [5.0 / 6, 1.1 / 3],
+    "size_circulation": [0.0, 14.0],
+    "size_inside": [16.0, 3.5],
+    "proportion_outside": [1.5, 50],
+    "proportion_circulation": [1.5, 0.5],
+    "proportion_inside": [1.5, 0.5],
+    "width_outside": [3.0, 0.3],
+    "width_circulation": [2.4, 0.2],
+    "width_inside": [4.0, 1.0],
+    "perpendicular_inside": 0.3,
+    "perpendicular_outside": 10.0,
+    "allow_sahn_circulation": 0,
+    "force_roof_garden": 1,
+    "evaluate_room_types": 1,
+}
+
+# Urb::Dom::Fitness::Base $COST
+COST_DEFAULTS: dict = {
+    "plot": 10.0,
+    "outside_covered_supported": 210.0,
+    "outside_covered": 110.0,
+    "outside_supported": 110.0,
+    "outside": 10.0,
+    "inside": 200.0,
+    "interior_wall": 200.0 / 3,
+    "exterior_wall": 100.0,
+    "boundary": 50.0 / 3,
+    "boundary_wall": 400.0 / 3,
+}
+
+# ProgrammeDriven::default_params ultimate fallbacks
+_PARAM_FALLBACKS = {
+    "size": [16.0, 3.5],
+    "width": [4.0, 1.0],
+    "proportion": [1.5, 0.5],
+}
+
+_E = 2.718281828  # Urb::Math::gaussian uses this truncated e, not math.e
+
+
+def gaussian(x: float, a: float, b: float, c: float) -> float:
+    """Bit-faithful port of ``Urb::Math::gaussian`` (note the truncated e)."""
+    return a * (_E ** (0 - ((x - b) ** 2 / (2 * c * c))))
+
+
+def load_config(directory: str | Path) -> tuple[dict, dict]:
+    """Load (patterns, costs) config for a corpus directory, mirroring
+    ``urb-fitness.pl``: project-level ``../<name>.config`` first, then the
+    local file's keys override it."""
+    directory = Path(directory)
+    conf: dict = {}
+    cost: dict = {}
+    for target, name in ((conf, "patterns.config"), (cost, "costs.config")):
+        for p in (directory.parent / name, directory / name):
+            if p.is_file():
+                with open(p) as fh:
+                    target.update(yaml.safe_load(fh) or {})
+    return conf, cost
+
+
+@dataclass
+class LeafEval:
+    level: int
+    id: str
+    type: str
+    area: float
+    rate: float
+    quality: float
+    factors: dict[str, float] = field(default_factory=dict)
+
+
+@dataclass
+class StoreyEval:
+    cost: float
+    value: float
+    leaves: list[LeafEval] = field(default_factory=list)
+
+
+def _t0(n: Node) -> str:
+    """First char of the type, lowercased ('' if untyped) — Urb's /^x/i tests."""
+    return n.type[0].lower() if n.type else ""
+
+
+def _height(n: Node) -> float:
+    """Floor-to-floor height of n's level; mirrors ``Urb::Quad::Height``."""
+    h = dom_mod._level_root(n).height
+    return h if h is not None else 3.0
+
+
+def _perimeter(n: Node) -> dict:
+    """Perimeter dict from the lowest level root (``Urb::Quad::Perimeter``)."""
+    lr = dom_mod._level_root(n)
+    while lr.below is not None:
+        lr = lr.below
+    return lr.perimeter or {}
+
+
+class Fitness:
+    """Programme-driven leaf quality + cost evaluation.
+
+    ``conf`` is the parsed patterns.config mapping (including ``spaces``);
+    ``cost`` the costs.config mapping.  Lookup falls back to the Base.pm
+    defaults, as ``Urb::Dom::Fitness::Base::Conf/Cost`` do.
+    """
+
+    def __init__(self, conf: dict | None = None, cost: dict | None = None):
+        self._conf = conf or {}
+        self._cost = cost or {}
+        self.spaces: dict = self._conf.get("spaces") or {}
+
+    def conf(self, key: str):
+        v = self._conf.get(key)
+        if v is not None:
+            return v
+        return CONF_DEFAULTS.get(key)
+
+    def cost(self, key: str) -> float:
+        v = self._cost.get(key)
+        if v is not None:
+            return v
+        return COST_DEFAULTS.get(key, 0.0)
+
+    def preprocess_building(self, root: Node) -> None:
+        """Sahn-to-Outside type conversion (``Building.pm::preprocess_building``).
+        Run BEFORE graph build and merge_divided — it changes merge outcomes."""
+        if self.conf("allow_sahn_circulation"):
+            return
+        for lvl in dom_mod.levels(root):
+            for leaf in lvl.leaves():
+                if _t0(leaf) == "s":
+                    leaf.type = "O"
+
+    # ------------------------------------------------------------------ #
+    # Programme-driven parameter lookup (ProgrammeDriven.pm:29-69)
+    # ------------------------------------------------------------------ #
+
+    def get_space_params(self, code: str, param: str) -> list[float]:
+        c0 = code[0].lower() if code else ""
+        if c0 == "c":
+            v = self.conf(f"{param}_circulation")
+            if v is not None:
+                return v
+        if c0 in ("o", "s"):
+            v = self.conf(f"{param}_outside")
+            if v is not None:
+                return v
+        sp = self.spaces.get(code)  # exact-key match, as in Perl
+        if sp is not None and param in sp:
+            return sp[param]
+        v = self.conf(f"{param}_inside")
+        if v is not None:
+            return v
+        return _PARAM_FALLBACKS.get(param)
+
+    # ------------------------------------------------------------------ #
+    # Quality terms (Leaf.pm)
+    # ------------------------------------------------------------------ #
+
+    def quality_perpendicular(self, leaf: Node) -> float:
+        sigma = self.conf(
+            "perpendicular_outside" if dom_mod.is_outside(leaf) else "perpendicular_inside"
+        )
+        score = 1.0
+        for i in range(4):
+            # 1.570796: Urb::Dom::Perpendicular hard-codes this, not pi/2
+            score *= gaussian(geometry.angle(leaf, i), 1.0, 1.570796, sigma)
+        return score
+
+    def quality_proportion(self, leaf: Node) -> float:
+        t0 = _t0(leaf)
+        if t0 in ("o", "s"):
+            params = self.conf("proportion_outside")
+        elif t0 == "c":
+            params = self.conf("proportion_circulation")
+        else:
+            params = self.get_space_params(leaf.type, "proportion")
+        aspect = geometry.aspect(leaf)
+        if aspect < params[0]:
+            return 1.0
+        return gaussian(aspect, 1.0, params[0], params[1])
+
+    def quality_size(self, leaf: Node) -> float:
+        t0 = _t0(leaf)
+        if t0 in ("o", "s"):
+            return 1.0
+        if t0 == "c":
+            params = self.conf("size_circulation")
+        else:
+            params = self.get_space_params(leaf.type, "size")
+        return gaussian(geometry.area(leaf), 1.0, params[0], params[1])
+
+    def quality_width(self, leaf: Node) -> float:
+        t0 = _t0(leaf)
+        if (
+            t0 in ("o", "s")
+            and not dom_mod.is_covered(leaf)
+            and not dom_mod.is_supported(leaf)
+            and dom_mod.level_of(leaf)
+        ):
+            return 1.0
+        if t0 in ("o", "s"):
+            params = self.conf("width_outside")
+        elif t0 == "c":
+            params = self.conf("width_circulation")
+        else:
+            params = self.get_space_params(leaf.type, "width")
+        width = geometry.length_narrowest(leaf)
+        if width > params[0]:
+            return 1.0
+        return gaussian(width, 1.0, params[0], params[1])
+
+    # --- simple crinkliness (URB_NO_OCCLUSION: illumination factor = 1) --- #
+
+    def area_outside(self, leaf: Node, G: nx.Graph, groups: dict) -> float:
+        """Illuminated external wall area; ``Urb::Dom::Area_Outside`` with the
+        CIEsky illumination factor pinned to 1 (simple crinkliness)."""
+        length = 0.0
+        for nb in G.neighbors(leaf):
+            if not dom_mod.is_outside(nb) or dom_mod.is_covered(nb):
+                continue
+            # Faithful loop over all internal boundaries: Overlap() is > 0
+            # only on a boundary both quads actually share an edge of.
+            for contributors in groups.values():
+                if geometry.boundary_pair_overlap(contributors, leaf, nb) > 0:
+                    length += G[leaf][nb]["width"]
+        perimeter = _perimeter(leaf)
+        for e in range(4):
+            bid = geometry.boundary_id(leaf, e)
+            if bid not in geometry._EXTERNAL:
+                continue
+            ptype = (perimeter.get(bid) or "").lower()
+            if ptype in ("private", "fortified"):
+                continue
+            length += geometry.edge_length(leaf, e)
+        return length * _height(leaf)
+
+    def crinkliness(self, leaf: Node, G: nx.Graph, groups: dict) -> float:
+        area = geometry.area(leaf)
+        if not area:
+            return 9999999999
+        return self.area_outside(leaf, G, groups) / area
+
+    def quality_uncrinkliness(self, leaf: Node, G: nx.Graph, groups: dict) -> float:
+        if dom_mod.is_outside(leaf) and not dom_mod.is_covered(leaf):
+            return 1.0
+        key = "uncrinkliness_circulation" if dom_mod.is_circulation(leaf) else "uncrinkliness"
+        distance, sigma = self.conf(key)
+        crink = self.crinkliness(leaf, G, groups)
+        if not crink:
+            return 0.0
+        return gaussian(1 / crink, 1.0, distance, sigma)
+
+    # --- access --- #
+
+    def neighbour_types(self, leaf: Node, G: nx.Graph) -> list[str]:
+        return sorted(nb.type or "" for nb in G.neighbors(leaf) if dom_mod.is_usable(nb))
+
+    def access(self, leaf: Node, G: nx.Graph) -> list[str]:
+        """Useful circulation/access neighbour types; ``Urb::Dom::Access``."""
+        types = self.neighbour_types(leaf, G)
+        if _t0(leaf) == "k":
+            return [t for t in types if t and t[0].lower() in ("l", "c", "s")]
+        if dom_mod.is_outside(leaf) or dom_mod.is_circulation(leaf):
+            return types
+        return [t for t in types if t and t[0].lower() in ("c", "s")]
+
+    # ------------------------------------------------------------------ #
+    # Leaf evaluation (Leaf.pm::evaluate_leaf)
+    # ------------------------------------------------------------------ #
+
+    def evaluate_leaf(
+        self, leaf: Node, G: nx.Graph, level_id: int, groups: dict, fail
+    ) -> tuple[float, dict[str, float]]:
+        """Return (quality, per-factor dict); appends failures via ``fail``.
+
+        Factor order and fail strings mirror ``evaluate_leaf`` exactly.
+        """
+        lid = leaf.id
+        factors: dict[str, float] = {}
+        quality = 1.0
+
+        f = self.quality_perpendicular(leaf)
+        if f < FAIL_THRESHOLD:
+            fail(f"{level_id}/{lid} perpendicular")
+        factors["perpendicular"] = f
+        quality *= f
+
+        f = self.quality_proportion(leaf)
+        if f < FAIL_THRESHOLD:
+            fail(f"{level_id}/{lid} proportion")
+        factors["proportion"] = f
+        quality *= f
+
+        f = self.quality_size(leaf)
+        if f < FAIL_THRESHOLD:
+            fail(f"{level_id}/{lid} size")
+        factors["size"] = f
+        quality *= f
+
+        f = self.quality_width(leaf)
+        if f < FAIL_THRESHOLD:
+            fail(f"{level_id}/{lid} width")
+        factors["width"] = f
+        quality *= f
+
+        f = self.quality_uncrinkliness(leaf, G, groups)
+        if f < FAIL_THRESHOLD:
+            fail(f"{level_id}/{lid} crinkliness")
+        factors["crinkliness"] = f
+        quality *= f
+
+        # Daylight pinned to 1 — URB_NO_OCCLUSION semantics (DESIGN.md §6).
+        factors["daylight"] = 1.0
+
+        if len(self.access(leaf, G)) > 0:
+            f = 1.0
+        elif not dom_mod.level_of(leaf) and dom_mod.is_outside(leaf):
+            f = 1.0
+        else:
+            f = 0.01
+            fail(f"{level_id}/{lid} access")
+        factors["access"] = f
+        quality *= f
+
+        return quality, factors
+
+    # ------------------------------------------------------------------ #
+    # Value rates and costs (Leaf.pm:146-251, Storey.pm:122-147)
+    # ------------------------------------------------------------------ #
+
+    def value_rate(self, leaf: Node) -> float:
+        t0 = _t0(leaf)
+        if t0 in ("o", "s") and dom_mod.level_of(leaf) == 0:
+            return self.conf("value_outside")
+        if t0 in ("o", "s"):
+            return self.conf("value_supported")
+        if t0 == "c":
+            return self.conf("value_circulation")
+        return self.conf("value_inside")
+
+    def leaf_cost(self, leaf: Node) -> float:
+        if dom_mod.is_outside(leaf):
+            covered = dom_mod.is_covered(leaf)
+            supported = dom_mod.is_supported(leaf)
+            if covered and supported:
+                rate = self.cost("outside_covered_supported")
+            elif covered:
+                rate = self.cost("outside_covered")
+            elif supported:
+                rate = self.cost("outside_supported")
+            else:
+                rate = self.cost("outside")
+        else:
+            rate = self.cost("inside")
+        return rate * geometry.area(leaf)
+
+    def edge_cost(self, G: nx.Graph, a: Node, b: Node, fail) -> float:
+        """Interior/exterior wall cost for one graph edge
+        (``Storey.pm::calculate_edge_cost``)."""
+        height = _height(a)
+        a_out, b_out = dom_mod.is_outside(a), dom_mod.is_outside(b)
+        if a_out and b_out:
+            rate = 0.0
+        elif not a_out and not b_out:
+            rate = self.cost("interior_wall")
+        else:
+            rate = self.cost("exterior_wall")
+        width = G[a][b]["width"]
+        if width > 8.0 and rate > 0.0:
+            fail(f"{dom_mod.level_of(a)}/{a.id} {b.id} edge too long")
+        return rate * width * height
+
+    def outside_edge_cost(self, leaf: Node, fail) -> float:
+        """Plot-boundary cost for a leaf's external edges
+        (``Leaf.pm::calculate_outside_edge_cost``)."""
+        rate = self.cost("boundary") if dom_mod.is_outside(leaf) else self.cost("boundary_wall")
+        length = 0.0
+        for e in range(4):
+            if geometry.boundary_id(leaf, e) not in geometry._EXTERNAL:
+                continue
+            edge_len = geometry.edge_length(leaf, e)
+            length += edge_len
+            if dom_mod.is_outside(leaf):
+                continue
+            if edge_len > 8.0:
+                fail(f"{dom_mod.level_of(leaf)}/{leaf.id} outside edge too long")
+        return rate * length * _height(leaf)
+
+    # ------------------------------------------------------------------ #
+    # Storey processing (Storey.pm::process_storey — cost/value/leaf scope)
+    # ------------------------------------------------------------------ #
+
+    def process_storey(self, level_root: Node, G: nx.Graph, level_id: int, fail) -> StoreyEval:
+        """Per-storey cost, value and leaf evaluations on the MERGED tree.
+
+        Covers the cost/value accumulation and per-leaf checks of
+        ``process_storey``; circulation connectivity, roof-garden, stair fit
+        and tracking-driven building checks are homemaker-py-hgg.
+        """
+        groups = geometry.boundary_groups(level_root)
+        cost = 0.0
+        value = 0.0
+        leaves_eval: list[LeafEval] = []
+
+        for leaf in level_root.leaves():
+            if dom_mod.is_outside(leaf) and dom_mod.is_covered(leaf) and level_id:
+                if not dom_mod.is_supported(leaf):
+                    fail(f"{level_id}/{leaf.id} unsupported covered outside")
+                fail(f"{level_id}/{leaf.id} covered outside above ground")
+
+            cost += self.leaf_cost(leaf)
+            if not dom_mod.is_usable(leaf):
+                continue
+
+            quality, factors = self.evaluate_leaf(leaf, G, level_id, groups, fail)
+            rate = self.value_rate(leaf)
+            value += quality * rate * geometry.area(leaf)
+            leaves_eval.append(
+                LeafEval(
+                    level=level_id,
+                    id=leaf.id,
+                    type=leaf.type or "",
+                    area=geometry.area(leaf),
+                    rate=rate,
+                    quality=quality,
+                    factors=factors,
+                )
+            )
+
+        for a, b in G.edges():
+            cost += self.edge_cost(G, a, b, fail)
+        for leaf in level_root.leaves():
+            cost += self.outside_edge_cost(leaf, fail)
+
+        return StoreyEval(cost=cost, value=value, leaves=leaves_eval)
+
+    def plot_cost(self, root: Node) -> float:
+        """The 'initial cost' term: plot rate x lowest-root area."""
+        return self.cost("plot") * geometry.area(root)
diff --git a/src/homemaker/geometry.py b/src/homemaker/geometry.py
index 2f4f761..2759279 100644
--- a/src/homemaker/geometry.py
+++ b/src/homemaker/geometry.py
@@ -86,7 +86,11 @@ def coord_b(n: Node) -> Point:
 
 
 def _dist(a: Point, b: Point) -> float:
-    return math.hypot(a[0] - b[0], a[1] - b[1])
+    # NOT math.hypot: Urb::Math::distance_2d is sqrt(dx**2 + dy**2) and the
+    # two differ in the last ULP. Boundary overlap tests feed the difference
+    # of near-equal lengths into a > 0 predicate (Urb::Boundary::Overlap), so
+    # a 1-ULP deviation flips adjacency decisions on exactly-touching quads.
+    return math.sqrt((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2)
 
 
 def _triangle_area(a: Point, b: Point, c: Point) -> float:
@@ -107,6 +111,29 @@ def edge_length(n: Node, idx: int) -> float:
     return _dist(coordinate(n, idx), coordinate(n, (idx + 1) % 4))
 
 
+def angle(n: Node, idx: int) -> float:
+    """Interior angle at corner ``idx`` in radians (``Urb::Quad::Angle``,
+    cosine rule). Clamped acos argument — Perl leaves it unclamped but only a
+    degenerate quad would push it out of [-1, 1]."""
+    a = edge_length(n, idx)
+    b = edge_length(n, (idx + 3) % 4)
+    c = _dist(coordinate(n, (idx + 1) % 4), coordinate(n, (idx + 3) % 4))
+    return math.acos(max(-1.0, min(1.0, (a * a + b * b - c * c) / (2 * a * b))))
+
+
+def aspect(n: Node) -> float:
+    """Plan aspect ratio, always >= 1 (``Urb::Quad::Aspect``)."""
+    asp = (edge_length(n, 0) + edge_length(n, 2)) / (edge_length(n, 1) + edge_length(n, 3))
+    if 0 < asp < 1:
+        asp = 1 / asp
+    return asp
+
+
+def length_narrowest(n: Node) -> float:
+    """Shortest of the four edge lengths (``Urb::Quad::Length_Narrowest``)."""
+    return min(edge_length(n, i) for i in range(4))
+
+
 # --------------------------------------------------------------------------- #
 # Plot wall inset (Urb::Quad::Coordinate_Offset, used on the root in Urb::Dom).
 # Positive offset moves a corner outward, negative inward. Computed per corner
@@ -210,6 +237,44 @@ def _edge_overlap(a: Node, edge_a: int, b: Node, edge_b: int) -> float:
     return max(0.0, len_a + len_b - max_dist)
 
 
+_EXTERNAL = frozenset("abcd")  # single-char external boundary ids
+
+
+def boundary_groups(level_root: Node) -> dict[str, list[tuple[Node, int]]]:
+    """Group (leaf, edge) pairs by shared internal boundary id; the data half
+    of ``Urb::Quad::Calc_Boundaries`` (external 'a'-'d' boundaries excluded —
+    Urb's ``Boundary::Overlap`` returns 0 for them anyway).
+
+    Membership test uses the frozenset, NOT ``bid not in "abcd"`` — the latter
+    is a substring check and silently drops the root-division boundary ('').
+    """
+    from collections import defaultdict
+
+    groups: dict[str, list[tuple[Node, int]]] = defaultdict(list)
+    for leaf in level_root.leaves():
+        for edge in range(4):
+            bid = boundary_id(leaf, edge)
+            if bid not in _EXTERNAL:
+                groups[bid].append((leaf, edge))
+    return groups
+
+
+def boundary_pair_overlap(contributors: list[tuple[Node, int]], a: Node, b: Node) -> float:
+    """Overlap of quads ``a`` and ``b`` on one boundary's contributor list;
+    mirrors ``Urb::Boundary::Overlap`` (last matching edge wins, as in the
+    Perl loop). Returns 0.0 if either quad has no edge on this boundary.
+    """
+    edge_a = edge_b = None
+    for leaf, edge in contributors:
+        if leaf is a:
+            edge_a = edge
+        if leaf is b:
+            edge_b = edge
+    if edge_a is None or edge_b is None:
+        return 0.0
+    return _edge_overlap(a, edge_a, b, edge_b)
+
+
 def leaf_graph(level_root: Node, door_width: float = 1.2):  # -> nx.Graph
     """Leaf-adjacency graph for one storey; mirrors ``Urb::Quad::Graph``.
 
@@ -220,19 +285,9 @@ def leaf_graph(level_root: Node, door_width: float = 1.2):  # -> nx.Graph
     never edges.  A single-leaf storey gets one isolated vertex.
     """
     import networkx as nx
-    from collections import defaultdict
 
-    _external = frozenset("abcd")  # single-char external boundary ids
     leaves = level_root.leaves()
-    # Group (leaf, edge) pairs by shared boundary id.
-    # Use frozenset membership, NOT 'bid not in "abcd"' — the latter is a
-    # substring check and silently drops the root-division boundary ('').
-    groups: dict[str, list[tuple[Node, int]]] = defaultdict(list)
-    for leaf in leaves:
-        for edge in range(4):
-            bid = boundary_id(leaf, edge)
-            if bid not in _external:
-                groups[bid].append((leaf, edge))
+    groups = boundary_groups(level_root)
 
     G: nx.Graph = nx.Graph()
     for leaf in leaves: