@JvmStatic
    fun main(args: Array<String>) {
        val keepOpen = args.contains("--keep")
        Options.enableVR = false
        Options.doVulkanSDK()

        outputDir.mkdirs()
        val window = Window("Field · Fabricate simulation chains")
        val renderer = Renderer(window)

        // Same reason as FabricatePassChainDemo: StandardRenderer.build allocates targets through a
        // bare `RenderToTexture { … }`, whose renderer defaults to `fromContext()`, so it needs a box
        // context. Without one the failure is a NullPointerException from inside the RTT constructor.
        val hostBox = Box("fabricate-simulation-chain-demo").apply { name = "fabricate-simulation-chain-demo" }
        hostBox.setPropertyRaw("renderer", renderer, doResync = false)
        hostBox.setPropertyRaw("window", window, doResync = false)

        try {
            Box.withBoxContextNow(hostBox) { run(renderer, window) }
            if (keepOpen) {
                println("-- --keep: leaving the window open, close it to exit --")
                while (!window.shouldClose()) pump(renderer, window, 1)
            }
        } catch (t: Throwable) {
            println("FAIL: ${t.message}")
            t.printStackTrace()
            failures++
        } finally {
            runCatching { renderer.cleanup() }
            runCatching { window.cleanup() }
        }

        if (failures == 0) {
            println("simulation chain tests PASSED — figures in ${outputDir.absolutePath}")
        } else {
            println("simulation chain tests FAILED ($failures)")
        }
        exitProcess(if (failures == 0) 0 else 1)
    }

    private fun run(renderer: Renderer, window: Window) {
        val bundle = StandardRenderer.build(
            renderer,
            StandardRenderer.Config(
                renderWidth = DIM,
                renderHeight = DIM,
                backgroundColor = Vector3f(0.05f, 0.06f, 0.08f)
            )
        )
        bundle.applyScenePlan()
        bundle.present(renderer = renderer)
        bundle.attachTo(renderer)
        pump(renderer, window, WARMUP_FRAMES)

        // Two fluids, identical but for the vorticity term, so "confinement does something" is an
        // A/B in one run rather than a number somebody picked.
        val fluid = bundle.addPassChain(
            FabricatePassLibrary.stableFluid(
                prefix = "fluid",
                resolution = FLUID_DIM,
                iterations = 20,
                dyeDissipation = DYE_DISSIPATION,
                vorticity = 24f
            ),
            renderer = renderer,
            printReport = true
        )
        val smooth = bundle.addPassChain(
            FabricatePassLibrary.stableFluid(
                prefix = "smooth",
                resolution = FLUID_DIM,
                iterations = 20,
                dyeDissipation = DYE_DISSIPATION,
                vorticity = 0f
            ),
            renderer = renderer
        )

        // The distance chains. Both synthesise their own mask, so neither depends on what the fluid
        // chains left in the post stack's output.
        val dot = bundle.addPassChain(dotMaskChain(), renderer = renderer, printReport = true)
        val disc = bundle.addPassChain(discMaskChain(), renderer = renderer)

        bundle.present(renderer = renderer)
        bundle.attachTo(renderer)

        guard("both fluids allocated a ping-pong pair per state field, and nothing else shares them") {
            listOf("fluid" to fluid, "smooth" to smooth).forEach { (name, chain) ->
                require(chain.stateNames() == listOf("$name-velocity", "$name-dye")) {
                    "$name state fields are ${chain.stateNames()}"
                }
                chain.stateNames().forEach { field ->
                    val state = chain.stateOf(field) ?: error("$field is not a ping-pong target")
                    require(state.frameCount == FabricatePassChain.PERSISTENT_IMAGE_COUNT) {
                        "$field has ${state.frameCount} frames"
                    }
                    require(state.width == FLUID_DIM) { "$field is ${state.width} wide" }
                }
                require(chain.plan.slots.count { it.persistent } == 2) {
                    "$name planned ${chain.plan.slots.count { it.persistent }} persistent slots"
                }
            }
            "two chains x two ${FLUID_DIM}² ping-pong fields, ${fluid.plan.passCount} passes each"
        }

        // --- paint, then stop painting ---

        splat(fluid, on = true)
        splat(smooth, on = true)
        pump(renderer, window, SPLAT_FRAMES)
        val painted = capture(renderer, fluid.output, "01-fluid-painted")

        splat(fluid, on = false)
        splat(smooth, on = false)
        pump(renderer, window, COAST_FRAMES)
        val early = capture(renderer, fluid.output, "02-fluid-coasting-early")
        val earlySmooth = capture(renderer, smooth.output, "03-fluid-coasting-early-no-vorticity")
        val velocity = capture(renderer, fluid.stateOf("fluid-velocity")!!, "04-fluid-velocity-field")

        pump(renderer, window, DRIFT_FRAMES)
        val late = capture(renderer, fluid.output, "05-fluid-coasting-late")
        val lateSmooth = capture(renderer, smooth.output, "06-fluid-coasting-late-no-vorticity")

        guard("the splat reached the dye field at all") {
            require(painted.meanLuma > 0.002) {
                "the dye field is effectively black (mean luma ${"%.5f".format(painted.meanLuma)}) — " +
                        "either the splat push constants never arrived or the persistent target is unbound"
            }
            require(painted.distinctLumaBuckets >= 6) {
                "the dye field has ${painted.distinctLumaBuckets} distinct levels — a uniform image is " +
                        "what a pass that never sampled its input renders"
            }
            // Checked before anything below draws a conclusion from a sum. A field over 1.0 reads
            // back clamped, and a clamped field's decay is invisible however fast it is decaying —
            // see SPLAT_MASS for the run where that cost an hour.
            require(painted.maxLuma < SATURATION_CEILING) {
                "the painted dye peaks at ${"%.4f".format(painted.maxLuma)}, at or over the readback's " +
                        "clamp — every energy measurement below would be measuring the *area* of the " +
                        "saturated region rather than the field. Lower SPLAT_MASS."
            }
            "painted dye mean luma ${"%.4f".format(painted.meanLuma)} over " +
                    "${painted.distinctLumaBuckets} levels, peak ${"%.3f".format(painted.maxLuma)} (unclamped)"
        }

        guard("the velocity field is being written, not left at its clear colour") {
            require(velocity.distinctLumaBuckets >= 4) {
                "the velocity field has ${velocity.distinctLumaBuckets} distinct levels"
            }
            "velocity field ${velocity.size}, ${velocity.distinctLumaBuckets} levels"
        }

        guard("with the input switched off the dye keeps moving — the feedback loop actually advances") {
            // The assertion this whole file exists for. A chain whose ping-pong halves were
            // snapshotted at build time plans, allocates and renders perfectly, and renders the
            // *same frame* forever. Both captures below are taken with every splat push constant at
            // zero, so anything that changes between them came from the previous frame's state.
            val moved = hypot(late.centroidU - early.centroidU, late.centroidV - early.centroidV)
            require(moved > DRIFT_FLOOR) {
                "the dye centroid moved only ${"%.4f".format(moved)} UV over $DRIFT_FRAMES frames with " +
                        "no input — under $DRIFT_FLOOR means the simulation is not advancing, which is " +
                        "what a frozen previous-frame binding looks like"
            }
            // And in the direction the splat pushed, which rules out the centroid drifting because
            // the dye is fading unevenly rather than being transported.
            require(late.centroidU > early.centroidU) {
                "the dye drifted to ${"%.4f".format(late.centroidU)} from " +
                        "${"%.4f".format(early.centroidU)} — the splat pushed along +u, so transport " +
                        "should carry it that way"
            }
            "centroid ${"%.4f".format(early.centroidU)} -> ${"%.4f".format(late.centroidU)} along u, " +
                    "${"%.4f".format(moved)} UV travelled over $DRIFT_FRAMES frames of no input"
        }

        guard("dissipation is compounding per frame, which only a persistent field can do") {
            // A field re-rendered from scratch every frame would land on the same total every time.
            // Geometric decay is evidence of accumulation, and the floor is set well under the
            // arithmetic so the assertion does not depend on the exact frame count.
            val ratio = late.totalEnergy / early.totalEnergy
            val expected = Math.pow(DYE_DISSIPATION.toDouble(), DRIFT_FRAMES.toDouble())
            require(ratio < 1.0) {
                "total dye rose from ${"%.1f".format(early.totalEnergy)} to ${"%.1f".format(late.totalEnergy)} " +
                        "with no input — with dissipation at $DYE_DISSIPATION it can only fall"
            }
            require(ratio > expected * 0.25) {
                "total dye fell to ${"%.3f".format(ratio)} of itself against an expected " +
                        "${"%.3f".format(expected)} — far more loss than dissipation accounts for, which " +
                        "is what advecting by a field with residual divergence does"
            }
            "total dye fell to ${"%.3f".format(ratio)} of itself over $DRIFT_FRAMES frames " +
                    "(dissipation alone predicts ${"%.3f".format(expected)})"
        }

        guard("vorticity confinement changes the fluid, measured against the same splat without it") {
            require(earlySmooth.meanLuma > 0.002) { "the no-vorticity fluid is black" }
            // Confinement puts back the curl a first-order advection dissipates, so the confined
            // fluid holds more structure at the top of the spectrum. Compared at the same instant
            // under the same splat schedule, so the only difference is the term.
            val confined = late.highFrequencyEnergy
            val plain = lateSmooth.highFrequencyEnergy
            require(abs(confined - plain) > plain * 0.05) {
                "the confined fluid's high-frequency energy (${"%.5f".format(confined)}) is within 5% of " +
                        "the unconfined one's (${"%.5f".format(plain)}) — the vorticity term is not " +
                        "reaching the shader, or the curl pass is bound to the wrong target"
            }
            "confined hf ${"%.5f".format(confined)} vs unconfined ${"%.5f".format(plain)} " +
                    "(${"%.2f".format(confined / plain)}x)"
        }

        // --- the distance field ---

        val dotSeed = capture(renderer, dot.targetOf("dot-seed")!!, "07-jfa-seed")
        val dotDistance = capture(renderer, dot.targetOf("dot-distance")!!, "08-jfa-distance")
        val dotBands = capture(renderer, dot.output, "09-jfa-bands")
        val discBands = capture(renderer, disc.output, "10-jfa-signed-bands")

        guard("the flood ran the full stride sequence and ended where the plan said") {
            val strides = dot.plan.passes.mapNotNull { it.spec.push["stride"] as? Int }
            require(strides.take(8) == listOf(128, 64, 32, 16, 8, 4, 2, 1)) { "strides were $strides" }
            require(dotDistance.width == JFA_DIM && dotDistance.height == JFA_DIM) {
                "the distance field came back ${dotDistance.size}, not ${JFA_DIM}²"
            }
            require(dotSeed.maxLuma > 0.05) { "nothing was seeded — the mask never crossed the threshold" }
            "${strides.size} flood passes over ${dotDistance.size}, seeded"
        }

        guard("the distance field agrees with the exact distance to the one seed it has") {
            // Ground truth, not a previous render and not a plausible band pattern. With a single
            // seed the answer is `length(p - seed) * JFA_SCALE`, so a flood one rung short, a
            // filtered read, or a stride sequence off by a factor of two all show up here as a
            // systematic error — and in nothing a person looking at the bands would notice.
            // Exactly one seed, checked rather than assumed: the ground truth below is the distance
            // to a single point, and a second seeded texel would make it quietly wrong in a wedge of
            // the field rather than everywhere.
            val zeros = dotDistance.countWhereRedCodeIs(0)
            require(zeros == 1) {
                "$zeros texels read as distance zero, not 1 — the mask seeded more than one texel, so " +
                        "nearest-seed distance is no longer distance to a point and the truth below is not " +
                        "the truth. Shrink SEED_RADIUS."
            }
            val seedX = SEED_TEXEL + 0.5
            val seedY = SEED_TEXEL + 0.5
            var worst = 0.0
            var worstAt = ""
            var samples = 0
            // Every eighth texel, skipping the outer ring, so the sample set spans the whole field
            // including the corners the largest stride has to reach.
            for (y in 8 until JFA_DIM - 8 step 8) {
                for (x in 8 until JFA_DIM - 8 step 8) {
                    val exact = hypot(x + 0.5 - seedX, y + 0.5 - seedY) * JFA_SCALE
                    // Both sides in 8-bit codes: the exact distance put through the same transfer
                    // function the readback applied, so the residual is quantisation and nothing else.
                    val exactCode = linearToSrgb(exact) * 255.0
                    val measuredCode = dotDistance.redCode(x, y)
                    val error = abs(measuredCode - exactCode)
                    samples++
                    if (error > worst) {
                        worst = error
                        worstAt = "($x, $y): measured code ${"%.0f".format(measuredCode)}, " +
                                "exact code ${"%.1f".format(exactCode)} (${"%.4f".format(exact)} linear)"
                    }
                }
            }
            require(samples > 200) { "only $samples usable samples" }
            require(worst <= DISTANCE_CODE_TOLERANCE) {
                "worst error ${"%.1f".format(worst)} codes at $worstAt, over a tolerance of " +
                        "$DISTANCE_CODE_TOLERANCE — a jump flood is near-exact for a single seed, so " +
                        "anything above the readback's own quantisation is a missing rung, a filtered " +
                        "read, or the wrong stride sequence"
            }
            "worst error ${"%.1f".format(worst)} of 255 codes over $samples samples spanning the whole " +
                    "field (tolerance $DISTANCE_CODE_TOLERANCE), worst at $worstAt"
        }

        guard("the bands are evenly spaced, which is the visual check that it is a distance") {
            require(dotBands.distinctLumaBuckets >= 8) {
                "the bands have ${dotBands.distinctLumaBuckets} levels — a distance field with no " +
                        "structure is one that never got a seed"
            }
            "bands ${dotBands.size}, ${dotBands.distinctLumaBuckets} levels"
        }

        guard("the signed field is negative inside the disc and positive outside it") {
            // The sign cannot be read out of an 8-bit PNG, so it is read out of the *colour* the
            // bands pass assigns each side. A signed chain that lost its mask binding produces an
            // everywhere-positive field, which is a perfectly plausible unsigned distance field and
            // is why this is asserted rather than eyeballed.
            val inside = discBands.meanChannelsIn(0.45, 0.55, 0.45, 0.55)
            val outside = discBands.meanChannelsIn(0.02, 0.12, 0.02, 0.12)
            // distanceBands' defaults: blue-dominant inside, orange (red-dominant) outside.
            require(inside.third > inside.first) {
                "the centre of the disc came back red-dominant (r ${"%.3f".format(inside.first)}, " +
                        "b ${"%.3f".format(inside.third)}) — it is on the *outside* of the field's sign"
            }
            require(outside.first > outside.third) {
                "the corner came back blue-dominant (r ${"%.3f".format(outside.first)}, " +
                        "b ${"%.3f".format(outside.third)}) — the field has no negative side, so either " +
                        "the mask binding or the boundary seeding is gone"
            }
            "inside rgb ${inside.fmt()} (blue side), outside rgb ${outside.fmt()} (orange side)"
        }

        guard("a live push retunes the fluid with no reallocation") {
            val before = fluid.plan.targetCount
            require(fluid.setPush("fluid-dye", "dissipation", 0.5f)) { "dissipation push was not accepted" }
            val was = capture(renderer, fluid.output, "11-fluid-before-fast-decay")
            pump(renderer, window, 12)
            val now = capture(renderer, fluid.output, "12-fluid-after-fast-decay")
            require(fluid.plan.targetCount == before) { "retuning changed the target count" }
            // 0.5^12 is about 0.0002, so twelve frames at half-per-frame is effectively a wipe.
            require(now.totalEnergy < was.totalEnergy * 0.05) {
                "dropping dissipation to 0.5 for 12 frames left ${"%.3f".format(now.totalEnergy / was.totalEnergy)} " +
                        "of the dye — the push constant is not reaching the shader"
            }
            "dissipation 0.985 -> 0.5 wiped the dye to " +
                    "${"%.5f".format(now.totalEnergy / was.totalEnergy)} of itself in 12 frames, no reallocation"
        }
    }

    private class Shot(val name: String, file: File) {
        private val image = ImageIO.read(file) ?: error("could not re-read $name from ${file.absolutePath}")
        val width = image.width
        val height = image.height
        val size get() = "${width}x$height"

        private val rgb = IntArray(width * height) { image.getRGB(it % width, it / width) }

        private val luma = DoubleArray(width * height) { index ->
            val c = rgb[index]
            0.2126 * (c shr 16 and 0xff) / 255.0 + 0.7152 * (c shr 8 and 0xff) / 255.0 + 0.0722 * (c and 0xff) / 255.0
        }

        /**
         * The same luminance in **linear light**.
         *
         * `TexturePreviewCache` blits a float target to `VK_FORMAT_R8G8B8A8_SRGB` on the way out, so
         * the PNG's bytes are display-encoded and *not* the numbers the shader wrote. Two of this
         * file's assertions failed on that before anything was wrong with a shader:
         *
         * - the total dye **rose** 1.7% over forty frames of pure dissipation, because a blob
         *   spreading out moves energy from high sRGB codes, where the curve compresses, to low ones
         *   where it expands — the encoded sum went up while the energy went down;
         * - the distance field read 0.3137 where the exact answer was 0.0804, a 290% error, which is
         *   `srgbToLinear(80/255) = 0.0802` and therefore a 0.25% error.
         *
         * Anything that is a *sum* or is compared against an *absolute* expectation has to be
         * decoded first. Ratios of the same statistic between two captures do not, which is why
         * [highFrequencyEnergy] is left encoded.
         */
        private val linearLuma = DoubleArray(width * height) { index ->
            val c = rgb[index]
            0.2126 * srgbToLinear((c shr 16 and 0xff) / 255.0) +
                    0.7152 * srgbToLinear((c shr 8 and 0xff) / 255.0) +
                    0.0722 * srgbToLinear((c and 0xff) / 255.0)
        }

        val meanLuma = luma.average()
        val maxLuma = luma.max()

        /** Total energy, in linear light — the statistic dissipation and mass conservation act on. */
        val totalEnergy = linearLuma.sum()
        val distinctLumaBuckets = luma.map { (it * 64).toInt() }.distinct().size

        /**
         * The red channel at a texel in linear light — the channel a distance pass writes its
         * distance into, decoded back out of the sRGB the readback encoded it into.
         */
        fun red(x: Int, y: Int): Double = srgbToLinear(redCode(x, y) / 255.0)

        /**
         * The red channel as the raw 0..255 code the PNG holds.
         *
         * The distance comparison is made here rather than in linear light, because sRGB is steep in
         * the dark and amplifies one code step into several percent of the linear value — see
         * [DISTANCE_CODE_TOLERANCE].
         */
        fun redCode(x: Int, y: Int): Double = (rgb[y * width + x] shr 16 and 0xff).toDouble()

        /** How many texels hold exactly [code] in red — used to count the seeds a mask produced. */
        fun countWhereRedCodeIs(code: Int): Int = rgb.count { (it shr 16 and 0xff) == code }

        /**
         * Luminance-weighted centroid, in UV.
         *
         * The statistic that distinguishes dye being *transported* from dye merely fading: a fading
         * blob keeps its centroid, a advected one moves it. Weighted by luminance rather than by a
         * threshold so a blob that is spreading as it moves does not jump as texels cross the line.
         */
        val centroidU: Double
        val centroidV: Double

        init {
            var wu = 0.0
            var wv = 0.0
            var total = 0.0
            for (y in 0 until height) {
                for (x in 0 until width) {
                    val w = luma[y * width + x]
                    wu += w * (x + 0.5) / width
                    wv += w * (y + 0.5) / height
                    total += w
                }
            }
            val safe = if (total > 1e-9) total else 1.0
            centroidU = wu / safe
            centroidV = wv / safe
        }

        /**
         * Mean absolute discrete Laplacian of luminance — energy at the top of the spectrum.
         *
         * Deliberately not a gradient magnitude: a gradient survives a blur reasonably well, which
         * would make the vorticity A/B far weaker than it looks.
         */
        val highFrequencyEnergy: Double by lazy {
            var total = 0.0
            var count = 0
            for (y in 1 until height - 1) {
                for (x in 1 until width - 1) {
                    val here = luma[y * width + x]
                    val lap = luma[y * width + x - 1] + luma[y * width + x + 1] +
                            luma[(y - 1) * width + x] + luma[(y + 1) * width + x] - 4 * here
                    total += abs(lap)
                    count++
                }
            }
            total / count
        }

        /** Mean (r, g, b) over a UV rectangle — how the sign of a signed field is read back. */
        fun meanChannelsIn(u0: Double, u1: Double, v0: Double, v1: Double): Triple<Double, Double, Double> {
            var r = 0.0
            var g = 0.0
            var b = 0.0
            var count = 0
            for (y in (v0 * height).toInt() until (v1 * height).toInt().coerceAtMost(height)) {
                for (x in (u0 * width).toInt() until (u1 * width).toInt().coerceAtMost(width)) {
                    val c = rgb[y * width + x]
                    r += (c shr 16 and 0xff) / 255.0
                    g += (c shr 8 and 0xff) / 255.0
                    b += (c and 0xff) / 255.0
                    count++
                }
            }
            require(count > 0) { "$name: empty region ($u0..$u1, $v0..$v1)" }
            return Triple(r / count, g / count, b / count)
        }
    }