untoldwind / KontrolSystem2

Autopilot scripting system for KSP2
Other
54 stars 14 forks source link

Feature request: Higher dimensional vector #167

Open versvisa opened 2 months ago

versvisa commented 2 months ago

For various numerical calculations it would be useful to have an N dimensional float vector, with standard operations like addition, scalar multiplication, dot product, etc.

It looks like right now the only option is to use arrays and functions, e.g. pub sync fn array_addition(a: float[], b: float[]) -> float[]

untoldwind commented 2 months ago

A potential alternative would be to support operator overloading in the language itself. ... I will have to play around a bit to see what impact that would have

untoldwind commented 2 months ago

The 0.5.9.0 pre-release (https://github.com/untoldwind/KontrolSystem2/releases/tag/v0.5.9.0) now has experimental support for operator overloading.

To quote the example from the release-description:

pub struct NVector(values: float[]) {
    dim: int = values.length
    values: float[] = values
}

impl NVector {
    sync fn add_scalar(self, scalar: float) -> NVector = NVector(self.values.map(fn(v) -> v + scalar))

    sync fn mul_scalar(self, scalar: float) -> NVector = NVector(self.values.map(fn(v) -> v * scalar))

    sync fn to_string(self) -> string = $"NVector({self.dim}, {self.values.to_string()})"
}

impl operators for NVector {
    sync fn neg(right: NVector) -> NVector = NVector(right.values.map(fn(r) -> -r))

    sync fn add(left: float, right: NVector) -> NVector = right.add_scalar(left)

    sync fn add(left: NVector, right: float) -> NVector = left.add_scalar(right)

    sync fn add(left: NVector, right: NVector) -> NVector = NVector((0..min_dim(left.dim, right.dim)).map(fn(i) -> left.values[i] + right.values[i]))

    sync fn mul(left: float, right: NVector) -> NVector = right.mul_scalar(left)

    sync fn mul(left: NVector, right: float) -> NVector = left.mul_scalar(right)

    sync fn div(left: NVector, right: float) -> NVector = NVector(left.values.map(fn(l) -> l / right))
}

pub sync fn zero_nvector(dim: int) -> NVector = NVector((0..dim).map(fn(i) -> 0.0))

sync fn min_dim(dim1: int, dim2: int) -> int = if(dim1 < dim2) dim1 else dim2

which allows:

use * from core::testing
use { NVector, zero_nvector } from nvector

test sync fn test_nvector_operators() -> Unit = {
    const vec1 = zero_nvector(4)

    assert_string("NVector(4, [0, 0, 0, 0])", vec1.to_string())

    const vec2 = 3.0 + vec1

    assert_string("NVector(4, [3, 3, 3, 3])", vec2.to_string())

    const vec3 = vec1 + 4.0

    assert_string("NVector(4, [4, 4, 4, 4])", vec3.to_string())

    const vec4 = NVector([1.0, 2.0, 3.0, 4.0])

    const vec5 = vec2 + vec4

    assert_string("NVector(4, [4, 5, 6, 7])", vec5.to_string())

    const vec6 = -vec5

    assert_string("NVector(4, [-4, -5, -6, -7])", vec6.to_string())

    const vec7 = 2 * vec4

    assert_string("NVector(4, [2, 4, 6, 8])", vec7.to_string())

    const vec8 = vec4 * 3

    assert_string("NVector(4, [3, 6, 9, 12])", vec8.to_string())

    const vec9 = vec8 / 2

    assert_string("NVector(4, [1.5, 3, 4.5, 6])", vec9.to_string())
}
versvisa commented 2 months ago

Awesome, will try it when I find some free time

versvisa commented 2 months ago

It does work for me. Here is a simplified version of the vector and a first use-case for it.

A square bracket operator for overloading indexing would be useful.

pub struct Vector(values: float[]) {
    values: float[] = values
}

impl Vector {
    sync fn to_string(self) -> string = $"Vector{self.values.to_string()}"
}

impl operators for Vector {
    sync fn add(left: Vector, right: Vector) -> Vector = Vector((0..left.values.length).map(fn(i) -> left.values[i] + right.values[i]))
    sync fn mul(scalar: float, vec: Vector) -> Vector = Vector(vec.values.map(fn(v) -> v * scalar))
    sync fn mul(vec: Vector, scalar: float) -> Vector = scalar * vec
}

pub sync fn RungeKuttaStep(ode: fn(Vector) -> Vector, x0: Vector, dt: float) -> Vector =
{
    const k1 = ode(x0)
    const k2 = ode(x0 + (0.5 * dt) * k1)
    const k3 = ode(x0 + (0.5 * dt) * k2)
    const k4 = ode(x0 + (      dt) * k3)
    x0 + (0.16666666666666666 * dt) * ((k1 + k4) + (2.0 * (k2 + k3)))
}

test sync fn test_RungeKuttaStep() -> Result<Unit> =
{
    const my_ode = fn(x:Vector) -> Vector([-x.values[1], x.values[0]]) // Simple harmonic oscillator
    for(i in 1..300)
    {
        const dt = i * 0.01
        const x1 = RungeKuttaStep(my_ode, Vector([2.0, 0.0]), dt)
        core::testing::assert_float(2.0 * core::math::cos(dt), x1.values[0], 0.02 * dt**(5.0))
        core::testing::assert_float(2.0 * core::math::sin(dt), x1.values[1], 0.02 * dt**(5.0))
    }
}
untoldwind commented 1 month ago

I will do a bit more testing over the weekend and make it part of the "official" release.