fft and deprecation

README.md

@@ -10,7 +10,7 @@ repositories {
 }
 
 dependencies {
-    implementation("com.martmists.ndarray-simd:ndarray-simd:1.2.3")
+    implementation("com.martmists.ndarray-simd:ndarray-simd:1.3.0")
 }
 ```
 

build.gradle.kts

@@ -14,7 +14,7 @@ plugins {
 }
 
 group = "com.martmists.ndarray-simd"
-version = "1.2.3"
+version = "1.3.0"
 val isProduction = (findProperty("production") ?: System.getProperty("production")) != null
 
 repositories {

src/commonMain/kotlin/com/martmists/ndarray/simd/F64Array.kt

@@ -1381,39 +1381,6 @@ interface F64Array {
      */
     fun hypot(other: F64Array): F64Array = copy().apply { hypotInPlace(other) }
 
-    /**
-     * Returns the diagonal of the array.
-     *
-     * @return the diagonal
-     */
-    @Deprecated("Will be moved to F64TwoAxisArray in the future")
-    fun diagonal(): F64FlatArray = unsupported()
-
-    /**
-     * Returns the determination of this matrix.
-     *
-     * @return the determinant
-     * @since 1.1.1
-     */
-    @Deprecated("Will be moved to F64TwoAxisArray in the future")
-    fun determinant(): Double = unsupported()
-
-    /**
-     * Returns the inverse matrix.
-     * Note: for MxN matrices, it assumes a full-rank matrix.
-     *
-     * @return the inverse matrix
-     * @since 1.1.1
-     */
-    @Deprecated("Will be moved to F64TwoAxisArray in the future")
-    fun inverse(): F64TwoAxisArray = unsupported()
-
-    /**
-     * Computes the matrix multiplication of this array with another array.
-     */
-    @Deprecated("Will be moved to F64TwoAxisArray in the future")
-    infix fun matmul(other: F64Array): F64TwoAxisArray = unsupported()
-
     /**
      * Returns the data from the array as a flat array.
      */
@@ -1698,6 +1665,40 @@ interface F64Array {
             return result
         }
 
+        /**
+         * Concatenates the given arrays along the specified axis.
+         *
+         * @param first the first array
+         * @param rest the rest of the arrays
+         * @param axis the axis to concatenate along
+         * @return the concatenated array
+         * @since 1.3.0
+         */
+        @JvmStatic
+        @JvmOverloads
+        fun concat(first: F64TwoAxisArray, vararg rest: F64TwoAxisArray, axis: Int = 0): F64TwoAxisArray {
+            for (other in rest) {
+                require(other.shape.remove(axis).contentEquals(first.shape.remove(axis))) {
+                    "input array shapes must be exactly equal for all dimensions except $axis"
+                }
+            }
+
+            val shape = first.shape.copyOf().apply {
+                this[axis] = first.shape[axis] + rest.sumOf { it.shape[axis] }
+            }
+
+            val result = invoke(shape[0], shape[1])
+            var offset = 0
+            for (a in arrayOf(first, *rest)) {
+                if (a.length > 0) {
+                    a.copyTo(result.slice(offset, offset + a.shape[axis], axis = axis))
+                    offset += a.shape[axis]
+                }
+            }
+
+            return result
+        }
+
         /**
          * Creates an array of the given shape with random values.
          * The values are uniformly distributed between 0 (inclusive) and 1 (exclusive).

src/commonMain/kotlin/com/martmists/ndarray/simd/F64FlatArray.kt

@@ -410,11 +410,13 @@ interface F64FlatArray : F64Array {
     override fun hypot(other: F64Array): F64FlatArray = copy().apply { hypotInPlace(other) }
 
     /**
-     * Unsupported on flat arrays.
+     * Performs a Real-to-Complex FFT on this array.
+     * The output array is NOT normalized.
      *
-     * @see F64Array.diagonal
+     * @return The FFT result, where the second axis is [real, imag]
+     * @since 1.3.0
      */
-    override fun diagonal(): F64FlatArray = unsupported()
+    fun fftR2C(): F64TwoAxisArray = F64Array(shape[0], 2) { i, x -> if (x == 1) 0.0 else this[i] }.apply { fftC2CInPlace() }
 
     companion object {
         /**

src/commonMain/kotlin/com/martmists/ndarray/simd/F64TwoAxisArray.kt

@@ -1,6 +1,7 @@
 package com.martmists.ndarray.simd
 
 import com.martmists.ndarray.simd.impl.create
+import com.martmists.ndarray.simd.impl.unsupported
 
 /**
  * A 2D specialization type for [F64Array].
@@ -342,6 +343,36 @@ interface F64TwoAxisArray : F64Array {
      */
     override fun hypot(other: F64Array): F64TwoAxisArray = copy().apply { hypotInPlace(other) }
 
+    /**
+     * Returns the diagonal of the array.
+     *
+     * @return the diagonal
+     */
+    fun diagonal(): F64FlatArray
+
+    /**
+     * Returns the determination of this matrix.
+     *
+     * @return the determinant
+     * @since 1.1.1
+     */
+    fun determinant(): Double
+
+    /**
+     * Returns the inverse matrix.
+     * Note: for MxN matrices, it assumes a full-rank matrix.
+     *
+     * @return the inverse matrix
+     * @since 1.1.1
+     */
+    fun inverse(): F64TwoAxisArray
+
+    /**
+     * Computes the matrix multiplication of this array with another array.
+     * This requires `this.shape[1] == other.shape[0]`
+     */
+    infix fun matmul(other: F64TwoAxisArray): F64TwoAxisArray
+
     /**
      * Calculates the eigenvalues and eigenvectors of this [F64TwoAxisArray] using the QR algorithm.
      *
@@ -350,6 +381,25 @@ interface F64TwoAxisArray : F64Array {
      */
     fun eigen(): Pair<F64FlatArray, F64TwoAxisArray>
 
+    /**
+     * Performs a Complex-to-Complex FFT on this array in-place.
+     * The second axis has to be [real, imag].
+     * The resulting array is not normalized.
+     *
+     * @since 1.3.0
+     */
+    fun fftC2CInPlace()
+
+    /**
+     * Performs a Complex-to-Complex FFT on this array.
+     * The second axis has to be [real, imag].
+     * The resulting array is not normalized.
+     *
+     * @return the FFT result
+     * @since 1.3.0
+     */
+    fun fftC2C() = copy().apply { fftC2CInPlace() }
+
     companion object {
         /**
          * Creates a [F64TwoAxisArray] from a [DoubleArray].

src/commonMain/kotlin/com/martmists/ndarray/simd/compat.kt

@@ -18,7 +18,7 @@ operator fun Double.times(arr: F64Array): F64Array = arr.times(this)
 /**
  * @see F64Array.div
  */
-operator fun Double.div(arr: F64Array): F64Array = arr.div(this)
+operator fun Double.div(arr: F64Array): F64Array = arr.div(1 / this)
 
 /**
  * @see F64Array.expBase

src/commonMain/kotlin/com/martmists/ndarray/simd/impl/F64ArrayImpl.kt

@@ -135,8 +135,8 @@ internal open class F64ArrayImpl internal constructor(
     }
 
     override fun slice(from: Int, to: Int, step: Int, axis: Int): F64Array {
+        checkAxis(axis)
         require(step > 0) { "slicing step must be positive, but was $step" }
-        require(axis in 0 until nDim) { "axis out of bounds: $axis" }
         require(from >= 0) { "slicing start index must be positive, but was $from" }
         val actualTo = if (to != -1) {
             require(to > from) { "slicing end index $to must be greater than start index $from" }
@@ -219,7 +219,7 @@ internal open class F64ArrayImpl internal constructor(
     }
 
     override fun reduce(axis: Int, operation: (Double, Double) -> Double): F64Array {
-        check(axis in 0 until nDim) { "axis out of bounds: $axis" }
+        checkAxis(axis)
 
         val items = along(axis)
         val base = items.first().copy()
@@ -247,7 +247,7 @@ internal open class F64ArrayImpl internal constructor(
     }
 
     override fun scan(axis: Int, operation: (Double, Double) -> Double) {
-        check(axis in 0 until nDim) { "axis out of bounds: $axis" }
+        checkAxis(axis)
         val items = along(axis)
         var tmp = items.first().copy()
         for (next in items) {
@@ -520,140 +520,6 @@ internal open class F64ArrayImpl internal constructor(
         commonUnrollToFlat(other, F64FlatArray::hypotInPlace)
     }
 
-    override fun matmul(other: F64Array): F64TwoAxisArray {
-        // FIXME: If both inputs are flattenable and above large dense threshold, use NativeSpeedup impl
-        check(nDim == 2) { "matmul is only supported for 2D arrays" }
-        check(other.nDim == 2) { "matmul is only supported for 2D arrays" }
-        check(shape[1] == other.shape[0]) {
-            "matmul dimensions do not match: ${shape[1]} != ${other.shape[0]}"
-        }
-        val result = F64Array.full(shape[0], other.shape[1], init=0.0)
-        for (i in 0 until shape[0]) {
-            for (j in 0 until other.shape[1]) {
-                for (k in 0 until shape[1]) {
-                    result[i, j] += this[i, k] * other[k, j]
-                }
-            }
-        }
-        return result
-    }
-
-    override fun diagonal(): F64FlatArray {
-        check(nDim == 2) { "Only 2D arrays are supported" }
-        check(shape[0] == shape[1]) { "Only square matrices are supported" }
-
-        return F64FlatArray.create(data, offset, strides[0] + strides[1], shape[0])
-    }
-
-    internal fun minor(row: Int, col: Int): F64Array {
-        val withoutRow = when (row) {
-            0 -> slice(1, shape[0])
-            shape[0] - 1 -> slice(0, row)
-            else -> {
-                val top = slice(0, row)
-                val bottom = slice(row + 1, shape[0])
-                F64Array.concat(top, bottom)
-            }
-        }
-        val withoutCol = when (col) {
-            0 -> withoutRow.slice(1, shape[1], axis = 1)
-            shape[1] - 1 -> withoutRow.slice(0, shape[1] - 1, axis = 1)
-            else -> {
-                val left = withoutRow.slice(0, col, axis = 1)
-                val right = withoutRow.slice(col + 1, shape[1], axis = 1)
-                F64Array.concat(left, right, axis = 1)
-            }
-        }
-        return withoutCol
-    }
-
-    override fun determinant(): Double {
-        check(shape.size == 2) { "Only 2D matrices are supported" }
-        check(shape[0] == shape[1]) { "Only square matrices are supported" }
-        check(shape[0] > 1) { "Matrix must be at least 2x2" }
-
-        // TODO: Consider cleaning up with an F64Array
-        val n = shape[0]
-        val lu = Array(n) { i -> DoubleArray(n) { this[i, it] } }
-        var det = 1.0
-
-        // LU Decomposition with partial pivoting
-        for (k in 0 until n) {
-            var max = 0.0
-            var maxRow = k
-            for (i in k until n) {
-                val abs = abs(lu[i][k])
-                if (abs > max) {
-                    max = abs
-                    maxRow = i
-                }
-            }
-
-            if (k != maxRow) {
-                val tmp = lu[k]
-                lu[k] = lu[maxRow]
-                lu[maxRow] = tmp
-
-                det *= -1
-            }
-
-            det *= lu[k][k]
-            if (lu[k][k] == 0.0) return 0.0
-
-            for (i in k + 1 until n) {
-                lu[i][k] /= lu[k][k]
-                for (j in k + 1 until n) {
-                    lu[i][j] -= lu[i][k] * lu[k][j]
-                }
-            }
-        }
-
-        return det
-    }
-
-    override fun inverse(): F64TwoAxisArray {
-        check(shape.size == 2) { "Only 2D matrices are supported" }
-
-        return when {
-            shape[0] == shape[1] -> {
-                if (shape[0] == 2) {
-                    val det = this[0, 0] * this[1, 1] - this[0, 1] * this[1, 0]
-                    val inverse = F64Array(shape[0], shape[1])
-                    inverse[0, 0] = this[1, 1]
-                    inverse[0, 1] = -this[0, 1]
-                    inverse[1, 0] = -this[1, 0]
-                    inverse[1, 1] = this[0, 0]
-                    inverse *= (1.0 / det)
-                    return inverse
-                }
-
-                val det = this.determinant()
-                val invDet = 1.0 / det
-                val res = F64Array.zeros(shape[0], shape[1])
-                for (i in 0 until shape[0]) {
-                    for (j in 0 until shape[1]) {
-                        if ((i + j) % 2 == 0) {
-                            val minorMatrix = this.minor(i, j)
-                            res[j, i] = invDet * minorMatrix.determinant()
-                        } else {
-                            val minorMatrix = this.minor(i, j)
-                            res[j, i] = -invDet * minorMatrix.determinant()
-                        }
-                    }
-                }
-                res
-            }
-            shape[0] < shape[1] -> {
-                val t = transpose()
-                return t.matmul(this.matmul(t).inverse())
-            }
-            else -> {
-                val t = transpose()
-                return t.matmul(this).inverse().matmul(t)
-            }
-        }
-    }
-
     override fun toString(
         maxDisplay: Int,
     ): String {

src/commonMain/kotlin/com/martmists/ndarray/simd/impl/F64FlatArrayImpl.kt

@@ -330,12 +330,6 @@ internal open class F64FlatArrayImpl internal constructor(
     override fun atanhInPlace() = transformInPlace(::atanh)
     override fun hypotInPlace(other: F64Array) = zipTransformInPlace(other, ::hypot)
 
-    override fun diagonal(): F64FlatArray = unsupported()
-
-    override fun determinant(): Double = unsupported()
-
-    override fun inverse(): F64TwoAxisArray = unsupported()
-
     override fun toDoubleArray() = DoubleArray(length) { unsafeGet(it) }
 
     private fun Double.format(digits: Int): String = when {