smile_interpolation_regression.kt

src/jvmMain/kotlin/afield/math/statistics/smile_interpolation_regression.kt · 168 lines

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

package afield.math.statistics

import afield.data.collections.invoke
import afield.data.collections.s
import afield.data.collections.times
import afield.document.graph.minusAssign
import afield.document.graph.plusAssign
import afield.document.properties.compareTo
import afield.document.workspace.file
import afield.geometry.mesh.normalized
import afield.graphics.render.rangeTo
import afield.graphics.shader.callable
import afield.math.autodiff.compareTo
import afield.math.autodiff.div
import afield.math.autodiff.minus
import afield.math.autodiff.plus
import afield.math.autodiff.times
import afield.math.linalg.invoke
import afield.math.linalg.remAssign
import afield.math.linalg.set
import afield.math.linalg.times
import afield.math.units.div
import afield.math.units.minus
import afield.math.units.plus
import afield.math.units.times
import afield.media.image.normalized

import afield.audio.synth.*
import afield.data.collections.get
import afield.math.autodiff.abs
import afield.math.autodiff.maxOf
import afield.math.linalg.i
import afield.math.linalg.x
import afield.math.gp.GpJoint
import afield.math.gp.GpPrediction
import smile.interpolation.CubicSplineInterpolation1D
import smile.math.kernel.GaussianKernel
import smile.regression.GaussianProcessRegression

Two 1D fitters that wrap the Smile library: smileCubicSpline (cubic spline interpolation through the data) and smileGaussianProcess (Gaussian-process regression with a Gaussian kernel, configured by SmileGpOptions).

Both take parallel xs/ys arrays, copy them, and require at least two finite samples with xs strictly increasing — a violation throws rather than producing a bad fit with no complaint. The returned SmileCubicSpline1D / SmileGaussianProcess1D are callable via invoke and also expose predict for a single x or for a whole array.

GP fitting failures propagate. No interpolator is substituted for a failed probabilistic model. For 1D/2D/3D inference, composed kernels, heteroscedastic noise and coherent samples see afield.math.gp.GaussianProcess.

data class SmileGpOptions

data class SmileGpOptions(
    val lengthScale: Double = 0.15,
    /** Smile uses noise VARIANCE in normalized output units when normalize is enabled. New GP APIs use noiseStd. */
    val noise: Double = 1e-3,
    val normalize: Boolean = true,
    val tol: Double = 1e-4,
    val maxIter: Int = 500
)

class SmileCubicSpline1D internal constructor

class SmileCubicSpline1D internal constructor(
    val xs: DoubleArray,
    val ys: DoubleArray,
    private val model: CubicSplineInterpolation1D
) {

fun predict(x: Double): Double = model.interpolate(x)

fun predict(xs: DoubleArray): DoubleArray

    fun predict(xs: DoubleArray): DoubleArray =
        DoubleArray(xs.size) { i -> predict(xs[i]) }

operator fun invoke(x: Double): Double = predict(x)

class SmileGaussianProcess1D internal constructor

class SmileGaussianProcess1D internal constructor(
    val xs: DoubleArray,
    val ys: DoubleArray,
    val options: SmileGpOptions,
    private val model: GaussianProcessRegression<DoubleArray>
) {

fun predict(x: Double): Double = model.predict(doubleArrayOf(x))

fun predictWithVariance(x: Double): GpPrediction

Latent posterior uncertainty in the original output units.

    fun predictWithVariance(x: Double): GpPrediction {
        require(x.isFinite())
        val estimate = DoubleArray(2)
        model.predict(doubleArrayOf(x), estimate)
        require(estimate.all { it.isFinite() }) { "Smile returned a non-finite prediction" }
        return GpPrediction(estimate[0], estimate[1] * estimate[1])
    }

fun posterior(xs: DoubleArray): GpJoint

Seedable coherent draws through Field's semidefinite-aware joint sampler.

    fun posterior(xs: DoubleArray): GpJoint {
        require(xs.isNotEmpty() && xs.size <= 4096 && xs.all { it.isFinite() })
        val prediction = model.query(Array(xs.size) { doubleArrayOf(xs[it]) })
        val covariance = Array(xs.size) { i -> DoubleArray(xs.size) { j -> prediction.cov.get(i, j) } }
        val scale = covariance.indices.maxOf { kotlin.math.abs(covariance[it][it]) }
        return GpJoint(prediction.mu.copyOf(), covariance, 1e-12, 1e-10 * scale)
    }

fun predict(xs: DoubleArray): DoubleArray

    fun predict(xs: DoubleArray): DoubleArray =
        DoubleArray(xs.size) { i -> predict(xs[i]) }

operator fun invoke(x: Double): Double = predict(x)

fun smileCubicSpline(xs: DoubleArray, ys: DoubleArray): SmileCubicSpline1D

fun smileCubicSpline(xs: DoubleArray, ys: DoubleArray): SmileCubicSpline1D {
    validateTrainingSeries(xs, ys, "smileCubicSpline")
    val xCopy = xs.copyOf()
    val yCopy = ys.copyOf()
    val spline = CubicSplineInterpolation1D(xCopy, yCopy)
    return SmileCubicSpline1D(xCopy, yCopy, spline)
}

fun smileGaussianProcess

fun smileGaussianProcess(
    xs: DoubleArray,
    ys: DoubleArray,
    options: SmileGpOptions = SmileGpOptions()
): SmileGaussianProcess1D {
    validateTrainingSeries(xs, ys, "smileGaussianProcess")
    require(options.lengthScale.isFinite() && options.lengthScale > 0.0) {
        "smileGaussianProcess: lengthScale must be > 0, found ${options.lengthScale}"
    }
    require(options.noise.isFinite() && options.noise >= 0.0) {
        "smileGaussianProcess: noise must be >= 0, found ${options.noise}"
    }
    require(options.tol.isFinite() && options.tol > 0.0) {
        "smileGaussianProcess: tol must be > 0, found ${options.tol}"
    }
    require(options.maxIter >= 0) {
        "smileGaussianProcess: maxIter must be >= 0, found ${options.maxIter}"
    }

    val xCopy = xs.copyOf()
    val yCopy = ys.copyOf()
    val xTrain = Array(xCopy.size) { i -> doubleArrayOf(xCopy[i]) }
    val kernel = GaussianKernel(options.lengthScale)
    // Smile's normalization divides by sample standard deviation; constant data has zero spread.
    val normalize = options.normalize && yCopy.any { it != yCopy[0] }
    val gpOptions = GaussianProcessRegression.Options(options.noise, normalize, options.tol, options.maxIter)
    val model = GaussianProcessRegression.fit(xTrain, yCopy, kernel, gpOptions)
    return SmileGaussianProcess1D(xCopy, yCopy, options, model)
}

private fun validateTrainingSeries(xs: DoubleArray, ys: DoubleArray, fn: String) {
    require(xs.size == ys.size) {
        "$fn: xs.size (${xs.size}) must match ys.size (${ys.size})"
    }
    require(xs.size >= 2) {
        "$fn: at least 2 samples required, got ${xs.size}"
    }

    for (i in xs.indices) {
        require(xs[i].isFinite()) { "$fn: xs[$i] must be finite, found ${xs[i]}" }
        require(ys[i].isFinite()) { "$fn: ys[$i] must be finite, found ${ys[i]}" }
    }

    for (i in 1 until xs.size) {
        require(xs[i] > xs[i - 1]) {
            "$fn: xs must be strictly increasing; xs[$i]=${xs[i]} <= xs[${i - 1}]=${xs[i - 1]}"
        }
    }
}