Skip to content

SolidInverse.kt

src/jvmMain/kotlin/afield/lang/solidsl/SolidInverse.kt · 634 lines

@file:afield.runtime.completion.CompletionCategory("Graphics")

package afield.lang.solidsl

import afield.math.autodiff.AD
import afield.math.autodiff.ADTape
import afield.math.autodiff.OptimizeMethod
import afield.math.autodiff.OptimizeOptions
import afield.math.autodiff.OptimizeRun
import afield.math.autodiff.clamp
import afield.math.autodiff.minimizeWithGradient
import afield.math.autodiff.valueAndGrad
import afield.math.autodiff.ADGradient
import kotlin.math.abs
import kotlin.math.max
import kotlin.math.sqrt
import kotlin.random.Random

The payoff of lowering (c): inverse design. Declare what a solid should look like and let a gradient find the numbers that make it so.

val truth  = (sphere(0.80) + cylinder(1.60, 0.30)).let { it }
val target = truth.sampleField(count = 3000, seed = 7)

val template = truth.parameterize {
    param("r") { it.node == "Sphere" && it.field == "radius" }
    param("h") { it.node == "Cylinder" && it.field == "height" }
}.orThrow()

val fit = template.fitParameters(target, start = doubleArrayOf(0.35, 2.60))
println(fit.report())
fit.fitted.march(96)      // an ordinary SolidExpr again

This is afield.lang.rectsl.Track.fitParameters in three dimensions, and it is the reason SolidGradient.kt exists. RectSLFn's KDoc: "a track that lowers to an Fn is a track you can take the derivative of, and therefore a track you can solve: declare where a rect should be at some instant and let a gradient find the numbers that put it there."

The target is a sampled field, not a mesh

A SolidSamples is a bag of points with a signed distance at each. Three things recommend it over fitting against a mesh:

  • It has a gradient. A mesh-to-mesh distance goes through a nearest-triangle search whose result is a discontinuous function of the parameters — the nearest triangle changes — and marching cubes has to be re-run at every step. A field residual is pred − target at a fixed point, and its derivative is exactly what SolidTemplate.adFieldWith hands over.
  • The target need not be a solid at all. sampleField builds one from another SolidExpr, which is what a recovery test wants; but a scan, a level set from somewhere else, or a hand- written list of "this point must be 2mm outside" constraints is the same type.
  • It is cheap and fixed. Sample once; every iteration reads the same array. A 3000-point target costs one march-free forward pass and one backward sweep per step.

The cost is that it fits the field and not the surface, and those are different objectives. Two shapes can have the same zero set and different fields — SolidExpr.Intersect's own SolidExpr.fieldNotes entry says it over-estimates distance outside a concave corner — so a field residual is a stricter target than "look the same". Strictly speaking that is a feature for recovery (there is one right answer and the loss knows it) and a hazard for design (you may not care about the field away from the surface). SolidFitOptions.band is the dial: clamp both sides into a narrow band around zero and the objective becomes almost purely about the surface.

What makes a parameter refuse to move

This is the failure mode of inverse design and it is worth naming before it happens, because it presents as a converged optimisation:

  • A hard SolidExpr.Union gives the whole derivative to whichever child is nearest. A child that is nowhere nearest — swallowed inside another — receives exactly zero from every sample, and its parameters have a zero gradient. Write the shape with smoothUnion and every child contributes within k of the seam.
  • A SolidExpr.Difference cutter entirely inside the body is the same story.
  • A SolidExpr.Empty subtree, or an Extrude sample past the end of the sweep, is a structural zero.

SolidFit.deadParameters reports exactly which parameters had a zero gradient at the start, so that reads as a diagnosis rather than as a mysteriously immediate convergence. And SolidExpr.gradientNotes says the same things about an expression before it is ever fitted.

@see SolidTemplate.checkGradient for the finite-difference second opinion on any single point @see SolidExpr.gradientAgreement for lowering (c)'s value measured against lowering (a)

The target

class SolidSamples

A bag of points with a signed distance at each — the thing a parametric solid is fitted to.

