build.gradle.kts

@@ -17,7 +17,7 @@ allprojects {
 subprojects {
     if (name.startsWith("kmath")) apply<MavenPublishPlugin>()
 
-    plugins.withId("org.jetbrains.dokka"){
+    plugins.withId("org.jetbrains.dokka") {
         tasks.withType<org.jetbrains.dokka.gradle.DokkaTaskPartial> {
             dependsOn(tasks["assemble"])
 
@@ -51,6 +51,18 @@ subprojects {
             }
         }
     }
+
+    plugins.withId("org.jetbrains.kotlin.multiplatform") {
+        configure<org.jetbrains.kotlin.gradle.dsl.KotlinMultiplatformExtension> {
+            sourceSets {
+                val commonTest by getting {
+                    dependencies {
+                        implementation(projects.testUtils)
+                    }
+                }
+            }
+        }
+    }
 }
 
 readme.readmeTemplate = file("docs/templates/README-TEMPLATE.md")

kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/Quaternion.kt

@@ -24,9 +24,12 @@ import kotlin.math.*
  * @property y The third component.
  * @property z The fourth component.
  */
-public data class Quaternion(
+public class Quaternion(
-    val w: Double, val x: Double, val y: Double, val z: Double,
+    public val w: Double,
-) {
+    public val x: Double,
+    public val y: Double,
+    public val z: Double,
+) : Buffer<Double> {
     init {
         require(!w.isNaN()) { "w-component of quaternion is not-a-number" }
         require(!x.isNaN()) { "x-component of quaternion is not-a-number" }
@@ -39,12 +42,49 @@ public data class Quaternion(
      */
     override fun toString(): String = "($w + $x * i + $y * j + $z * k)"
 
-    public companion object : MemorySpec<Quaternion> {
+    override val size: Int get() = 4
-        override val objectSize: Int
-            get() = 32
 
-        override fun MemoryReader.read(offset: Int): Quaternion =
+    override fun get(index: Int): Double = when (index) {
-            Quaternion(readDouble(offset), readDouble(offset + 8), readDouble(offset + 16), readDouble(offset + 24))
+        0 -> w
+        1 -> x
+        2 -> y
+        3 -> z
+        else -> error("Index $index out of bounds [0,3]")
+    }
+
+    override fun iterator(): Iterator<Double> = listOf(w, x, y, z).iterator()
+
+    override fun equals(other: Any?): Boolean {
+        if (this === other) return true
+        if (other == null || this::class != other::class) return false
+
+        other as Quaternion
+
+        if (w != other.w) return false
+        if (x != other.x) return false
+        if (y != other.y) return false
+        if (z != other.z) return false
+
+        return true
+    }
+
+    override fun hashCode(): Int {
+        var result = w.hashCode()
+        result = 31 * result + x.hashCode()
+        result = 31 * result + y.hashCode()
+        result = 31 * result + z.hashCode()
+        return result
+    }
+
+    public companion object : MemorySpec<Quaternion> {
+        override val objectSize: Int get() = 32
+
+        override fun MemoryReader.read(offset: Int): Quaternion = Quaternion(
+            readDouble(offset),
+            readDouble(offset + 8),
+            readDouble(offset + 16),
+            readDouble(offset + 24)
+        )
 
         override fun MemoryWriter.write(offset: Int, value: Quaternion) {
             writeDouble(offset, value.w)
@@ -66,30 +106,22 @@ public fun Quaternion(w: Number, x: Number = 0.0, y: Number = 0.0, z: Number = 0
  * This quaternion's conjugate.
  */
 public val Quaternion.conjugate: Quaternion
-    get() = QuaternionField { z - x * i - y * j - z * k }
+    get() = Quaternion(w, -x, -y, -z)
 
 /**
  * This quaternion's reciprocal.
  */
 public val Quaternion.reciprocal: Quaternion
     get() {
-        QuaternionField {
+        val norm2 = (w * w + x * x + y * y + z * z)
-            val n = norm(this@reciprocal)
+        return Quaternion(w / norm2, -x / norm2, -y / norm2, -z / norm2)
-            return conjugate / (n * n)
     }
-    }
-
-/**
- * Absolute value of the quaternion.
- */
-public val Quaternion.r: Double
-    get() = sqrt(w * w + x * x + y * y + z * z)
 
 /**
  * A field of [Quaternion].
  */
 @OptIn(UnstableKMathAPI::class)
-public object QuaternionField : Field<Quaternion>, Norm<Quaternion, Quaternion>, PowerOperations<Quaternion>,
+public object QuaternionField : Field<Quaternion>, Norm<Quaternion, Double>, PowerOperations<Quaternion>,
     ExponentialOperations<Quaternion>, NumbersAddOps<Quaternion>, ScaleOperations<Quaternion> {
     override val zero: Quaternion = Quaternion(0.0)
     override val one: Quaternion = Quaternion(1.0)
@@ -217,7 +249,12 @@ public object QuaternionField : Field<Quaternion>, Norm<Quaternion, Quaternion>,
         Quaternion(toDouble() * arg.w, toDouble() * arg.x, toDouble() * arg.y, toDouble() * arg.z)
 
     override fun Quaternion.unaryMinus(): Quaternion = Quaternion(-w, -x, -y, -z)
-    override fun norm(arg: Quaternion): Quaternion = sqrt(arg.conjugate * arg)
+    override fun norm(arg: Quaternion): Double = sqrt(
+        arg.w.pow(2) +
+                arg.x.pow(2) +
+                arg.y.pow(2) +
+                arg.z.pow(2)
+    )
 
     override fun bindSymbolOrNull(value: String): Quaternion? = when (value) {
         "i" -> i

kmath-complex/src/commonTest/kotlin/space/kscience/kmath/complex/QuaternionTest.kt

@@ -6,10 +6,23 @@
 package space.kscience.kmath.complex
 
 import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.testutils.assertBufferEquals
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
-internal class QuaternionFieldTest {
+internal class QuaternionTest {
+
+    @Test
+    fun testNorm() {
+        assertEquals(2.0, QuaternionField.norm(Quaternion(1.0, 1.0, 1.0, 1.0)))
+    }
+
+    @Test
+    fun testInverse() = QuaternionField {
+        val q = Quaternion(1.0, 2.0, -3.0, 4.0)
+        assertBufferEquals(one, q * q.reciprocal, 1e-4)
+    }
+
     @Test
     fun testAddition() {
         assertEquals(Quaternion(42, 42), QuaternionField { Quaternion(16, 16) + Quaternion(26, 26) })

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/Euclidean3DSpace.kt

@@ -8,6 +8,7 @@ package space.kscience.kmath.geometry
 import space.kscience.kmath.linear.Point
 import space.kscience.kmath.operations.ScaleOperations
 import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.Buffer
 import kotlin.math.sqrt
 
 public interface Vector3D : Point<Double>, Vector {
@@ -29,6 +30,19 @@ public interface Vector3D : Point<Double>, Vector {
 @Suppress("FunctionName")
 public fun Vector3D(x: Double, y: Double, z: Double): Vector3D = Vector3DImpl(x, y, z)
 
+public fun Buffer<Double>.asVector3D(): Vector3D = object : Vector3D {
+    init {
+        require(this@asVector3D.size == 3) { "Buffer of size 3 is required for Vector3D" }
+    }
+
+    override val x: Double get() = this@asVector3D[0]
+    override val y: Double get() = this@asVector3D[1]
+    override val z: Double get() = this@asVector3D[2]
+
+    override fun toString(): String = this@asVector3D.toString()
+
+}
+
 public val Vector3D.r: Double get() = Euclidean3DSpace { norm() }
 
 private data class Vector3DImpl(

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/projections.kt

@@ -5,6 +5,8 @@
 
 package space.kscience.kmath.geometry
 
+//TODO move vector to receiver
+
 /**
  * Project vector onto a line.
  * @param vector to project

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/rotations3D.kt

@@ -0,0 +1,100 @@
+/*
+ * Copyright 2018-2021 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.geometry
+
+import space.kscience.kmath.complex.Quaternion
+import space.kscience.kmath.complex.QuaternionField
+import space.kscience.kmath.complex.reciprocal
+import space.kscience.kmath.linear.LinearSpace
+import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.linear.linearSpace
+import space.kscience.kmath.linear.matrix
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.operations.invoke
+import kotlin.math.pow
+import kotlin.math.sqrt
+
+internal fun Vector3D.toQuaternion(): Quaternion = Quaternion(0.0, x, y, z)
+
+/**
+ * Angle in radians denoted by this quaternion rotation
+ */
+public val Quaternion.theta: Double get() = kotlin.math.acos(w) * 2
+
+/**
+ * An axis of quaternion rotation
+ */
+public val Quaternion.vector: Vector3D
+    get() {
+        val sint2 = sqrt(1 - w * w)
+
+        return object : Vector3D {
+            override val x: Double get() = this@vector.x/sint2
+            override val y: Double get() = this@vector.y/sint2
+            override val z: Double get() = this@vector.z/sint2
+            override fun toString(): String = listOf(x, y, z).toString()
+        }
+    }
+
+/**
+ * Rotate a vector in a [Euclidean3DSpace]
+ */
+public fun Euclidean3DSpace.rotate(vector: Vector3D, q: Quaternion): Vector3D = with(QuaternionField) {
+    val p = vector.toQuaternion()
+    (q * p * q.reciprocal).vector
+}
+
+/**
+ * Use a composition of quaternions to create a rotation
+ */
+public fun Euclidean3DSpace.rotate(vector: Vector3D, composition: QuaternionField.() -> Quaternion): Vector3D =
+    rotate(vector, QuaternionField.composition())
+
+public fun Euclidean3DSpace.rotate(vector: Vector3D, matrix: Matrix<Double>): Vector3D {
+    require(matrix.colNum == 3 && matrix.rowNum == 3) { "Square 3x3 rotation matrix is required" }
+    return with(DoubleField.linearSpace) { matrix.dot(vector).asVector3D() }
+}
+
+/**
+ * Convert a [Quaternion] to a rotation matrix
+ */
+@OptIn(UnstableKMathAPI::class)
+public fun Quaternion.toRotationMatrix(
+    linearSpace: LinearSpace<Double, *> = DoubleField.linearSpace,
+): Matrix<Double> {
+    val s = QuaternionField.norm(this).pow(-2)
+    return linearSpace.matrix(3, 3)(
+        1.0 - 2 * s * (y * y + z * z), 2 * s * (x * y - z * w), 2 * s * (x * z + y * w),
+        2 * s * (x * y + z * w), 1.0 - 2 * s * (x * x + z * z), 2 * s * (y * z - x * w),
+        2 * s * (x * z - y * w), 2 * s * (y * z + x * w), 1.0 - 2 * s * (x * x + y * y)
+    )
+}
+
+/**
+ * taken from https://d3cw3dd2w32x2b.cloudfront.net/wp-content/uploads/2015/01/matrix-to-quat.pdf
+ */
+public fun Quaternion.Companion.fromRotationMatrix(matrix: Matrix<Double>): Quaternion {
+    val t: Double
+    val q = if (matrix[2, 2] < 0) {
+        if (matrix[0, 0] > matrix[1, 1]) {
+            t = 1 + matrix[0, 0] - matrix[1, 1] - matrix[2, 2]
+            Quaternion(t, matrix[0, 1] + matrix[1, 0], matrix[2, 0] + matrix[0, 2], matrix[1, 2] - matrix[2, 1])
+        } else {
+            t = 1 - matrix[0, 0] + matrix[1, 1] - matrix[2, 2]
+            Quaternion(matrix[0, 1] + matrix[1, 0], t, matrix[1, 2] + matrix[2, 1], matrix[2, 0] - matrix[0, 2])
+        }
+    } else {
+        if (matrix[0, 0] < -matrix[1, 1]) {
+            t = 1 - matrix[0, 0] - matrix[1, 1] + matrix[2, 2]
+            Quaternion(matrix[2, 0] + matrix[0, 2], matrix[1, 2] + matrix[2, 1], t, matrix[0, 1] - matrix[1, 0])
+        } else {
+            t = 1 + matrix[0, 0] + matrix[1, 1] + matrix[2, 2]
+            Quaternion(matrix[1, 2] - matrix[2, 1], matrix[2, 0] - matrix[0, 2], matrix[0, 1] - matrix[1, 0], t)
+        }
+    }
+    return QuaternionField.invoke { q * (0.5 / sqrt(t)) }
+}

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/RotationTest.kt

@@ -0,0 +1,35 @@
+/*
+ * Copyright 2018-2021 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.geometry
+
+import space.kscience.kmath.complex.Quaternion
+import space.kscience.kmath.testutils.assertBufferEquals
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class RotationTest {
+
+    @Test
+    fun rotations() = with(Euclidean3DSpace) {
+        val vector = Vector3D(1.0, 1.0, 1.0)
+        val q = Quaternion(1.0, 2.0, -3.0, 4.0)
+        val rotatedByQ = rotate(vector, q)
+        val matrix = q.toRotationMatrix()
+        val rotatedByM = rotate(vector,matrix)
+
+        assertBufferEquals(rotatedByQ, rotatedByM, 1e-4)
+    }
+
+    @Test
+    fun rotationConversion() {
+
+        val q = Quaternion(1.0, 2.0, -3.0, 4.0)
+
+        val matrix = q.toRotationMatrix()
+
+        assertEquals(q, Quaternion.fromRotationMatrix(matrix))
+    }
+}

kmath-tensorflow/src/main/kotlin/space/kscience/kmath/tensorflow/DoubleTensorFlowAlgebra.kt

@@ -7,6 +7,7 @@ import org.tensorflow.op.core.Constant
 import org.tensorflow.types.TFloat64
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.misc.PerformancePitfall
+import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.DefaultStrides
 import space.kscience.kmath.nd.Shape
 import space.kscience.kmath.nd.StructureND
@@ -74,6 +75,7 @@ public class DoubleTensorFlowAlgebra internal constructor(
  *
  * The resulting tensor is available outside of scope
  */
+@UnstableKMathAPI
 public fun DoubleField.produceWithTF(
     block: DoubleTensorFlowAlgebra.() -> StructureND<Double>,
 ): StructureND<Double> = Graph().use { graph ->

settings.gradle.kts

@@ -20,6 +20,7 @@ dependencyResolutionManagement {
 }
 
 include(
+    ":test-utils",
     ":kmath-memory",
     ":kmath-complex",
     ":kmath-core",

test-utils/build.gradle.kts

@@ -0,0 +1,13 @@
+plugins {
+    id("ru.mipt.npm.gradle.mpp")
+    id("ru.mipt.npm.gradle.native")
+}
+
+kotlin.sourceSets {
+    commonMain {
+        dependencies {
+            api(projects.kmath.kmathCore)
+            api(kotlin("test"))
+        }
+    }
+}

test-utils/src/commonMain/kotlin/FieldVerifier.kt

@@ -10,7 +10,7 @@ import space.kscience.kmath.operations.invoke
 import kotlin.test.assertEquals
 import kotlin.test.assertNotEquals
 
-internal class FieldVerifier<T, out A : Field<T>>(
+public class FieldVerifier<T, out A : Field<T>>(
     algebra: A, a: T, b: T, c: T, x: Number,
 ) : RingVerifier<T, A>(algebra, a, b, c, x) {
 

test-utils/src/commonMain/kotlin/RingVerifier.kt

@@ -10,7 +10,7 @@ import space.kscience.kmath.operations.ScaleOperations
 import space.kscience.kmath.operations.invoke
 import kotlin.test.assertEquals
 
-internal open class RingVerifier<T, out A>(algebra: A, a: T, b: T, c: T, x: Number) :
+public open class RingVerifier<T, out A>(algebra: A, a: T, b: T, c: T, x: Number) :
     SpaceVerifier<T, A>(algebra, a, b, c, x) where A : Ring<T>, A : ScaleOperations<T> {
 
     override fun verify() {

test-utils/src/commonMain/kotlin/SpaceVerifier.kt

@@ -11,12 +11,12 @@ import space.kscience.kmath.operations.invoke
 import kotlin.test.assertEquals
 import kotlin.test.assertNotEquals
 
-internal open class SpaceVerifier<T, out S>(
+public open class SpaceVerifier<T, out S>(
     override val algebra: S,
-    val a: T,
+    public val a: T,
-    val b: T,
+    public val b: T,
-    val c: T,
+    public val c: T,
-    val x: Number,
+    public val x: Number,
 ) : AlgebraicVerifier<T, Ring<T>> where S : Ring<T>, S : ScaleOperations<T> {
     override fun verify() {
         algebra {

test-utils/src/commonMain/kotlin/asserts.kt

@@ -0,0 +1,20 @@
+/*
+ * Copyright 2018-2021 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.testutils
+
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.indices
+import kotlin.test.assertEquals
+import kotlin.test.fail
+
+public fun assertBufferEquals(expected: Buffer<Double>, result: Buffer<Double>, tolerance: Double = 1e-4) {
+    if (expected.size != result.size) {
+        fail("Expected size is ${expected.size}, but the result size is ${result.size}")
+    }
+    expected.indices.forEach {
+        assertEquals(expected[it], result[it], tolerance)
+    }
+}