diff --git a/README.md b/README.md index 5678fc530..34761e838 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,12 @@ -Bintray: [ ![Download](https://api.bintray.com/packages/mipt-npm/scientifik/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/scientifik/kmath-core/_latestVersion) - -Bintray-dev: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/scientifik/kmath-core/_latestVersion) - +[![JetBrains Research](https://jb.gg/badges/research.svg)](https://confluence.jetbrains.com/display/ALL/JetBrains+on+GitHub) [![DOI](https://zenodo.org/badge/129486382.svg)](https://zenodo.org/badge/latestdoi/129486382) +![Gradle build](https://github.com/mipt-npm/kmath/workflows/Gradle%20build/badge.svg) + +Bintray: [ ![Download](https://api.bintray.com/packages/mipt-npm/scientifik/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/scientifik/kmath-core/_latestVersion) + +Bintray-dev: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/scientifik/kmath-core/_latestVersion) + # KMath Could be pronounced as `key-math`. The Kotlin MATHematics library is intended as a Kotlin-based analog to Python's `numpy` library. In contrast to `numpy` and `scipy` it is modular and has a lightweight core. @@ -40,6 +43,8 @@ can be used for a wide variety of purposes from high performance calculations to * **Streaming** Streaming operations on mathematical objects and objects buffers. +* **Type-safe dimensions** Type-safe dimensions for matrix operations. + * **Commons-math wrapper** It is planned to gradually wrap most parts of [Apache commons-math](http://commons.apache.org/proper/commons-math/) library in Kotlin code and maybe rewrite some parts to better suit the Kotlin programming paradigm, however there is no fixed roadmap for that. Feel free to submit a feature request if you want something to be done first. diff --git a/build.gradle.kts b/build.gradle.kts index a63966dfe..6ab33d31c 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -1,8 +1,8 @@ plugins { - id("scientifik.publish") version "0.2.6" apply false + id("scientifik.publish") version "0.4.2" apply false } -val kmathVersion by extra("0.1.4-dev-1") +val kmathVersion by extra("0.1.4-dev-4") val bintrayRepo by extra("scientifik") val githubProject by extra("kmath") diff --git a/examples/build.gradle.kts b/examples/build.gradle.kts index ad59afc5c..8853b78a5 100644 --- a/examples/build.gradle.kts +++ b/examples/build.gradle.kts @@ -4,8 +4,8 @@ import org.jetbrains.kotlin.gradle.tasks.KotlinCompile plugins { java kotlin("jvm") - kotlin("plugin.allopen") version "1.3.60" - id("kotlinx.benchmark") version "0.2.0-dev-5" + kotlin("plugin.allopen") version "1.3.71" + id("kotlinx.benchmark") version "0.2.0-dev-7" } configure { @@ -13,10 +13,9 @@ configure { } repositories { - maven("https://dl.bintray.com/kotlin/kotlin-eap") maven("http://dl.bintray.com/kyonifer/maven") - maven ("https://dl.bintray.com/orangy/maven") - + maven("https://dl.bintray.com/mipt-npm/scientifik") + maven("https://dl.bintray.com/mipt-npm/dev") mavenCentral() } @@ -29,12 +28,11 @@ dependencies { implementation(project(":kmath-coroutines")) implementation(project(":kmath-commons")) implementation(project(":kmath-koma")) + implementation(project(":kmath-viktor")) + implementation(project(":kmath-dimensions")) implementation("com.kyonifer:koma-core-ejml:0.12") - implementation("org.jetbrains.kotlinx:kotlinx-io-jvm:${Scientifik.ioVersion}") - - implementation("org.jetbrains.kotlinx:kotlinx.benchmark.runtime:0.2.0-dev-2") - - + implementation("org.jetbrains.kotlinx:kotlinx-io-jvm:0.2.0-npm-dev-6") + implementation("org.jetbrains.kotlinx:kotlinx.benchmark.runtime:0.2.0-dev-7") "benchmarksCompile"(sourceSets.main.get().compileClasspath) } diff --git a/examples/src/benchmarks/kotlin/scientifik/kmath/structures/ViktorBenchmark.kt b/examples/src/benchmarks/kotlin/scientifik/kmath/structures/ViktorBenchmark.kt new file mode 100644 index 000000000..be4115d81 --- /dev/null +++ b/examples/src/benchmarks/kotlin/scientifik/kmath/structures/ViktorBenchmark.kt @@ -0,0 +1,71 @@ +package scientifik.kmath.structures + +import org.jetbrains.bio.viktor.F64Array +import org.openjdk.jmh.annotations.Benchmark +import org.openjdk.jmh.annotations.Scope +import org.openjdk.jmh.annotations.State +import scientifik.kmath.operations.RealField +import scientifik.kmath.viktor.ViktorNDField + + +@State(Scope.Benchmark) +class ViktorBenchmark { + final val dim = 1000 + final val n = 100 + + // automatically build context most suited for given type. + final val autoField = NDField.auto(RealField, dim, dim) + final val realField = NDField.real(dim, dim) + + final val viktorField = ViktorNDField(intArrayOf(dim, dim)) + + @Benchmark + fun `Automatic field addition`() { + autoField.run { + var res = one + repeat(n) { + res += 1.0 + } + } + } + + @Benchmark + fun `Viktor field addition`() { + viktorField.run { + var res = one + repeat(n) { + res += one + } + } + } + + @Benchmark + fun `Raw Viktor`() { + val one = F64Array.full(init = 1.0, shape = *intArrayOf(dim, dim)) + var res = one + repeat(n) { + res = res + one + } + } + + @Benchmark + fun `Real field log`() { + realField.run { + val fortyTwo = produce { 42.0 } + var res = one + + repeat(n) { + res = ln(fortyTwo) + } + } + } + + @Benchmark + fun `Raw Viktor log`() { + val fortyTwo = F64Array.full(dim, dim, init = 42.0) + var res: F64Array + repeat(n) { + res = fortyTwo.log() + } + } +} \ No newline at end of file diff --git a/examples/src/benchmarks/kotlin/scientifik/kmath/utils/utils.kt b/examples/src/benchmarks/kotlin/scientifik/kmath/utils/utils.kt new file mode 100644 index 000000000..6ec9e9c17 --- /dev/null +++ b/examples/src/benchmarks/kotlin/scientifik/kmath/utils/utils.kt @@ -0,0 +1,8 @@ +package scientifik.kmath.utils + +import kotlin.system.measureTimeMillis + +internal inline fun measureAndPrint(title: String, block: () -> Unit) { + val time = measureTimeMillis(block) + println("$title completed in $time millis") +} \ No newline at end of file diff --git a/examples/src/main/kotlin/scientifik/kmath/structures/NDField.kt b/examples/src/main/kotlin/scientifik/kmath/structures/NDField.kt index 7e8c433c0..cfd1206ff 100644 --- a/examples/src/main/kotlin/scientifik/kmath/structures/NDField.kt +++ b/examples/src/main/kotlin/scientifik/kmath/structures/NDField.kt @@ -4,7 +4,13 @@ import kotlinx.coroutines.GlobalScope import scientifik.kmath.operations.RealField import kotlin.system.measureTimeMillis -fun main(args: Array) { +internal inline fun measureAndPrint(title: String, block: () -> Unit) { + val time = measureTimeMillis(block) + println("$title completed in $time millis") +} + + +fun main() { val dim = 1000 val n = 1000 @@ -15,8 +21,7 @@ fun main(args: Array) { //A generic boxing field. It should be used for objects, not primitives. val genericField = NDField.boxing(RealField, dim, dim) - - val autoTime = measureTimeMillis { + measureAndPrint("Automatic field addition") { autoField.run { var res = one repeat(n) { @@ -25,18 +30,14 @@ fun main(args: Array) { } } - println("Automatic field addition completed in $autoTime millis") - - val elementTime = measureTimeMillis { + measureAndPrint("Element addition"){ var res = genericField.one repeat(n) { res += 1.0 } } - println("Element addition completed in $elementTime millis") - - val specializedTime = measureTimeMillis { + measureAndPrint("Specialized addition") { specializedField.run { var res: NDBuffer = one repeat(n) { @@ -45,10 +46,7 @@ fun main(args: Array) { } } - println("Specialized addition completed in $specializedTime millis") - - - val lazyTime = measureTimeMillis { + measureAndPrint("Lazy addition") { val res = specializedField.one.mapAsync(GlobalScope) { var c = 0.0 repeat(n) { @@ -60,9 +58,7 @@ fun main(args: Array) { res.elements().forEach { it.second } } - println("Lazy addition completed in $lazyTime millis") - - val genericTime = measureTimeMillis { + measureAndPrint("Generic addition") { //genericField.run(action) genericField.run { var res: NDBuffer = one @@ -72,6 +68,4 @@ fun main(args: Array) { } } - println("Generic addition completed in $genericTime millis") - } \ No newline at end of file diff --git a/examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt b/examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt new file mode 100644 index 000000000..fdc09ed5d --- /dev/null +++ b/examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt @@ -0,0 +1,35 @@ +package scientifik.kmath.structures + +import scientifik.kmath.dimensions.D2 +import scientifik.kmath.dimensions.D3 +import scientifik.kmath.dimensions.DMatrixContext +import scientifik.kmath.dimensions.Dimension +import scientifik.kmath.operations.RealField + +fun DMatrixContext.simple() { + val m1 = produce { i, j -> (i + j).toDouble() } + val m2 = produce { i, j -> (i + j).toDouble() } + + //Dimension-safe addition + m1.transpose() + m2 +} + + +object D5 : Dimension { + override val dim: UInt = 5u +} + +fun DMatrixContext.custom() { + val m1 = produce { i, j -> (i + j).toDouble() } + val m2 = produce { i, j -> (i - j).toDouble() } + val m3 = produce { i, j -> (i - j).toDouble() } + + (m1 dot m2) + m3 +} + +fun main() { + DMatrixContext.real.run { + simple() + custom() + } +} \ No newline at end of file diff --git a/gradle/wrapper/gradle-wrapper.jar b/gradle/wrapper/gradle-wrapper.jar index cc4fdc293..490fda857 100644 Binary files a/gradle/wrapper/gradle-wrapper.jar and b/gradle/wrapper/gradle-wrapper.jar differ diff --git a/gradle/wrapper/gradle-wrapper.properties b/gradle/wrapper/gradle-wrapper.properties index 6ce793f21..a4b442974 100644 --- a/gradle/wrapper/gradle-wrapper.properties +++ b/gradle/wrapper/gradle-wrapper.properties @@ -1,5 +1,5 @@ distributionBase=GRADLE_USER_HOME distributionPath=wrapper/dists -distributionUrl=https\://services.gradle.org/distributions/gradle-6.0-bin.zip +distributionUrl=https\://services.gradle.org/distributions/gradle-6.3-bin.zip zipStoreBase=GRADLE_USER_HOME zipStorePath=wrapper/dists diff --git a/gradlew.bat b/gradlew.bat index 9618d8d96..62bd9b9cc 100644 --- a/gradlew.bat +++ b/gradlew.bat @@ -29,6 +29,9 @@ if "%DIRNAME%" == "" set DIRNAME=. set APP_BASE_NAME=%~n0 set APP_HOME=%DIRNAME% +@rem Resolve any "." and ".." in APP_HOME to make it shorter. +for %%i in ("%APP_HOME%") do set APP_HOME=%%~fi + @rem Add default JVM options here. You can also use JAVA_OPTS and GRADLE_OPTS to pass JVM options to this script. set DEFAULT_JVM_OPTS="-Xmx64m" "-Xms64m" diff --git a/kmath-commons/build.gradle.kts b/kmath-commons/build.gradle.kts index 9b446f872..5ce1b935a 100644 --- a/kmath-commons/build.gradle.kts +++ b/kmath-commons/build.gradle.kts @@ -8,7 +8,6 @@ dependencies { api(project(":kmath-core")) api(project(":kmath-coroutines")) api(project(":kmath-prob")) + api(project(":kmath-functions")) api("org.apache.commons:commons-math3:3.6.1") - testImplementation("org.jetbrains.kotlin:kotlin-test") - testImplementation("org.jetbrains.kotlin:kotlin-test-junit") } \ No newline at end of file diff --git a/kmath-commons/src/test/kotlin/scientifik/kmath/commons/expressions/AutoDiffTest.kt b/kmath-commons/src/test/kotlin/scientifik/kmath/commons/expressions/AutoDiffTest.kt index 100f49948..c77f30eb7 100644 --- a/kmath-commons/src/test/kotlin/scientifik/kmath/commons/expressions/AutoDiffTest.kt +++ b/kmath-commons/src/test/kotlin/scientifik/kmath/commons/expressions/AutoDiffTest.kt @@ -1,7 +1,7 @@ package scientifik.kmath.commons.expressions -import org.junit.Test import scientifik.kmath.expressions.invoke +import kotlin.test.Test import kotlin.test.assertEquals inline fun diff(order: Int, vararg parameters: Pair, block: DerivativeStructureField.() -> R) = diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/FeaturedMatrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/FeaturedMatrix.kt index 9f6b9f600..297586fc6 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/FeaturedMatrix.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/FeaturedMatrix.kt @@ -41,16 +41,6 @@ fun Structure2D.Companion.square(vararg elements: T): FeaturedMatrix Structure2D.Companion.build(rows: Int, columns: Int): MatrixBuilder = MatrixBuilder(rows, columns) - -class MatrixBuilder(val rows: Int, val columns: Int) { - operator fun invoke(vararg elements: T): FeaturedMatrix { - if (rows * columns != elements.size) error("The number of elements ${elements.size} is not equal $rows * $columns") - val buffer = elements.asBuffer() - return BufferMatrix(rows, columns, buffer) - } -} - val Matrix<*>.features get() = (this as? FeaturedMatrix)?.features?: emptySet() /** diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgebra.kt similarity index 100% rename from kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt rename to kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgebra.kt diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixBuilder.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixBuilder.kt new file mode 100644 index 000000000..466dbea6e --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixBuilder.kt @@ -0,0 +1,14 @@ +package scientifik.kmath.linear + +import scientifik.kmath.structures.Structure2D +import scientifik.kmath.structures.asBuffer + +class MatrixBuilder(val rows: Int, val columns: Int) { + operator fun invoke(vararg elements: T): FeaturedMatrix { + if (rows * columns != elements.size) error("The number of elements ${elements.size} is not equal $rows * $columns") + val buffer = elements.asBuffer() + return BufferMatrix(rows, columns, buffer) + } +} + +fun Structure2D.Companion.build(rows: Int, columns: Int): MatrixBuilder = MatrixBuilder(rows, columns) \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/functions.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/functions.kt deleted file mode 100644 index cef8209ed..000000000 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/functions.kt +++ /dev/null @@ -1,37 +0,0 @@ -package scientifik.kmath.coroutines - -import scientifik.kmath.operations.RealField -import scientifik.kmath.operations.SpaceOperations -import kotlin.jvm.JvmName - -/** - * A suspendable univariate function defined in algebraic context - */ -interface UFunction> { - suspend operator fun C.invoke(arg: T): T -} - -suspend fun UFunction.invoke(arg: Double) = RealField.invoke(arg) - -/** - * A suspendable multivariate (N->1) function defined on algebraic context - */ -interface MFunction> { - /** - * The input dimension of the function - */ - val dimension: UInt - - suspend operator fun C.invoke(vararg args: T): T -} - -suspend fun MFunction.invoke(args: DoubleArray) = RealField.invoke(*args.toTypedArray()) -@JvmName("varargInvoke") -suspend fun MFunction.invoke(vararg args: Double) = RealField.invoke(*args.toTypedArray()) - -/** - * A suspendable univariate function with parameter - */ -interface ParametricUFunction> { - suspend operator fun C.invoke(arg: T, parameter: P): T -} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt index 0ed769db8..485185526 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt @@ -6,14 +6,14 @@ annotation class KMathContext /** * Marker interface for any algebra */ -interface Algebra +interface Algebra -inline operator fun T.invoke(block: T.() -> R): R = run(block) +inline operator fun , R> T.invoke(block: T.() -> R): R = run(block) /** * Space-like operations without neutral element */ -interface SpaceOperations : Algebra { +interface SpaceOperations : Algebra { /** * Addition operation for two context elements */ diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt index 4e8dbf36d..bfb4199a3 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt @@ -1,4 +1,15 @@ package scientifik.kmath.operations -fun Space.sum(data : Iterable): T = data.fold(zero) { left, right -> add(left,right) } -fun Space.sum(data : Sequence): T = data.fold(zero) { left, right -> add(left, right) } \ No newline at end of file +fun Space.sum(data: Iterable): T = data.fold(zero) { left, right -> add(left, right) } +fun Space.sum(data: Sequence): T = data.fold(zero) { left, right -> add(left, right) } + +fun > Iterable.sumWith(space: S): T = space.sum(this) + +//TODO optimized power operation +fun RingOperations.power(arg: T, power: Int): T { + var res = arg + repeat(power - 1) { + res *= arg + } + return res +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/BigInt.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/BigInt.kt new file mode 100644 index 000000000..1661170d3 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/BigInt.kt @@ -0,0 +1,484 @@ +package scientifik.kmath.operations + +import scientifik.kmath.operations.BigInt.Companion.BASE +import scientifik.kmath.operations.BigInt.Companion.BASE_SIZE +import kotlin.math.log2 +import kotlin.math.max +import kotlin.math.min +import kotlin.math.sign + + +typealias Magnitude = UIntArray +typealias TBase = ULong + +/** + * Kotlin Multiplatform implementation of Big Integer numbers (KBigInteger). + * + * @author Robert Drynkin (https://github.com/robdrynkin) and Peter Klimai (https://github.com/pklimai) + */ +object BigIntField : Field { + override val zero: BigInt = BigInt.ZERO + override val one: BigInt = BigInt.ONE + + override fun add(a: BigInt, b: BigInt): BigInt = a.plus(b) + + override fun multiply(a: BigInt, k: Number): BigInt = a.times(k.toLong()) + + override fun multiply(a: BigInt, b: BigInt): BigInt = a.times(b) + + operator fun String.unaryPlus(): BigInt = this.parseBigInteger() ?: error("Can't parse $this as big integer") + + operator fun String.unaryMinus(): BigInt = + -(this.parseBigInteger() ?: error("Can't parse $this as big integer")) + + override fun divide(a: BigInt, b: BigInt): BigInt = a.div(b) +} + +class BigInt internal constructor( + private val sign: Byte, + private val magnitude: Magnitude +) : Comparable { + + override fun compareTo(other: BigInt): Int { + return when { + (this.sign == 0.toByte()) and (other.sign == 0.toByte()) -> 0 + this.sign < other.sign -> -1 + this.sign > other.sign -> 1 + else -> this.sign * compareMagnitudes(this.magnitude, other.magnitude) + } + } + + override fun equals(other: Any?): Boolean { + if (other is BigInt) { + return this.compareTo(other) == 0 + } else error("Can't compare KBigInteger to a different type") + } + + override fun hashCode(): Int { + return magnitude.hashCode() + this.sign + } + + fun abs(): BigInt = if (sign == 0.toByte()) this else BigInt(1, magnitude) + + operator fun unaryMinus(): BigInt { + return if (this.sign == 0.toByte()) this else BigInt((-this.sign).toByte(), this.magnitude) + } + + operator fun plus(b: BigInt): BigInt { + return when { + b.sign == 0.toByte() -> this + this.sign == 0.toByte() -> b + this == -b -> ZERO + this.sign == b.sign -> BigInt(this.sign, addMagnitudes(this.magnitude, b.magnitude)) + else -> { + val comp: Int = compareMagnitudes(this.magnitude, b.magnitude) + + if (comp == 1) { + BigInt(this.sign, subtractMagnitudes(this.magnitude, b.magnitude)) + } else { + BigInt((-this.sign).toByte(), subtractMagnitudes(b.magnitude, this.magnitude)) + } + } + } + } + + operator fun minus(b: BigInt): BigInt { + return this + (-b) + } + + operator fun times(b: BigInt): BigInt { + return when { + this.sign == 0.toByte() -> ZERO + b.sign == 0.toByte() -> ZERO +// TODO: Karatsuba + else -> BigInt((this.sign * b.sign).toByte(), multiplyMagnitudes(this.magnitude, b.magnitude)) + } + } + + operator fun times(other: UInt): BigInt { + return when { + this.sign == 0.toByte() -> ZERO + other == 0U -> ZERO + else -> BigInt(this.sign, multiplyMagnitudeByUInt(this.magnitude, other)) + } + } + + operator fun times(other: Int): BigInt { + return if (other > 0) + this * kotlin.math.abs(other).toUInt() + else + -this * kotlin.math.abs(other).toUInt() + } + + operator fun div(other: UInt): BigInt { + return BigInt(this.sign, divideMagnitudeByUInt(this.magnitude, other)) + } + + operator fun div(other: Int): BigInt { + return BigInt( + (this.sign * other.sign).toByte(), + divideMagnitudeByUInt(this.magnitude, kotlin.math.abs(other).toUInt()) + ) + } + + private fun division(other: BigInt): Pair { + // Long division algorithm: + // https://en.wikipedia.org/wiki/Division_algorithm#Integer_division_(unsigned)_with_remainder + // TODO: Implement more effective algorithm + var q: BigInt = ZERO + var r: BigInt = ZERO + + val bitSize = + (BASE_SIZE * (this.magnitude.size - 1) + log2(this.magnitude.lastOrNull()?.toFloat() ?: 0f + 1)).toInt() + for (i in bitSize downTo 0) { + r = r shl 1 + r = r or ((abs(this) shr i) and ONE) + if (r >= abs(other)) { + r -= abs(other) + q += (ONE shl i) + } + } + + return Pair(BigInt((this.sign * other.sign).toByte(), q.magnitude), r) + } + + operator fun div(other: BigInt): BigInt { + return this.division(other).first + } + + infix fun shl(i: Int): BigInt { + if (this == ZERO) return ZERO + if (i == 0) return this + + val fullShifts = i / BASE_SIZE + 1 + val relShift = i % BASE_SIZE + val shiftLeft = { x: UInt -> if (relShift >= 32) 0U else x shl relShift } + val shiftRight = { x: UInt -> if (BASE_SIZE - relShift >= 32) 0U else x shr (BASE_SIZE - relShift) } + + val newMagnitude: Magnitude = Magnitude(this.magnitude.size + fullShifts) + + for (j in this.magnitude.indices) { + newMagnitude[j + fullShifts - 1] = shiftLeft(this.magnitude[j]) + if (j != 0) { + newMagnitude[j + fullShifts - 1] = newMagnitude[j + fullShifts - 1] or shiftRight(this.magnitude[j - 1]) + } + } + + newMagnitude[this.magnitude.size + fullShifts - 1] = shiftRight(this.magnitude.last()) + + return BigInt(this.sign, stripLeadingZeros(newMagnitude)) + } + + infix fun shr(i: Int): BigInt { + if (this == ZERO) return ZERO + if (i == 0) return this + + val fullShifts = i / BASE_SIZE + val relShift = i % BASE_SIZE + val shiftRight = { x: UInt -> if (relShift >= 32) 0U else x shr relShift } + val shiftLeft = { x: UInt -> if (BASE_SIZE - relShift >= 32) 0U else x shl (BASE_SIZE - relShift) } + if (this.magnitude.size - fullShifts <= 0) { + return ZERO + } + val newMagnitude: Magnitude = Magnitude(this.magnitude.size - fullShifts) + + for (j in fullShifts until this.magnitude.size) { + newMagnitude[j - fullShifts] = shiftRight(this.magnitude[j]) + if (j != this.magnitude.size - 1) { + newMagnitude[j - fullShifts] = newMagnitude[j - fullShifts] or shiftLeft(this.magnitude[j + 1]) + } + } + + return BigInt(this.sign, stripLeadingZeros(newMagnitude)) + } + + infix fun or(other: BigInt): BigInt { + if (this == ZERO) return other; + if (other == ZERO) return this; + val resSize = max(this.magnitude.size, other.magnitude.size) + val newMagnitude: Magnitude = Magnitude(resSize) + for (i in 0 until resSize) { + if (i < this.magnitude.size) { + newMagnitude[i] = newMagnitude[i] or this.magnitude[i] + } + if (i < other.magnitude.size) { + newMagnitude[i] = newMagnitude[i] or other.magnitude[i] + } + } + return BigInt(1, stripLeadingZeros(newMagnitude)) + } + + infix fun and(other: BigInt): BigInt { + if ((this == ZERO) or (other == ZERO)) return ZERO; + val resSize = min(this.magnitude.size, other.magnitude.size) + val newMagnitude: Magnitude = Magnitude(resSize) + for (i in 0 until resSize) { + newMagnitude[i] = this.magnitude[i] and other.magnitude[i] + } + return BigInt(1, stripLeadingZeros(newMagnitude)) + } + + operator fun rem(other: Int): Int { + val res = this - (this / other) * other + return if (res == ZERO) 0 else res.sign * res.magnitude[0].toInt() + } + + operator fun rem(other: BigInt): BigInt { + return this - (this / other) * other + } + + fun modPow(exponent: BigInt, m: BigInt): BigInt { + return when { + exponent == ZERO -> ONE + exponent % 2 == 1 -> (this * modPow(exponent - ONE, m)) % m + else -> { + val sqRoot = modPow(exponent / 2, m) + (sqRoot * sqRoot) % m + } + } + } + + override fun toString(): String { + if (this.sign == 0.toByte()) { + return "0x0" + } + var res: String = if (this.sign == (-1).toByte()) "-0x" else "0x" + var numberStarted = false + + for (i in this.magnitude.size - 1 downTo 0) { + for (j in BASE_SIZE / 4 - 1 downTo 0) { + val curByte = (this.magnitude[i] shr 4 * j) and 0xfU + if (numberStarted or (curByte != 0U)) { + numberStarted = true + res += hexMapping[curByte] + } + } + } + + return res + } + + companion object { + const val BASE = 0xffffffffUL + const val BASE_SIZE: Int = 32 + val ZERO: BigInt = BigInt(0, uintArrayOf()) + val ONE: BigInt = BigInt(1, uintArrayOf(1u)) + + private val hexMapping: HashMap = hashMapOf( + 0U to "0", 1U to "1", 2U to "2", 3U to "3", + 4U to "4", 5U to "5", 6U to "6", 7U to "7", + 8U to "8", 9U to "9", 10U to "a", 11U to "b", + 12U to "c", 13U to "d", 14U to "e", 15U to "f" + ) + + private fun compareMagnitudes(mag1: Magnitude, mag2: Magnitude): Int { + when { + mag1.size > mag2.size -> return 1 + mag1.size < mag2.size -> return -1 + else -> { + for (i in mag1.size - 1 downTo 0) { + if (mag1[i] > mag2[i]) { + return 1 + } else if (mag1[i] < mag2[i]) { + return -1 + } + } + return 0 + } + } + } + + private fun addMagnitudes(mag1: Magnitude, mag2: Magnitude): Magnitude { + val resultLength: Int = max(mag1.size, mag2.size) + 1 + val result = Magnitude(resultLength) + var carry: TBase = 0UL + + for (i in 0 until resultLength - 1) { + val res = when { + i >= mag1.size -> mag2[i].toULong() + carry + i >= mag2.size -> mag1[i].toULong() + carry + else -> mag1[i].toULong() + mag2[i].toULong() + carry + } + result[i] = (res and BASE).toUInt() + carry = (res shr BASE_SIZE) + } + result[resultLength - 1] = carry.toUInt() + return stripLeadingZeros(result) + } + + private fun subtractMagnitudes(mag1: Magnitude, mag2: Magnitude): Magnitude { + val resultLength: Int = mag1.size + val result = Magnitude(resultLength) + var carry = 0L + + for (i in 0 until resultLength) { + var res: Long = + if (i < mag2.size) mag1[i].toLong() - mag2[i].toLong() - carry + else mag1[i].toLong() - carry + + carry = if (res < 0) 1 else 0 + res += carry * (BASE + 1UL).toLong() + + result[i] = res.toUInt() + } + + return stripLeadingZeros(result) + } + + private fun multiplyMagnitudeByUInt(mag: Magnitude, x: UInt): Magnitude { + val resultLength: Int = mag.size + 1 + val result = Magnitude(resultLength) + var carry: ULong = 0UL + + for (i in mag.indices) { + val cur: ULong = carry + mag[i].toULong() * x.toULong() + result[i] = (cur and BASE.toULong()).toUInt() + carry = cur shr BASE_SIZE + } + result[resultLength - 1] = (carry and BASE).toUInt() + + return stripLeadingZeros(result) + } + + private fun multiplyMagnitudes(mag1: Magnitude, mag2: Magnitude): Magnitude { + val resultLength: Int = mag1.size + mag2.size + val result = Magnitude(resultLength) + + for (i in mag1.indices) { + var carry: ULong = 0UL + for (j in mag2.indices) { + val cur: ULong = result[i + j].toULong() + mag1[i].toULong() * mag2[j].toULong() + carry + result[i + j] = (cur and BASE.toULong()).toUInt() + carry = cur shr BASE_SIZE + } + result[i + mag2.size] = (carry and BASE).toUInt() + } + + return stripLeadingZeros(result) + } + + private fun divideMagnitudeByUInt(mag: Magnitude, x: UInt): Magnitude { + val resultLength: Int = mag.size + val result = Magnitude(resultLength) + var carry: ULong = 0UL + + for (i in mag.size - 1 downTo 0) { + val cur: ULong = mag[i].toULong() + (carry shl BASE_SIZE) + result[i] = (cur / x).toUInt() + carry = cur % x + } + return stripLeadingZeros(result) + } + + } + +} + + +private fun stripLeadingZeros(mag: Magnitude): Magnitude { + if (mag.isEmpty() || mag.last() != 0U) { + return mag + } + var resSize: Int = mag.size - 1 + while (mag[resSize] == 0U) { + if (resSize == 0) + break + resSize -= 1 + } + return mag.sliceArray(IntRange(0, resSize)) +} + +fun abs(x: BigInt): BigInt = x.abs() + +/** + * Convert this [Int] to [BigInt] + */ +fun Int.toBigInt() = BigInt(sign.toByte(), uintArrayOf(kotlin.math.abs(this).toUInt())) + +/** + * Convert this [Long] to [BigInt] + */ +fun Long.toBigInt() = BigInt( + sign.toByte(), stripLeadingZeros( + uintArrayOf( + (kotlin.math.abs(this).toULong() and BASE).toUInt(), + ((kotlin.math.abs(this).toULong() shr BASE_SIZE) and BASE).toUInt() + ) + ) +) + +/** + * Convert UInt to [BigInt] + */ +fun UInt.toBigInt() = BigInt(1, uintArrayOf(this)) + +/** + * Convert ULong to [BigInt] + */ +fun ULong.toBigInt() = BigInt( + 1, + stripLeadingZeros( + uintArrayOf( + (this and BigInt.BASE).toUInt(), + ((this shr BigInt.BASE_SIZE) and BigInt.BASE).toUInt() + ) + ) +) + +/** + * Create a [BigInt] with this array of magnitudes with protective copy + */ +fun UIntArray.toBigInt(sign: Byte): BigInt { + if (sign == 0.toByte() && isNotEmpty()) error("") + return BigInt(sign, this.copyOf()) +} + +val hexChToInt = hashMapOf( + '0' to 0, '1' to 1, '2' to 2, '3' to 3, + '4' to 4, '5' to 5, '6' to 6, '7' to 7, + '8' to 8, '9' to 9, 'A' to 10, 'B' to 11, + 'C' to 12, 'D' to 13, 'E' to 14, 'F' to 15 +) + +/** + * Returns null if a valid number can not be read from a string + */ +fun String.parseBigInteger(): BigInt? { + val sign: Int + val sPositive: String + when { + this[0] == '+' -> { + sign = +1 + sPositive = this.substring(1) + } + this[0] == '-' -> { + sign = -1 + sPositive = this.substring(1) + } + else -> { + sPositive = this + sign = +1 + } + } + var res = BigInt.ZERO + var digitValue = BigInt.ONE + val sPositiveUpper = sPositive.toUpperCase() + if (sPositiveUpper.startsWith("0X")) { // hex representation + val sHex = sPositiveUpper.substring(2) + for (ch in sHex.reversed()) { + if (ch == '_') continue + res += digitValue * (hexChToInt[ch] ?: return null) + digitValue *= 16.toBigInt() + } + } else { // decimal representation + for (ch in sPositiveUpper.reversed()) { + if (ch == '_') continue + if (ch !in '0'..'9') { + return null + } + res += digitValue * (ch.toInt() - '0'.toInt()) + digitValue *= 10.toBigInt() + } + } + return res * sign +} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt index 2a77e36ef..bd83932e7 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt @@ -29,7 +29,7 @@ fun >> ctg(arg: T): T = arg.conte /** * A context extension to include power operations like square roots, etc */ -interface PowerOperations { +interface PowerOperations : Algebra { fun power(arg: T, pow: Number): T fun sqrt(arg: T) = power(arg, 0.5) @@ -42,7 +42,7 @@ fun >> sqr(arg: T): T = arg pow 2.0 /* Exponential */ -interface ExponentialOperations { +interface ExponentialOperations: Algebra { fun exp(arg: T): T fun ln(arg: T): T } diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt index c63384fb9..f02fd8dd0 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt @@ -73,6 +73,8 @@ fun Buffer.asSequence(): Sequence = Sequence(::iterator) fun Buffer.asIterable(): Iterable = asSequence().asIterable() +val Buffer<*>.indices: IntRange get() = IntRange(0, size - 1) + interface MutableBuffer : Buffer { operator fun set(index: Int, value: T) diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt index 370dc646e..3437644ff 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt @@ -2,7 +2,6 @@ package scientifik.kmath.structures import scientifik.kmath.operations.* - interface ExtendedNDField> : NDField, TrigonometricOperations, diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt index de13c9888..7d1209963 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt @@ -1,26 +1,13 @@ package scientifik.kmath.linear import scientifik.kmath.structures.Matrix +import scientifik.kmath.structures.NDStructure +import scientifik.kmath.structures.as2D import kotlin.test.Test import kotlin.test.assertEquals class MatrixTest { - @Test - fun testSum() { - val vector1 = RealVector(5) { it.toDouble() } - val vector2 = RealVector(5) { 5 - it.toDouble() } - val sum = vector1 + vector2 - assertEquals(5.0, sum[2]) - } - - @Test - fun testVectorToMatrix() { - val vector = RealVector(5) { it.toDouble() } - val matrix = vector.asMatrix() - assertEquals(4.0, matrix[4, 0]) - } - @Test fun testTranspose() { val matrix = MatrixContext.real.one(3, 3) @@ -28,21 +15,6 @@ class MatrixTest { assertEquals(matrix, transposed) } - - @Test - fun testDot() { - val vector1 = RealVector(5) { it.toDouble() } - val vector2 = RealVector(5) { 5 - it.toDouble() } - - val matrix1 = vector1.asMatrix() - val matrix2 = vector2.asMatrix().transpose() - val product = MatrixContext.real.run { matrix1 dot matrix2 } - - - assertEquals(5.0, product[1, 0]) - assertEquals(6.0, product[2, 2]) - } - @Test fun testBuilder() { val matrix = Matrix.build(2, 3)( @@ -74,4 +46,20 @@ class MatrixTest { val toTenthPower = transitionMatrix pow 10 } + + @Test + fun test2DDot() { + val firstMatrix = NDStructure.auto(2,3){ (i, j) -> (i + j).toDouble() }.as2D() + val secondMatrix = NDStructure.auto(3,2){ (i, j) -> (i + j).toDouble() }.as2D() + MatrixContext.real.run { +// val firstMatrix = produce(2, 3) { i, j -> (i + j).toDouble() } +// val secondMatrix = produce(3, 2) { i, j -> (i + j).toDouble() } + val result = firstMatrix dot secondMatrix + assertEquals(2, result.rowNum) + assertEquals(2, result.colNum) + assertEquals(8.0, result[0,1]) + assertEquals(8.0, result[1,0]) + assertEquals(14.0, result[1,1]) + } + } } \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntAlgebraTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntAlgebraTest.kt new file mode 100644 index 000000000..5ae977196 --- /dev/null +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntAlgebraTest.kt @@ -0,0 +1,50 @@ +package scientifik.kmath.operations + +import kotlin.test.Test +import kotlin.test.assertEquals + +class BigIntAlgebraTest { + @Test + fun testKBigIntegerRingSum() { + val res = BigIntField { + 1_000L.toBigInt() * 1_000L.toBigInt() + } + assertEquals(res, 1_000_000.toBigInt()) + } + + @Test + fun testKBigIntegerRingSum_100_000_000__100_000_000() { + BigIntField { + val sum = +"100_000_000" + +"100_000_000" + assertEquals(sum, "200_000_000".parseBigInteger()) + } + } + + @Test + fun test_mul_3__4() { + BigIntField { + val prod = +"0x3000_0000_0000" * +"0x4000_0000_0000_0000_0000" + assertEquals(prod, "0xc00_0000_0000_0000_0000_0000_0000_0000".parseBigInteger()) + } + } + + @Test + fun test_div_big_1() { + BigIntField { + val res = +"1_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000" / + +"555_000_444_000_333_000_222_000_111_000_999_001" + assertEquals(res, +"1801800360360432432518919022699") + } + } + + @Test + fun test_rem_big_1() { + BigIntField { + val res = +"1_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000" % + +"555_000_444_000_333_000_222_000_111_000_999_001" + assertEquals(res, +"324121220440768000291647788404676301") + } + } + +} + diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntConstructorTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntConstructorTest.kt new file mode 100644 index 000000000..2af1b7e50 --- /dev/null +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntConstructorTest.kt @@ -0,0 +1,26 @@ +package scientifik.kmath.operations + +import kotlin.test.Test +import kotlin.test.assertEquals + +class BigIntConstructorTest { + @Test + fun testConstructorZero() { + assertEquals(0.toBigInt(), uintArrayOf().toBigInt(0)) + } + + @Test + fun testConstructor8() { + assertEquals( + 8.toBigInt(), + uintArrayOf(8U).toBigInt(1) + ) + } + + @Test + fun testConstructor_0xffffffffaL() { + val x = -0xffffffffaL.toBigInt() + val y = uintArrayOf(0xfffffffaU, 0xfU).toBigInt(-1) + assertEquals(x, y) + } +} \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntConversionsTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntConversionsTest.kt new file mode 100644 index 000000000..51b9509e0 --- /dev/null +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntConversionsTest.kt @@ -0,0 +1,43 @@ +package scientifik.kmath.operations + +import kotlin.test.Test +import kotlin.test.assertEquals + +@kotlin.ExperimentalUnsignedTypes +class BigIntConversionsTest { + @Test + fun testToString0x10() { + val x = 0x10.toBigInt() + assertEquals("0x10", x.toString()) + } + + @Test + fun testToString0x17ffffffd() { + val x = 0x17ffffffdL.toBigInt() + assertEquals("0x17ffffffd", x.toString()) + } + + @Test + fun testToString_0x17ead2ffffd() { + val x = -0x17ead2ffffdL.toBigInt() + assertEquals("-0x17ead2ffffd", x.toString()) + } + + @Test + fun testToString_0x17ead2ffffd11223344() { + val x = uintArrayOf(0x11223344U, 0xad2ffffdU, 0x17eU).toBigInt(-1) + assertEquals("-0x17ead2ffffd11223344", x.toString()) + } + + @Test + fun testFromString_0x17ead2ffffd11223344() { + val x = "0x17ead2ffffd11223344".parseBigInteger() + assertEquals("0x17ead2ffffd11223344", x.toString()) + } + + @Test + fun testFromString_7059135710711894913860() { + val x = "-7059135710711894913860".parseBigInteger() + assertEquals("-0x17ead2ffffd11223344", x.toString()) + } +} \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntOperationsTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntOperationsTest.kt new file mode 100644 index 000000000..72ac9f229 --- /dev/null +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntOperationsTest.kt @@ -0,0 +1,381 @@ +package scientifik.kmath.operations + +import kotlin.test.Test +import kotlin.test.assertEquals + +@kotlin.ExperimentalUnsignedTypes +class BigIntOperationsTest { + @Test + fun testPlus_1_1() { + val x = 1.toBigInt() + val y = 1.toBigInt() + + val res = x + y + val sum = 2.toBigInt() + + assertEquals(sum, res) + } + + @Test + fun testPlusBigNumbers() { + val x = 0x7fffffff.toBigInt() + val y = 0x7fffffff.toBigInt() + val z = 0x7fffffff.toBigInt() + + val res = x + y + z + val sum = uintArrayOf(0x7ffffffdU, 0x1U).toBigInt(1) + + assertEquals(sum, res) + } + + @Test + fun testUnaryMinus() { + val x = 1234.toBigInt() + val y = -1234.toBigInt() + assertEquals(-x, y) + } + + @Test + fun testMinus_2_1() { + val x = 2.toBigInt() + val y = 1.toBigInt() + + val res = x - y + val sum = 1.toBigInt() + + assertEquals(sum, res) + } + + @Test + fun testMinus__2_1() { + val x = -2.toBigInt() + val y = 1.toBigInt() + + val res = x - y + val sum = -3.toBigInt() + + assertEquals(sum, res) + } + + @Test + fun testMinus___2_1() { + val x = -2.toBigInt() + val y = 1.toBigInt() + + val res = -x - y + val sum = 1.toBigInt() + + assertEquals(sum, res) + } + + @Test + fun testMinusBigNumbers() { + val x = 12345.toBigInt() + val y = 0xffffffffaL.toBigInt() + + val res = x - y + val sum = -0xfffffcfc1L.toBigInt() + + assertEquals(sum, res) + } + + @Test + fun testMultiply_2_3() { + val x = 2.toBigInt() + val y = 3.toBigInt() + + val res = x * y + val prod = 6.toBigInt() + + assertEquals(prod, res) + } + + @Test + fun testMultiply__2_3() { + val x = -2.toBigInt() + val y = 3.toBigInt() + + val res = x * y + val prod = -6.toBigInt() + + assertEquals(prod, res) + } + + @Test + fun testMultiply_0xfff123_0xfff456() { + val x = 0xfff123.toBigInt() + val y = 0xfff456.toBigInt() + + val res = x * y + val prod = 0xffe579ad5dc2L.toBigInt() + + assertEquals(prod, res) + } + + @Test + fun testMultiplyUInt_0xfff123_0xfff456() { + val x = 0xfff123.toBigInt() + val y = 0xfff456U + + val res = x * y + val prod = 0xffe579ad5dc2L.toBigInt() + + assertEquals(prod, res) + } + + @Test + fun testMultiplyInt_0xfff123__0xfff456() { + val x = 0xfff123.toBigInt() + val y = -0xfff456 + + val res = x * y + val prod = -0xffe579ad5dc2L.toBigInt() + + assertEquals(prod, res) + } + + @Test + fun testMultiply_0xffffffff_0xffffffff() { + val x = 0xffffffffL.toBigInt() + val y = 0xffffffffL.toBigInt() + + val res = x * y + val prod = 0xfffffffe00000001UL.toBigInt() + + assertEquals(prod, res) + } + + @Test + fun test_shr_20() { + val x = 20.toBigInt() + assertEquals(10.toBigInt(), x shr 1) + } + + @Test + fun test_shl_20() { + val x = 20.toBigInt() + assertEquals(40.toBigInt(), x shl 1) + } + + @Test + fun test_shl_1_0() { + assertEquals( + BigInt.ONE, + BigInt.ONE shl 0 + ) + } + + @Test + fun test_shl_1_32() { + assertEquals( + 0x100000000UL.toBigInt(), + BigInt.ONE shl 32 + ) + } + + @Test + fun test_shl_1_33() { + assertEquals( + 0x200000000UL.toBigInt(), + BigInt.ONE shl 33 + ) + } + + @Test + fun test_shr_1_33_33() { + assertEquals( + BigInt.ONE, + (BigInt.ONE shl 33) shr 33 + ) + } + + @Test + fun test_shr_1_32() { + assertEquals( + BigInt.ZERO, + BigInt.ONE shr 32 + ) + } + + @Test + fun test_and_123_456() { + val x = 123.toBigInt() + val y = 456.toBigInt() + assertEquals(72.toBigInt(), x and y) + } + + @Test + fun test_or_123_456() { + val x = 123.toBigInt() + val y = 456.toBigInt() + assertEquals(507.toBigInt(), x or y) + } + + @Test + fun test_asd() { + assertEquals( + BigInt.ONE, + BigInt.ZERO or ((20.toBigInt() shr 4) and BigInt.ONE) + ) + } + + @Test + fun test_square_0x11223344U_0xad2ffffdU_0x17eU() { + val num = + uintArrayOf(0x11223344U, 0xad2ffffdU, 0x17eU).toBigInt(-1) + println(num) + val res = num * num + assertEquals( + res, + uintArrayOf(0xb0542a10U, 0xbbd85bc8U, 0x2a1fa515U, 0x5069e03bU, 0x23c09U).toBigInt(1) + ) + } + + @Test + fun testDivision_6_3() { + val x = 6.toBigInt() + val y = 3U + + val res = x / y + val div = 2.toBigInt() + + assertEquals(div, res) + } + + @Test + fun testBigDivision_6_3() { + val x = 6.toBigInt() + val y = 3.toBigInt() + + val res = x / y + val div = 2.toBigInt() + + assertEquals(div, res) + } + + @Test + fun testDivision_20__3() { + val x = 20.toBigInt() + val y = -3 + + val res = x / y + val div = -6.toBigInt() + + assertEquals(div, res) + } + + @Test + fun testBigDivision_20__3() { + val x = 20.toBigInt() + val y = -3.toBigInt() + + val res = x / y + val div = -6.toBigInt() + + assertEquals(div, res) + } + + @Test + fun testDivision_0xfffffffe00000001_0xffffffff() { + val x = 0xfffffffe00000001UL.toBigInt() + val y = 0xffffffffU + + val res = x / y + val div = 0xffffffffL.toBigInt() + + assertEquals(div, res) + } + + @Test + fun testBigDivision_0xfffffffeabcdef01UL_0xfffffffeabc() { + val res = 0xfffffffeabcdef01UL.toBigInt() / 0xfffffffeabc.toBigInt() + assertEquals(res, 0x100000.toBigInt()) + } + + @Test + fun testBigDivision_0xfffffffe00000001_0xffffffff() { + val x = 0xfffffffe00000001UL.toBigInt() + val y = 0xffffffffU.toBigInt() + + val res = x / y + val div = 0xffffffffL.toBigInt() + + assertEquals(div, res) + } + + @Test + fun testMod_20_3() { + val x = 20.toBigInt() + val y = 3 + + val res = x % y + val mod = 2 + + assertEquals(mod, res) + } + + @Test + fun testBigMod_20_3() { + val x = 20.toBigInt() + val y = 3.toBigInt() + + val res = x % y + val mod = 2.toBigInt() + + assertEquals(mod, res) + } + + @Test + fun testMod_0xfffffffe00000001_12345() { + val x = 0xfffffffe00000001UL.toBigInt() + val y = 12345 + + val res = x % y + val mod = 1980 + + assertEquals(mod, res) + } + + @Test + fun testBigMod_0xfffffffe00000001_12345() { + val x = 0xfffffffe00000001UL.toBigInt() + val y = 12345.toBigInt() + + val res = x % y + val mod = 1980.toBigInt() + + assertEquals(mod, res) + } + + @Test + fun testModPow_3_10_17() { + val x = 3.toBigInt() + val exp = 10.toBigInt() + val mod = 17.toBigInt() + + val res = 8.toBigInt() + + return assertEquals(res, x.modPow(exp, mod)) + } + + @Test + fun testModPowBigNumbers() { + val x = 0xfffffffeabcdef01UL.toBigInt() + val exp = 2.toBigInt() + val mod = 0xfffffffeabcUL.toBigInt() + + val res = 0xc2253cde01.toBigInt() + + return assertEquals(res, x.modPow(exp, mod)) + } + + @Test + fun testModBigNumbers() { + val x = 0xfffffffeabcdef01UL.toBigInt() + val mod = 0xfffffffeabcUL.toBigInt() + + val res = 0xdef01.toBigInt() + + return assertEquals(res, x % mod) + } +} \ No newline at end of file diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumers.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumbers.kt similarity index 100% rename from kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumers.kt rename to kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumbers.kt diff --git a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt index 8dba9e215..1c2872d17 100644 --- a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt +++ b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt @@ -16,7 +16,9 @@ package scientifik.kmath.chains +import kotlinx.coroutines.InternalCoroutinesApi import kotlinx.coroutines.flow.Flow +import kotlinx.coroutines.flow.FlowCollector import kotlinx.coroutines.sync.Mutex import kotlinx.coroutines.sync.withLock @@ -25,7 +27,7 @@ import kotlinx.coroutines.sync.withLock * A not-necessary-Markov chain of some type * @param R - the chain element type */ -interface Chain { +interface Chain: Flow { /** * Generate next value, changing state if needed */ @@ -36,14 +38,15 @@ interface Chain { */ fun fork(): Chain + @InternalCoroutinesApi + override suspend fun collect(collector: FlowCollector) { + kotlinx.coroutines.flow.flow { while (true) emit(next()) }.collect(collector) + } + companion object } -/** - * Chain as a coroutine flow. The flow emit affects chain state and vice versa - */ -fun Chain.flow(): Flow = kotlinx.coroutines.flow.flow { while (true) emit(next()) } fun Iterator.asChain(): Chain = SimpleChain { next() } fun Sequence.asChain(): Chain = iterator().asChain() @@ -127,6 +130,21 @@ fun Chain.map(func: suspend (T) -> R): Chain = object : Chain { override fun fork(): Chain = this@map.fork().map(func) } +/** + * [block] must be a pure function or at least not use external random variables, otherwise fork could be broken + */ +fun Chain.filter(block: (T) -> Boolean): Chain = object : Chain { + override suspend fun next(): T { + var next: T + do { + next = this@filter.next() + } while (!block(next)) + return next + } + + override fun fork(): Chain = this@filter.fork().filter(block) +} + /** * Map the whole chain */ diff --git a/kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/BufferFlowTest.kt b/kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/BufferFlowTest.kt index af2d6cd43..147f687f0 100644 --- a/kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/BufferFlowTest.kt +++ b/kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/BufferFlowTest.kt @@ -3,11 +3,12 @@ package scientifik.kmath.streaming import kotlinx.coroutines.* import kotlinx.coroutines.flow.asFlow import kotlinx.coroutines.flow.collect -import org.junit.Test +import org.junit.jupiter.api.Timeout import scientifik.kmath.coroutines.async import scientifik.kmath.coroutines.collect import scientifik.kmath.coroutines.mapParallel import java.util.concurrent.Executors +import kotlin.test.Test @ExperimentalCoroutinesApi @@ -17,7 +18,8 @@ class BufferFlowTest { val dispatcher = Executors.newFixedThreadPool(4).asCoroutineDispatcher() - @Test(timeout = 2000) + @Test + @Timeout(2000) fun map() { runBlocking { (1..20).asFlow().mapParallel( dispatcher) { @@ -31,7 +33,8 @@ class BufferFlowTest { } } - @Test(timeout = 2000) + @Test + @Timeout(2000) fun async() { runBlocking { (1..20).asFlow().async(dispatcher) { diff --git a/kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/RingBufferTest.kt b/kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/RingBufferTest.kt index 25af1f589..c14d1a26c 100644 --- a/kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/RingBufferTest.kt +++ b/kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/RingBufferTest.kt @@ -2,8 +2,8 @@ package scientifik.kmath.streaming import kotlinx.coroutines.flow.* import kotlinx.coroutines.runBlocking -import org.junit.Test import scientifik.kmath.structures.asSequence +import kotlin.test.Test import kotlin.test.assertEquals class RingBufferTest { diff --git a/kmath-dimensions/build.gradle.kts b/kmath-dimensions/build.gradle.kts new file mode 100644 index 000000000..dda6cd2f0 --- /dev/null +++ b/kmath-dimensions/build.gradle.kts @@ -0,0 +1,19 @@ +plugins { + id("scientifik.mpp") +} + +description = "A proof of concept module for adding typ-safe dimensions to structures" + +kotlin.sourceSets { + commonMain { + dependencies { + api(project(":kmath-core")) + } + } + + jvmMain{ + dependencies{ + api(kotlin("reflect")) + } + } +} \ No newline at end of file diff --git a/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt new file mode 100644 index 000000000..37e89c111 --- /dev/null +++ b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt @@ -0,0 +1,35 @@ +package scientifik.kmath.dimensions + +import kotlin.reflect.KClass + +/** + * An abstract class which is not used in runtime. Designates a size of some structure. + * Could be replaced later by fully inline constructs + */ +interface Dimension { + + val dim: UInt + companion object { + + } +} + +fun KClass.dim(): UInt = Dimension.resolve(this).dim + +expect fun Dimension.Companion.resolve(type: KClass): D + +expect fun Dimension.Companion.of(dim: UInt): Dimension + +inline fun Dimension.Companion.dim(): UInt = D::class.dim() + +object D1 : Dimension { + override val dim: UInt get() = 1U +} + +object D2 : Dimension { + override val dim: UInt get() = 2U +} + +object D3 : Dimension { + override val dim: UInt get() = 31U +} \ No newline at end of file diff --git a/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt new file mode 100644 index 000000000..f447866c0 --- /dev/null +++ b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt @@ -0,0 +1,161 @@ +package scientifik.kmath.dimensions + +import scientifik.kmath.linear.GenericMatrixContext +import scientifik.kmath.linear.MatrixContext +import scientifik.kmath.linear.Point +import scientifik.kmath.linear.transpose +import scientifik.kmath.operations.Ring +import scientifik.kmath.structures.Matrix +import scientifik.kmath.structures.Structure2D + +/** + * A matrix with compile-time controlled dimension + */ +interface DMatrix : Structure2D { + companion object { + /** + * Coerces a regular matrix to a matrix with type-safe dimensions and throws a error if coercion failed + */ + inline fun coerce(structure: Structure2D): DMatrix { + if (structure.rowNum != Dimension.dim().toInt()) { + error("Row number mismatch: expected ${Dimension.dim()} but found ${structure.rowNum}") + } + if (structure.colNum != Dimension.dim().toInt()) { + error("Column number mismatch: expected ${Dimension.dim()} but found ${structure.colNum}") + } + return DMatrixWrapper(structure) + } + + /** + * The same as [coerce] but without dimension checks. Use with caution + */ + fun coerceUnsafe(structure: Structure2D): DMatrix { + return DMatrixWrapper(structure) + } + } +} + +/** + * An inline wrapper for a Matrix + */ +inline class DMatrixWrapper( + val structure: Structure2D +) : DMatrix { + override val shape: IntArray get() = structure.shape + override fun get(i: Int, j: Int): T = structure[i, j] +} + +/** + * Dimension-safe point + */ +interface DPoint : Point { + companion object { + inline fun coerce(point: Point): DPoint { + if (point.size != Dimension.dim().toInt()) { + error("Vector dimension mismatch: expected ${Dimension.dim()}, but found ${point.size}") + } + return DPointWrapper(point) + } + + fun coerceUnsafe(point: Point): DPoint { + return DPointWrapper(point) + } + } +} + +/** + * Dimension-safe point wrapper + */ +inline class DPointWrapper(val point: Point) : + DPoint { + override val size: Int get() = point.size + + override fun get(index: Int): T = point[index] + + override fun iterator(): Iterator = point.iterator() +} + + +/** + * Basic operations on dimension-safe matrices. Operates on [Matrix] + */ +inline class DMatrixContext>(val context: GenericMatrixContext) { + + inline fun Matrix.coerce(): DMatrix { + if (rowNum != Dimension.dim().toInt()) { + error("Row number mismatch: expected ${Dimension.dim()} but found $rowNum") + } + if (colNum != Dimension.dim().toInt()) { + error("Column number mismatch: expected ${Dimension.dim()} but found $colNum") + } + return DMatrix.coerceUnsafe(this) + } + + /** + * Produce a matrix with this context and given dimensions + */ + inline fun produce(noinline initializer: (i: Int, j: Int) -> T): DMatrix { + val rows = Dimension.dim() + val cols = Dimension.dim() + return context.produce(rows.toInt(), cols.toInt(), initializer).coerce() + } + + inline fun point(noinline initializer: (Int) -> T): DPoint { + val size = Dimension.dim() + return DPoint.coerceUnsafe( + context.point( + size.toInt(), + initializer + ) + ) + } + + inline infix fun DMatrix.dot( + other: DMatrix + ): DMatrix { + return context.run { this@dot dot other }.coerce() + } + + inline infix fun DMatrix.dot(vector: DPoint): DPoint { + return DPoint.coerceUnsafe(context.run { this@dot dot vector }) + } + + inline operator fun DMatrix.times(value: T): DMatrix { + return context.run { this@times.times(value) }.coerce() + } + + inline operator fun T.times(m: DMatrix): DMatrix = + m * this + + + inline operator fun DMatrix.plus(other: DMatrix): DMatrix { + return context.run { this@plus + other }.coerce() + } + + inline operator fun DMatrix.minus(other: DMatrix): DMatrix { + return context.run { this@minus + other }.coerce() + } + + inline operator fun DMatrix.unaryMinus(): DMatrix { + return context.run { this@unaryMinus.unaryMinus() }.coerce() + } + + inline fun DMatrix.transpose(): DMatrix { + return context.run { (this@transpose as Matrix).transpose() }.coerce() + } + + /** + * A square unit matrix + */ + inline fun one(): DMatrix = produce { i, j -> + if (i == j) context.elementContext.one else context.elementContext.zero + } + + inline fun zero(): DMatrix = produce { _, _ -> + context.elementContext.zero + } + + companion object { + val real = DMatrixContext(MatrixContext.real) + } +} \ No newline at end of file diff --git a/kmath-dimensions/src/commonTest/kotlin/scientifik/dimensions/DMatrixContextTest.kt b/kmath-dimensions/src/commonTest/kotlin/scientifik/dimensions/DMatrixContextTest.kt new file mode 100644 index 000000000..74d20205c --- /dev/null +++ b/kmath-dimensions/src/commonTest/kotlin/scientifik/dimensions/DMatrixContextTest.kt @@ -0,0 +1,30 @@ +package scientifik.dimensions + +import scientifik.kmath.dimensions.D2 +import scientifik.kmath.dimensions.D3 +import scientifik.kmath.dimensions.DMatrixContext +import kotlin.test.Test + + +class DMatrixContextTest { + @Test + fun testDimensionSafeMatrix() { + val res = DMatrixContext.real.run { + val m = produce { i, j -> (i + j).toDouble() } + + //The dimension of `one()` is inferred from type + (m + one()) + } + } + + @Test + fun testTypeCheck() { + val res = DMatrixContext.real.run { + val m1 = produce { i, j -> (i + j).toDouble() } + val m2 = produce { i, j -> (i + j).toDouble() } + + //Dimension-safe addition + m1.transpose() + m2 + } + } +} \ No newline at end of file diff --git a/kmath-dimensions/src/jsMain/kotlin/scientifik/kmath/dimensions/dim.kt b/kmath-dimensions/src/jsMain/kotlin/scientifik/kmath/dimensions/dim.kt new file mode 100644 index 000000000..bbd580629 --- /dev/null +++ b/kmath-dimensions/src/jsMain/kotlin/scientifik/kmath/dimensions/dim.kt @@ -0,0 +1,22 @@ +package scientifik.kmath.dimensions + +import kotlin.reflect.KClass + +private val dimensionMap = hashMapOf( + 1u to D1, + 2u to D2, + 3u to D3 +) + +@Suppress("UNCHECKED_CAST") +actual fun Dimension.Companion.resolve(type: KClass): D { + return dimensionMap.entries.find { it.value::class == type }?.value as? D ?: error("Can't resolve dimension $type") +} + +actual fun Dimension.Companion.of(dim: UInt): Dimension { + return dimensionMap.getOrPut(dim) { + object : Dimension { + override val dim: UInt get() = dim + } + } +} \ No newline at end of file diff --git a/kmath-dimensions/src/jvmMain/kotlin/scientifik/kmath/dimensions/dim.kt b/kmath-dimensions/src/jvmMain/kotlin/scientifik/kmath/dimensions/dim.kt new file mode 100644 index 000000000..e8fe8f59b --- /dev/null +++ b/kmath-dimensions/src/jvmMain/kotlin/scientifik/kmath/dimensions/dim.kt @@ -0,0 +1,18 @@ +package scientifik.kmath.dimensions + +import kotlin.reflect.KClass + +actual fun Dimension.Companion.resolve(type: KClass): D{ + return type.objectInstance ?: error("No object instance for dimension class") +} + +actual fun Dimension.Companion.of(dim: UInt): Dimension{ + return when(dim){ + 1u -> D1 + 2u -> D2 + 3u -> D3 + else -> object : Dimension { + override val dim: UInt get() = dim + } + } +} \ No newline at end of file diff --git a/kmath-for-real/build.gradle.kts b/kmath-for-real/build.gradle.kts new file mode 100644 index 000000000..a8a8975bc --- /dev/null +++ b/kmath-for-real/build.gradle.kts @@ -0,0 +1,11 @@ +plugins { + id("scientifik.mpp") +} + +kotlin.sourceSets { + commonMain { + dependencies { + api(project(":kmath-core")) + } + } +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt similarity index 100% rename from kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt rename to kmath-for-real/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt diff --git a/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt new file mode 100644 index 000000000..d9c6b5eab --- /dev/null +++ b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt @@ -0,0 +1,146 @@ +package scientifik.kmath.real + +import scientifik.kmath.linear.MatrixContext +import scientifik.kmath.linear.RealMatrixContext.elementContext +import scientifik.kmath.linear.VirtualMatrix +import scientifik.kmath.operations.sum +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.Matrix +import scientifik.kmath.structures.asSequence +import kotlin.math.pow + +/* + * Functions for convenient "numpy-like" operations with Double matrices. + * + * Initial implementation of these functions is taken from: + * https://github.com/thomasnield/numky/blob/master/src/main/kotlin/org/nield/numky/linear/DoubleOperators.kt + * + */ + +/* + * Functions that help create a real (Double) matrix + */ + +fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double) = + MatrixContext.real.produce(rowNum, colNum, initializer) + +fun Sequence.toMatrix() = toList().let { + MatrixContext.real.produce(it.size, it[0].size) { row, col -> it[row][col] } +} + +fun Matrix.repeatStackVertical(n: Int) = VirtualMatrix(rowNum*n, colNum) { + row, col -> get(if (row == 0) 0 else row % rowNum, col) +} + +/* + * Operations for matrix and real number + */ + +operator fun Matrix.times(double: Double) = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row, col] * double +} + +operator fun Matrix.plus(double: Double) = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row, col] + double +} + +operator fun Matrix.minus(double: Double) = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row, col] - double +} + +operator fun Matrix.div(double: Double) = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row, col] / double +} + +operator fun Double.times(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { + row, col -> this * matrix[row, col] +} + +operator fun Double.plus(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { + row, col -> this * matrix[row, col] +} + +operator fun Double.minus(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { + row, col -> this - matrix[row, col] +} + +// TODO: does this operation make sense? Should it be 'this/matrix[row, col]'? +//operator fun Double.div(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { +// row, col -> matrix[row, col] / this +//} + +/* + * Per-element (!) square and power operations + */ + +fun Matrix.square() = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row, col].pow(2) +} + +fun Matrix.pow(n: Int) = MatrixContext.real.produce(rowNum, colNum) { + i, j -> this[i, j].pow(n) +} + +/* + * Operations on two matrices (per-element!) + */ + +operator fun Matrix.times(other: Matrix) = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row, col] * other[row, col] +} + +operator fun Matrix.plus(other: Matrix) = MatrixContext.real.add(this, other) + +operator fun Matrix.minus(other: Matrix) = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row,col] - other[row,col] +} + +/* + * Operations on columns + */ + +inline fun Matrix.appendColumn(crossinline mapper: (Buffer) -> Double) = + MatrixContext.real.produce(rowNum,colNum+1) { + row, col -> + if (col < colNum) + this[row, col] + else + mapper(rows[row]) + } + +fun Matrix.extractColumns(columnRange: IntRange) = MatrixContext.real.produce(rowNum, columnRange.count()) { + row, col -> this[row, columnRange.first + col] +} + +fun Matrix.extractColumn(columnIndex: Int) = extractColumns(columnIndex..columnIndex) + +fun Matrix.sumByColumn() = MatrixContext.real.produce(1, colNum) { _, j -> + val column = columns[j] + with(elementContext) { + sum(column.asSequence()) + } +} + +fun Matrix.minByColumn() = MatrixContext.real.produce(1, colNum) { + _, j -> columns[j].asSequence().min() ?: throw Exception("Cannot produce min on empty column") +} + +fun Matrix.maxByColumn() = MatrixContext.real.produce(1, colNum) { + _, j -> columns[j].asSequence().max() ?: throw Exception("Cannot produce min on empty column") +} + +fun Matrix.averageByColumn() = MatrixContext.real.produce(1, colNum) { + _, j -> columns[j].asSequence().average() +} + +/* + * Operations processing all elements + */ + +fun Matrix.sum() = elements().map { (_, value) -> value }.sum() + +fun Matrix.min() = elements().map { (_, value) -> value }.min() + +fun Matrix.max() = elements().map { (_, value) -> value }.max() + +fun Matrix.average() = elements().map { (_, value) -> value }.average() diff --git a/kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt b/kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt new file mode 100644 index 000000000..808794442 --- /dev/null +++ b/kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt @@ -0,0 +1,16 @@ +package scientific.kmath.real + +import scientifik.kmath.real.average +import scientifik.kmath.real.realMatrix +import scientifik.kmath.real.sum +import kotlin.test.Test +import kotlin.test.assertEquals + +class RealMatrixTest { + @Test + fun testSum() { + val m = realMatrix(10, 10) { i, j -> (i + j).toDouble() } + assertEquals(m.sum(), 900.0) + assertEquals(m.average(), 9.0) + } +} diff --git a/kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt b/kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt new file mode 100644 index 000000000..0e73ee4a5 --- /dev/null +++ b/kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt @@ -0,0 +1,36 @@ +package scientifik.kmath.linear + +import kotlin.test.Test +import kotlin.test.assertEquals + +class VectorTest { + @Test + fun testSum() { + val vector1 = RealVector(5) { it.toDouble() } + val vector2 = RealVector(5) { 5 - it.toDouble() } + val sum = vector1 + vector2 + assertEquals(5.0, sum[2]) + } + + @Test + fun testVectorToMatrix() { + val vector = RealVector(5) { it.toDouble() } + val matrix = vector.asMatrix() + assertEquals(4.0, matrix[4, 0]) + } + + @Test + fun testDot() { + val vector1 = RealVector(5) { it.toDouble() } + val vector2 = RealVector(5) { 5 - it.toDouble() } + + val matrix1 = vector1.asMatrix() + val matrix2 = vector2.asMatrix().transpose() + val product = MatrixContext.real.run { matrix1 dot matrix2 } + + + assertEquals(5.0, product[1, 0]) + assertEquals(6.0, product[2, 2]) + } + +} \ No newline at end of file diff --git a/kmath-for-real/src/jvmMain/kotlin/.gitkeep b/kmath-for-real/src/jvmMain/kotlin/.gitkeep new file mode 100644 index 000000000..e69de29bb diff --git a/kmath-functions/build.gradle.kts b/kmath-functions/build.gradle.kts new file mode 100644 index 000000000..4c158a32e --- /dev/null +++ b/kmath-functions/build.gradle.kts @@ -0,0 +1,11 @@ +plugins { + id("scientifik.mpp") +} + +kotlin.sourceSets { + commonMain { + dependencies { + api(project(":kmath-core")) + } + } +} diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Piecewise.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Piecewise.kt new file mode 100644 index 000000000..16f8aa12b --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Piecewise.kt @@ -0,0 +1,57 @@ +package scientifik.kmath.functions + +import scientifik.kmath.operations.Ring + +interface Piecewise { + fun findPiece(arg: T): R? +} + +interface PiecewisePolynomial : + Piecewise> + +/** + * Ordered list of pieces in piecewise function + */ +class OrderedPiecewisePolynomial>(delimeter: T) : + PiecewisePolynomial { + + private val delimiters: ArrayList = arrayListOf(delimeter) + private val pieces: ArrayList> = ArrayList() + + /** + * Dynamically add a piece to the "right" side (beyond maximum argument value of previous piece) + * @param right new rightmost position. If is less then current rightmost position, a error is thrown. + */ + fun putRight(right: T, piece: Polynomial) { + require(right > delimiters.last()) { "New delimiter should be to the right of old one" } + delimiters.add(right) + pieces.add(piece) + } + + fun putLeft(left: T, piece: Polynomial) { + require(left < delimiters.first()) { "New delimiter should be to the left of old one" } + delimiters.add(0, left) + pieces.add(0, piece) + } + + override fun findPiece(arg: T): Polynomial? { + if (arg < delimiters.first() || arg >= delimiters.last()) { + return null + } else { + for (index in 1 until delimiters.size) { + if (arg < delimiters[index]) { + return pieces[index - 1] + } + } + error("Piece not found") + } + } +} + +/** + * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise definition. + */ +fun , C : Ring> PiecewisePolynomial.value(ring: C, arg: T): T? = + findPiece(arg)?.value(ring, arg) + +fun , C : Ring> PiecewisePolynomial.asFunction(ring: C): (T) -> T? = { value(ring, it) } \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Polynomial.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Polynomial.kt new file mode 100644 index 000000000..b747b521d --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Polynomial.kt @@ -0,0 +1,73 @@ +package scientifik.kmath.functions + +import scientifik.kmath.operations.Ring +import scientifik.kmath.operations.Space +import kotlin.math.max +import kotlin.math.pow + +/** + * Polynomial coefficients without fixation on specific context they are applied to + * @param coefficients constant is the leftmost coefficient + */ +inline class Polynomial(val coefficients: List) { + constructor(vararg coefficients: T) : this(coefficients.toList()) +} + +fun Polynomial.value() = + coefficients.reduceIndexed { index: Int, acc: Double, d: Double -> acc + d.pow(index) } + + +fun > Polynomial.value(ring: C, arg: T): T = ring.run { + if (coefficients.isEmpty()) return@run zero + var res = coefficients.first() + var powerArg = arg + for (index in 1 until coefficients.size) { + res += coefficients[index] * powerArg + //recalculating power on each step to avoid power costs on long polynomials + powerArg *= arg + } + return@run res +} + +/** + * Represent a polynomial as a context-dependent function + */ +fun > Polynomial.asMathFunction(): MathFunction = object : + MathFunction { + override fun C.invoke(arg: T): T = value(this, arg) +} + +/** + * Represent the polynomial as a regular context-less function + */ +fun > Polynomial.asFunction(ring: C): (T) -> T = { value(ring, it) } + +/** + * An algebra for polynomials + */ +class PolynomialSpace>(val ring: C) : Space> { + + override fun add(a: Polynomial, b: Polynomial): Polynomial { + val dim = max(a.coefficients.size, b.coefficients.size) + ring.run { + return Polynomial(List(dim) { index -> + a.coefficients.getOrElse(index) { zero } + b.coefficients.getOrElse(index) { zero } + }) + } + } + + override fun multiply(a: Polynomial, k: Number): Polynomial { + ring.run { + return Polynomial(List(a.coefficients.size) { index -> a.coefficients[index] * k }) + } + } + + override val zero: Polynomial = + Polynomial(emptyList()) + + operator fun Polynomial.invoke(arg: T): T = value(ring, arg) +} + +fun , R> C.polynomial(block: PolynomialSpace.() -> R): R { + return PolynomialSpace(this).run(block) +} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/functions.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/functions.kt new file mode 100644 index 000000000..2b822b3ba --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/functions.kt @@ -0,0 +1,33 @@ +package scientifik.kmath.functions + +import scientifik.kmath.operations.Algebra +import scientifik.kmath.operations.RealField + +/** + * A regular function that could be called only inside specific algebra context + * @param T source type + * @param C source algebra constraint + * @param R result type + */ +interface MathFunction, R> { + operator fun C.invoke(arg: T): R +} + +fun MathFunction.invoke(arg: Double): R = RealField.invoke(arg) + +/** + * A suspendable function defined in algebraic context + */ +interface SuspendableMathFunction, R> { + suspend operator fun C.invoke(arg: T): R +} + +suspend fun SuspendableMathFunction.invoke(arg: Double) = RealField.invoke(arg) + + +/** + * A parametric function with parameter + */ +interface ParametricFunction> { + operator fun C.invoke(arg: T, parameter: P): T +} diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt new file mode 100644 index 000000000..8d83e4198 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt @@ -0,0 +1,45 @@ +package scientifik.kmath.interpolation + +import scientifik.kmath.functions.PiecewisePolynomial +import scientifik.kmath.functions.value +import scientifik.kmath.operations.Ring +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.asBuffer + +interface Interpolator { + fun interpolate(points: XYPointSet): (X) -> Y +} + +interface PolynomialInterpolator> : Interpolator { + val algebra: Ring + + fun getDefaultValue(): T = error("Out of bounds") + + fun interpolatePolynomials(points: XYPointSet): PiecewisePolynomial + + override fun interpolate(points: XYPointSet): (T) -> T = { x -> + interpolatePolynomials(points).value(algebra, x) ?: getDefaultValue() + } +} + +fun > PolynomialInterpolator.interpolatePolynomials( + x: Buffer, + y: Buffer +): PiecewisePolynomial { + val pointSet = BufferXYPointSet(x, y) + return interpolatePolynomials(pointSet) +} + +fun > PolynomialInterpolator.interpolatePolynomials( + data: Map +): PiecewisePolynomial { + val pointSet = BufferXYPointSet(data.keys.toList().asBuffer(), data.values.toList().asBuffer()) + return interpolatePolynomials(pointSet) +} + +fun > PolynomialInterpolator.interpolatePolynomials( + data: List> +): PiecewisePolynomial { + val pointSet = BufferXYPointSet(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer()) + return interpolatePolynomials(pointSet) +} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt new file mode 100644 index 000000000..98beb4391 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt @@ -0,0 +1,26 @@ +package scientifik.kmath.interpolation + +import scientifik.kmath.functions.OrderedPiecewisePolynomial +import scientifik.kmath.functions.PiecewisePolynomial +import scientifik.kmath.functions.Polynomial +import scientifik.kmath.operations.Field + +/** + * Reference JVM implementation: https://github.com/apache/commons-math/blob/master/src/main/java/org/apache/commons/math4/analysis/interpolation/LinearInterpolator.java + */ +class LinearInterpolator>(override val algebra: Field) : PolynomialInterpolator { + + override fun interpolatePolynomials(points: XYPointSet): PiecewisePolynomial = algebra.run { + require(points.size > 0) { "Point array should not be empty" } + insureSorted(points) + + OrderedPiecewisePolynomial(points.x[0]).apply { + for (i in 0 until points.size - 1) { + val slope = (points.y[i + 1] - points.y[i]) / (points.x[i + 1] - points.x[i]) + val const = points.y[i] - slope * points.x[i] + val polynomial = Polynomial(const, slope) + putRight(points.x[i + 1], polynomial) + } + } + } +} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt new file mode 100644 index 000000000..6707bd8bc --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt @@ -0,0 +1,296 @@ +package scientifik.kmath.interpolation +// +//import scientifik.kmath.functions.PiecewisePolynomial +//import scientifik.kmath.operations.Ring +//import scientifik.kmath.structures.Buffer +//import kotlin.math.abs +//import kotlin.math.sqrt +// +// +///** +// * Original code: https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/LoessInterpolator.java +// */ +//class LoessInterpolator>(override val algebra: Ring) : PolynomialInterpolator { +// /** +// * The bandwidth parameter: when computing the loess fit at +// * a particular point, this fraction of source points closest +// * to the current point is taken into account for computing +// * a least-squares regression. +// * +// * +// * A sensible value is usually 0.25 to 0.5. +// */ +// private var bandwidth = 0.0 +// +// /** +// * The number of robustness iterations parameter: this many +// * robustness iterations are done. +// * +// * +// * A sensible value is usually 0 (just the initial fit without any +// * robustness iterations) to 4. +// */ +// private var robustnessIters = 0 +// +// /** +// * If the median residual at a certain robustness iteration +// * is less than this amount, no more iterations are done. +// */ +// private var accuracy = 0.0 +// +// /** +// * Constructs a new [LoessInterpolator] +// * with a bandwidth of [.DEFAULT_BANDWIDTH], +// * [.DEFAULT_ROBUSTNESS_ITERS] robustness iterations +// * and an accuracy of {#link #DEFAULT_ACCURACY}. +// * See [.LoessInterpolator] for an explanation of +// * the parameters. +// */ +// fun LoessInterpolator() { +// bandwidth = DEFAULT_BANDWIDTH +// robustnessIters = DEFAULT_ROBUSTNESS_ITERS +// accuracy = DEFAULT_ACCURACY +// } +// +// fun LoessInterpolator(bandwidth: Double, robustnessIters: Int) { +// this(bandwidth, robustnessIters, DEFAULT_ACCURACY) +// } +// +// fun LoessInterpolator(bandwidth: Double, robustnessIters: Int, accuracy: Double) { +// if (bandwidth < 0 || +// bandwidth > 1 +// ) { +// throw OutOfRangeException(LocalizedFormats.BANDWIDTH, bandwidth, 0, 1) +// } +// this.bandwidth = bandwidth +// if (robustnessIters < 0) { +// throw NotPositiveException(LocalizedFormats.ROBUSTNESS_ITERATIONS, robustnessIters) +// } +// this.robustnessIters = robustnessIters +// this.accuracy = accuracy +// } +// +// fun interpolate( +// xval: DoubleArray, +// yval: DoubleArray +// ): PolynomialSplineFunction { +// return SplineInterpolator().interpolate(xval, smooth(xval, yval)) +// } +// +// fun XYZPointSet.smooth(): XYPointSet { +// checkAllFiniteReal(x) +// checkAllFiniteReal(y) +// checkAllFiniteReal(z) +// MathArrays.checkOrder(xval) +// if (size == 1) { +// return doubleArrayOf(y[0]) +// } +// if (size == 2) { +// return doubleArrayOf(y[0], y[1]) +// } +// val bandwidthInPoints = (bandwidth * size).toInt() +// if (bandwidthInPoints < 2) { +// throw NumberIsTooSmallException( +// LocalizedFormats.BANDWIDTH, +// bandwidthInPoints, 2, true +// ) +// } +// val res = DoubleArray(size) +// val residuals = DoubleArray(size) +// val sortedResiduals = DoubleArray(size) +// val robustnessWeights = DoubleArray(size) +// // Do an initial fit and 'robustnessIters' robustness iterations. +// // This is equivalent to doing 'robustnessIters+1' robustness iterations +// // starting with all robustness weights set to 1. +// Arrays.fill(robustnessWeights, 1.0) +// for (iter in 0..robustnessIters) { +// val bandwidthInterval = intArrayOf(0, bandwidthInPoints - 1) +// // At each x, compute a local weighted linear regression +// for (i in 0 until size) { +//// val x = x[i] +// // Find out the interval of source points on which +// // a regression is to be made. +// if (i > 0) { +// updateBandwidthInterval(x, z, i, bandwidthInterval) +// } +// val ileft = bandwidthInterval[0] +// val iright = bandwidthInterval[1] +// // Compute the point of the bandwidth interval that is +// // farthest from x +// val edge: Int +// edge = if (x[i] - x[ileft] > x[iright] - x[i]) { +// ileft +// } else { +// iright +// } +// // Compute a least-squares linear fit weighted by +// // the product of robustness weights and the tricube +// // weight function. +// // See http://en.wikipedia.org/wiki/Linear_regression +// // (section "Univariate linear case") +// // and http://en.wikipedia.org/wiki/Weighted_least_squares +// // (section "Weighted least squares") +// var sumWeights = 0.0 +// var sumX = 0.0 +// var sumXSquared = 0.0 +// var sumY = 0.0 +// var sumXY = 0.0 +// val denom: Double = abs(1.0 / (x[edge] - x[i])) +// for (k in ileft..iright) { +// val xk = x[k] +// val yk = y[k] +// val dist = if (k < i) x - xk else xk - x[i] +// val w = tricube(dist * denom) * robustnessWeights[k] * z[k] +// val xkw = xk * w +// sumWeights += w +// sumX += xkw +// sumXSquared += xk * xkw +// sumY += yk * w +// sumXY += yk * xkw +// } +// val meanX = sumX / sumWeights +// val meanY = sumY / sumWeights +// val meanXY = sumXY / sumWeights +// val meanXSquared = sumXSquared / sumWeights +// val beta: Double +// beta = if (sqrt(abs(meanXSquared - meanX * meanX)) < accuracy) { +// 0.0 +// } else { +// (meanXY - meanX * meanY) / (meanXSquared - meanX * meanX) +// } +// val alpha = meanY - beta * meanX +// res[i] = beta * x[i] + alpha +// residuals[i] = abs(y[i] - res[i]) +// } +// // No need to recompute the robustness weights at the last +// // iteration, they won't be needed anymore +// if (iter == robustnessIters) { +// break +// } +// // Recompute the robustness weights. +// // Find the median residual. +// // An arraycopy and a sort are completely tractable here, +// // because the preceding loop is a lot more expensive +// java.lang.System.arraycopy(residuals, 0, sortedResiduals, 0, size) +// Arrays.sort(sortedResiduals) +// val medianResidual = sortedResiduals[size / 2] +// if (abs(medianResidual) < accuracy) { +// break +// } +// for (i in 0 until size) { +// val arg = residuals[i] / (6 * medianResidual) +// if (arg >= 1) { +// robustnessWeights[i] = 0.0 +// } else { +// val w = 1 - arg * arg +// robustnessWeights[i] = w * w +// } +// } +// } +// return res +// } +// +// fun smooth(xval: DoubleArray, yval: DoubleArray): DoubleArray { +// if (xval.size != yval.size) { +// throw DimensionMismatchException(xval.size, yval.size) +// } +// val unitWeights = DoubleArray(xval.size) +// Arrays.fill(unitWeights, 1.0) +// return smooth(xval, yval, unitWeights) +// } +// +// /** +// * Given an index interval into xval that embraces a certain number of +// * points closest to `xval[i-1]`, update the interval so that it +// * embraces the same number of points closest to `xval[i]`, +// * ignoring zero weights. +// * +// * @param xval Arguments array. +// * @param weights Weights array. +// * @param i Index around which the new interval should be computed. +// * @param bandwidthInterval a two-element array {left, right} such that: +// * `(left==0 or xval[i] - xval[left-1] > xval[right] - xval[i])` +// * and +// * `(right==xval.length-1 or xval[right+1] - xval[i] > xval[i] - xval[left])`. +// * The array will be updated. +// */ +// private fun updateBandwidthInterval( +// xval: Buffer, weights: Buffer, +// i: Int, +// bandwidthInterval: IntArray +// ) { +// val left = bandwidthInterval[0] +// val right = bandwidthInterval[1] +// // The right edge should be adjusted if the next point to the right +// // is closer to xval[i] than the leftmost point of the current interval +// val nextRight = nextNonzero(weights, right) +// if (nextRight < xval.size && xval[nextRight] - xval[i] < xval[i] - xval[left]) { +// val nextLeft = nextNonzero(weights, bandwidthInterval[0]) +// bandwidthInterval[0] = nextLeft +// bandwidthInterval[1] = nextRight +// } +// } +// +// /** +// * Return the smallest index `j` such that +// * `j > i && (j == weights.length || weights[j] != 0)`. +// * +// * @param weights Weights array. +// * @param i Index from which to start search. +// * @return the smallest compliant index. +// */ +// private fun nextNonzero(weights: Buffer, i: Int): Int { +// var j = i + 1 +// while (j < weights.size && weights[j] == 0.0) { +// ++j +// } +// return j +// } +// +// /** +// * Compute the +// * [tricube](http://en.wikipedia.org/wiki/Local_regression#Weight_function) +// * weight function +// * +// * @param x Argument. +// * @return `(1 - |x|3)3` for |x| < 1, 0 otherwise. +// */ +// private fun tricube(x: Double): Double { +// val absX: Double = FastMath.abs(x) +// if (absX >= 1.0) { +// return 0.0 +// } +// val tmp = 1 - absX * absX * absX +// return tmp * tmp * tmp +// } +// +// /** +// * Check that all elements of an array are finite real numbers. +// * +// * @param values Values array. +// * @throws org.apache.commons.math4.exception.NotFiniteNumberException +// * if one of the values is not a finite real number. +// */ +// private fun checkAllFiniteReal(values: DoubleArray) { +// for (i in values.indices) { +// MathUtils.checkFinite(values[i]) +// } +// } +// +// override fun interpolatePolynomials(points: Collection>): PiecewisePolynomial { +// TODO("not implemented") //To change body of created functions use File | Settings | File Templates. +// } +// +// companion object { +// /** Default value of the bandwidth parameter. */ +// const val DEFAULT_BANDWIDTH = 0.3 +// +// /** Default value of the number of robustness iterations. */ +// const val DEFAULT_ROBUSTNESS_ITERS = 2 +// +// /** +// * Default value for accuracy. +// */ +// const val DEFAULT_ACCURACY = 1e-12 +// } +//} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/SplineInterpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/SplineInterpolator.kt new file mode 100644 index 000000000..e1af0c1a2 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/SplineInterpolator.kt @@ -0,0 +1,58 @@ +package scientifik.kmath.interpolation + +import scientifik.kmath.functions.OrderedPiecewisePolynomial +import scientifik.kmath.functions.PiecewisePolynomial +import scientifik.kmath.functions.Polynomial +import scientifik.kmath.operations.Field +import scientifik.kmath.structures.MutableBufferFactory + +/** + * Generic spline interpolator. Not recommended for performance critical places, use platform-specific and type specific ones. + * Based on https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/SplineInterpolator.java + */ +class SplineInterpolator>( + override val algebra: Field, + val bufferFactory: MutableBufferFactory +) : PolynomialInterpolator { + + //TODO possibly optimize zeroed buffers + + override fun interpolatePolynomials(points: XYPointSet): PiecewisePolynomial = algebra.run { + if (points.size < 3) { + error("Can't use spline interpolator with less than 3 points") + } + insureSorted(points) + + // Number of intervals. The number of data points is n + 1. + val n = points.size - 1 + // Differences between knot points + val h = bufferFactory(points.size) { i -> points.x[i + 1] - points.x[i] } + val mu = bufferFactory(points.size - 1) { zero } + val z = bufferFactory(points.size) { zero } + + for (i in 1 until n) { + val g = 2.0 * (points.x[i + 1] - points.x[i - 1]) - h[i - 1] * mu[i - 1] + mu[i] = h[i] / g + z[i] = + (3.0 * (points.y[i + 1] * h[i - 1] - points.x[i] * (points.x[i + 1] - points.x[i - 1]) + points.y[i - 1] * h[i]) / (h[i - 1] * h[i]) + - h[i - 1] * z[i - 1]) / g + } + + // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants) + + OrderedPiecewisePolynomial(points.x[points.size - 1]).apply { + var cOld = zero + for (j in n - 1 downTo 0) { + val c = z[j] - mu[j] * cOld + val a = points.y[j] + val b = (points.y[j + 1] - points.y[j]) / h[j] - h[j] * (cOld + 2.0 * c) / 3.0 + val d = (cOld - c) / (3.0 * h[j]) + val polynomial = Polynomial(a, b, c, d) + cOld = c + putLeft(points.x[j], polynomial) + } + } + + } + +} diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/XYPointSet.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/XYPointSet.kt new file mode 100644 index 000000000..d8e10b880 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/XYPointSet.kt @@ -0,0 +1,54 @@ +package scientifik.kmath.interpolation + +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.Structure2D + +interface XYPointSet { + val size: Int + val x: Buffer + val y: Buffer +} + +interface XYZPointSet : XYPointSet { + val z: Buffer +} + +internal fun > insureSorted(points: XYPointSet) { + for (i in 0 until points.size - 1) { + if (points.x[i + 1] <= points.x[i]) error("Input data is not sorted at index $i") + } +} + +class NDStructureColumn(val structure: Structure2D, val column: Int) : Buffer { + init { + require(column < structure.colNum) { "Column index is outside of structure column range" } + } + + override val size: Int get() = structure.rowNum + + override fun get(index: Int): T = structure[index, column] + + override fun iterator(): Iterator = sequence { + repeat(size) { + yield(get(it)) + } + }.iterator() +} + +class BufferXYPointSet(override val x: Buffer, override val y: Buffer) : XYPointSet { + init { + require(x.size == y.size) { "Sizes of x and y buffers should be the same" } + } + + override val size: Int + get() = x.size +} + +fun Structure2D.asXYPointSet(): XYPointSet { + require(shape[1] == 2) { "Structure second dimension should be of size 2" } + return object : XYPointSet { + override val size: Int get() = this@asXYPointSet.shape[0] + override val x: Buffer get() = NDStructureColumn(this@asXYPointSet, 0) + override val y: Buffer get() = NDStructureColumn(this@asXYPointSet, 1) + } +} \ No newline at end of file diff --git a/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt b/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt new file mode 100644 index 000000000..23acd835c --- /dev/null +++ b/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt @@ -0,0 +1,27 @@ +package scientifik.kmath.interpolation + +import scientifik.kmath.functions.PiecewisePolynomial +import scientifik.kmath.functions.asFunction +import scientifik.kmath.operations.RealField +import kotlin.test.Test +import kotlin.test.assertEquals + + +class LinearInterpolatorTest { + @Test + fun testInterpolation() { + val data = listOf( + 0.0 to 0.0, + 1.0 to 1.0, + 2.0 to 3.0, + 3.0 to 4.0 + ) + val polynomial: PiecewisePolynomial = LinearInterpolator(RealField).interpolatePolynomials(data) + val function = polynomial.asFunction(RealField) + + assertEquals(null, function(-1.0)) + assertEquals(0.5, function(0.5)) + assertEquals(2.0, function(1.5)) + assertEquals(3.0, function(2.0)) + } +} \ No newline at end of file diff --git a/kmath-geometry/build.gradle.kts b/kmath-geometry/build.gradle.kts new file mode 100644 index 000000000..39aa833ad --- /dev/null +++ b/kmath-geometry/build.gradle.kts @@ -0,0 +1,9 @@ +plugins { + id("scientifik.mpp") +} + +kotlin.sourceSets.commonMain { + dependencies { + api(project(":kmath-core")) + } +} \ No newline at end of file diff --git a/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Euclidean2DSpace.kt b/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Euclidean2DSpace.kt new file mode 100644 index 000000000..2313b2170 --- /dev/null +++ b/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Euclidean2DSpace.kt @@ -0,0 +1,58 @@ +package scientifik.kmath.geometry + +import scientifik.kmath.linear.Point +import scientifik.kmath.operations.SpaceElement +import scientifik.kmath.operations.invoke +import kotlin.math.sqrt + + +interface Vector2D : Point, Vector, SpaceElement { + val x: Double + val y: Double + + override val size: Int get() = 2 + + override fun get(index: Int): Double = when (index) { + 1 -> x + 2 -> y + else -> error("Accessing outside of point bounds") + } + + override fun iterator(): Iterator = listOf(x, y).iterator() + + override val context: Euclidean2DSpace get() = Euclidean2DSpace + + override fun unwrap(): Vector2D = this + + override fun Vector2D.wrap(): Vector2D = this +} + +val Vector2D.r: Double get() = Euclidean2DSpace.run { sqrt(norm()) } + +@Suppress("FunctionName") +fun Vector2D(x: Double, y: Double): Vector2D = Vector2DImpl(x, y) + +private data class Vector2DImpl( + override val x: Double, + override val y: Double +) : Vector2D + +/** + * 2D Euclidean space + */ +object Euclidean2DSpace : GeometrySpace { + fun Vector2D.norm(): Double = sqrt(x * x + y * y) + + override fun Vector2D.distanceTo(other: Vector2D): Double = (this - other).norm() + + override fun add(a: Vector2D, b: Vector2D): Vector2D = + Vector2D(a.x + b.x, a.y + b.y) + + override fun multiply(a: Vector2D, k: Number): Vector2D = + Vector2D(a.x * k.toDouble(), a.y * k.toDouble()) + + override val zero: Vector2D = Vector2D(0.0, 0.0) + + override fun Vector2D.dot(other: Vector2D): Double = + x * other.x + y * other.y +} \ No newline at end of file diff --git a/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Euclidean3DSpace.kt b/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Euclidean3DSpace.kt new file mode 100644 index 000000000..dd1776342 --- /dev/null +++ b/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Euclidean3DSpace.kt @@ -0,0 +1,57 @@ +package scientifik.kmath.geometry + +import scientifik.kmath.linear.Point +import scientifik.kmath.operations.SpaceElement +import kotlin.math.sqrt + + +interface Vector3D : Point, Vector, SpaceElement { + val x: Double + val y: Double + val z: Double + + override val size: Int get() = 3 + + override fun get(index: Int): Double = when (index) { + 1 -> x + 2 -> y + 3 -> z + else -> error("Accessing outside of point bounds") + } + + override fun iterator(): Iterator = listOf(x, y, z).iterator() + + override val context: Euclidean3DSpace get() = Euclidean3DSpace + + override fun unwrap(): Vector3D = this + + override fun Vector3D.wrap(): Vector3D = this +} + +@Suppress("FunctionName") +fun Vector3D(x: Double, y: Double, z: Double): Vector3D = Vector3DImpl(x, y, z) + +val Vector3D.r: Double get() = Euclidean3DSpace.run { sqrt(norm()) } + +private data class Vector3DImpl( + override val x: Double, + override val y: Double, + override val z: Double +) : Vector3D + +object Euclidean3DSpace : GeometrySpace { + override val zero: Vector3D = Vector3D(0.0, 0.0, 0.0) + + fun Vector3D.norm(): Double = sqrt(x * x + y * y + z * z) + + override fun Vector3D.distanceTo(other: Vector3D): Double = (this - other).norm() + + override fun add(a: Vector3D, b: Vector3D): Vector3D = + Vector3D(a.x + b.x, a.y + b.y, a.z + b.z) + + override fun multiply(a: Vector3D, k: Number): Vector3D = + Vector3D(a.x * k.toDouble(), a.y * k.toDouble(), a.z * k.toDouble()) + + override fun Vector3D.dot(other: Vector3D): Double = + x * other.x + y * other.y + z * other.z +} \ No newline at end of file diff --git a/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/GeometrySpace.kt b/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/GeometrySpace.kt new file mode 100644 index 000000000..b65a8dd3a --- /dev/null +++ b/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/GeometrySpace.kt @@ -0,0 +1,17 @@ +package scientifik.kmath.geometry + +import scientifik.kmath.operations.Space + +interface Vector + +interface GeometrySpace: Space { + /** + * L2 distance + */ + fun V.distanceTo(other: V): Double + + /** + * Scalar product + */ + infix fun V.dot(other: V): Double +} \ No newline at end of file diff --git a/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Line.kt b/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Line.kt new file mode 100644 index 000000000..d802a103f --- /dev/null +++ b/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Line.kt @@ -0,0 +1,6 @@ +package scientifik.kmath.geometry + +data class Line(val base: V, val direction: V) + +typealias Line2D = Line +typealias Line3D = Line diff --git a/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/ReferenceFrame.kt b/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/ReferenceFrame.kt new file mode 100644 index 000000000..420e38ce2 --- /dev/null +++ b/kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/ReferenceFrame.kt @@ -0,0 +1,4 @@ +package scientifik.kmath.geometry + +interface ReferenceFrame { +} \ No newline at end of file diff --git a/kmath-histograms/build.gradle.kts b/kmath-histograms/build.gradle.kts index 39aa833ad..993bfed8e 100644 --- a/kmath-histograms/build.gradle.kts +++ b/kmath-histograms/build.gradle.kts @@ -5,5 +5,6 @@ plugins { kotlin.sourceSets.commonMain { dependencies { api(project(":kmath-core")) + api(project(":kmath-for-real")) } } \ No newline at end of file diff --git a/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt b/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt index a2f0fc516..85f078fda 100644 --- a/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt +++ b/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt @@ -45,8 +45,7 @@ class RealHistogram( private val values: NDStructure = NDStructure.auto(strides) { LongCounter() } - //private val weight: NDStructure = ndStructure(strides){null} - + private val weights: NDStructure = NDStructure.auto(strides) { DoubleCounter() } override val dimension: Int get() = lower.size @@ -102,21 +101,33 @@ class RealHistogram( return MultivariateBin(getDef(index), getValue(index)) } +// fun put(point: Point){ +// val index = getIndex(point) +// values[index].increment() +// } + override fun putWithWeight(point: Buffer, weight: Double) { - if (weight != 1.0) TODO("Implement weighting") val index = getIndex(point) values[index].increment() + weights[index].add(weight) } - override fun iterator(): Iterator> = values.elements().map { (index, value) -> + override fun iterator(): Iterator> = weights.elements().map { (index, value) -> MultivariateBin(getDef(index), value.sum()) }.iterator() /** * Convert this histogram into NDStructure containing bin values but not bin descriptions */ - fun asNDStructure(): NDStructure { - return NDStructure.auto(this.values.shape) { values[it].sum() } + fun values(): NDStructure { + return NDStructure.auto(values.shape) { values[it].sum() } + } + + /** + * Sum of weights + */ + fun weights():NDStructure{ + return NDStructure.auto(weights.shape) { weights[it].sum() } } companion object { diff --git a/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt b/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt index 30a87228a..df2e9847a 100644 --- a/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt +++ b/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt @@ -1,6 +1,10 @@ package scientifik.memory import java.nio.ByteBuffer +import java.nio.channels.FileChannel +import java.nio.file.Files +import java.nio.file.Path +import java.nio.file.StandardOpenOption /** @@ -11,7 +15,7 @@ actual fun Memory.Companion.allocate(length: Int): Memory { return ByteBufferMemory(buffer) } -class ByteBufferMemory( +private class ByteBufferMemory( val buffer: ByteBuffer, val startOffset: Int = 0, override val size: Int = buffer.limit() @@ -90,4 +94,13 @@ class ByteBufferMemory( } override fun writer(): MemoryWriter = writer +} + +/** + * Use direct memory-mapped buffer from file to read something and close it afterwards. + */ +fun Path.readAsMemory(position: Long = 0, size: Long = Files.size(this), block: Memory.() -> R): R { + return FileChannel.open(this, StandardOpenOption.READ).use { + ByteBufferMemory(it.map(FileChannel.MapMode.READ_ONLY, position, size)).block() + } } \ No newline at end of file diff --git a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Distribution.kt b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Distribution.kt index 44af84d5d..3b874adaa 100644 --- a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Distribution.kt +++ b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Distribution.kt @@ -58,9 +58,13 @@ fun Sampler.sampleBuffer( //creating temporary storage once val tmp = ArrayList(size) return sample(generator).collect { chain -> - for (i in tmp.indices) { - tmp[i] = chain.next() + //clear list from previous run + tmp.clear() + //Fill list + repeat(size){ + tmp.add(chain.next()) } + //return new buffer with elements from tmp bufferFactory(size) { tmp[it] } } } diff --git a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/RandomChain.kt b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/RandomChain.kt index 9b581afd7..47fc6e4c5 100644 --- a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/RandomChain.kt +++ b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/RandomChain.kt @@ -11,5 +11,4 @@ class RandomChain(val generator: RandomGenerator, private val gen: suspen override fun fork(): Chain = RandomChain(generator.fork(), gen) } -fun RandomGenerator.chain(gen: suspend RandomGenerator.() -> R): RandomChain = RandomChain(this, gen) -fun RandomGenerator.flow(gen: suspend RandomGenerator.() -> R) = chain(gen).fork() \ No newline at end of file +fun RandomGenerator.chain(gen: suspend RandomGenerator.() -> R): RandomChain = RandomChain(this, gen) \ No newline at end of file diff --git a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Statistic.kt b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Statistic.kt index 1af2570b0..804aed089 100644 --- a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Statistic.kt +++ b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Statistic.kt @@ -11,6 +11,7 @@ import scientifik.kmath.coroutines.mapParallel import scientifik.kmath.operations.* import scientifik.kmath.structures.Buffer import scientifik.kmath.structures.asIterable +import scientifik.kmath.structures.asSequence /** * A function, that transforms a buffer of random quantities to some resulting value @@ -83,9 +84,9 @@ class Mean(val space: Space) : ComposableStatistic, T> { /** * Non-composable median */ -class Median(comparator: Comparator) : Statistic { +class Median(private val comparator: Comparator) : Statistic { override suspend fun invoke(data: Buffer): T { - return data.asIterable().toList()[data.size / 2] //TODO check if this is correct + return data.asSequence().sortedWith(comparator).toList()[data.size / 2] //TODO check if this is correct } companion object { diff --git a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/distributions.kt b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/UniformDistribution.kt similarity index 100% rename from kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/distributions.kt rename to kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/UniformDistribution.kt diff --git a/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/SamplerTest.kt b/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/SamplerTest.kt new file mode 100644 index 000000000..1152f3057 --- /dev/null +++ b/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/SamplerTest.kt @@ -0,0 +1,17 @@ +package scientifik.kmath.prob + +import kotlinx.coroutines.runBlocking +import kotlin.test.Test + +class SamplerTest { + + @Test + fun bufferSamplerTest(){ + val sampler: Sampler = + BasicSampler { it.chain { nextDouble() } } + val data = sampler.sampleBuffer(RandomGenerator.default, 100) + runBlocking { + println(data.next()) + } + } +} \ No newline at end of file diff --git a/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/StatisticTest.kt b/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/StatisticTest.kt index 8b46cac52..af069810f 100644 --- a/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/StatisticTest.kt +++ b/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/StatisticTest.kt @@ -3,7 +3,7 @@ package scientifik.kmath.prob import kotlinx.coroutines.flow.drop import kotlinx.coroutines.flow.first import kotlinx.coroutines.runBlocking -import scientifik.kmath.chains.flow + import scientifik.kmath.streaming.chunked import kotlin.test.Test @@ -13,7 +13,7 @@ class StatisticTest { //Create a stateless chain from generator. val data = generator.chain { nextDouble() } //Convert a chaint to Flow and break it into chunks. - val chunked = data.flow().chunked(1000) + val chunked = data.chunked(1000) @Test fun testParallelMean() { diff --git a/kmath-viktor/build.gradle.kts b/kmath-viktor/build.gradle.kts new file mode 100644 index 000000000..52ee7c497 --- /dev/null +++ b/kmath-viktor/build.gradle.kts @@ -0,0 +1,10 @@ +plugins { + id("scientifik.jvm") +} + +description = "Binding for https://github.com/JetBrains-Research/viktor" + +dependencies { + api(project(":kmath-core")) + api("org.jetbrains.bio:viktor:1.0.1") +} \ No newline at end of file diff --git a/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorBuffer.kt b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorBuffer.kt new file mode 100644 index 000000000..040eee951 --- /dev/null +++ b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorBuffer.kt @@ -0,0 +1,20 @@ +package scientifik.kmath.viktor + +import org.jetbrains.bio.viktor.F64FlatArray +import scientifik.kmath.structures.MutableBuffer + +@Suppress("NOTHING_TO_INLINE", "OVERRIDE_BY_INLINE") +inline class ViktorBuffer(val flatArray: F64FlatArray) : MutableBuffer { + override val size: Int get() = flatArray.size + + override inline fun get(index: Int): Double = flatArray[index] + override inline fun set(index: Int, value: Double) { + flatArray[index] = value + } + + override fun copy(): MutableBuffer { + return ViktorBuffer(flatArray.copy().flatten()) + } + + override fun iterator(): Iterator = flatArray.data.iterator() +} \ No newline at end of file diff --git a/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt new file mode 100644 index 000000000..84e927721 --- /dev/null +++ b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt @@ -0,0 +1,86 @@ +package scientifik.kmath.viktor + +import org.jetbrains.bio.viktor.F64Array +import scientifik.kmath.operations.RealField +import scientifik.kmath.structures.DefaultStrides +import scientifik.kmath.structures.MutableNDStructure +import scientifik.kmath.structures.NDField + +@Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE") +inline class ViktorNDStructure(val f64Buffer: F64Array) : MutableNDStructure { + + override val shape: IntArray get() = f64Buffer.shape + + override inline fun get(index: IntArray): Double = f64Buffer.get(*index) + + override inline fun set(index: IntArray, value: Double) { + f64Buffer.set(*index, value = value) + } + + override fun elements(): Sequence> { + return DefaultStrides(shape).indices().map { it to get(it) } + } +} + +fun F64Array.asStructure(): ViktorNDStructure = ViktorNDStructure(this) + +@Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE") +class ViktorNDField(override val shape: IntArray) : NDField { + override val zero: ViktorNDStructure + get() = F64Array.full(init = 0.0, shape = *shape).asStructure() + override val one: ViktorNDStructure + get() = F64Array.full(init = 1.0, shape = *shape).asStructure() + + val strides = DefaultStrides(shape) + + override val elementContext: RealField get() = RealField + + override fun produce(initializer: RealField.(IntArray) -> Double): ViktorNDStructure = F64Array(*shape).apply { + this@ViktorNDField.strides.indices().forEach { index -> + set(value = RealField.initializer(index), indices = *index) + } + }.asStructure() + + override fun map(arg: ViktorNDStructure, transform: RealField.(Double) -> Double): ViktorNDStructure = + F64Array(*shape).apply { + this@ViktorNDField.strides.indices().forEach { index -> + set(value = RealField.transform(arg[index]), indices = *index) + } + }.asStructure() + + override fun mapIndexed( + arg: ViktorNDStructure, + transform: RealField.(index: IntArray, Double) -> Double + ): ViktorNDStructure = F64Array(*shape).apply { + this@ViktorNDField.strides.indices().forEach { index -> + set(value = RealField.transform(index, arg[index]), indices = *index) + } + }.asStructure() + + override fun combine( + a: ViktorNDStructure, + b: ViktorNDStructure, + transform: RealField.(Double, Double) -> Double + ): ViktorNDStructure = F64Array(*shape).apply { + this@ViktorNDField.strides.indices().forEach { index -> + set(value = RealField.transform(a[index], b[index]), indices = *index) + } + }.asStructure() + + override fun add(a: ViktorNDStructure, b: ViktorNDStructure): ViktorNDStructure { + return (a.f64Buffer + b.f64Buffer).asStructure() + } + + override fun multiply(a: ViktorNDStructure, k: Number): ViktorNDStructure = + (a.f64Buffer * k.toDouble()).asStructure() + + override inline fun ViktorNDStructure.plus(b: ViktorNDStructure): ViktorNDStructure = + (f64Buffer + b.f64Buffer).asStructure() + + override inline fun ViktorNDStructure.minus(b: ViktorNDStructure): ViktorNDStructure = + (f64Buffer - b.f64Buffer).asStructure() + + override inline fun ViktorNDStructure.times(k: Number): ViktorNDStructure = (f64Buffer * k.toDouble()).asStructure() + + override inline fun ViktorNDStructure.plus(arg: Double): ViktorNDStructure = (f64Buffer.plus(arg)).asStructure() +} \ No newline at end of file diff --git a/settings.gradle.kts b/settings.gradle.kts index f52edfac5..a08d5f7ee 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -1,10 +1,10 @@ pluginManagement { plugins { - id("scientifik.mpp") version "0.2.5" - id("scientifik.jvm") version "0.2.5" - id("scientifik.atomic") version "0.2.5" - id("scientifik.publish") version "0.2.5" + id("scientifik.mpp") version "0.4.1" + id("scientifik.jvm") version "0.4.1" + id("scientifik.atomic") version "0.4.1" + id("scientifik.publish") version "0.4.1" } repositories { @@ -13,13 +13,14 @@ pluginManagement { gradlePluginPortal() maven("https://dl.bintray.com/kotlin/kotlin-eap") maven("https://dl.bintray.com/mipt-npm/scientifik") + maven("https://dl.bintray.com/mipt-npm/dev") maven("https://dl.bintray.com/kotlin/kotlinx") } resolutionStrategy { eachPlugin { when (requested.id.id) { - "scientifik.mpp", "scientifik.publish" -> useModule("scientifik:gradle-tools:${requested.version}") + "scientifik.mpp", "scientifik.jvm", "scientifik.publish" -> useModule("scientifik:gradle-tools:${requested.version}") } } } @@ -29,12 +30,17 @@ rootProject.name = "kmath" include( ":kmath-memory", ":kmath-core", + ":kmath-functions", // ":kmath-io", ":kmath-coroutines", ":kmath-histograms", ":kmath-commons", + ":kmath-viktor", ":kmath-koma", ":kmath-prob", ":kmath-io", + ":kmath-dimensions", + ":kmath-for-real", + ":kmath-geometry", ":examples" )