Four parallel DoubleArrays rather than a List<Pair<Vec3, Double>> because the inner loop of a fit walks all of them once per iteration and per line-search probe, and because Vec3 is Vector3f: Float, and mutable, which is the objection P2 was introduced to answer.

class SolidSamples(
    val xs: DoubleArray,
    val ys: DoubleArray,
    val zs: DoubleArray,
    val distances: DoubleArray,
) {

init

    init {
        require(xs.size == ys.size && ys.size == zs.size && zs.size == distances.size) {
            "a SolidSamples needs four arrays of one length, got ${xs.size}, ${ys.size}, ${zs.size}, ${distances.size}"
        }
        require(xs.isNotEmpty()) { "a SolidSamples needs at least one point" }
    }

val size: Int get() = xs.size

val insideCount: Int get() = distances.count { it < 0.0 }

How many of these points are inside the shape they were sampled from — a sanity number.

val rms: Double get() = sqrt(distances.sumOf { it * it } / size)

The root-mean-square distance in this target — the scale a loss should be read against.

A loss of 0.01 means nothing on its own; 0.01 against an rms of 0.5 is a 2% residual.

override fun toString(): String

    override fun toString(): String = "SolidSamples($size points, $insideCount inside, rms ${"%.4g".format(rms)})"
}

fun SolidExpr.sampleField

Sample this shape's field at count pseudo-random points in its own bounds, expanded by margin as a fraction of the longest axis.

Uniform in the box rather than concentrated near the surface: a near-surface sampling is a better surface objective but it has to know where the surface is, which for a target that is being recovered is circular. Uniform-in-a-box needs only the bounds, and SolidFitOptions.band recovers the near-surface behaviour by clamping the residual instead of by moving the points.

The sequence is a plain Random seeded by seed, so a fit is reproducible to the last bit and a test that depends on one is not depending on the wall clock.

@param bounds overrides the shape's own SolidExpr.bounds — required for an unbounded shape (one containing a SolidExpr.Gyroid), which has no window of its own.

fun SolidExpr.sampleField(
    count: Int = 2000,
    seed: Int = 0,
    bounds: afield.geometry.spatial.AxisAlignedBoundingBox3D? = null,
    margin: Double = 0.25,
): SolidSamples {
    require(count > 0) { "count must be positive, got $count" }
    val raw = bounds ?: this.bounds() ?: error(
        "cannot sample $this: it is unbounded (it contains a Gyroid, or a Union with one) so there " +
            "is no window to sample in — pass `bounds` explicitly, or intersect it with something bounded",
    )
    val s = raw.size
    val longest = maxOf(s.x, s.y, s.z).toDouble()
    val box = raw.expandedBy(longest * margin)
    val lo = box.min
    val hi = box.max
    val f = toField3()
    val rng = Random(seed)
    val xs = DoubleArray(count); val ys = DoubleArray(count)
    val zs = DoubleArray(count); val ds = DoubleArray(count)
    for (i in 0 until count) {
        val x = lo.x + rng.nextDouble() * (hi.x - lo.x)
        val y = lo.y + rng.nextDouble() * (hi.y - lo.y)
        val z = lo.z + rng.nextDouble() * (hi.z - lo.z)
        xs[i] = x; ys[i] = y; zs[i] = z
        ds[i] = f.at(xs[i], ys[i], zs[i])
    }
    return SolidSamples(xs, ys, zs, ds)
}

The measurement: lowering (c) against lowering (a)

data class SolidFieldAgreement

Lowering (c)'s value against lowering (a)'s, over one shape.

SolidAgreement's argument, applied to a third column: "Two implementations that share no code agreeing to five figures is a much better argument that either one is right than any number of assertions against remembered constants."

It is a separate measurement from crossCheck rather than a fourth set of fields on SolidAgreement, and the reason is that they measure different objects. crossCheck compares surfaces — volume, area, genus, Hausdorff deviation — and lowering (c) does not produce a surface; it produces a field, exactly as lowering (a) does. So the comparison is field-against-field, pointwise, and folding it into a surface report would have meant inventing a "volume" for something that has none. What the two measurements share is the argument, and this KDoc and SolidAgreement's cross-reference each other so a reader finds both.

