city.js
544 lines · in natura-src-1.0.0.zip as natura/city.js · a module · its entry in the module reference · its test: city-test.js
// city.js -- Natura CITY (the SOCIETY ontology's BUILT-WORLD rung; the first HUMAN-WORLD domain the
// framework has touched -- not a nature system but the artefact a society raises on the land it lives on).
// A city is EMERGENT: nobody plans the core, yet a monocentric downtown rises out of a flat, uniform plain
// because every parcel answers the same local question -- "what is this land worth, and what should stand on
// it?" -- and the answers COUPLE. Value depends on ACCESS; access runs along ROADS whose cost rises as they
// CONGEST; development follows value; density near you RAISES your value (AGGLOMERATION); a raised core draws
// population that commutes back along the roads that then congest harder. Out of that loop, with no central
// hand, comes the shape of a city: a dense downtown, a gradient of density falling with distance, ribbon
// development crawling out along a new road, and -- once the centre chokes on its own traffic -- SUB-CENTRES
// budding off it (monocentric -> polycentric). Buildings AGE; a neglected, under-valued district DECAYS and
// its BLIGHT spreads to its neighbours; and when access and demand lift the land's worth back above what
// stands on it, the blight is CLEARED and REDEVELOPED denser than before -- the renewal cycle that turns a
// city over from within. This is the Alonso/von-Thunen bid-rent gradient, Krugman's agglomeration, a gravity
// traffic assignment, and the filtering/renewal cycle, all as ONE coupled grid.
//
// SELF-CONTAINED: it seeds its own plain, its own founding road + settlement, and its own exogenous regional
// growth (a slow inflow of households + firms looking for somewhere to land), so it RUNS and TESTS standalone.
// It reads no peer system, but exposes rich fields() (use / density / landValue / accessibility / blight /
// congestion / population) that a renderer or a coupled world can read. Deterministic: all randomness comes
// from init(rng) and ctx.rng; time is dt seconds. UMD (node require / window.City).
(function (root, factory) {
if (typeof module !== 'undefined' && module.exports) module.exports = factory();
else root.City = factory();
})(typeof self !== 'undefined' ? self : this, function () {
'use strict';
const clamp = (v, a, b) => v < a ? a : v > b ? b : v;
// land uses
const UND = 0, RES = 1, COM = 2, IND = 3;
function makeCity(config) {
config = config || {};
const N = config.N || 24; // grid is N x N parcels
const P = Object.assign({
// --- exogenous regional growth: households + firms arrive from "outside" looking to locate ------
pop0: 40, // founding settlement population (seeded at centre)
jobs0: 12,
popGrowth: 0.9, // households arriving per sim-DAY (regional inflow)
jobGrowth: 0.02, // rate firms open/close toward the labour-market target (per day) -- slow, so jobs don't thrash
jobsPerWorker: 0.85, // balanced labour market: regional jobs track this ratio of housed workers
popCap: 6000, jobCap: 100000, // regional household pool ceiling (a finite region, so the city reaches a finite size)
// --- accessibility (gravity + distance decay over the road network) ----------------------------
beta: 0.35, // distance-decay: access = sum_j opportunity_j * exp(-beta * cost_ij)
relaxIters: 26, // shortest-path relaxation sweeps (enough to cross an N=24 grid)
offRoadPenalty: 6.0, // moving a cell WITHOUT a road under it is this many times costlier than a road hop
roadBaseCost: 1.0, // free-flow cost of traversing one road cell
// --- congestion: a road's cost rises with the flow it carries (BPR-like) -----------------------
congestAlpha: 1.4, congestBeta: 2.4, // reported congestion ratio = 1 + alpha*(flow/cap)^beta
congCostCap: 2.5, // the congestion PENALTY that enters the access solve saturates near this (a jam ~3.5x cost, not infinite -> no core collapse)
roadCap: 16, // free-flow capacity of a road cell (flow beyond this bites)
flowSmooth: 0.2, // congestion adapts toward the new load at this rate/step (slow -> hysteresis, no thrash)
// --- land value (bid-rent): access + agglomeration, minus a distaste for industry nearby -------
valAccessRes: 0.8, // households bid MODERATELY + broadly on access to jobs (win the ring)
valAccessCom: 2.6, // firms bid STEEPLY on centrality (access^1.5) -> a tight core
valAccessInd: 0.5, // industry bids shallowly, mostly on labour access
indEdgeBonus: 0.02, // industry tolerates cheap peripheral land (a small floor so it wins SOME edge)
agglomGain: 0.85, // COMMERCIAL value rises with nearby developed density (the core-forming feedback)
agglomRadius: 3,
industryNuisance: 0.35, // residential value falls near industry (a repulsion that segregates uses)
demandGain: 0.5, // how strongly unmet regional demand lifts rents (market tightness -> value)
valDecay: 0.15, // land value relaxes toward its target at this rate/step (sticky, not instant)
// --- development + density ----------------------------------------------------------------------
buildThresh: 0.05, // land value above which an undeveloped parcel is worth building on
buildRate: 0.10, // per-step probability a profitable parcel breaks ground (construction is not instant)
maxDensity: 12, // storeys / intensity cap
densRate: 0.06, // how fast density climbs toward what value supports (construction time)
densSupport: 9, // a unit of land value supports this much density (value*densSupport = target floors)
capacityPerDens: 9, // residents (or jobs) a unit of density houses
moveRate: 0.12, // how fast households/firms migrate toward their target parcel (finite, not instant)
resBase: 0.35, // base willingness of a household to occupy a built parcel (breaks the ghost-town cold start)
jobBase: 0.3, // base willingness of a firm to occupy a built parcel
leaveStick: 0.35, // households leave at this fraction of the move-in rate (moving is costly -> hysteresis)
// --- aging / blight / renewal -------------------------------------------------------------------
ageRate: 1.0, // building age advances at 1x real years
blightAge: 28, // years before age alone starts to erode a building
blightRate: 0.05, // /step blight accrual when under-valued or old
blightSpread: 0.10, // blight bleeds to neighbours at this fraction (contagion)
blightValueDrag: 0.6, // blight this severe drags realised land value down by up to this fraction
abandonBlight: 0.85, // blight past this and the parcel empties (goes vacant/undeveloped) unless renewed
renewValueRatio: 1.5, // renewal fires when land value exceeds what STANDS there by this ratio (worth more cleared)
renewMinBlight: 0.35, // ...and only once the parcel is at least this blighted/aged (ripe for turnover)
// --- roads --------------------------------------------------------------------------------------
roadSeed: 'cross' // founding road layout: 'cross' = one central N-S + E-W artery through the seed
}, config.params || config);
// allow flat overrides too (makeCity({N, popGrowth}) as well as makeCity({params:{...}}))
for (const k in config) if (k in P && k !== 'params') P[k] = config[k];
// OPTIONAL coupling: read a traffic microsim's analytic fundamental diagram (vdCost) as the volume-delay
// curve, grounding congestion in real car-following instead of the assumed BPR exponent. Default null ->
// the BPR path, byte-identical to the standalone city (so the existing tests are untouched).
const congestFrom = config.congestFrom || (config.params && config.params.congestFrom) || null;
const idx = (i, j) => j * N + i;
const inB = (i, j) => i >= 0 && i < N && j >= 0 && j < N;
const mid = (N - 1) / 2;
// ---- accessibility solver: shortest travel COST from every cell, by relaxation ---------------------
// We compute, for each cell, its network distance to the two opportunity fields (jobs, population)
// using a multi-source distance-decay. Rather than an all-pairs solve (N^4), we run a small number of
// "potential" relaxations that spread an opportunity's reach along cheap edges -- a discrete analogue of
// solving access = sum_j O_j exp(-beta d_ij) by diffusion on the road graph. Congestion makes edge cost
// depend on flow, so accessibility is ENDOGENOUS (it falls where the roads it rides are jammed).
function edgeCost(st, i, j, ni, nj) {
// cost to step from (i,j) to neighbour (ni,nj): cheap if BOTH ends carry a road, else off-road penalty
const onRoad = st.road[idx(i, j)] > 0 && st.road[idx(ni, nj)] > 0;
if (!onRoad) return P.roadBaseCost * P.offRoadPenalty;
return st.roadCost[idx(ni, nj)]; // congested per-cell traversal cost (updated each step)
}
// GRAVITY ACCESSIBILITY as a discounted POTENTIAL solved by relaxation. access[n] = the opportunity
// (jobs, or people) reachable from n, each unit discounted by exp(-beta * network-cost) to reach it:
// the classic gravity/Hansen accessibility A_n = sum_m O_m exp(-beta d_nm). We solve it as a fixed point
// instead of an O(N^4) all-pairs sum: a cell's potential is its OWN opportunity plus the best discounted
// potential flowing in from a road neighbour, pot[n] = opp[n] + max_nbr( pot[nbr] * exp(-beta*edgeCost) ).
// Using the max-plus (best corridor) keeps it bounded and mass-AWARE -- a tall potential radiates from a
// big cluster and decays with network distance, so being near the MASS matters, not merely near any one
// job. Because edgeCost carries CONGESTION, the accessibility a jam produces is endogenous: a choked
// artery raises the decay along it, lowering the access of everything downstream. This is what gives a
// real monocentric gradient (peak over the cluster, falling outward) that a nearest-distance proxy loses
// the moment opportunity spreads out.
function solvePotential(st, opp, pot) {
for (let n = 0; n < N * N; n++) pot[n] = opp[n];
for (let it = 0; it < P.relaxIters; it++) {
let maxChange = 0;
for (let pass = 0; pass < 2; pass++) {
const lo = pass === 0 ? 0 : N * N - 1, hi = pass === 0 ? N * N : -1, dir = pass === 0 ? 1 : -1;
for (let n = lo; n !== hi; n += dir) {
const i = n % N, j = (n / N) | 0;
let best = 0;
if (i > 0) { const v = pot[idx(i - 1, j)] * Math.exp(-P.beta * edgeCost(st, i - 1, j, i, j)); if (v > best) best = v; }
if (i < N - 1) { const v = pot[idx(i + 1, j)] * Math.exp(-P.beta * edgeCost(st, i + 1, j, i, j)); if (v > best) best = v; }
if (j > 0) { const v = pot[idx(i, j - 1)] * Math.exp(-P.beta * edgeCost(st, i, j - 1, i, j)); if (v > best) best = v; }
if (j < N - 1) { const v = pot[idx(i, j + 1)] * Math.exp(-P.beta * edgeCost(st, i, j + 1, i, j)); if (v > best) best = v; }
const nv = opp[n] + best;
const ch = Math.abs(nv - pot[n]); if (ch > maxChange) maxChange = ch;
pot[n] = nv;
}
}
if (maxChange < 1e-4) break;
}
}
// ---- helpers to house population/jobs on the grid --------------------------------------------------
const capacityAt = (st, n) => st.use[n] === RES || st.use[n] === COM || st.use[n] === IND ? st.density[n] * P.capacityPerDens : 0;
return {
name: config.name || 'city',
domain: 'city',
scale: 'built',
params: P,
init(rng) {
const st = {
N,
use: new Uint8Array(N * N),
density: new Float32Array(N * N),
age: new Float32Array(N * N),
landValue: new Float32Array(N * N),
blight: new Float32Array(N * N),
road: new Uint8Array(N * N),
roadCost: new Float32Array(N * N), // congested per-cell road traversal cost
flow: new Float32Array(N * N), // traffic load on each road cell (smoothed)
access: new Float32Array(N * N), // accessibility to jobs (what households bid on)
accessPop: new Float32Array(N * N), // accessibility to population/customers (what firms bid on)
pop: new Float32Array(N * N), // residents housed per cell
jobs: new Float32Array(N * N), // jobs located per cell
// scratch (dedicated, no fragile aliasing): _potJobs/_potPop hold the gravity potentials, _use holds
// the winning-bid use per parcel, _load holds the traffic load, _wt a per-cell weight, _nb the
// neighbour-blight accumulation
_potJobs: new Float32Array(N * N), _potPop: new Float32Array(N * N),
_use: new Uint8Array(N * N), _load: new Float32Array(N * N),
_wt: new Float32Array(N * N), _nb: new Float32Array(N * N),
// aggregates / observables
totalPop: 0, totalJobs: 0, poolPop: P.pop0, poolJobs: P.jobs0, demandPop: 1, demandJob: 1, laborBase: P.pop0,
peakDensity: 0, developed: 0, meanBlight: 0, maxCongestion: 1, meanAccess: 0,
age_yr: 0, rng: null, evs: []
};
st.rng = rng;
// free-flow cost everywhere to start
for (let n = 0; n < N * N; n++) st.roadCost[n] = P.roadBaseCost;
// --- found the city: a central seed settlement + a founding road ------------------------------
const c = Math.round(mid);
if (P.roadSeed === 'cross') {
for (let k = 0; k < N; k++) { st.road[idx(c, k)] = 1; st.road[idx(k, c)] = 1; } // N-S + E-W arteries
}
// seed a small residential settlement + a few founding jobs at the crossroads
const seedR = 2;
let rc = 0;
for (let dj = -seedR; dj <= seedR; dj++) for (let di = -seedR; di <= seedR; di++) {
const i = c + di, j = c + dj; if (!inB(i, j)) continue; rc++;
st.use[idx(i, j)] = (di === 0 && dj === 0) ? COM : RES;
st.density[idx(i, j)] = (di === 0 && dj === 0) ? 2.0 : 1.0;
st.age[idx(i, j)] = 1 + rng() * 3;
}
st.jobs[idx(c, c)] = P.jobs0;
// scatter tiny founding population over the seed
let pr = P.pop0;
for (let dj = -seedR; dj <= seedR; dj++) for (let di = -seedR; di <= seedR; di++) {
const i = c + di, j = c + dj; if (!inB(i, j) || st.use[idx(i, j)] !== RES) continue;
const share = pr / Math.max(1, rc); st.pop[idx(i, j)] = share;
}
return st;
},
step(st, ctx, dt) {
st.evs.length = 0;
const rng = (ctx && ctx.rng) || st.rng;
const dd = clamp((dt || 0) / 86400, 0, 400); // step length in sim-days
const yrs = (dt || 0) / (365 * 86400);
st.age_yr += yrs;
// === 1. EXOGENOUS REGIONAL GROWTH ============================================================
// households + firms trickle in from outside, from a finite regional pool (logistic-capped so the
// inflow is a founding pressure, not a runaway). These are the demand the city must house/employ.
// The two pools are a LABOUR MARKET, not independent: firms want workers and workers want jobs, so
// the job pool is tethered to the housed population (jobs chase workers within a band, jobsPerWorker),
// and household inflow eases as the built housing stock fills -- so neither side runs away from the
// other and the city reaches a finite, balanced size instead of exploding.
st.poolPop += P.popGrowth * dd * (1 - st.poolPop / P.popCap);
// desired regional jobs = a fixed ratio of the workers actually housed (a balanced labour market),
// approached slowly (firms open/close toward that target, they do not teleport).
// the labour base is a SLOW average of the housed population (firms plan on the trend, not the
// month-to-month wobble), so a transient dip in occupancy does not evict a wave of jobs and touch
// off a boom-bust. This is the main damper on the coupled pop<->jobs oscillation.
st.laborBase += (st.totalPop - st.laborBase) * clamp(0.02 * dd, 0, 1);
const jobTarget = Math.max(P.jobs0, st.laborBase * P.jobsPerWorker);
st.poolJobs += (jobTarget - st.poolJobs) * clamp(P.jobGrowth * dd, 0, 1);
// DEMAND PRESSURE (market clearing): households/firms in the region that the city has not yet housed
// bid RENTS UP. A tight market (pool >> housed) lifts land value everywhere, pushing the build
// threshold outward and intensifying density -- so the city GROWS to absorb its people instead of
// stalling at its first cluster. Slack (housed catches up to the pool) lets pressure fall back to 1.
const dpTarget = 1 + P.demandGain * clamp((st.poolPop - st.totalPop) / Math.max(1, st.totalPop), 0, 4);
const djTarget = 1 + P.demandGain * clamp((st.poolJobs - st.totalJobs) / Math.max(1, st.totalJobs), 0, 4);
st.demandPop += (dpTarget - st.demandPop) * 0.25; // smoothed, so the market does not thrash the whole grid each step
st.demandJob += (djTarget - st.demandJob) * 0.25;
// === 2. ACCESSIBILITY over the (congested) road network ======================================
// opportunity fields: jobs (households commute TO jobs), population (firms sell TO people). The
// gravity potential radiates from wherever the mass concentrates and decays with network cost, so
// being near the CLUSTER is what scores -- the gradient survives the opportunity spreading out.
solvePotential(st, st.jobs, st._potJobs); // gravity potential of JOB access
solvePotential(st, st.pop, st._potPop); // gravity potential of POPULATION access
let totJobs = 0, totPop = 0, maxAj = 1e-6, maxAp = 1e-6;
for (let n = 0; n < N * N; n++) {
totJobs += st.jobs[n]; totPop += st.pop[n];
if (st._potJobs[n] > maxAj) maxAj = st._potJobs[n];
if (st._potPop[n] > maxAp) maxAp = st._potPop[n];
}
// NORMALISE each potential to its own peak -> a 0..1 accessibility SHAPE (where the core is), then
// scale by a bounded log market-size factor so a bigger city's centre is worth a bit more overall
// without the field exploding. Shape * bounded level = a gradient that sorts the grid into core+ring
// instead of every cell clearing the build threshold at once.
const mkt = Math.log10(10 + Math.max(totJobs, totPop)) / 2.2; // ~0.45 tiny .. ~2.3 huge, bounded
for (let n = 0; n < N * N; n++) {
st.access[n] = mkt * (st._potJobs[n] / maxAj);
st.accessPop[n] = mkt * (st._potPop[n] / maxAp);
}
// === 3. TRAFFIC ASSIGNMENT -> CONGESTION =====================================================
// a gravity traffic model: commuters flow from where they live to where the jobs are, along the
// cheapest roads. We load each road cell with a flow proportional to how much through-travel it
// carries, approximated by (local population + jobs it connects) attenuated by access, then update
// each road cell's cost with a BPR-like curve so busy roads get dearer -> access becomes endogenous.
// total commuters on the network = the workers who travel from home to a job each day (min of the two
// sides -- an unemployed worker or an unfilled job generates no commute). They ride the roads, and a
// road cell carries more of them the more CENTRAL it is (higher access = more shortest paths pass
// through it). So flow = commuters * centrality-share, a bounded gravity assignment: central arteries
// jam first, exactly where a real city does.
const commuters = Math.min(st.totalPop, st.totalJobs);
const load = st._load; load.fill(0);
let accSum = 0;
for (let n = 0; n < N * N; n++) if (st.road[n]) accSum += st.access[n] + st.accessPop[n];
accSum = Math.max(1e-6, accSum);
for (let n = 0; n < N * N; n++) {
if (st.road[n] === 0) continue;
// this road cell's share of all through-travel, by its centrality; plus the local trips it originates
const centrality = (st.access[n] + st.accessPop[n]) / accSum;
load[n] = commuters * centrality + 0.2 * (st.pop[n] + st.jobs[n]);
}
// smear load one ring along roads so arteries accumulate their corridor's traffic
for (let n = 0; n < N * N; n++) {
if (st.road[n] === 0) continue; const i = n % N, j = (n / N) | 0; let s = load[n], cnt = 1;
if (inB(i - 1, j) && st.road[idx(i - 1, j)]) { s += load[idx(i - 1, j)]; cnt++; }
if (inB(i + 1, j) && st.road[idx(i + 1, j)]) { s += load[idx(i + 1, j)]; cnt++; }
if (inB(i, j - 1) && st.road[idx(i, j - 1)]) { s += load[idx(i, j - 1)]; cnt++; }
if (inB(i, j + 1) && st.road[idx(i, j + 1)]) { s += load[idx(i, j + 1)]; cnt++; }
st.flow[n] += (s / cnt - st.flow[n]) * P.flowSmooth; // hysteresis: congestion builds/eases, not instant
}
// volume-delay: by default an assumed BPR power law. If COUPLED to a traffic microsim (congestFrom), read
// its analytic fundamental diagram (vdCost, indexed by utilization x = flow/capacity over [0, vdXmax])
// and use THAT as the congestion curve -- the city's macroscopic delay is then grounded in emergent
// car-following, exactly the layer traffic.js says the city lacks. Same input x, same saturation wrapper.
const tvd = congestFrom && ctx && ctx.fields[congestFrom] && ctx.fields[congestFrom].vdCost;
const tvdMax = tvd ? (ctx.fields[congestFrom].vdXmax || 2) : 2;
let maxCong = 1;
for (let n = 0; n < N * N; n++) {
if (st.road[n] === 0) { st.roadCost[n] = P.roadBaseCost * P.offRoadPenalty; continue; }
const x = st.flow[n] / P.roadCap;
let raw;
if (tvd) {
const xf = clamp(x, 0, tvdMax) / tvdMax * (tvd.length - 1);
const i0 = Math.floor(xf), i1 = Math.min(tvd.length - 1, i0 + 1);
raw = tvd[i0] + (xf - i0) * (tvd[i1] - tvd[i0]); // IDM-grounded volume-delay
} else {
raw = 1 + P.congestAlpha * Math.pow(Math.max(0, x), P.congestBeta); // BPR: the assumed congestion ratio
}
if (raw > maxCong) maxCong = raw;
// the cost that ENTERS the accessibility solve is SATURATED: congestion makes central travel dearer
// and pushes activity outward (dispersal), but a jam raises commute cost by a bounded factor rather
// than annihilating access -- so the core is discouraged, not destroyed (no runaway collapse).
const effRatio = 1 + P.congCostCap * (1 - Math.exp(-(raw - 1) / P.congCostCap)); // -> saturates at ~1+congCostCap
st.roadCost[n] = P.roadBaseCost * effRatio;
}
st.maxCongestion = maxCong;
// === 4. LAND VALUE (bid-rent) + AGGLOMERATION ================================================
// the target land value of each parcel is the highest bid across uses. Households bid on access to
// jobs; firms (commercial) bid HARDEST on access to jobs+customers AND on nearby developed density
// (AGGLOMERATION -- being near other activity is itself valuable), which is the positive feedback
// that concentrates a core. Industry bids modestly and repels nearby residential value (nuisance).
// Realised value relaxes toward the target (sticky) and is dragged down by local blight.
const R = P.agglomRadius;
for (let n = 0; n < N * N; n++) {
const i = n % N, j = (n / N) | 0;
// agglomeration term: developed density in the neighbourhood (self-reinforcing concentration)
let agg = 0, an = 0;
for (let dj = -R; dj <= R; dj++) for (let di = -R; di <= R; di++) {
const ii = i + di, jj = j + dj; if (!inB(ii, jj)) continue;
const m = idx(ii, jj); const w = 1 / (1 + di * di + dj * dj);
if (st.use[m] !== UND) { agg += st.density[m] * w; an += w; }
}
agg = an > 0 ? agg / an : 0;
// industry nuisance: how much industry sits nearby (drags residential bids down)
let nui = 0;
for (let dj = -1; dj <= 1; dj++) for (let di = -1; di <= 1; di++) {
const ii = i + di, jj = j + dj; if (!inB(ii, jj)) continue; if (st.use[idx(ii, jj)] === IND) nui += 1;
}
nui *= P.industryNuisance / 9;
// the three BID-RENT curves. What makes a city SORT into rings (von Thunen / Alonso) is not the
// level of each bid but its SLOPE against access: the use that bids most STEEPLY on centrality wins
// the centre, and flatter bidders are pushed outward to where they still win. So:
// - COMMERCIAL bids STEEPLY (access^1.5) and gains from AGGLOMERATION -- it will pay almost anything
// to be at the peak and nothing at the edge, so it wins a tight core and the agglomeration
// feedback makes that core self-reinforce.
// - RESIDENTIAL bids MODERATELY and broadly on access to jobs -- outbid at the very peak, it wins
// the wide RING around the core.
// - INDUSTRY bids SHALLOWLY, cares more about labour access (accessPop) than the prime centre, and
// tolerates cheap peripheral land -- so it takes the outer edge, away from homes it would blight.
const a = st.access[n], ap = st.accessPop[n];
// each bid is lifted by MARKET DEMAND PRESSURE (tight market -> higher rents -> development pushes
// outward). Commercial + industrial ride job demand, residential rides household demand.
const bidCom = (P.valAccessCom * Math.pow(a, 1.5) * (1 + P.agglomGain * agg / P.maxDensity) + 0.25 * ap) * st.demandJob;
const bidRes = (P.valAccessRes * a) * st.demandPop - nui;
const bidInd = (P.valAccessInd * (0.5 * a + 0.8 * ap) + P.indEdgeBonus) * st.demandJob;
// winning use = highest bidder
let bid = bidRes, wUse = RES;
if (bidCom > bid) { bid = bidCom; wUse = COM; }
if (bidInd > bid) { bid = bidInd; wUse = IND; }
const target = Math.max(0, bid);
st.landValue[n] += (target - st.landValue[n]) * P.valDecay;
// stash the winning use as a hint for development (only applied when (re)building)
st._use[n] = wUse;
}
// === 5. DEVELOPMENT + DENSITY ================================================================
// an undeveloped parcel whose land value clears the build threshold breaks ground (stochastically --
// construction is lumpy and takes time). A developed parcel intensifies: density climbs toward what
// its value supports, at a finite construction rate. Blight drags realised value, so a blighted
// parcel's supported density falls (it can't attract investment).
for (let n = 0; n < N * N; n++) {
const i = n % N, j = (n / N) | 0;
const valEff = st.landValue[n] * (1 - P.blightValueDrag * st.blight[n]);
if (st.use[n] === UND) {
if (valEff > P.buildThresh && rng() < P.buildRate * (0.5 + valEff)) {
// break ground as the use that bids highest here
st.use[n] = st._use[n] || RES;
st.density[n] = 0.6;
st.age[n] = 0;
st.evs.push({ kind: 'develop', i, j, use: st.use[n], value: +valEff.toFixed(3) });
}
} else {
const targetDens = clamp(valEff * P.densSupport, 0, P.maxDensity);
st.density[n] += (targetDens - st.density[n]) * P.densRate;
if (st.density[n] < 0.05) st.density[n] = 0.05;
}
}
// === 6. AGING / BLIGHT / RENEWAL =============================================================
// buildings age. An old or under-valued building accrues BLIGHT; blight spreads to neighbours
// (contagion); a badly blighted parcel is abandoned (emptied). When land value rises back above
// what stands there (worth more cleared than kept) AND the parcel is ripe, it is REDEVELOPED --
// cleared and rebuilt fresh + denser. This is the filtering/renewal cycle that churns a city.
const nb = st._nb; // neighbour-blight accumulation
for (let n = 0; n < N * N; n++) nb[n] = 0;
for (let n = 0; n < N * N; n++) {
if (st.use[n] === UND) continue;
st.age[n] += yrs * P.ageRate;
// under-valuation: the building is denser than the (blight-dragged) value now supports, or simply old
const valEff = st.landValue[n] * (1 - P.blightValueDrag * st.blight[n]);
const supported = clamp(valEff * P.densSupport, 0, P.maxDensity);
const overbuilt = st.density[n] > supported + 0.5;
const old = st.age[n] > P.blightAge;
if (overbuilt || old) {
st.blight[n] = clamp(st.blight[n] + P.blightRate * (0.5 + (overbuilt ? 0.5 : 0) + (old ? 0.4 : 0)), 0, 1);
} else {
st.blight[n] = clamp(st.blight[n] - P.blightRate * 0.7, 0, 1); // upkeep/reinvestment heals mild blight
}
}
// blight contagion: a blighted parcel bleeds decay onto its neighbours
for (let n = 0; n < N * N; n++) {
if (st.blight[n] < 0.2) continue; const i = n % N, j = (n / N) | 0;
for (let dj = -1; dj <= 1; dj++) for (let di = -1; di <= 1; di++) {
if (di === 0 && dj === 0) continue; const ii = i + di, jj = j + dj; if (!inB(ii, jj)) continue;
if (st.use[idx(ii, jj)] !== UND) nb[idx(ii, jj)] += st.blight[n];
}
}
for (let n = 0; n < N * N; n++) {
if (nb[n] > 0 && st.use[n] !== UND) st.blight[n] = clamp(st.blight[n] + P.blightSpread * (nb[n] / 8), 0, 1);
// abandonment + RENEWAL
if (st.use[n] === UND) continue;
const valEff = st.landValue[n] * (1 - P.blightValueDrag * st.blight[n]);
const standing = st.density[n] / P.maxDensity; // "what stands there" as a fraction
// renewal: land is worth clearly more than the tired thing on it, and the thing is ripe -> rebuild
if (st.landValue[n] > standing * P.renewValueRatio && st.blight[n] > P.renewMinBlight && st.landValue[n] > P.buildThresh) {
const i = n % N, j = (n / N) | 0;
st.use[n] = st._use[n] || st.use[n]; // rebuild as today's highest bidder
st.density[n] = clamp(st.density[n] * 0.5 + 1.0, 0.6, P.maxDensity); // clear + start fresh, denser floor
st.age[n] = 0; st.blight[n] = 0;
st.evs.push({ kind: 'renew', i, j, use: st.use[n], value: +st.landValue[n].toFixed(3) });
} else if (st.blight[n] > P.abandonBlight) {
const i = n % N, j = (n / N) | 0;
st.use[n] = UND; st.density[n] = 0; st.pop[n] = 0; st.jobs[n] = 0; st.age[n] = 0;
st.evs.push({ kind: 'abandon', i, j });
}
}
// === 7. POPULATION + JOBS LOCATE by utility ==================================================
// households and firms already in the region settle into the parcels that maximise their utility,
// subject to each parcel's housing/floor CAPACITY (density * capacityPerDens). Household utility ~
// access to jobs minus rent (land value) minus a crowding cost; firm utility ~ access to customers +
// agglomeration. We do a soft assignment: fill parcels in proportion to utility*free-capacity.
// Population/jobs already housed decay slightly when a parcel blights/empties (people leave decay).
// -- residents
{
let wsum = 0; const wbuf = st._wt; wbuf.fill(0);
for (let n = 0; n < N * N; n++) {
if (st.use[n] !== RES) continue;
const cap = capacityAt(st, n); const free = Math.max(0, cap - st.pop[n]);
// household desirability: mostly access to jobs, plus a base willingness to live in a built,
// un-blighted parcel at all (so a home near work is preferred but the city is not a ghost town
// waiting on perfect access), gently penalised by rent as a tie-breaker between good parcels.
const util = (P.resBase + st.access[n] * 1.2 - st.landValue[n] * 0.15) * (1 - 0.8 * st.blight[n]);
const w = Math.max(0, util) * free; wbuf[n] = w; wsum += w;
}
st.totalPop = 0;
// total households to house = the smaller of the regional pool and the housing capacity that
// exists. Distribute that total across RES parcels IN PROPORTION to utility*free-capacity, and
// relax each parcel toward its target share (moving IN takes time -- people do not teleport). This
// is a proper allocation of a fixed demand, not an unbounded per-step accretion.
if (wsum > 1e-9) {
let totCap = 0; for (let n = 0; n < N * N; n++) if (st.use[n] === RES) totCap += capacityAt(st, n);
const want = Math.min(st.poolPop, totCap);
for (let n = 0; n < N * N; n++) {
if (st.use[n] !== RES) { continue; }
const cap = capacityAt(st, n);
const target = clamp(want * (wbuf[n] / wsum), 0, cap);
// asymmetric migration: households move IN readily but are STICKY to leave (moving is costly), so
// a transient jam or dip does not evacuate a district overnight -- damps the boom-bust into a
// gentle cycle while preserving the real congestion-driven dispersal.
const rate = target > st.pop[n] ? P.moveRate : P.moveRate * P.leaveStick;
st.pop[n] += (target - st.pop[n]) * rate;
}
}
// shed residents where capacity shrank or blight rose
for (let n = 0; n < N * N; n++) {
if (st.use[n] !== RES) { st.pop[n] = 0; continue; }
st.pop[n] = clamp(st.pop[n], 0, capacityAt(st, n));
st.totalPop += st.pop[n];
}
}
// -- jobs (commercial + industrial)
{
let wsum = 0; const wbuf = st._wt; wbuf.fill(0);
for (let n = 0; n < N * N; n++) {
if (st.use[n] !== COM && st.use[n] !== IND) continue;
const cap = capacityAt(st, n); const free = Math.max(0, cap - st.jobs[n]);
// firm desirability: access to customers/labour, plus a base willingness to operate at all, minus
// rent as a tie-breaker -- so firms cluster toward the centre but the whole commercial stock fills.
const util = Math.max(0, P.jobBase + st.accessPop[n] * 1.0 + st.access[n] * 0.4 - st.landValue[n] * 0.15) * (1 - 0.8 * st.blight[n]);
const w = util * free; wbuf[n] = w; wsum += w;
}
st.totalJobs = 0;
if (wsum > 1e-9) {
let totCap = 0; for (let n = 0; n < N * N; n++) if (st.use[n] === COM || st.use[n] === IND) totCap += capacityAt(st, n);
const want = Math.min(st.poolJobs, totCap);
for (let n = 0; n < N * N; n++) {
if (st.use[n] !== COM && st.use[n] !== IND) continue;
const cap = capacityAt(st, n);
const target = clamp(want * (wbuf[n] / wsum), 0, cap);
st.jobs[n] += (target - st.jobs[n]) * P.moveRate;
}
}
for (let n = 0; n < N * N; n++) {
if (st.use[n] !== COM && st.use[n] !== IND) { st.jobs[n] = 0; continue; }
st.jobs[n] = clamp(st.jobs[n], 0, capacityAt(st, n));
st.totalJobs += st.jobs[n];
}
}
// === 8. OBSERVABLES ==========================================================================
let peak = 0, developed = 0, blk = 0, blkN = 0, commuteW = 0, commuteAcc = 0;
for (let n = 0; n < N * N; n++) {
if (st.density[n] > peak) peak = st.density[n];
if (st.use[n] !== UND) { developed++; blk += st.blight[n]; blkN++; }
// mean commute proxy: population-weighted job-ACCESS (higher access = shorter effective commute).
// We report its inverse-ish as an accessibility, and separately the population-weighted access.
if (st.pop[n] > 0) { commuteW += st.pop[n]; commuteAcc += st.pop[n] * st.access[n]; }
}
st.peakDensity = peak; st.developed = developed;
st.meanBlight = blkN > 0 ? blk / blkN : 0;
st.meanAccess = commuteW > 0 ? commuteAcc / commuteW : 0; // how good, on average, a resident's job access is
st.totalPop = st.totalPop || 0;
return st;
},
fields(st) {
return {
N: st.N,
use: st.use, density: st.density, landValue: st.landValue,
access: st.access, accessPop: st.accessPop, blight: st.blight,
road: st.road, flow: st.flow, roadCost: st.roadCost,
pop: st.pop, jobs: st.jobs,
totalPop: st.totalPop, totalJobs: st.totalJobs,
peakDensity: st.peakDensity, developed: st.developed,
meanBlight: st.meanBlight, maxCongestion: st.maxCongestion, meanAccess: st.meanAccess,
poolPop: st.poolPop, poolJobs: st.poolJobs
};
},
events(st) { return st.evs; }
};
}
return { makeCity };
});