The expected residual is not zero and not the faceting deficit either: it is Float. Lowering (a) reaches a caller through toScalarField, which is (Vec3) -> Float, so about seven decimal digits survive. maxRelativeError against 1e-6 is therefore the tight assertion, and anything larger is a real disagreement.

data class SolidFieldAgreement(
    val expr: SolidExpr,
    val samples: Int,
    /** Largest `|c − a|` over the samples, in model units. */
    val maxAbsoluteError: Double,
    /** Largest `|c − a| / max(1, |a|)` — relative, but floored so a near-zero value is not a blow-up. */
    val maxRelativeError: Double,
    /** Where the worst disagreement was, for a `println` that can be re-evaluated by hand. */
    val worstAt: Triple<Double, Double, Double>,
    /** How many samples produced a non-finite value on either side — must be 0. */
    val nonFinite: Int,
) {

fun report(): String

    fun report(): String =
        "field agreement over $samples samples: max abs ${"%.4e".format(maxAbsoluteError)}, " +
            "max rel ${"%.4e".format(maxRelativeError)} at " +
            "(${"%.4g".format(worstAt.first)}, ${"%.4g".format(worstAt.second)}, ${"%.4g".format(worstAt.third)})" +
            (if (nonFinite > 0) ", $nonFinite NON-FINITE" else "")
}

fun SolidExpr.gradientAgreement

Evaluate this shape both through toADField and through toScalarField at the same points, and report how far apart they land.

One tape for the whole sweep, and no backward pass: this measures the value half of the differentiable reading, which is the half that has an independent second opinion. The gradient half is checked against central differences instead, by SolidTemplate.checkGradient — there is no second implementation of the derivative to compare against, which is precisely why the numerical oracle is used there and not here.

fun SolidExpr.gradientAgreement(
    count: Int = 4000,
    seed: Int = 0,
    bounds: afield.geometry.spatial.AxisAlignedBoundingBox3D? = null,
    margin: Double = 0.25,
): SolidFieldAgreement {
    val target = sampleField(count, seed, bounds, margin)
    val tape = ADTape()
    val ad = toADField(SolidNumerics(tape))
    var maxAbs = 0.0
    var maxRel = 0.0
    var worst = Triple(0.0, 0.0, 0.0)
    var bad = 0
    val approx = toScalarField()
    for (i in 0 until target.size) {
        val x = target.xs[i]; val y = target.ys[i]; val z = target.zs[i]
        val a = approx(afield.math.linalg.Vec3(x.toFloat(), y.toFloat(), z.toFloat())).toDouble()
        val c = ad.at(tape.constant(x), tape.constant(y), tape.constant(z)).value
        if (!a.isFinite() || !c.isFinite()) { bad++; continue }
        val e = abs(c - a)
        val rel = e / max(1.0, abs(a))
        if (rel > maxRel) { maxRel = rel; worst = Triple(x, y, z) }
        if (e > maxAbs) maxAbs = e
    }
    return SolidFieldAgreement(this, target.size, maxAbs, maxRel, worst, bad)
}

The finite-difference second opinion

data class SolidGradientCheck

An analytic gradient beside a central-difference one, component by component.

The oracle is central rather than forward differences for two reasons, and the second is the one that matters. Central differences are O(h²) accurate instead of O(h), which buys three digits; and — as autodiff.kt's kink convention says explicitly — a central difference at a kink reports the average of the two one-sided slopes, which is exactly the 0.5/0.5 split the analytic gradient takes there. So the two agree at kinks as well as away from them, and a CSG seam does not have to be excluded from a census.

The verdict is per component, and there is a verdict about the verdict

usableComponents marks which entries the oracle is willing to speak about. A central difference is O(h²), so halving h should barely move it; if it moves a lot, the third derivative is enormous or the function is not differentiable at this point at all, and the numerical answer is not evidence about anything. That happens at a sphere's own centre, on a polygon's outline, and within h of any kink — a real property of distance fields, not a weakness of the analytic side.

Per component rather than per point because the alternative loses most of the census: a twisted star extrusion has twenty-six parameters, and one of them being near its own kink is not a reason to discard the verdict on the other twenty-five.

And a census that counts the unusable ones is the only kind worth running. A check with nothing usable has no verdict, and a suite of those passes silently — which is exactly the failure RectSLFnBackendTest documents guarding against. usableCount against size is the number to print and then floor.

data class SolidGradientCheck(
    val analytic: DoubleArray,
    val numerical: DoubleArray,
    /** The same central difference at half the step — what [usableComponents] is decided from. */
    val numericalHalfStep: DoubleArray,
    val step: Double,
) {

val size: Int get() = analytic.size

How many components there are.

val usableComponents: BooleanArray

Which components the numerical oracle is prepared to have an opinion about.

A statement about the oracle only. It says nothing about whether the analytic side produced a number, and keeping those two apart is what makes a census readable: a predicate that also required analytic[i](../../../../../../reference/source/jvmMain/afield/audio/synth/Subharmonic.md#L240).isFinite() would file a NaN from the thing under test as "no verdict here", in the same bucket as a genuine singularity, and report high coverage over a derivative that was NaN throughout. See nonFiniteAnalytic, and adBox3's KDoc.

    val usableComponents: BooleanArray = BooleanArray(analytic.size) { i ->
        numerical[i].isFinite() && numericalHalfStep[i].isFinite() &&
            abs(numerical[i] - numericalHalfStep[i]) /
            max(1.0, max(abs(numerical[i]), abs(numericalHalfStep[i]))) < 1e-4
    }

val nonFiniteAnalytic: List<Int>

Components where the oracle is trustworthy and the analytic gradient is not a number.

Always a bug, never a property of the geometry: the numerical side having a confident opinion means the function really is differentiable there, so a NaN or an Infinity from the lowering is the lowering's fault. Must be empty; SolidSLGradientTest asserts it over the whole census.

    val nonFiniteAnalytic: List<Int>
        get() = analytic.indices.filter { usableComponents[it] && !analytic[it].isFinite() }

val usableCount: Int get() = usableComponents.count { it }

How many of them, i.e. how much evidence this check actually carries.

val usable: Boolean get() = usableCount > 0

Whether there is any verdict here at all.

val maxError: Double

The disagreement over the usable components, measured against a floored scale — max(1, |analytic|, |numerical|).

A relative error against a derivative that happens to be 1e-9 is meaningless and fails any tolerance; a floor of 1 makes it an absolute test exactly where an absolute test is the right one, and a relative one where the derivative is large (a Scale factor's 1/s² reaches the thousands).

0.0 when nothing is usable — which is why a caller must consult usableCount and not read a zero here as a pass.

    val maxError: Double
        get() = analytic.indices.filter { usableComponents[it] }.maxOfOrNull { i ->
            abs(analytic[i] - numerical[i]) / max(1.0, max(abs(analytic[i]), abs(numerical[i])))
        } ?: 0.0

val worstComponent: Int

Which component disagreed the most, or -1 when there is no verdict.

    val worstComponent: Int
        get() = analytic.indices.filter { usableComponents[it] }.maxByOrNull { i ->
            abs(analytic[i] - numerical[i]) / max(1.0, max(abs(analytic[i]), abs(numerical[i])))
        } ?: -1

override fun equals(other: Any?): Boolean

    override fun equals(other: Any?): Boolean =
        other is SolidGradientCheck && analytic.contentEquals(other.analytic) &&
            numerical.contentEquals(other.numerical) &&
            numericalHalfStep.contentEquals(other.numericalHalfStep) && step == other.step

override fun hashCode(): Int

    override fun hashCode(): Int =
        31 * (31 * analytic.contentHashCode() + numerical.contentHashCode()) + step.hashCode()

override fun toString(): String

    override fun toString(): String =
        "$usableCount/$size usable" +
            (if (nonFiniteAnalytic.isNotEmpty()) ", NON-FINITE ANALYTIC at ${nonFiniteAnalytic}" else "") +
            ", max floored relative error ${"%.3e".format(maxError)} at " +
            "component $worstComponent  analytic [${analytic.joinToString(", ") { "%.6g".format(it) }}]  " +
            "numerical [${numerical.joinToString(", ") { "%.6g".format(it) }}]"
}

fun SolidTemplate.checkGradient

The parameter gradient at one point, beside its central-difference oracle.

Both directions matter and they are checked separately, because they fail differently: a wrong position gradient is a wrong normal and shows up in a render, while a wrong parameter gradient is a fit that walks the wrong way and shows up as a plausible-looking non-answer.

@param h the difference step. 1e-5 is the sweet spot for a Double field of unit scale — big enough that the O(h²) truncation is around 1e-10 and small enough that cancellation in (f(x+h) − f(x−h)) leaves ten digits.

fun SolidTemplate.checkGradient(
    x: Double,
    y: Double,
    z: Double,
    values: DoubleArray = seed,
    h: Double = 1e-5,
): SolidGradientCheck {
    val analytic = gradientAt(x, y, z, values).dParameters
    val numerical = DoubleArray(size) { central(x, y, z, values, it, h) }
    val half = DoubleArray(size) { central(x, y, z, values, it, h / 2.0) }
    return SolidGradientCheck(analytic, numerical, half, h)
}

fun SolidTemplate.checkPositionGradient

checkGradient for the position gradient — the surface normal.

fun SolidTemplate.checkPositionGradient(
    x: Double,
    y: Double,
    z: Double,
    values: DoubleArray = seed,
    h: Double = 1e-5,
): SolidGradientCheck {
    val analytic = gradientAt(x, y, z, values).dPosition
    val field = withValues(values).toField3()
    fun at(step: Double) = doubleArrayOf(
        (field.at(x + step, y, z) - field.at(x - step, y, z)) / (2.0 * step),
        (field.at(x, y + step, z) - field.at(x, y - step, z)) / (2.0 * step),
        (field.at(x, y, z + step) - field.at(x, y, z - step)) / (2.0 * step),
    )
    return SolidGradientCheck(analytic, at(h), at(h / 2.0), h)
}

One component of the central difference.

The tree is rebuilt and re-lowered per probe, which is O(nodes) and looks wasteful. The waste is the method: substituting into a SolidExpr and lowering it afresh is exactly what the analytic side is being compared against, so a shared cache between the two would be a shared assumption, and a shared assumption is the one thing a second opinion must not have.

private fun SolidTemplate.central(
    x: Double,
    y: Double,
    z: Double,
    values: DoubleArray,
    i: Int,
    h: Double,
): Double {
    val up = values.copyOf(); up[i] += h
    val down = values.copyOf(); down[i] -= h
    val fu = withValues(up).toField3().at(x, y, z)
    val fd = withValues(down).toField3().at(x, y, z)
    return (fu - fd) / (2.0 * h)
}

The fit

data class SolidFitOptions

The dials on SolidTemplate.fitParameters.

@param band clamps both the prediction and the target into [-band, +band] before differencing, turning a whole-field objective into a near-surface one. null leaves the residual alone. This is the standard truncated-SDF trick. It stops a sample far outside the shape, where every candidate is roughly equally wrong, from outvoting the samples that decide what the surface looks like — and it does so by killing the gradient of every sample whose prediction is outside the band, whatever the target says there. So only samples near the model's current surface steer the descent: a sphere seeded at half its target's radius with band = 0.1 still walks all the way there, but a target sampled only where the model reads beyond the band leaves every parameter in SolidFit.deadParameters. On an approximation the band trades the far field for the surface (three spheres against a rounded rod, band = 0.05: near-surface rms 0.0595 → 0.0484, whole-field rms 0.0506 → 0.0620). The default is null because on a recovery problem the far field is informative. SolidSLFitBandTest pins all of this. @param method OptimizeMethod.LBFGS by default, and measured to be the right default on every fit problem solidsl has (SolidSLOptimizerComparisonTest, four starts each, after adBox3's interior NaN was fixed): an exact plate recovery, three spheres approximating a rod, a cube and bar approximating a sphere and bar, and the point-cloud lobes. On every one LBFGS reached the same minimum Adam reached — or, on the exact recovery, a lower one — in 6 to 38 iterations against Adam's 340 to 6000, and 4 to 130 times less wall time. The line search absorbs the bad relative scaling SolidExpr.gradientNotes describes (an angle in degrees beside a length), which is exactly what costs Adam's single learningRate thousands of iterations. Earlier advice to prefer OptimizeMethod.ADAM on "rough" landscapes was measured while the gradient was NaN inside every box, and does not survive the fix. Adam remains available; reach for it when an LBFGS run stops on line search after very few iterations — but check SolidGradientCheck.nonFiniteAnalytic first, because that is far more often a broken gradient than a hard landscape. @param iterations bounded, like everything else in autodiff_optimize.kt, so a descent cannot disappear into an animation loop.

data class SolidFitOptions(
    val band: Double? = null,
    val method: OptimizeMethod = OptimizeMethod.LBFGS,
    val iterations: Int = 200,
    val gradientTolerance: Double = 1e-10,
    val learningRate: Double = 0.05,
    val memory: Int = 8,
) {

internal fun toOptimizeOptions(): OptimizeOptions

    internal fun toOptimizeOptions(): OptimizeOptions = OptimizeOptions(
        method = method,
        iterations = iterations,
        gradientTolerance = gradientTolerance,
        learningRate = learningRate,
        memory = memory,
        // A fit is a library call, not an animation: a progress reporter on it would print from
        // inside every test and every headless docs capture.
        showProgress = false,
    )
}

data class SolidFit

What a fit did, and what it started from.

initialLoss is not decoration. "The fit improved things" is unfalsifiable without the counterfactual in the same object, measured in the same run, by the same code: a remembered before-number is a number that can go stale while the after-number stays green. improvementFactor is the ratio, and it is the thing to assert a floor on.

data class SolidFit(
    val template: SolidTemplate,
    /** Where the descent started. */
    val start: DoubleArray,
    /** Where it finished. */
    val values: DoubleArray,
    /** The loss at [start], measured now. */
    val initialLoss: Double,
    /** The loss at [values], measured now. */
    val loss: Double,
    /** The fitted parameters substituted back in — an ordinary [SolidExpr]. */
    val fitted: SolidExpr,
    /** The whole descent, for a loss curve or a convergence reason. */
    val run: OptimizeRun,
    /**
     * Parameters whose gradient was **exactly zero** at [start], or that the whole descent left
     * exactly where they began — the ones that could not, or did not, move. The second catches a
     * child that is never the nearest surface anywhere: a cube wholly inside a bar is invisible to
     * the loss, and its parameter comes back bit-identical to the seed.
     *
     * See the note on [afield.lang.solidsl.SolidFit] about what causes it. This is reported rather than
     * thrown on, because a legitimately-dead parameter is a fine thing to have in a vector (a fit
     * over six holes where one is genuinely irrelevant), but an *unexpectedly* dead one is the
     * failure that looks like a success.
     */
    val deadParameters: List<String>,
) {

val improvementFactor: Double get()

initialLoss / loss — how many times better the answer is than the seed.

    val improvementFactor: Double get() = if (loss <= 0.0) Double.POSITIVE_INFINITY else initialLoss / loss

fun shapeAt(index: Int): SolidExpr

The shape the descent was holding at iteration index — an ordinary SolidExpr, marchable like any other.

Indexed like afield.math.autodiff.OptimizeRun.trajectory and therefore like afield.math.autodiff.OptimizeRun.losses: shapeAt(0) is the seed and shapeAt(run.losses.lastIndex) is fitted. This is what makes a convergence filmstrip a recording rather than a set of re-runs truncated at different iteration caps — the frames are the points the optimizer actually passed through, not points a second descent would have reached.

    fun shapeAt(index: Int): SolidExpr {
        val t = run.trajectory
        require(index in t.indices) { "iteration $index is outside the recorded trajectory 0..${t.lastIndex}" }
        return template.withValues(t[index])
    }

fun shapesAt(indices: List<Int>): List<SolidExpr> = indices.map(::shapeAt)

shapeAt for a whole strip, in the order given.

val rmsError: Double get() = sqrt(max(0.0, loss))

The residual as an RMS distance, which is comparable with SolidSamples.rms.

fun report(): String

The whole story, for a println in a box.

    fun report(): String = buildString {
        appendLine("fit ${template.names.size} parameter(s) in ${run.iterations} iterations (${run.reason})")
        template.names.forEachIndexed { i, name ->
            appendLine("  %-14s %12.6f  ->  %12.6f".format(name, start[i], values[i]))
        }
        appendLine("  loss           ${"%.6e".format(initialLoss)}  ->  ${"%.6e".format(loss)}")
        appendLine("  rms distance   ${"%.6e".format(sqrt(max(0.0, initialLoss)))}  ->  ${"%.6e".format(rmsError)}")
        appendLine("  improvement    ${"%.4g".format(improvementFactor)}x")
        if (deadParameters.isNotEmpty()) {
            appendLine("  DEAD (zero gradient at the start, or never moved): ${deadParameters.joinToString(", ")}")
        }
        append("  fitted         $fitted")
    }

override fun equals(other: Any?): Boolean = this === other

override fun hashCode(): Int = System.identityHashCode(this)

fun SolidTemplate.fitParameters

Fit this template's named parameters so its field matches target, by gradient descent on the mean squared distance residual.

The loss is mean over samples of (predicted − target)², and every gradient in it is exact: one ADTape per evaluation carries all SolidSamples.size residuals and all SolidTemplate.size parameters, and one backward sweep reads every derivative off it. That is the whole reason to have built lowering (c) — the alternative, central-differencing the loss, costs 2n extra full sweeps per step and needs an fdStep tuned per problem, which is the complaint autodiff.kt opens with.

val fit = template.fitParameters(target, start = doubleArrayOf(0.35, 2.6))
fit.improvementFactor          // assert a floor on this
fit.fitted.crossCheck()        // and then check the answer with the other two lowerings

@param start where to begin; the template's own SolidTemplate.seed by default. A recovery test passes something measurably far away and asserts how far back it came.

body — 49 lines, from line 580
fun SolidTemplate.fitParameters(
    target: SolidSamples,
    start: DoubleArray = seed,
    options: SolidFitOptions = SolidFitOptions(),
): SolidFit {
    require(start.size == size) {
        "this template has $size parameter(s) (${names.joinToString(", ")}), got ${start.size} start values"
    }
    val band = options.band
    require(band == null || band > 0.0) { "band must be positive when given, got $band" }

    val valueAndGradient: (DoubleArray) -> ADGradient = { p ->
        valueAndGrad(p) { parameters ->
            val tape = parameters[0].tape
            val field = adFieldWith(parameters)
            var loss = tape.constant(0.0)
            for (i in 0 until target.size) {
                val raw = field.at(
                    tape.constant(target.xs[i]),
                    tape.constant(target.ys[i]),
                    tape.constant(target.zs[i]),
                )
                val pred = if (band == null) raw else raw.clamp(-band, band)
                val want = if (band == null) target.distances[i] else target.distances[i].coerceIn(-band, band)
                val e = pred - want
                loss = loss + e * e
            }
            loss / target.size.toDouble()
        }
    }

    val at0 = valueAndGradient(start)
    val dead = names.indices.filter { at0.gradient[it] == 0.0 }.map { names[it] }

    val run = minimizeWithGradient(start, options.toOptimizeOptions(), valueAndGradient)
    val end = valueAndGradient(run.x)
    val stuck = neverMoved(start, run).map { names[it] }

    return SolidFit(
        template = this,
        start = start.copyOf(),
        values = run.x.copyOf(),
        initialLoss = at0.value,
        loss = end.value,
        fitted = withValues(run.x),
        run = run,
        deadParameters = names.filter { it in dead || it in stuck },
    )
}

internal fun neverMoved(start: DoubleArray, run: OptimizeRun): List<Int>

Parameters the descent never moved off their seed: every recorded iterate holds start exactly.

internal fun neverMoved(start: DoubleArray, run: OptimizeRun): List<Int> =
    if (run.iterations == 0) emptyList() else start.indices.filter { i -> run.trajectory.all { it[i] == start[i] } && run.x[i] == start[i] }